file mbe9code.txt  has Maxima input code
   for Maxima by Example, Chapter 9,
   Bigfloats and High Accuracy Quadrature.

    Edwin L Woollett, 2011
    woollett@charter.net
    http://www.csulb.edu/~woollett
-------------------------------------------------

1 Introduction
------------------------------------------
2 Bigfloat Numbers in Maxima

(%i1) load(bfloat);
(%i2) [fpprec,fpprintprec,display2d];
(%i3) bfloat(1);
(%i4) bfdigits(%);
(%i5) bfloat(10^5);
(%i6) bfloat(10^-5);
(%i7) bfdigits(%);
(%i8) bftrunc:false$
(%i9) bfloat(10^-5);
(%i10) bfdigits(%);
(%i11) bftrunc:true$
(%i12) fpprec;
(%i13) q : 1.0b-5;
(%i14) bfdigits(q);
(%i15) bftrunc:false$
(%i16) q;
(%i17) bftrunc:true$

(%i1) load(bfloat);
(%i2) [fpprec,fpprintprec];
(%i3) q1 : bfloat(1/10);
(%i4) q2 : 1b-1;
(%i5) is(equal(q2,q1));
(%i6) q3 : 10b0^(-1);
(%i7) is(equal(q3,q1));
(%i8) 10b0;

(%i1) fpprec;
(%i2) ?fpprec;
(%i3) :lisp $fpprec
(%i3) :lisp fpprec
---------------------------------------------------------
2.1 Experimenting with Bigfloats

(%i1) sin(1);
(%i2) sin(1.0);
(%i3) sin(1),numer;
(%i4) float(sin(1));

(%i5) fpprec;
(%i6) q : 1.0b0;
(%i7) bfdigits(q);
(%i8) q1 : sin(q);
(%i9) q2 : bfloat(sin(1));
(%i10) map ('bfdigits,[q1,q2]);
(%i11) is (equal (q1,q2));
(%i12) q3 : bfloat(sin(1)), fpprec:20;
(%i13) q4 : sin(1.0b0), fpprec : 20;
(%i14) map ('bfdigits,[q3,q4]);
(%i15) is (equal (q3,q4));
(%i16) fpprec;

(%i17) s20 : bfloat (sin(1)), fpprec:20;
(%i18) sfp : sin(1.0);
(%i19) sfp_20 : bfloat(sfp),fpprec:20;
(%i20) abs(sfp_20 - s20);
(%i21) abs(sfp_20 - s20), fpprec:20;
(%i22) abs(sfp - s20), fpprec:20;

(%i23) fpprec;
(%i24) fpprec:20$
(%i25) sfp : sin(1.0);
(%i26) sfp_20 : bfloat(sfp);
(%i27) s20 : bfloat (sin(1));
(%i28) abs(sfp_20 - s20);
(%i29) abs(sfp - s20);

(%i30) fpprec;
(%i31) erf(1.0);
(%i32) erf(1.0b0);
(%i33) bessel_j(0,1.0);
(%i34) bessel_j(0,1.0b0);

(%i1) b1 : 1.57b0;
(%i2) load(bfloat);
(%i3) bfdigits(b1);
(%i4) bfloatp(b1);
(%i5) bftrunc;
(%i6) bftrunc : false$
(%i7) b1;
(%i8) fpprec;
(%i9) fpprintprec;
(%i10) floatnump(b1);
(%i11) bfdigits(1.000000000000000b0);
(%i12) float(exp(%i*%pi/sin(1)));
(%i13) float(exp(%i*%pi/sin(1))),numer;
(%i14) float( rectform (exp(%i*%pi/sin(1))));

(%i15) fpprec;
(%i16) q : bfloat(exp(%i*%pi/sin(1)));
(%i17) bfloatp(q);
(%i18) realpart(q);
(%i19) imagpart(q);
(%i20) imagpart(1.2b0);

(%i1) load(bfloat);
(%i2) q : bfloat(exp(%i*%pi/sin(1)));
(%i3) bfdigits(q);
(%i4) q : bfloat(exp(%i*%pi/sin(1))), fpprec:20;
(%i5) bfdigits(q);

(%i6) [fpprec,fpprintprec,bftrunc];
(%i7) q : 1.0b0;
(%i8) bfdigits(q);
(%i9) bftrunc:false$
(%i10) q;
(%i11) bftrunc:true$
(%i12) q : bfloat(1),fpprec:20;
(%i13) bfdigits(q);

(%i14) fpprec;
(%i15) q1 : -q;
(%i16) bfdigits(q1);
(%i17) q2 : abs(q);
(%i18) bfdigits(q2);

(%i19) q5 : -q, fpprec:20;
(%i20) bfdigits(q5);
(%i21) q6 : abs(q), fpprec:20;
(%i22) bfdigits(q6);

(%i23) q3 : sin(q);
(%i24) bfdigits(q3);
(%i25) q4 : sin(q), fpprec:20;
(%i26) bfdigits(q4);

--------------------------------------------
2.2 Bigfloat Numbers and fpprintprec

(%i1) load(bfloat);
(%i2) [fpprec,fpprintprec];
(%i3) q : integrate(exp(x),x,-1,1);
(%i4) q : bfloat(q),fpprec : 45;
(%i5) bfdigits(q);
(%i6) print (q),fpprintprec : 12$
(%i7) fpprec;
(%i8) print (q),fpprintprec : 15$
(%i9) print (q),fpprintprec : 16$
(%i10) print (q),fpprintprec : 17$
(%i11) print (q),fpprec:18, fpprintprec : 17$
(%i12) print (q),fpprec:19, fpprintprec : 18$
(%i13) fpprintprec;
(%i14) print(q)$

(%i15) [fpprec,fpprintprec];
(%i16) piby2 : block([fpprintprec,fpprec:30,val],
                val:bfloat(%pi/2),
                fpprintprec:8,
                disp (val),
                display (val),
                print("  val =   ",val),
                print("  bfdigits(val) = ", bfdigits(val) ),
                val);
(%i17) bfdigits(piby2);
(%i18) [fpprec,fpprintprec];

(%i19) piby2 : (val:bfloat(%pi/2),
                disp (val),
                display (val),
                print("  val =   ",val),
                print("  bfdigits(val) = ", bfdigits(val) ),val ), fpprintprec:8, fpprec:30;
(%i20) piby2;
(%i21) bfdigits(piby2);
(%i22) val;

(%i23) f2(w) := block([v2 ],
                  print ("   in function f2 "),
                  display([w,fpprec,fpprintprec]),
                  print ("  bfdigits(w) = ", bfdigits(w)),
                  v2 : sin(w),
                  print (" v2 = sin(w) =  ", v2),
                  print ("     bfdigits(v2) = ", bfdigits(v2) ),
                  print("  "),
                  v2 )$
(%i24) f1(x, fpp, fpprt) := block([fpprintprec,fpprec:fpp,v1],
                  fpprintprec:fpprt,             
                  display([x,fpprec,fpprintprec]),
                  print("  in function f1, call f2(x) "),
                  v1 : f2(bfloat(x))^2,			  
                  print(" back  in f1, v1 = sin(x)^2 =  ",v1),
                  print(" bfdigits(v1) = ",bfdigits(v1)),
                  v1 )$
				  
(%i25) q : f1(0.5,30,8);
(%i26) q;
(%i27) bfdigits(q);
(%i28) [fpprec,fpprintprec];

(%i29) fpprec:20$
(%i30) q : bfloat(2);
(%i31) bfdigits(q);
(%i32) q1 : 1 + q;
(%i33) bfdigits(q1);
(%i34) q2 : 1.0 + q;
(%i35) bfdigits(q2);
(%i36) q3 : sin(q2);
(%i37) bfdigits(q3);

(%i38) a : 2.38b-1$
(%i39) a;
(%i40) bfdigits(a);
(%i41) q4 : exp(%pi*a);
(%i42) q4 : bfloat(q4);
(%i43) bfdigits(q4);

(%i44) fpprec;
(%i45) q1 : bfloat(1);
(%i46) bfdigits(q1);
(%i47) q2 : bfloat(10^(-25));
(%i48) bfdigits(q2);
(%i49) q3 : q1 + q2;
(%i50) q3, fpprec:25;
(%i51) q3, fpprec:26;
(%i52) q3 : q1 + q2,fpprec:26;
(%i53) bfdigits(q3);
----------------------------------------------
2.3 Using print with Bigfloats

(%i1) load(bfloat);
(%i2) [fpprec,fpprintprec];
(%i3) integrate(exp(x),x,-1,1);
(%i4) q : bfloat(%),fpprec : 45;
(%i5) bfdigits(q);

(%i6) print(q),fpprintprec:5$
(%i7) print(q),fpprintprec:15$
(%i8) print(q),fpprintprec:16$
(%i9) print(q),fpprintprec:0$
(%i10) print(q),fpprec:21,fpprintprec:20$
(%i11) print(q),fpprec:17,fpprintprec:16$
(%i12) [fpprec,fpprintprec];

(%i13) bfprint(" q16 = ",q,16)$
(%i14) bfprint(" q20 = ",q,20)$
------------------------------------------------------
2.4 Using printf with Bigfloats

(%i15) q : bfloat(exp(-20)),fpprec : 30;
(%i16) bfdigits(q);
(%i17) printf(true,"~d~a",3,string(q))$
(%i18) printf(true,"  ~d    ~a",3,string(q))$
(%i19) (printf(true,"  ~d    ~a~%",3,string(q)),
                  printf(true,"  ~d    ~a",3,string(q)))$

(%i20) printf(true,"  ~d    ~a",3,string(q)), fpprintprec:8$
(%i21) printf(true,"  ~d    ~f",3,q)$
(%i22) printf(true,"  ~d    ~e",3,q)$
(%i23) printf(true,"  ~d    ~h",3,q)$
(%i24) printf(true,"  ~d    ~h",3,q),fpprintprec : 12$
-----------------------------------------------------------------------
2.5 A Table of Bigfloats using printf

(%i25) print_test(fp) :=
        block([fpprec,fpprintprec,val],
         fpprec : fp,
         fpprintprec : 8,
         display(fpprec),
         print("  k          value "),
         print("  "),
         for k thru 4 do
          ( val : bfloat(exp(k^2)), 
            printf(true,"  ~d       ~a  ~%",k,string(val) ) ))$
(%i26) print_test(30)$

(%i27) print_test2(fp) :=
        block([fpprec,fpprintprec,val],
         fpprec : fp,
         fpprintprec : 8,
         display(fpprec),
         printf(true,"~%  ~a          ~a ~%~%",k,value), 
         for k thru 4 do
          ( val : bfloat(exp(k^2)), 
            printf(true,"  ~d       ~a  ~%",k,string(val) ) ))$
(%i28) print_test2(30)$
----------------------------------------------------------------
2.6  Adding Bigfloats Having Differing Accuracy

(%i1) fpprintprec:8$
(%i2) fpprec;
(%i3) pi100 : bfloat(%pi),fpprec:100;
(%i4) pi50 : bfloat(%pi),fpprec:50;
(%i5) abs(pi50 - pi100),fpprec:101;
(%i6) twopi : bfloat(2*%pi),fpprec:100;
(%i7) pisum : pi50 + pi100,fpprec:100;
(%i8) pisum - twopi,fpprec:101;
-------------------------------------------------------------------

2.7  Polynomial Roots Using bfallroots

(%i1) load(bfloat);
(%i2) fpprec;
(%i3) pi50 : bfloat(%pi),fpprec:50;
(%i4) bfdigits(pi50);
(%i5) e : expand ( (x-%pi)^3);
(%i6) e16 : float(e);
(%i7) sar16 : map ('rhs, allroots (%i*e16));
(%i8) for s in sar16 do disp( expand (subst (s, x, e16)))$
(%i9) for s in sar16 do disp( expand (subst (s, x, %i*e16)))$
(%i10) for s in sar16 do disp (pi50 - realpart(s))$
(%i11) sbfar16 : map ('rhs, bfallroots (%i*e16));
(%i12) for s in sbfar16 do  disp( expand (subst (s,x,e16)))$
(%i13) for s in sbfar16 do disp (pi50 - realpart(s))$
(%i14) e;
(%i15) fpprec:40$
(%i16) e40 : bfloat(e);
(%i17) bftorat;
(%i18) ratprint:true$
(%i19) ratepsilon;
(%i20) sbfar40 : map ('rhs, bfallroots (%i*e40));
(%i21) for s in sbfar40 do disp (expand  (subst (s,x,e40)))$
(%i22) for s in sbfar40 do disp (pi50 - realpart(s))$
(%i23) 9.424777960769379715387930149838508652592b0 -
                    9.424777960769379480849831552612202452135b0;
(%i24) bfloat (9.424777960769379715387930149838508652592b0 -
                   9.424777960769379480849831552612202452135b0);
(%i25) ratepsilon : 1.0e-41$
(%i26) sbfar40 : map ('rhs, bfallroots (%i*e40));
(%i27) for s in sbfar40 do disp (expand  (subst (s,x,e40)))$
(%i28) for s in sbfar40 do disp (pi50 - realpart(s))$
(%i29) fpprec;
(%i30) ratepsilon:2.0e-15;
(%i31) bftorat:true$
(%i32) sbfar40 : map ('rhs, bfallroots (%i*e40));
(%i33) for s in sbfar40 do disp (expand  (subst (s,x,e40)))$
(%i34) for s in sbfar40 do disp (pi50 - realpart(s))$
----------------------------------------------------------------------------
2.8 Finding Highly Accurate Roots using bf_find_root

(%i1) p4 : legendre_p(4,x);
(%i2) solve(p4,x);
(%i3) float(%);
(%i4) fpprintprec:8$
(%i5) [x1,x2,x3,x4] : map ('rhs,%o3);
(%i6) for z in % do print (z, ev(p4,x = z))$
(%i7) fpprec : 40$
(%i8) rt : bf_find_root(p4,x,0.3,0.4);
(%i9) p4,x = rt;
(%i10) fpprintprec : 0$
(%i11) rt;
(%i10) plot2d([1-2*x/3 - sin(x),[discrete,[[0,0],[1,0]]]],[x,0,1],
                  [style,[lines,2]],[legend,false])$
(%i1) g(x) := 1 - 2*x/3 - sin(x)$
(%i2) g(0.62);
(%i3) fpprec : 40$
(%i4) rt : bf_find_root(g,0.6,0.64);
(%i5) g(rt);
------------------------------------------------------------------------------------
2.9   Bigfloat Number Gaps and Binary Arithmetic

(%i1) fpprec:2$
(%i2) ?fpprec;
(%i3) :lisp $fpprec
(%i3) :lisp fpprec
(%i3) x : 0.0b0$
(%i4) :lisp $x;
(%i4) x : 1.0b0;
(%i5) ?print(x);
(%i6) :lisp $x
(%i6) :lisp (msetq $x '((BIGFLOAT SIMP 9) 0 0))
(%i6) x;				  
(%i7) :lisp (msetq $u '((BIGFLOAT SIMP 9) 1 0))
(%i7) u;

(%i1) fpprec;
(%i2) fpprec:4$
(%i3) ?fpprec;
(%i4) x :bfloat(2/3);
(%i5) u : bfloat(2^(-18));
(%i6) x1 : x + u;
(%i7) x1 - x;
(%i8) x2 : x + 3.814b-6;
(%i9) x2 - x;
(%i10) ulp : bfloat(2^(-17));

(%i11) fpprec:16$
(%i12) ?fpprec;
(%i13) x :bfloat(2/3);
(%i14) u : bfloat(2^(-58));
(%i15) x1 : x + u;
(%i16) x1 - x;
(%i17) x2 : x + 3.469446951953613b-18;
(%i18) x2 - x;
(%i19) ulp : bfloat(2^(-57));
--------------------------------------------------------------
2.10   Effect of Floating Point Precision on Function Evaluation

(%i1) load(bfloat);
(%i2) fpprintprec:8$
(%i3) g(x) := sin(x/2)$
(%i4) fpprec;
(%i5) fdf(g,1,10);
(%i6) fdf(g,1,10),fpprec:30;
(%i7) fpprec;

====================================
3.   High Accuracy Quadrature with Maxima

3.1  Using bromberg for High Accuracy Quadrature

(%i1) load(bfloat);
(%i2) fpprintprec:8$
(%i3) load(brmbrg);
(%i4) [brombergtol,brombergabs,brombergit,brombergmin,fpprec,fpprintprec];
(%i5) integrate(exp(x),x,-1,1);
(%i6) tval : bfloat(%),fpprec:42;
(%i7) fpprec;
(%i8) (brombergtol:0.0b0,brombergit:100)$
(%i9) b15:(brombergabs:1.0b-15, bromberg (exp(x),x,-1,1) ), fpprec:30;
(%i10) abs (b15 - tval),fpprec:42;
(%i11) b20:(brombergabs:1.0b-20, bromberg (exp(x),x,-1,1) ), fpprec:30;
(%i12) abs (b20 - tval),fpprec:42;

(%i1) load(bfloat);
(%i2) fpprintprec:8$
(%i3) load(brmbrg);
(%i4) [brombergtol,brombergabs,brombergit,brombergmin,fpprec,fpprintprec];
(%i5) qbr20 : qbromberg(exp,-1,1,20,40,100);
(%i6) abs(qbr20 - tval);
(%i7) abs(qbr20 - tval),fpprec:40;
(%i8) qbrlist(exp,-1,1,[10,15,17 ],20,100)$
(%i9) qbrlist(exp,-1,1,[10,20,27 ],30,100)$
(%i10) qbrlist(exp,-1,1,[10,20,30,35],40,100)$

(%i1) load(bfloat);
(%i2) g(x):= sqrt(x)*log(x)$
(%i3) i1 : integrate(g(t),t,0,1);
(%i4) fpprec : 20$
(%i5) tval : bfloat(i1);
(%i6) i18 : bromberg_abs(g,0,1,1.0b-18);
------------------------------------------------------------------
3.2  A Double Exponential Quadrature Method for a <= x < inf

(%i1) fpprec;
(%i2) fpprintprec:8$
(%i3) g(x):= exp(-x)$
(%i4) ival : integrate(g(x),x,0,inf);
(%i5) tval : bfloat(ival),fpprec:45;
(%i6) load(quad_de);
(%i7) quad_de(g,0,30,40);
(%i8) abs(first(%) - tval),fpprec:45;
(%i9) idek(g,0,4,40);
(%i10) abs(% - tval),fpprec:45;
(%i11) idek_e(g,0,4,40);
(%i12) abs(first(%) - tval),fpprec:45;
(%i13) ide(g,0,30,40)$
(%i14) ide_test(g,0,30,40)$

Test Integral 1

(%i15) g(x):= exp(-x)/sqrt(x)$
(%i16) integrate(g(t),t,0,inf);
(%i17) tval : bfloat(%),fpprec:45;
(%i18) quad_de(g,0,30,40);
(%i19) abs(first(%) - tval),fpprec:45;
(%i20) idek_e(g,0,4,40);

Test Integral 2

(%i21) g(x) := exp(-x^2/2)$
(%i22) tval : bfloat(sqrt(%pi/2)),fpprec:45$
(%i23) quad_de(g,0,30,40);
(%i24) abs(first(%) - tval),fpprec:45;
(%i25) idek_e(g,0,5,40);

Test Integral 3

(%i26) g(x) := exp(-x)*cos(x)$
(%i27) integrate(g(x),x,0,inf);
(%i28) tval : bfloat(%),fpprec:45$
(%i29) quad_de(g,0,30,40);
(%i30) abs(first(%) - tval),fpprec:45;
(%i31) idek_e(g,0,5,40);

-------------------------------------------------------------
3.3   The tanh-sinh Quadrature Method for a <= x  <= b

(%i1) fpprec;
(%i2) fpprintprec:8$
(%i3) tval : bfloat( integrate( exp(x),x,-1,1 )),fpprec:45;
(%i4) load(quad_ts);
(%i5) bfprint(tval,45)$
(%i6) quad_ts(exp,-1,1,30,40);
(%i7) abs(first(%) - tval),fpprec:45;
(%i8) qtsk(exp,-1,1,5,40);
(%i9) abs(% - tval),fpprec:45;
(%i10) qtsk_e(exp,-1,1,5,40);
(%i11) abs(first(%) - tval),fpprec:45;
(%i12) last(_yw%[8,40]);
(%i13) fpxy(40)$
(%i14) qts(exp,-1,1,30,40)$
(%i15) qts_test(exp,-1,1,30,40)$   

   Test Integral 1
   
(%i16) g(x):= sqrt(x)*log(x)$
(%i17) tval : bfloat(integrate(g(t),t,0,1)),fpprec:45;
(%i18) quad_ts(g,0,1,30,40);
(%i19) abs(first(%) - tval),fpprec:45;
(%i20) qtsk_e(g,0,1,5,40);   

  Test Integral 2
  
(%i21) g(x):= atan(sqrt(2+x^2))/(sqrt(2+x^2)*(1+x^2))$
(%i22) integrate(g(t),t,0,1);
(%i23) quad_qags(g(t),t,0,1);
(%i24) float(5*%pi^2/96);
(%i25) tval: bfloat(5*%pi^2/96),fpprec:45;
(%i26) quad_ts(g,0,1,30,40);
(%i27) abs(first(%) - tval),fpprec:45;
(%i28) qtsk_e(g,0,1,5,40);

  Test Integral 3

(%i29) g(x):= sqrt(x)/sqrt(1 - x^2)$
(%i30) quad_qags(g(t),t,0,1);
(%i31) integrate(g(t),t,0,1);
(%i32) tval : bfloat(%),fpprec:45;
(%i33) quad_ts(g,0,1,30,40);
(%i34) abs(first(%) - tval),fpprec:45;
(%i35) qtsk_e(g,0,1,5,40);
(%i36) integrate(g(t),t,0,1);
(%i37) makegamma(%);
(%i38) bfloat(%),fpprec:45;

  Test Integral 4
  
(%i39) g(x) := log(x)^2$
(%i40) integrate(g(t),t,0,1);
(%i41) tval : bfloat(%), fpprec:45;
(%i42) quad_ts(g,0,1,30,40);
(%i43) abs( first(%) - tval ),fpprec:45;
(%i44) qtsk_e(g,0,1,5,40);

   Test Integral 5

(%i45) g(x) := log( cos(x) )$
(%i46) quad_qags(g(t),t,0,%pi/2);
(%i47) integrate(g(t),t,0,%pi/2);
(%i48) tval : bfloat(%),fpprec:45;
(%i49) rectform(%);
(%i50) float(-%pi*log(2)/2);
(%i51) tval : bfloat(-%pi*log(2)/2),fpprec:45;
(%i52) quad_ts(g,0,%pi/2,30,40);
(%i53) ans: realpart( first(%) ),fpprec:45;
(%i54) abs(ans - tval),fpprec:45;
(%i55) qtsk_e(g,0,%pi/2,5,40);
------------------------------------------------------------------
3.3   The Gauss-Legendre Quadrature Method for a <= x <= b

(%i1) load("quad-maxima.lisp");
(%i2) fpprintprec:8$
(%i3) tval : bfloat(integrate(exp(x),x,-1,1)),fpprec:45;
(%i4) load(quad_gs);
(%i5) arrays;
(%i6) arrayinfo(ab_and_wts);
(%i7) gaussunit(exp,4);
(%i8) abs(% - tval),fpprec:45;
(%i9) fpprec;
(%i10) arrayinfo(ab_and_wts);
(%i11) first( ab_and_wts[4, 16] );
(%i12) second( ab_and_wts[4, 16] );
(%i13) lp4 : legenp(4,x);
(%i14) float(solve(lp4));
(%i15) ab_and_wts[4,16];
(%i16) lfp4(x) := legenp(4,x)$
(%i17) map('lfp4,%o11);
(%i18) map( lambda([z], legenp (4, z)),%o11 );
(%i19) gaussunit_e(exp,4);
(%i20) arrayinfo(ab_and_wts);
(%i21) gaussab(exp,-1,1,4);
(%i22) abs(% - tval),fpprec:45;
(%i23) gaussab_e(exp,-1,1,4);
(%i24) arrayinfo(ab_and_wts);
(%i25) quad_gs(exp,-1,1,10);
(%i26) abs(first(%) -tval),fpprec:45;
(%i27) arrayinfo(ab_and_wts);
(%i28) fpprec:30$
(%i29) quad_gs(exp,-1,1,20);
(%i30) abs(first(%) -tval),fpprec:45;
(%i31) arrayinfo(ab_and_wts);
(%i32) fpprec:40$
(%i33) quad_gs(exp,-1,1,30);
(%i34) abs(first(%) -tval),fpprec:45;
(%i35) arrayinfo(ab_and_wts);
(%i36) gaussab_e(exp,-1,1,40);
(%i37) arrayinfo(ab_and_wts);
(%i38) quad_gs_table(exp,-1,1,30)$
(%i39) arrayinfo(ab_and_wts);