# Missing bignum sqrt

Bug #180170 reported by Jens Axel Søgaard on 2008-01-03
Affects Status Importance Assigned to Milestone
High

### Bug Description

I stumbled upon this error while factorizing 5673483433473648736487364.

Unhandled exception:
Condition components:
1. &assertion
2. &who: sqrt
3. &message: "BUG: bignum sqrt not implemented"

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-03: #1

Will fix. Thanks.

 Changed in ikarus: assignee: nobody → aghuloum importance: Undecided → High status: New → Confirmed
 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-03: Re: [Bug 180170] Missing bignum sqrt #2

Is it possible that you actually wanted exact-integer-sqrt?

> (exact-integer-sqrt 5673483433473648736487364)
2381907519924
3148279521588
> (+ (* 2381907519924 2381907519924) 3148279521588)
5673483433473648736487364
>

 Jens Axel Søgaard (soegaard) wrote on 2008-01-03: #3

Abdulaziz Ghuloum wrote:
> Is it possible that you actually wanted exact-integer-sqrt?

In the case of the factorization function sqrt was used
used to break out of a loop early, so it doesn't
matter which one of sqrt and integer-sqrt is used here.

Grepping for sqrt reveals the following use, where
it is not possible to replace sqrt with integer-sqrt.

(define (second-degree-solutions a b c)
; return list of solutions to a x^2 + b x + c = 0
(let ([d (- (* b b) (* 4 a c))])
(cond
[(< d 0)
'()]
[(= d 0)
(list (/ b (* -2 a)))]
[else
(let ([sqrt-d (sqrt d)])
(list (/ (- (- b) sqrt-d) (* 2 a))
(/ (+ (- b) sqrt-d) (* 2 a))))])))

It is used like this to examine whether a number is
triangular, pentagonal, etc.

(define (natural? n)
(and (integer? n) (positive? n)))

(define (second-degree-natural-solutions a b c)
; return list of integer solutions to a x^2 + b x + c = 0
(filter natural? (second-degree-solutions a b c)))

(define (triangle? n)
(not (null? (second-degree-natural-solutions 1/2 1/2 (- n)))))

(define (pentagonal? n)
(not (null? (second-degree-natural-solutions 3/2 -1/2 (- n)))))

Hmm. I suppose one could use the following definition of s instead
of sqrt.

(define (square n)
(* n n))

(define (s n)
(let ([sqrt-n (integer-sqrt n)])
(if (= (square sqrt-n) n)
sqrt-n
(sqrt (exact->inexact n)))))

/Jens Axel

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-03: #4

On Jan 3, 2008, at 5:23 PM, Jens Axel Søgaard wrote:

> (define (square n)
> (* n n))
>
> (define (s n)
> (let ([sqrt-n (integer-sqrt n)])
> (if (= (square sqrt-n) n)
> sqrt-n
> (sqrt (exact->inexact n)))))

Unfortunately, (sqrt (inexact n)) != (sqrt n) for large values of n
since
(sqrt (inexact N)) ;;; N is big
== (sqrt +inf.0)
== +inf.0

And even for not-so-large values of N, it might be possible to obtain
more accurate results than what that formula yields.

Doesn't sound so easy now, so, it might actually take some time :-(

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-05: #5

Fixed (kind of) in revision 1330.

> (* 248957894759872489574895 248957894759872489574895)
61980033363267746003611829355008910457814261025
> (sqrt 61980033363267746003611829355008910457814261025)
248957894759872489574895
> (sqrt 619800333632677460036118293550089104578142610251)
7.872739889216952e23
> (* 7.872739889216952e23 7.872739889216952e23)
6.198003336326774e47
> (inexact 619800333632677460036118293550089104578142610251)
6.198003336326774e47
> (sqrt (expt 2 2000))
10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

This is an interesting case:
> (sqrt (+ (expt 2 2000) 1))
1.0715086071862673e301
> (inexact (expt 2 1000))
1.0715086071862673e301

 Changed in ikarus: status: Confirmed → Fix Committed
 Jens Axel Søgaard (soegaard) wrote on 2008-01-05: Re: [Bug 180170] Re: Missing bignum sqrt #6

Abdulaziz Ghuloum wrote:
> Fixed (kind of) in revision 1330.

Thanks!

/Jens Axel

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-05: #7

On Jan 5, 2008, at 6:48 AM, Jens Axel Søgaard wrote:

> Thanks!

You're most welcome. Let me know if you ever come across a case
where you notice it's losing precision. I tried my best, failed, and
did not bother proving its optimal correctness. :-)

 Jens Axel Søgaard (soegaard) wrote on 2008-01-05: #8

Abdulaziz Ghuloum wrote:
> On Jan 5, 2008, at 6:48 AM, Jens Axel Søgaard wrote:
>
>> Thanks!
>
> You're most welcome. Let me know if you ever come across a case
> where you notice it's losing precision. I tried my best, failed, and
> did not bother proving its optimal correctness. :-)

Do you handle the case:

x a bignum, where (exact->inexact x) is +inf.0,
but where (exact->inexact (sqrt x)) is representable as a
floating point

Yes you do!

(let ()
(define (first-unrepresentable)
(let loop ([i 1] [x 2])
(display i) (newline)
(if (= (exact->inexact x) +inf.0)
x
(loop (+ i 1) (* 2 x)))))

(let ([x (first-unrepresentable)])
(display (list (sqrt x)
(sqrt (+ x 1))))))

Ikarus:

(13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096
1.3407807929942597e154)

But the others don't ...

MzScheme:

(13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096
+inf.0)

Gambit:

13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096
+inf.0

Larceny (0.95)
(+inf.0 +inf.0)

Except Chez Petite:

(13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084096
1.3407807929942597e154)

/Jens Axel

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-05: #9

On Jan 5, 2008, at 7:56 AM, Jens Axel Søgaard wrote:

> Except Chez Petite:

Didn't I tell you before that I would bet my life on Chez? Now you
should believe me.

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-05: #10

On Jan 5, 2008, at 7:56 AM, Jens Axel Søgaard wrote:

> (+inf.0 +inf.0)

And that's not even funny!

PS. No, I'm not sending a bug report due to the sad state of my
previous attempts at reporting bugs there.

 Jens Axel Søgaard (soegaard) wrote on 2008-01-05: #11

Abdulaziz Ghuloum wrote:
> On Jan 5, 2008, at 7:56 AM, Jens Axel Søgaard wrote:
>
>> Except Chez Petite:
>
> Didn't I tell you before that I would bet my life on Chez? Now you
> should believe me.

Beginning to :-)

80:~/Scheme soegaard\$ ikarus
Ikarus Scheme version 0.0.2patched+ (revision 1330, build 2008-01-05)

> (log (expt 2 1024))
+inf.0

80:~/Scheme soegaard\$ petite
Petite Chez Scheme Version 7.3

> (log (expt 2 1024))
709.782712893384

80:~/Scheme soegaard\$ gsi
Gambit v4.1.2

> (log (expt 2 1024))
+inf.0

80:~/Scheme soegaard\$ MzScheme
Welcome to MzScheme v371.3 [3m], Copyright (c) 2004-2007 PLT Scheme Inc.
> (log (expt 2 1024))
709.782712893384

/Jens Axel

 Jens Axel Søgaard (soegaard) wrote on 2008-01-05: #12

Abdulaziz Ghuloum wrote:
> On Jan 5, 2008, at 7:56 AM, Jens Axel Søgaard wrote:
>
>> Except Chez Petite:
>
> Didn't I tell you before that I would bet my life on Chez? Now you
> should believe me.
>

Brad Lucier found the following example:

80:~/Scheme soegaard\$ petite
Petite Chez Scheme Version 7.3

> (sqrt (+ 1 (expt 5/2 1400)))
3.614149143438584e278

80:~/Scheme soegaard\$ ikarus
Ikarus Scheme version 0.0.2patched+ (revision 1330, build 2008-01-05)

> (sqrt (+ 1 (expt 5/2 1400)))
+inf.0

Tricky.

/Jens Axel

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-05: #13

On Jan 5, 2008, at 11:34 AM, Jens Axel Søgaard wrote:
>
>> (sqrt (+ 1 (expt 5/2 1400)))
> +inf.0

:-)

Will fix some day.

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-06: #14

On Jan 5, 2008, at 11:27 AM, Jens Axel Søgaard wrote:

>
>> (log (expt 2 1024))
> +inf.0

This just got fixed in rev 1331.

log is broken for rationals (updated bug 162334 with this regard)

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-06: #15

On Jan 5, 2008, at 11:34 AM, Jens Axel Søgaard wrote:

>> (sqrt (+ 1 (expt 5/2 1400)))
> +inf.0

Tried fixing that today and failed. I updated bug 162334 with this.
I must say, producing good inexact results from exact arguments is
very hard.

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-13: #16

Fixed a round-off error in bignum sqrt in the overflow case in revision 1341. The following code is derived from bug #24 in gambit's tracker.

(import (ikarus))

(define (square x) (* x x))

(define first-unrepresentable
(let loop ((x 2))
(if (= (exact->inexact x) +inf.0)
x
(loop (* 2 x)))))

(define (exact-int.sqrt x)
(call-with-values (lambda () (exact-integer-sqrt x)) cons))

(let* ((s-f-u
(sqrt first-unrepresentable))
(i-s-f-u
(exact->inexact s-f-u))
(n-a-i-s-f-u ;; next-after-...
(let loop ((increment 1.))
(if (= (+ i-s-f-u increment)
i-s-f-u)
(loop (* increment 2.))
(+ i-s-f-u increment))))
(a-s-f-u ;; exact average
(quotient (+ s-f-u
(inexact->exact n-a-i-s-f-u))
2))
(y (+ (square a-s-f-u) 1))
(s-y (exact-int.sqrt y)))
(pretty-print (list i-s-f-u
n-a-i-s-f-u
(sqrt y)
(number->string (inexact->exact n-a-i-s-f-u) 2)
(number->string a-s-f-u 2)
(number->string (car s-y) 2)
(number->string (cdr s-y) 2))))

 Abdulaziz Ghuloum (aghuloum) wrote on 2008-01-29: #17

This bug report is about to be closed as the fix comitted
previously will be incorporated in the next 0.0.3 release of
Ikarus Scheme, scheduled for January 31, 2008. A release