This is the mail archive of the guile@cygnus.com mailing list for the guile project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Jim Blandy <jimb@red-bean.com> writes: > > Here's the kind of thing I was thinking of. This is just off the top > > of my head, so there are probably much better ways of doing these > > things that maybe a numerical analyst on this list could suggest, but > > here's a start: > > Yeesh. The algorithm I posted runs in a constant number of bignum > ops, and lets you specify an exact lower bound on the number of bits > of precision you get. At least it works and is runnable... But, I didn't mean to make a proposal, just to give an example of how these things can be done. BTW, the algorithms I posted also operate in a constant number of bignum ops. It's just a large constant :). You said: > Guile could do better. To compute B / C, where B and C are > bignums, as a floating-point value with a mantissa at most p bits > long, Guile should find the closest integer to (B * 2^p) / C == > (B/C) * 2^p, using only exact bignum arithmetic, and then use the > IEEE 754 scalb function to simultaneously divide by 2^p and convert > to double, without losing bits. The integer you refer to is (quotient (* B (expt 2 mantissa_bits)) C) or this number + 1 (depending on whether the remainder is > or < C/2). But this fails when this number is > the largest double, which is when the result is within 2^p of max_double. This also fails when C is large relative to B. In particular, if C > B*2^p, then the closest integer is zero, even though the quotient is easily representable in a double. Also, I can't find scalb on my system (linux) either in the header files or in the libs. How does it differ from ldexp? I guess you could get around using scalb & get somewhat better accuracy by doing as you suggest, but in 2 pieces - right-shift the bignum p bits to get the integer part, convert it to a double, then add in the bits that were dropped/2^p (as a float). The GNU MP package has code for right-shifting bignums, but however guile implements them this shouldn't be too hard. > What I could really use is a good bignum sqrt algorithm. There's also one in GNU MP. -- Harvey J. Stein BFM Financial Research hjstein@bfr.co.il