Big float Round off errors



Try:

(%i1) fpprec : 15;
(%o1) 15


(%i2) :lisp(defun $bigfloatepsilon () ($bfloat (div 1 (expt 2 fpprec)))));
$BIGFLOATEPSILON

(%i2) e : bigfloatepsilon();
(%o2) 2.22044604925031b-16

Look at 1.0b0, 1.0b0 + e, and 1.0b0 + 2 * e:

(%i3) ?print([1.0b0, 1.0b0 + e, 1.0b0 + 2 * e]);
((MLIST SIMP) ((BIGFLOAT SIMP 52) 2251799813685248 1)
              ((BIGFLOAT SIMP 52) 2251799813685248 1)
              ((BIGFLOAT SIMP 52) 2251799813685249 1)) 

Notice that 1.0b0 = 1.0b0 + e and 1.0b0 # 1.0b0 + 2 * e.

Barton

maxima-bounces at math.utexas.edu wrote on 05/06/2008 10:37:38 AM:


> Try something like this: http://netlib.org/blas/machar.f