#4148 isqrt 'make check' failure

obsolete: 8.6a2
closed-accepted
5
2008-10-05
2008-10-02
Anonymous
No

'make check' produces -

==== expr-47.12 isqrt of various sizes of integer FAILED
==== Contents of test case:

set faults 0
for {set i 0} {$faults < 10 && $i <= 1024} {incr i} {
set root [expr {1 << $i}]
set rm1 [expr {$root - 1}]
set arg [expr {1 << (2 * $i)}]
set tval [expr {isqrt($arg-1)}]
if {$tval != $rm1} {
append trouble "i = " $i ": isqrt(" $arg "-1) == " $tval " != " $rm1 "\n"
incr faults
}
set tval [expr {isqrt($arg)}]
if {$tval != $root} {
append trouble "i = " $i ": isqrt(" $arg ") == " $tval " != " $root "\n"
incr faults
}
set tval [expr {isqrt($arg+1)}]
if {$tval != $root} {
append trouble "i = " $i ": isqrt(" $arg "+1) == " $tval " != " $root "\n"
incr faults
}
}
set trouble

---- Result was:
i = 28: isqrt(72057594037927936-1) == 279620267 != 268435455
i = 28: isqrt(72057594037927936+1) == 279620267 != 268435456
i = 29: isqrt(288230376151711744-1) == 545818760 != 536870911
i = 29: isqrt(288230376151711744+1) == 545818760 != 536870912
i = 30: isqrt(1152921504606846976-1) == 1079707056 != 1073741823
i = 30: isqrt(1152921504606846976+1) == 1079707056 != 1073741824
i = 31: isqrt(4611686018427387904-1) == 2150992608 != 2147483647
i = 31: isqrt(4611686018427387904+1) == 2150992608 != 2147483648
i = 32: isqrt(18446744073709551616-1) == 4296881274 != 4294967295
i = 41: isqrt(4835703278458516698824704-1) == 2199023259647 != 2199023255551

---- Result should have been (exact matching):

==== expr-47.12 FAILED

This is with gcc-4.2.1, -O2. Without -O2, the test passes.

I believe the cause is in bn_mp_sqrt.c; the following patch fixes it.

--- tcl8.6a2/libtommath/bn_mp_sqrt.c 2006-12-01 05:48:23.000000000 +0000
+++ tcl8.6a2-patched/libtommath/bn_mp_sqrt.c 2008-10-01 00:02:55.000000000 +0100
@@ -69,7 +69,7 @@
if (dig) {
t1.used = i+2;
d -= ldexp((double) dig, DIGIT_BIT);
- if (d != 0.0) {
+ if (d >= 1.0) {
t1.dp[i+1] = dig;
t1.dp[i] = ((mp_digit) d) - 1;
} else {

The reason is that d can become small, much less than 1.0 but not exactly 0.0 - in that case t1 is not set to a valid representation.
t1.dp[i] is -1, outside the range of a 28bit mp_digit.
Without -O2, d is exactly 0.0 for the first failing case, hence t1 is valid and the final result is correct.
This could be due to intermediate results being rounded to 64 bit double, but with -O2 everything is kept as 80 bit long double.

With attached code extracted from bn_mp_isqrt.c, I get
% make CC=/usr/local/gcc-4.2.1/bin/gcc CFLAGS=-O2 LDFLAGS=-lm isqrt && ./isqrt
/usr/local/gcc-4.2.1/bin/gcc -O2 -lm isqrt.c -o isqrt
d=-1.86265e-09=-0x1p-29 dig=1

System is
Linux piglet 2.4.36 #1 Tue Jan 22 01:12:23 GMT 2008 i686 unknown unknown GNU/Linux

Malc.
malcolm.boffey@virgin.net

Discussion

  • Nobody/Anonymous

    extract from bn_mp_isqrt.c

     
  • Malcolm Boffey

    Malcolm Boffey - 2008-10-03

    I see that this issue has come up before - there are several closed tickets that look like exactly the
    same issue (2030480 1609936 (which in turn mentions 1735583)).
    Compiling without optimisation (or debugging prints) gets expr.test to pass, however there are
    still values that fail, for example (1<<56) + (1<<5) -

    malc@piglet:~/src/tcl/tcl8.6a2/unix> ./tcltest ../tests/expr.test
    expr.test: Total 1749 Passed 1743 Skipped 6 Failed 0
    Number of tests skipped for each constraint:
    5 !ieeeFloatingPoint
    1 longIs64bit
    malc@piglet:~/src/tcl/tcl8.6a2/unix> ./tcltest ~/isqrt.tcl
    d=100000000000020 q=10aaaaab qq=115c71c7ce38e39
    malc@piglet:~/src/tcl/tcl8.6a2/unix> tclsh8.6 ~/isqrt.tcl
    d=100000000000020 q=10000000 qq=100000000000000

    (where ./tcltest is unmodified compiled with -O0, tclsh8.6 has the above path applied).
    isqrt.tcl -

    set d [expr (1<<56) + (1<<5)]
    set q [expr isqrt($d)]
    set qq [expr $q*$q]
    puts [format "d=%llx q=%llx qq=%llx" $d $q $qq]

     
  • Kevin B KENNY

    Kevin B KENNY - 2008-10-05
    • status: open --> closed-accepted
     
  • Kevin B KENNY

    Kevin B KENNY - 2008-10-05

    Your analysis appears correct. Thanks for the patch!

     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks