From: Paulette L. <Pau...@rs...> - 2008-01-18 00:27:18
|
OS: debian VXL: CVS Nov 2007 compiler: gcc Dear all, It looks like there is a pb with bignum dividion. I have verified the output with magma (computational algebra package). Any other software will do. Let's suppose I want to compute (n choose k). This is the code I am using : vnl_bignum fact(vnl_bignum n) { if (n > 1) return n * fact(n-1); return vnl_bignum(1); } vnl_bignum C(vnl_bignum n, vnl_bignum k) // (n choose k) { vnl_bignum CC(1); if (n == k || k == 0) return CC; /* speed up */ else { for (vnl_bignum i(n); i >= n - k + 1; i--) CC = CC * i; vcl_cout << "\nN " << CC << "\n";); vcl_cout << "fact(k) " << fact(k) << "\n";); CC = CC / fact(k); vcl_cout << "\nCC " << CC << "\n";); return CC; } } So, (n choose k) can be written as N / D where N = CC *= i for i = n-k+1 ... n and D = fact(k) Let n = 100, k = 44. One verifies that computation of N and D is correct: N 131260760632562500136146347534839927958052506820465444563169450468678041600000000000 fact(k) 2658271574788448768043625811014615890319638528000000000 CC 48656022762127140676080092503 but that the result (CC) is not. In magma (or anything suitable), we have: > n := 100; k := 44; > N := 1; > for i := n-k+1 to n do for> N *:= i; for> end for; > N; 131260760632562500136146347534839927958052506820465444563169450468678041600000000000 > Factorial(k); 2658271574788448768043625811014615890319638528000000000 > N / Factorial(k); 49378235797073715747364762200 which is equivalent to saying > Binomial(n, k); 49378235797073715747364762200 In the meantine, if we implement (n choose k) as follows, we get the right answer: vnl_bignum C(vnl_bignum n, vnl_bignum k) // (n choose k) { vnl_bignum CC(1); if (n == k || k == 0) return CC; /* speed up */ else { for (vnl_bignum i(n); i >= n - k + 1; i--) CC *= i; IF_DEB(vcl_cout << "\nN " << CC << "\n";); for (vnl_bignum i(k); i >= 1; i--) CC /= i; IF_DEB(vcl_cout << "\nCC " << CC << "\n";); return CC; } } N 131260760632562500136146347534839927958052506820465444563169450468678041600000000000 CC 49378235797073715747364762200 That is, instead of divide by a bignum, we actually divide by an int. I am sorry, I don't have the time to look into the inplementation of divide, but I thought I should notify this. Thanks! Paulette. -- ============================== Paulette Lieby NICTA Tower A 7 London Circuit Canberra ACT 2601 T: +61 2 6267 6232 F: +61 2 6267 6210 |
From: Peter V. <pet...@ya...> - 2008-01-20 16:02:29
|
Thanks for the bug report, Paulette! > CC 48656022762127140676080092503 > 49378235797073715747364762200 This must indeed be an implementation error in vnl_bignum. I'll have a look at it when I find some time ... -- Peter. -- __________________________________________________________ Går det långsamt? Skaffa dig en snabbare bredbandsuppkoppling. Sök och jämför hos Yahoo! Shopping. http://shopping.yahoo.se/c-100015813-bredband.html?partnerId=96914325 |
From: Peter V. <pet...@ya...> - 2008-04-28 20:56:25
|
Paulette Lieby wrote: > It looks like there is a problem with bignum division. The bug even shows up when trying to divide e.g. 1e200 by 1e100. Tests have been added in core/vnl/tests; see the vxl Dashboard: http://www.cs.rpi.edu/research/groups/vxl/Testing/Sites/lems.brown.edu/Linux-2.6.18_gcc-4.1.2/20080428-0300-Nightly/Results/__core_vnl_tests_vnl_test_bignum.html Strangly enough, some dashboard builds don't have any error, see e.g. http://www.cs.rpi.edu/research/groups/vxl/Testing/Sites/imorphics.com/Linux-2.6.22_gcc-4.1.2_RelWithDebInfo/20080428-0300-Nightly/Results/__core_vnl_tests_vnl_test_bignum.html The only difference I can see are the "-g" and/or "-O" compile switches; seems like adding "-g" to the compilation of vnl_bignum.cxx solves the problem, at least with the Dashboard gcc compilers. I was not able to locate the error spot in the source code; maybe someone more familiar with compiler behaviour, esp. with -g, could have a look? -- Peter. -- __________________________________________________________ Går det långsamt? Skaffa dig en snabbare bredbandsuppkoppling. Sök och jämför hos Yahoo! Shopping. http://shopping.yahoo.se/c-100015813-bredband.html?partnerId=96914325 |