Missing bignum sqrt
Affects  Status  Importance  Assigned to  Milestone  

 Ikarus Scheme 
High

Abdulaziz Ghuloum 
Bug Description
I stumbled upon this error while factorizing 567348343347364
Unhandled exception:
Condition components:
1. &assertion
2. &who: sqrt
3. &message: "BUG: bignum sqrt not implemented"
Abdulaziz Ghuloum (aghuloum) wrote :  #1 
Changed in ikarus:  
assignee:  nobody → aghuloum 
importance:  Undecided → High 
status:  New → Confirmed 
Is it possible that you actually wanted exactintegersqrt?
> (exactintegersqrt 567348343347364
2381907519924
3148279521588
> (+ (* 2381907519924 2381907519924) 3148279521588)
567348343347364
>
Jens Axel Søgaard (soegaard) wrote :  #3 
Abdulaziz Ghuloum wrote:
> Is it possible that you actually wanted exactintegersqrt?
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 integersqrt is used here.
Grepping for sqrt reveals the following use, where
it is not possible to replace sqrt with integersqrt.
(define (second
; 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 ([sqrtd (sqrt d)])
(list (/ ( ( b) sqrtd) (* 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
; return list of integer solutions to a x^2 + b x + c = 0
(filter natural? (second
(define (triangle? n)
(not (null? (second
(define (pentagonal? n)
(not (null? (second
Hmm. I suppose one could use the following definition of s instead
of sqrt.
(define (square n)
(* n n))
(define (s n)
(let ([sqrtn (integersqrt n)])
(if (= (square sqrtn) n)
sqrtn
(sqrt (exact>inexact n)))))
/Jens Axel
Abdulaziz Ghuloum (aghuloum) wrote :  #4 
On Jan 3, 2008, at 5:23 PM, Jens Axel Søgaard wrote:
> (define (square n)
> (* n n))
>
> (define (s n)
> (let ([sqrtn (integersqrt n)])
> (if (= (square sqrtn) n)
> sqrtn
> (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 notsolarge 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 :  #5 
Fixed (kind of) in revision 1330.
> (* 248957894759872
619800333632677
> (sqrt 619800333632677
248957894759872
> (sqrt 619800333632677
7.8727398892169
> (* 7.8727398892169
6.1980033363267
> (inexact 619800333632677
6.1980033363267
> (sqrt (expt 2 2000))
107150860718626
This is an interesting case:
> (sqrt (+ (expt 2 2000) 1))
1.0715086071862
> (inexact (expt 2 1000))
1.0715086071862
Changed in ikarus:  
status:  Confirmed → Fix Committed 
Abdulaziz Ghuloum wrote:
> Fixed (kind of) in revision 1330.
Thanks!
/Jens Axel
Abdulaziz Ghuloum (aghuloum) wrote :  #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 :  #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
(let loop ([i 1] [x 2])
(display i) (newline)
(if (= (exact>inexact x) +inf.0)
x
(loop (+ i 1) (* 2 x)))))
(let ([x (first
(display (list (sqrt x)
Ikarus:
(13407807929942
1.3407807929942
But the others don't ...
MzScheme:
(13407807929942
+inf.0)
Gambit:
134078079299425
+inf.0
Larceny (0.95)
(+inf.0 +inf.0)
Except Chez Petite:
(13407807929942
1.3407807929942
/Jens Axel
Abdulaziz Ghuloum (aghuloum) wrote :  #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 :  #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 :  #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 20080105)
Copyright (c) 20062007 Abdulaziz Ghuloum
> (log (expt 2 1024))
+inf.0
80:~/Scheme soegaard$ petite
Petite Chez Scheme Version 7.3
Copyright (c) 19852006 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) 20042007 PLT Scheme Inc.
> (log (expt 2 1024))
709.782712893384
/Jens Axel
Jens Axel Søgaard (soegaard) wrote :  #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
Copyright (c) 19852006 Cadence Research Systems
> (sqrt (+ 1 (expt 5/2 1400)))
3.6141491434385
80:~/Scheme soegaard$ ikarus
Ikarus Scheme version 0.0.2patched+ (revision 1330, build 20080105)
Copyright (c) 20062007 Abdulaziz Ghuloum
> (sqrt (+ 1 (expt 5/2 1400)))
+inf.0
Tricky.
/Jens Axel
Abdulaziz Ghuloum (aghuloum) wrote :  #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 :  #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 :  #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 :  #16 
Fixed a roundoff 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 firstunreprese
(let loop ((x 2))
(if (= (exact>inexact x) +inf.0)
x
(loop (* 2 x)))))
(define (exactint.sqrt x)
(callwithvalues (lambda () (exactintegersqrt x)) cons))
(let* ((sfu
(sqrt firstunreprese
(isfu
(naisfu ;; nextafter...
(let loop ((increment 1.))
(if (= (+ isfu increment)
(loop (* increment 2.))
(+ isfu increment))))
(asfu ;; exact average
(quotient (+ sfu
(y (+ (square asfu) 1))
(sy (exactint.sqrt y)))
(prettyprint (list isfu
Abdulaziz Ghuloum (aghuloum) wrote :  #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
candidate tarball is available for download from:
http://
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.)
Changed in ikarus:  
milestone:  none → 0.0.3 
Changed in ikarus:  
status:  Fix Committed → Fix Released 
Will fix. Thanks.