(
/* 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 ceil(logb(2,fpprec)) do /* should be enough for full prec. */
       ([val,deriv]:legenpd(n,v), 
        v:v-val/deriv),
    push(v,abscissae),
    push(legenwt(n,v),weights)),
    if oddp(n) then (push(0,abscissae),
                     push(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)) */
/*
legenwt(k,x):=  
 block([t0:1,t1:x,ans:0,pkp1,dpk],
   for i:2 thru k do
    (ans: ((2*i-1)*x*t1 -(i-1)*t0)/i,
     t0:t1,
     t1:ans),
   pkp1:((2*k+1)*x*t1 -(k+1)*t0)/(k+1), dpk: k*(x*t1-t0)/(x^2-1),
 -2/((k+1)*dpk*pkp1)), */
/*Here's a correct formula .. */

legenwt(k,x):=2/(1-x^2)/legendrediff[k](x)^2,
legendrediff[i](z):=diff(legendre_p[i](z),z),
/* another formula is 
legenwt2(k,x):= 2*(1-x^2)/k^2/(legendre_p[k-1](x)-x*legendre_p[k](x))^2,
*/

/*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 */
  gq1(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)),

/*Tanh/Sinh method, at least the basics. See the BJL paper and its code for
arranging phases so that the previous phase's result is updated, going from
n points to 2n to 4n etc points. */

quadtsN(%%g,n):= /* Naive Tanh/Sinh method, integral of g from -1 to 1, 2n+1 points. */
 block([sum:0, piby2: bfloat(%pi/2), h:4/(n+2), hj, u1,u2, weight],
  for j from -n thru n do 
  (hj: bfloat(j*h),
   u1: piby2*cosh(hj),
   u2: piby2*sinh(hj),
   weight:u1/cosh(u2)^2,
   sum: sum + weight*%%g(bfloat( tanh(u2)))),
   h*sum),

 quadts(%%g,n):= /* Tanh/Sinh method, integral of g from -1 to 1 */
 block([sum:0, piby2: bfloat(%pi/2), h:4/n, hj, u1,u2, weight, ab],
   sum:piby2*%%g(0),
  for j from 1 thru n do 
  (hj: bfloat(j*h),
   t2:exp(hj),
   u1: piby2*(t2+1/t2)/2,   /*cosh */
   u2: piby2*(t2-1/t2)/2,   /*sinh */
   t3: exp(u2), 
   t4:(t3+1/t3)/2,         /* cosh(u2) */
   t3:(t3-1/t3)/2,         /* sinh(u2) */
   weight:u1/(t4*t4),
   ab: t3/t4,              /*tanh(u2) */
   sum: sum + weight*(%%g(ab)+%%g(-ab))),
   h*sum),

/* return a list of the TS approximation and the fp error */

 quadts_e(%ggg,n):=
  quadts(lambda([%z],uncert(%ggg,[%z])),n)
)
$