re LargeIntegers and sqrt,etc

The method for square roots does not work with LargeIntegers, it gives a
Float overflow error, a similar error occurs with ln as well.

        I have written a work round for sqrt that will work with
LargeIntegers, but I have not got one for ln.

    Is there a solution to this problem?


    Ron

Ronald Hallam wrote:

>         I have written a work round for sqrt that will work with
> LargeIntegers, but I have not got one for ln.
>
>     Is there a solution to this problem?

The appended package overrides the #log and #ln methods of LargeInteger and
Fraction to avoid creating unecessarily overflowing intermediate Floats.

A similar approach could be used for #sqrt.

NOTE: this has NOT been extensively tested, I just hacked it up to
demonstrate the idea.

Try:
    (2 raisedToInteger: 5000) log: 2.

>     Ron

    -- chris


========= CU Safer Maths.pac ==========
| package |
package := Package name: 'CU Safer LargeInteger'.
package paxVersion: 0;
 basicComment: 'This adds a few methods to LargeInteger and Fraction which
eliminate some avoidable floating point overflows when working with very
high precision Integers and Fractions.

Note that *most* such overflows are *not* avoidable!!

 -- chris ([hidden email])'.

package basicPackageVersion: ''.

"Add the package scripts"

"Add the class names, loose method names, global names, resource names"
package classNames
 yourself.

package methodNames
 add: #Fraction -> #ln;
 add: #Fraction -> #log;
 add: #LargeInteger -> #ln;
 add: #LargeInteger -> #log;
 add: #LargeInteger -> #safeLog2;
 yourself.

package globalNames
 yourself.

package resourceNames
 yourself.

"Binary Global Names"
package binaryGlobalNames: (Set new
 yourself).
"Resource Names"
package allResourceNames: (Set new
 yourself).

"Add the prerequisite names"
package setPrerequisites: (IdentitySet new
 add: 'Dolphin';
 yourself).

package!

"Class Definitions"!

"Loose Methods"!

!Fraction methodsFor!

ln
 "Answer a <Float> which is the natural logarithm of the receiver."
#CUadded.

 "overridden to avoid overflowing intermediate Floats"
 ^ numerator ln - denominator ln.
!

log
 "Answer the log to the base 10 of the receiver"

#CUadded.

 "overridden to avoid overflowing intermediate Floats"
 ^ numerator log - denominator log.
! !
!Fraction categoriesFor: #ln!*-CU-added!mathematical!public! !
!Fraction categoriesFor: #log!*-CU-added!mathematical!public! !

!LargeInteger methodsFor!

ln
 "Answer a <Float> which is the natural logarithm of the receiver."
#CUadded.

 "overridden to correct for overflowing intermediate Floats"
 [^ super ln]
  on: FloatingPointException
  do: [:error | error isOverflow
    ifTrue: [^ ##(2 ln) * self safeLog2]
    ifFalse: [error pass]].!

log
 "Answer the log to the base 10 of the receiver"
#CUadded.

 "overridden to correct for overflowing intermediate Floats"
 [^ super log]
  on: FloatingPointException
  do: [:error | error isOverflow
    ifTrue: [^ ##(2 log) * self safeLog2]
    ifFalse: [error pass]].
!

safeLog2
 "private - Answer a <Float> which is the log to the base 2 of the
receiver."
 "This uses an algorithm which avoids floating point overflow by
bit-twiddling
 with the underlying representation.  Note that we assume that we
 are already known to be > 0"
 | source digits ignore approx scale |
#CUadded.

 source := self normalize.
 digits := source digitSize.

 "we use at most 64bits of precision and will scale by the number of digits
we ignored"
 ignore := digits - 8 max: 1.
 scale := ignore - 1 * 8.

 "construct approximation"
 approx := 0.0.
 digits to: ignore by: -1 do: [ :i | approx := approx * 256.0 + (source
byteAt: i) asFloat].

 "now we can do the scaling by *addition* after taking the log"
 ^ (approx log: 2) + scale.
! !
!LargeInteger categoriesFor: #ln!*-CU-added!mathematical!public! !
!LargeInteger categoriesFor: #log!*-CU-added!mathematical!public! !
!LargeInteger categoriesFor: #safeLog2!*-CU-added!mathematical!private! !

"End of package definition"!


"Binary Globals"!

"Resources"!

Chris
    Thanks for the attached, I have got  a the sqrt working all  right as
all I needed for that was the nearest integer.

sqrVal1: aNumber raisedTo: aPower
 |sNum p w num numHold diff |
 sNum := aNumber sqrt.
 p := aPower//2.
 num := LargeInteger new64.
 1 to: (((sNum printString)size) - 1) do: [:x | (p > 0) ifTrue:[w := (sNum *
10)//10.
      sNum := (sNum*10) - (w*10).
      num := num + (w * (10 raisedTo: p)).
      p := p -1]].
 diff := self - (num raisedTo: 2).
 [(diff > ((2 * num) + 1)) &  (p > 0)] whileTrue:[w := 0.
       numHold:= num.
      [self > (numHold raisedTo:2)] whileTrue:[
      w := w + 1.
      numHold := num + (w * (10 raisedTo: p))].
      w := w -1.
      num := num + (w * (10 raisedTo: p)).
      diff := self - (num raisedTo: 2).
      p := p -1].
 diff := self - (num raisedTo: 2).
 [diff > ((2 * num) + 1)] whileTrue:[num := num + 1.
     diff := self - (num raisedTo: 2)].
 ^num

sqrVal: aNumber raisedTo: aPower
 | a p|
 ((aPower\\2)=0) ifTrue:[^self sqrVal1:aNumber raisedTo: aPower]
    ifFalse:[^self sqrVal1:(aNumber*10) raisedTo: aPower - 1]

lsqrt
 | a aPower |
 aPower := (self printString) size.
 a := self // (1* (10 raisedTo: aPower -1)).
 ^self sqrVal: a raisedTo: aPower-1


Ron

 I have added a similar method to lsqrt to smallInteger, to return an
integer instead of a Float.
Chris Uppal wrote in message
<[hidden email]>...

Ron,

>     Thanks for the attached, I have got  a the sqrt working all  right as
> all I needed for that was the nearest integer.

[implementation snipped]

I'm afraid that your implementation has problems.  It can be *very*
inefficient, for instance on this machine it takes about 45 seconds to
evaluate

    (2 raisedToInteger: 48) sqrt.

I does depend on the number, though.  In this case the numbers on either
side of 2^48 also take most of a minute.  I think too, that there's a bad
relation between the size of the number and the execution time, it seems to
get slow roughly in proportion to the size of the number, rather than the
number of bits in its representation.

Another problem is that it gives off-by-one answers for some powers of two.
E.g. 2^34 sqrt comes out to be one less than 2^17

One way to avoid this is to use my earlier package and then override
LargeInteger>>sqrt:

========================
sqrt
 "Answer a <number> which is the positive square root of the receiver."
#CUadded.

 "overridden to correct for overflowing intermediate Floats"
 [^ super sqrt]
  on: FloatingPointException
  do: [:error | error isOverflow
    ifTrue: [^ 2 raisedTo: (self safeLog2 / 2)]
    ifFalse: [error pass]].
========================


You may as well do Fraction too:
========================
sqrt
 "Answer a <number> which is the positive square root of the receiver."
#CUadded.

 "overridden to avoid overflowing intermediate Floats"
 ^ numerator sqrt / denominator sqrt.
========================

That runs quickly (a handful of microseconds).

However, that has the problem that it can be inaccurate for LargeIntegers
even when an exact answer is possible.

Another problem is that it is attempting to return a Float, so it won't work
if the square root is too large to be representable.   E.g. try:

     (10 raisedTo: 1000) sqrt.

(Blair, if you're reading this, that expression answers "1.#INF" instead of
throwing a FloatingPointError.  I think that's an error and I'll post a
separate bug report).

So I think what's really needed is a separate method which calculates the
integer square root (the largest integer such that its square is <= the
input number).

There's standard iterative formula for calculating square roots, I think
it's the called "Newton-Ralphson" (but the days when I could remember such
things are long gone).

The implementation doesn't depend on any properties of LargeInteger, so I've
defined it against Number:

========================
isqrt
 "Answer the <integer> square root of the receiver"
 | x try  new |
#CUadded.

 self < 0 ifTrue: [ArithmeticError signal: 'Square root of negative
number'].
 self < 1 ifTrue: [^ 0].
 self < 3 ifTrue: [^ 1]. "iteration doesn't work on numbers < 2"
 x := self.
 try := self // 2.
 [new := (try + (x // try)) // 2. new < try] whileTrue: [try := new].
 ^ try.
========================

It runs pretty fast.  It'll do 2^48 in about 1/2 a millisecond, and even
handles (16 * 10^1000) in less than a second.

The problem with that is that I'm not mathematician enough to feel happy
with the convergence of the iteration, nor am I confident that the integer
rounding (the uses of #//) aren't making it inaccurate.

So, the last formulation is the one I like best.  It's a bit slower than the
Newton-Ralphson (by a constant factor), but I can understand it better and I
feel *much* more confident of it.  It also can be defined against Number:

========================
isqrt2
 "Answer the <integer> square root of the receiver"
 | digits result |
#CUadded.

 "this amounts to a binary chop across the space of all positive Integers
 with at most half as many bits as the receiver"

 self < 0 ifTrue: [ArithmeticError signal: 'Square root of negative
number'].

 "how many binary digits can root have ?"
 digits := 0.
 [(1 bitShift: digits) < self] whileTrue:  [digits := digits + 1].
 digits := digits // 2.

 "try setting each bit of the root in turn, starting at the high bit, if
doing so doesn't make the square too high, then it should be set
in the final result"
 result := 0.
 digits to: 0 by: -1 do:
  [:place || try |
  try := result bitOr: (1 bitShift: place).
  try * try > self ifFalse: [result := try]].

 ^ result.
========================

I'm sure there are even faster ways of doing these things, but the above
will do for me for now.

I hope this helped.  If not then I've had fun trying things out anyway...

> Ron

    -- chris