R Fateman wrote:
> If a,b,c,d are all bfloats, or even if only one of a,c and one of
> b,d are bfloats, then [a+c,b+d] is all you have to do. If they are
> not, why would you round the endpoints of an interval, potentially
> losing some of it?
>
> Nevertheless.. If a, b,c,d are all double-floats or rationals or
> integers or ....
>
> then all you have to do is [bfloat(a+c), bfloat(b+d)]
>
> But bfloats are not all the same, since they depend on fpprec.
>
> Also what you are doing is actually "wrong" because of rounding issues,
> if you want rigorous intervals.
> Are you aware of the work by Chee Yap (at NYU) or other work on
> geometric calculations, e.g.
> at Berkeley, Jonathan Shewchuk?
>
>
> And yes, iterated rational computations often make numerator and
> denominator grow in size.
>
>
> Finding a "small" approximation to a rational number is often possible,
> though.
>
> Robert Dodier wrote:
>
>> On 1/26/08, S. Newhouse <sen1 at math.msu.edu> wrote:
>>
>>
>>
>>> For instance, if one defines the sum of intervals [a,b], [c,d] (with 0
>>> < a < b and 0 < c < d) as
>>> Iplus([a,b],[c,d])=[a+c, b+d],
>>> then I would define
>>> bIplus([a,b],[c,d]) = [bfloat(a)+bfloat(c), bfloat(b) + bfloat(d)]
>>>
>>>
>> (On a tangent: this shows a "destructuring bind" style of function definition
>> that could be generally useful. Maxima can do something similar with
>> pattern matching but expressed in Maxima's pattern matching language,
>> it is quite a bit more obscure.)
>>
>>
>>
>>> It seems that I have to put in 'bfloat' in every representation of a
>>> real number.
>>>
>>> Is there a way to declare all variables to be evaluated using 'bfloat'
>>> without simply applying 'bfloat' to every occurrence?
>>>
>>>
>> Not that I know of. Another way to approach this is to define an
>> operator, let's say "bf+", which applies bigfloat conversion to its
>> arguments, e.g.
>>
>> nary ("bf+");
>> "bf+" ([L]) := apply ("+", bfloat (L));
>>
>> 1 bf+ %pi bf+ %e;
>> => 6.859874482048838b0
>>
>> bIplus (I1, I2) := [I1[1] bf+ I2[1], I1[2] bf+ I2[2]];
>> bIplus ([%e, 10], [1, %pi]);
>> => [3.718281828459045b0, 1.314159265358979b1]
>>
>> Writing "bf+" has the advantage that one can tell from reading
>> an isolated expression (e.g. I1[2] bf+ I2[2]) that bigfloat conversion
>> is applied to the arguments. If the conversion were controlled by
>> some global variable, you would have to read the whole program
>> and keep the global state in mind when reading each expression.
>> (This is a general argument against global variables btw.)
>>
>> I wonder if you really need to enforce bigfloat arithmetic.
>> Users might expect, extrapolating from Maxima's general treatment
>> of integers and rationals, that an operation on exact intervals
>> yields another exact interval.
>>
>> HTH
>>
>> Robert
>> _______________________________________________
>> Maxima mailing list
>> Maxima at math.utexas.edu
>> http://www.math.utexas.edu/mailman/listinfo/maxima
>>
>>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
> .
>
>
Thanks for the comments. I was not aware of the work of Chee Yap and
Shewchuk. I'll look at their papers.
But let me way that the example I used before was only to simplify the
code to make the question about declaring variables easy to understand.
My current defintions of Iplus and bIplus are actually the following.
Iplus(x,y):=[ x[1]+y[1]-r_err*(abs(x[1])+abs(y[1])),
x[2]+y[2]+r_err*(abs(x[1])+abs(y[1]))];
bIplus(x,y):= [ bf(x[1])+bf(y[1])-bf(r_err)*( bf(abs(x[1]))+
bf(abs(y[1]))), bf(x[2])+bf(y[2])+bf(r_err)*(bf(abs(x[1]))+bf(abs(y[1])))];
Of course, I have them for Iminus, Iprod, Idiv, etc and evaluations of
polynomials as well.
Here r_err is set to give a particular round-off error:
e.g. fpprec: 16, r_err: 1e-16 or
fpprec: 50, r_err: 1e-40 (less to avoid the dropping off of the last few
bits and get a generous enclosure)
Do these still look wrong to you? I am hoping that I get rigorous
enclosures of the actual floating point (and big floating point)
calculations, but I haven't checked closely yet. I am using that, if op
represents a standard arithmetic operation and fop represents its
floating point counterpart with precistion eps, then
x fop y = (1 + eps)*(x op y)
as long as x and y are not too big or too small (I was going to look at
this more closely when I get to actually claim something rigorous, but
maybe I should do this now). I got this with a quick reading of Kahan's
IEEE standard paper. Again, am I missing something here?
I did not try yet to go through all the arguments needed to try to prove
the statements in Kahan's paper, and I hope I did not misunderstand his
paper.
For rational data a, b, since bfloat(a+b) only does one conversion, it
seems that bfloat(a+b) is closer to the actual sum a+b than bfloat(a) +
bfloat(b), so I'll change the definitions accordingly. However, by
taking a few extra bits off as in the (fpprec: 50, r_err: 1e-40), isn't
the difference between (bfloat(a+b) and bfloat(a) + bfloat(b)) not
important?
Incidentally, I tried Robert's suggestions about 'bf+', etc. and they
work fine, but in doing many operations, the code is about 10% slower
than putting the bfloats in directly. I guess this is because the bfloat
conversions are done with every operation in maxima whereas using the
bfloats for addition directly does many operations in the underlying lisp.
I find these kinds of discussions, very useful, but if they are more
appropriately done by private email--off the maxima list-- I am happy to
do that.
Comments?
TIA,
-sen