( /* Integrate function over n points from -1 to 1. Compute abscissas and weights if not already available. Precision is that of the (global) setting for fpprec. */ gaussunit(%ggg,n):= gaussunit1(%ggg,n,ab_and_wts[n,fpprec]), gaussunit1(%%g,n,aw):= /* the function summing the terms */ block([sum:0], map(lambda([%%a,%%w],sum:sum+%%w*(%%g(%%a)+%%g(-%%a))),aw[1],aw[2]), sum+ (if oddp(n) then -%%g(0)*aw[2][1] else 0 )), /* Integrate function from x=lo to x=hi using n points and fpprec */ gaussab(%hh,lo,hi,n):= block([a:(hi-lo)/2, b:(hi+lo)/2], local(%zz), define (%zz(x),%hh(a*x+b)), /* transformed function */ a* gaussunit(%zz,n)), /* Background subroutines needed for Gauss integration */ legenpd(k,x):= /* return P[k](x) and P'[k](x) */ if (k=0) then [1,0] else if (k=1) then [x,1] else block([t0:1,t1:x,ans:0], for i:2 thru k do (ans: ((2*i-1)*x*t1 -(i-1)*t0)/i, t0:t1, t1:ans), [t1, k*(x*t1-t0)/(x^2-1)]), legenp(k,x):= /* return P[k](x) */ if (k=0) then 1 else if (k=1) then x else block([t0:1,t1:x,ans:0], for i:2 thru k do (ans: ((2*i-1)*x*t1 -(i-1)*t0)/i, t0:t1, t1:ans), t1), /* return a pair of lists of abscissae and weights */ ab_and_wts[n, fpprec]:= block([a:0,v:0.0d0,np1:n+1,nph:1.0d0/(n+1/2),halfn:floor(n/2), val,deriv,abscissae:[],weights:[],float2bf:true], for i:0 thru halfn-1 do (v:cos(?pi*((i+3/4)*nph)), /*an approx zero of legendre-p[n] */ /*refine it by Newton iteration */ v:bfloat(v), for k:0 thru ceiling(log(fpprec)/log(2)) do /* should be enough for full prec. */ ( /*[val,deriv]:legenpd(n,v), ..if syntax supported... else*/ deriv:legenpd(n,v), val:deriv[1], deriv:deriv[2], v:v-val/deriv), abscissae:cons(v,abscissae), weights:cons (legenwt(n+1,v),weights)), if oddp(n) then (abscissae:cons(0,abscissae), weights:cons(2-2*apply("+",weights), weights)), [abscissae,weights]), /* The weight at x=root of _P[k]: w[k]:= -2/ ( (k+1)* P'[k](x)*P[k+1](x)). call on k+1, and compute compute -2/(k*P'[k-1]*P[k]) */ legenwt(k,x):= block([t0:1,t1:x,ans:0,pkp1,dpk], for i:2 thru k-1 do (ans: ((2*i-1)*x*t1 -(i-1)*t0)/i, t0:t1, t1:ans), pk:((2*k-1)*x*t1 -(k-1)*t0)/k, /*P[k](x)*/ dpkm1: (k-1)*(x*t1-t0)/(x^2-1), /*P'[k-1](x)*/ -2/(k*dpkm1*pk)), /*UNCERT: a Macsyma program to provide value of function f and uncertainty, heuristic. uncert(f,v) evaluates function f at (in general, vector) point v at some floating-point precision somewhat in excess of current setting of fpprec. It returns a list of two items, [y,u], where y is approximate value of f(v), and u = (nonnegative) uncertainly in the provided value y. Both y and u are bigfloats. Some functions may be devious enough as to mislead this calculation, but this should be exceedingly rare. Example. [q(x):=1/(asin(atan(x))-atan(asin(x))), uncert(q,[1/10000]), bfloat(q(1/10000))]; Try the above variously with fpprec:20, fpprec:100 */ bfapply(%fun,%args,fpprec):= apply(%fun,map(bfloat,%args)), uncert(%fun,%args):= block([%ll, %hh,%dd,oldprec:fpprec], %ll: bfapply(%fun,%args,fpprec), fpprec: fpprec+10, %hh: bfapply(%fun,%args,fpprec), %dd: abs(%hh-%ll), fpprec: oldprec, [bfloat(%hh), bfloat(%dd)]), /* putting this together with quadrature */ gaussunit_e(%ggg,n):= /*with error */ gaussunit1(lambda([%z],uncert(%ggg,[%z])), n,ab_and_wts[n,fpprec]), gaussab_e(%hh,lo,hi,n):= block([a:(hi-lo)/2, b:(hi+lo)/2], local(%zz), define (%zz(x),%hh(a*x+b)), a* gaussunit_e(%zz,n)), quadts(%%g,n):= /* Tanh/Sinh method, integral of g from -1 to 1 */ block([sum:0, piby2: bfloat(%pi/2), h:4/n, he, t2:1, t3,t4, ab, correction, oldsum:0], sum:piby2*%%g(0), he:bfloat(exp(h)), for j from 1 thru 2*n do (t2:he*t2, /*exp(h*j)*/ t3: exp( piby2*(t2-1/t2)/2), t4:(t3+1/t3)/2, /* cosh */ ab: (t3-1/t3)/2/t4, /* tanh */ correction: ((piby2*(t2+1/t2)/2)/(t4*t4))*(%%g(ab)+%%g(-ab)), oldsum:sum, sum:sum+correction, if abs((sum-oldsum)/oldsum)<10^(-fpprec) then return(sum*h)), h*sum), /* return a list of the TS approximation and the fp error */ quadts_e(%ggg,n):= quadts(lambda([%z],uncert(%ggg,[%z])),n) )\$