On Sun, Oct 4, 2009 at 6:57 AM, rjf <fate...@gmail.com> wrote:
> On Oct 3, 5:11 am, Fredrik Johansson <fredrik.johans...@gmail.com>
> wrote:
>
> My guess is that you have not talked this over with a numerical
> analyst.
No, and I suppose a might if a compelling case were presented to me
that I'm doing something wrong. So far you have failed to present such
a case.
>> The purpose of this code is *not* to add a list of binary
>> fractions accurately.
>
> It is unlikely that the best way to add a list of any numbers that you
> are given is
> to start by throwing out information that provides the detailed values
> of those numbers. Which is what you appear to be doing.
Most of the time the terms are generated within the evalf code, and
are known to be approximations. The evalf code *does* allow you to add
exact binary numbers exactly, if that's what you really want to do --
it will do this by creating an intermediate mantissa that is large
enough to hold the exact sum at any time. This just happens to be
inefficient when adding terms that are several hundred orders of
magnitude apart, which is rare. (If we're talking about differences
less than 100 orders of magnitude or so, the overhead of doing extra
arbitrary-precision operations is much higher than the cost of
increasing the precision. The simple algorithm is faster in the common
case.)
It would be possible to check for the case where all terms are exact
binary fractions of widely different magnitudes and detect exact
cancellations. But exact addition of huge and tiny binary fractions is
a relatively uncommon thing to do and optimizing for it isn't
necessarily worth the added complexity. I've never needed faster code
for such a thing in SymPy, and no one has asked for it.
Generally, the evalf code is written subject to the assumption that
rational parts of symbolic expressions have already been simplified.
SymPy has separate facilities for exact rational arithmetic, and those
are generally better for exact arithmetic.
>>It is to approximate the sum of symbolically
>> defined terms, which with a high probability don't have exact binary
>> representations.
>
> So you are saying that I am lying, or don't know what I'm doing when I
> type in s=1.234501000, (maybe even 1234501.0/1000000.0)
>
> You are saying that you are so sure that I don't mean that number that
> you are going to throw out some of the digits.
>
> And let us say that I didn't type in that number, but I typed in
> r = 5559698243588505001/4503599627370496000
> and then coerced it to a float, call it rf.
> Are you still so sure that I meant some irrational number that you are
> going to throw out some of the digits that I supplied?
>
> would you then compute s-rf by rounding each of them to 5 digits, and
> computing it as zero? Or would you compute the difference between s,
> and r exactly, about 7.07e-17 ?
I'm saying that the code *isn't optimized* for handling extremely
special cases involving exact binary numbers. evalf deals *correctly*
with them. It definitely assumes that user-input floating-point
numbers are what they are (that the user knows what he's doing), and
does not pretend that they are magical approximations of something
else. If you input the Python float (=double) 0.1 it will consider it
to be the exact fraction 3602879701896397/36028797018963968 because
there is no reason to assume anything else.
For example:
>>> Add(Rational(-1,10), x).evalf(subs={x:0.1})
5.55111512312578e-18
>>> Add(Rational(-1,10), x).evalf(30, subs={x:0.1})
5.55111512312578270211815834045e-18
>
>> In this context, including much smaller terms is a
>> waste of time.
>
> I think people would expect Sage to get the arithmetically closest
> representable answer, and would not be so concerned about the
> possibility of wasting a little time.
>
>>
>> I'll illustrate with a simple example (in decimal). Consider three
>> symbolic terms (with high probability irrational or at least not exact
>> binary numbers)
>
> How could you possibly know that the terms you are viewing are not
> exact binary numbers?
> Your best bet is to believe that the numbers you are asked to add
> together are exactly those binary numbers you are asked to add,
> because any other bet is less accurate.
Answered above.
>> which to 10 digit accuracy are:
>>
>> 1.234501000
>> 1.000000000e-10
>> -1.234500000
>>
>> Suppose we want to evaluate their sum to a precision of 5 digits.
>
> If you want to do arithmetic on data that is 10 decimal digits
> precisely in decimal,
That's not what I said. I said that those were 10 digit accurate
representations of some numbers (possibly exact decimals and possibly
not).
> it makes sense to do that in 10 (or more) decimal digit arithmetic. If
> you want the answer to be correct to 10 decimals, you ordinarily have
> to do a little extra work to get it right. If you want to evaluate the
> sum of a new set of data that is the rounded value (to 5 digits) of a
> set of numbers using 5 digit arithmetic, that is a different task.
> The first two numbers add up to
> 1.2345010001 if you did it accurately, but maybe you are not... for 5
> digits, you are adding 1.2345 and 1.0000e-10
> which rounds, in 5 digit arithmetic, to
> 1.2345
> Next, if you add -1.2345 the sum is zero.
>
>
>> The
>> sum should be about 1.0001e-6.
>
> No, you are not adding in 5 digits at all. You apparently are using 10
> digits.
The exact sum of the original numbers, rounded to 5 digits, is
1.0001e-6. This is what evalf computes when called with a precision
argument of 5 digits.
> First we evaluate each term to 5
>> digits:
>>
>> 1.2345
>> 1.0000e-10
>> -1.2345
>>
>> Now, adding these terms in order naively yields 0.0 which is wrong.
>
> Actually, 0.0 is perfectly correct, in 5 digit arithmetic.
Yes, that's exactly what I said. 0.0 is the sum using 5 digit
arithmetic, but the wrong value for the exact sum, and the exact sum
is what evalf is trying to approximate.
>> Adding them using compensated (or exact) summation yields 1.0000e-10,
>> which is the exact sum of the approximated terms but the wrong value
>> for the sum of the exact terms.
>
> Well, it is a far more accurate sum of the three numbers you provided,
> than you
> would expect with a 5-digit adder.
>
>> The only way to compute the sum of the
>> exact terms accurately is to restart with higher precision
>> approximations for the terms.
>
> You can compute the sum perfectly accurately. If you restart with
> different values, you can compute THAT sum perfectly accurately as
> well. You are confusing summation with
> some idea about better approximation of terms.
I am definitely not.
> You are missing a basic fact, and that is you think there are some
> "other" numbers, perhaps irrational (why not transcendental? why not
> complex with itsy bitsy imaginary parts?) that are the "accurate"
> values, and that your summation program is required to somehow read
> the minds of the values to get "higher precision".
That's right, the main purpose of the evalf code is to "read the
minds" of the "other" numbers. It's to be able to convert an exact
expression like exp(1/pi^1000)-1 into a floating-point approximation
of the real (or complex) number it represents. Since "reading the
minds" is generally impossible (determining whether general
expressions are zero is undecidable), it uses a heuristic (tracking
the error and adaptively increasing the precision) that produces the
right result most of the time, and otherwise reports that it failed to
reach the accuracy goal (in the form of a low precision result),
allowing the user to decide what to do (e.g. repeating the calculation
with even higher precision, try to rewrite it in a more numerically
stable form, or simplifying it symbolically...). For example:
>>> (exp(1/pi**1000)-1).evalf()
.0e-125
>>> (exp(1/pi**1000)-1).evalf(maxprec=1000)
7.08153336783938e-498
>>> (2*exp(1/pi**1000/2)*sinh(1/pi**1000/2)).evalf()
7.08153336783938e-498
The best way to handle sub-evaluations that may have exact results is
to try to do them symbolically / with exact rationals before calling
evalf. Ultimately the user needs to be aware of these issues and
choose the appropriate sequence of operations.
Fredrik
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel-unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
