[Maxima-bugs] [ maxima-Bugs-2910001 ] optimization of bfloat input

[SourceForge.net <no...@so...>]

Bugs item #2910001, was opened at 2009-12-07 08:16
Message generated for change (Comment added) made by rtoy
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2910001&group_id=4933

Please note that this message will contain a full copy of the comment thread,
including the initial issue submission, for this request,
not just the latest update.
Category: Lisp Core - Floating point
Group: None
>Status: Closed
>Resolution: Fixed
Priority: 5
Private: No
Submitted By: Clarphimous (clarphimous)
Assigned to: Nobody/Anonymous (nobody)
Summary: optimization of bfloat input

Initial Comment:
When using the input method that goes like 2.53b10, it tends to get very slow with large exponents.  This applies for large negative exponents as well.  I believe the reason is that it is converting all the digits up to the decimal point into binary, where it really only needs to convert them up to the number of bits of precision.  For example:

1b1000000;
takes quite a while to perform.  Not far above this, Maxima crashes on my system (another bug I was originally investigating, which applies to rational numbers as well).

On the other hand, if we do a calculation like
x:1b1000;
x^40000000;
you can get a result in less than a second, although it does have some rounding errors at the end.

----------------------------------------------------------------------

>Comment By: Raymond Toy (rtoy)
Date: 2009-12-14 09:39

Message:
Some fixes have been checked in.  The fast algorithm is used if the
exponent exceeds 100000, but that is user-changeable.  Further tests
indicate that the fast algorithm matches the exact algorithm most of the
time (.2% difference) if the number of digits is less than fpprec.  If the
number of digits exceeds fpprec, the fast method differs about 5% of the
time.

Perhaps this is good enough, since it's easy to change the bfloat
precision.

Closing report.


----------------------------------------------------------------------

Comment By: Raymond Toy (rtoy)
Date: 2009-12-10 19:30

Message:
Extra digits are already used, based on the size of the exponent.  (Hmm. 
Should include  the exponent and the number in case someone does 12345<many
digits>90.123<many more>90 b <bigexp>)

And with your example 1.12233445566...b0, maxima computes
11223344556677889900112233/10^25, exactly, and then rounds that (once) to a
bfloat with 10 digits.

The fast method always rounds at least 3 times, perhaps many more when
computing 10b0^n.

All of the differences occur with negative exponents.  Don't know why that
is.

----------------------------------------------------------------------

Comment By: Clarphimous (clarphimous)
Date: 2009-12-10 17:12

Message:
Alrighty.  I think what would prevent roundoff errors is if you give it a
few extra bits of precision while doing the calculations.  In other words,
guard digits/bits.  If the input is shorter than the number of digits of
precision, add some extra zeroes.  The calculators I've used usually have 3
guard digits (although they were using binary coded decimal).  For
example:

fpprec: 10;
x:1.1223344556677889900112233b0

For this, cut it off 3 digits after the 10th digit.  That would be
1.122334455668, with rounding.  Then, convert that to binary with some
extra bits of precision.  Finally, round that answer off to the number of
bits you usually have.

Another example:

a:1.23456789b0;

For this I would assume that the actual input is 1.234567890000.  Then
convert that to binary, and cut off the excess bits when you're done.

It might take some experimentation to figure out how how many extra digits
and bits you need, but I'm pretty sure this will work.  Here is a quick
experiment I did where I manually rounded the numbers:

fpprec: 10;
x:1.1223344556677889900112233b0;
x-1.122334455b0;
y:1.1223344557b0;
y-1.122334455b0;
z:1.12233445567b0;
z-1.122334455b0;

http://www.disillusionary.com/files/exp.png

----------------------------------------------------------------------

Comment By: Raymond Toy (rtoy)
Date: 2009-12-10 12:40

Message:
I've done some work on this, and bfloat input for large inputs is now much
faster. 2.53b10 is converted by computing bfloat(253)*10b0^(10-2).   Except
the exponent is much larger.  Small exponents use the current algorithm.

The new algorithm has roundoff issues and about 5% of the time differs
from the old algorithm.  Perhaps this is acceptable for large exponents.

----------------------------------------------------------------------

Comment By: Clarphimous (clarphimous)
Date: 2009-12-08 02:22

Message:
Okay, I don't know exactly how you'd implement it, as I'm not very
knowledgeable about base conversion algorithms (although I do know a
little).  However, I can tell you that the Sage CAS performs the exact same
conversion without taking nearly as long for really large numbers.  (On the
other hand, the default Sage interface seems to only accept floating-point
numbers up to around 323,000,000 digits before marking them as
"infinity.")

First, a 1000 digit bfloat in Maxima:

set_display('ascii)$
fpprec:35;
x:3.4000000000000000000000000000000001b1000;
fpprec:1001;
(x+1)-1;

result: http://www.disillusionary.com/files/maxima_big_a.png


Next, the same calculation in Sage.

x = 3.4000000000000000000000000000000001e1000
x
parent(x)
y = ZZ(x)
parent(y)
y

result: http://www.disillusionary.com/files/sage_big_a.png

Apparently with this calculation the same number of bits of precision are
used.  For the next one, though, I had to pad Sage's input with an extra
zero to get the same result.

In Maxima:

fpprec:35;
x:3.4000000000000000000000000000000001b5000000;
fpprec:1001;
x+1b4999000;

result: http://www.disillusionary.com/files/maxima_big_b.png


And now for Sage

x = 3.40000000000000000000000000000000010e5000000
x
parent(x)
y = ZZ(x)
parent(y)
y

result: http://www.disillusionary.com/files/sage_big_b.png

Here is where the results in speed and memory usage come into play.  Sage
doesn't take any longer than it usually does to store the value of x.

So... the point is that there is a way, but I do not know how to do it,
myself.  I was thinking maybe you could use the same process you use to
output the decimal value of the bfloat each time it's displayed... but I
don't know if it works both ways.

----------------------------------------------------------------------

Comment By: Raymond Toy (rtoy)
Date: 2009-12-07 14:48

Message:
The algorithm for reading bfloats converts the bfloat to an exact rational.
 Thus 2.53b10 is 253/100*10^10.  This is then converted to a bfloat.  So
1b1000000 is converted to the exact rational 10^1000000 and converted to a
bfloat.  This should incur exactly one rounding operation.

I think doing it in any other way will incur additional roundoff errors,
and we would very much like that the internal representation be as close as
possible to the input number.

This is also complicated by the fact that bfloats are internally
represented using base 2 exponents, not base 10.

Marking as pending/wontfix.

----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=104933&aid=2910001&group_id=4933