Accuracy and error analysis (was Re: [Maxima] primes)

>>>>> "Doug" == Doug Stewart  writes:

    Doug> here is one way to do the Numerical error estimation of measurement
    Doug> error.

    Doug> This is Octave code

    Doug> r1=75.1
    Doug> r2=99.3
    Doug> r3=300

    Doug> y=[ 0 0  0
    Doug>       0 0 1
    Doug> 0 1 0
    Doug> 0 1 1
    Doug> 1 0 0
    Doug> 1 0 1
    Doug> 1 1 0
    Doug> 1 1 1
    Doug> ];

    Doug> % the .05 is the   +-5%  error  in the measurement
    Doug> % If your measurement is accurate to 1% then use 0.01
    Doug> for k=1:8
    Doug> r11=r1+r1*.05*(-1)^y(k,1);
    Doug> r22=r2+r2*.05*(-1)^y(k,2);
    Doug> r33=r3+r3*.05*(-1)^y(k,3);
    Doug> %r(k)=1/(1/(r11)+1/(r22)+1/(r33));
    Doug> r(k)=r11*r22*r33;
    Doug> end
    Doug> range=max(r)-min(r)
    Doug> rt=r1*r2*r3
    Doug> range/2
    Doug> rt*0.15


    Doug> This code adds and subtracts the % error to each input and calculates 
    Doug> all possible combinations and then takes the largest and smallest 
    Doug> answers and gets a range.

    Doug> This example is trivial because we know the total will be the sum of the

    Doug> 3 input ranges  5%+5%+5%=15%.
    Doug> But for some complex formulas this would is the way that I get a handle 
    Doug> on the range of significant digits.

FWIW, With cmucl (and sbcl):

        (defun zot (r1 r2 r3)
	   (declare (type (double-float 71.345d0 78.75d0) r1)
		    (type (double-float 94.335d0 104.625d0) r2)
		    (type (double-float 285d0 315d0) r3))
	   (* r1 r2 r3))
       (describe 'zot)
       ZOT is an internal symbol in the COMMON-LISP-USER package.
       Function: #
       Function arguments:
         (r1 r2 r3)
       Its defined argument types are:
         ((DOUBLE-FLOAT 71.345d0 78.75d0) (DOUBLE-FLOAT 94.335d0 104.625d0)
          (DOUBLE-FLOAT 285.0d0 315.0d0))
       Its result type is:
         (DOUBLE-FLOAT 1918144.2138749997d0 2595353.90625d0)

So that if r1, r2, and r3 are +/- 5% of the nominal values (as given
in the declarations for r1, r2, r3), CMUCL thinks that the final
product is in the range [1918144.2138749997d0 2595353.90625d0].

Some Lisp compilers can do interval arithmetic in the compiler.

Ray