Missing bignum sqrt

[Jens Axel Søgaard]

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)]

Will fix. Thanks.

[Abdulaziz Ghuloum (aghuloum)]

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)]

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)]

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)]

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

[Jens Axel Søgaard (soegaard)]

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

Thanks!

/Jens Axel

[Abdulaziz Ghuloum (aghuloum)]

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)]

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)]

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)]

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)]

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)
Copyright (c) 2006-2007 Abdulaziz Ghuloum

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

80:~/Scheme soegaard$ petite
Petite Chez Scheme Version 7.3
Copyright (c) 1985-2006 Cadence Research Systems

 > (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)]

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
Copyright (c) 1985-2006 Cadence Research Systems

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

80:~/Scheme soegaard$ ikarus
Ikarus Scheme version 0.0.2patched+ (revision 1330, build 2008-01-05)
Copyright (c) 2006-2007 Abdulaziz Ghuloum

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

Tricky.

/Jens Axel

[Abdulaziz Ghuloum (aghuloum)]

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)]

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)]

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)]

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)]

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
candidate tarball is available for download from:
http://www.cs.indiana.edu/~aghuloum/ikarus/ikarus-0.0.3-rc1.tar.gz
Please do test it if you have the time and report any issues
you might encounter. Thank you very much for your support.
(Sorry for the duplicates; I'm updating every open bug.)