/************************************************************************
    quad_de.mac is a software file which accompanies Chapter 9 of
    Maxima by Example, Bigfloats and High Accuracy Quadrature.
    
    Copyright (C) 2009, 2016
    Edwin L Woollett,  woollett@charter.net
    http://www.csulb.edu/~woollett
    
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU GENERAL PUBLIC LICENSE, Version 2, June 1991,
    as published by the Free Software Foundation.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details. You should have received a
    copy of the GNU General Public License along with this program.
    If not, see http://www.fsf.org/licensing/.
************************************************************************/





/* file quad_de.mac integrates over [a,inf] using a
    [a,inf] domain double exponential method. After transforming
    to [0,inf] the transformation from x to u is
    x = exp(u - exp(-u) ) suitable for integrands which
    already have some kind of exponential decay in their
    integrands.
    functions defined in this file
     1. de_sum(%%g,kk)

     2. idek(%fg,a,kk,fp)

     3. fdf (%ff, %xx, %dfp)
     
     4. idek_e(%fg,a,kk,fp)

     5. ide(%fg,a,ra,fp)

     6. quad_de(%%fg,a,ra,fp)

     7. ide_test(%fg,a,ra,fp)
     
     8. bfprint(bf,fpp)

    
     */

/********************************************/

 /* de_sum(f,k) does sums needed for integration
    of f over [0,inf] using h:1/k. 
     de_sum(f,k ) should be called by a function which
       has set fpprec.              
  de_sum adds up contributions for a given
     level k, skipping every other grid point
     unless k = 1,  and stops when abs(newterm) <= eps1 
        */
    

de_sum(%%g,kk ) :=
 block([fp,h,jstep,eps1,jmax:500,u,a,
          x,w,asum,nt,ant ],
    fp:fpprec,
   
   if kk = 1 then jstep:1 else jstep:2,   
   
   h : 1/2^kk,
   eps1: bfloat(10^(-fpprec)),
   
   asum : 0.0b0,
   
       /* sum first over negative u values */
    
    for j:1 step jstep thru (jmax+4) do 
     (  block([fpprec],
          fpprec : fp + 40,    
          u: bfloat(-h*j),    
          a:exp(-u),
          x:exp(u-a)),
          
       w:exp(-a) + x,
       nt :bfloat( %%g(x)) * w ,        
       asum: asum + nt ,        
       
       if listp(nt) then ant:abs(first(nt)) else ant:abs(nt),       
       
       if j = jmax then 
         (print("j = jmax limit reached"),return(done)),
             
        if ant < eps1  then return(done)  ), 
        

     
     /* now sum over positive values of u  */
    

    
    
    for j:1 step jstep thru (jmax+4) do 
     ( u: bfloat(h*j),    
       a:exp(-u),
       x:exp(u-a),
       w:exp(-a)+x,
       nt :bfloat( %%g(x)) * w ,        
       asum: asum + nt ,
       if listp(nt) then ant:abs(first(nt)) else ant:abs(nt),                    
              
       if j = jmax then 
         (print("j = jmax limit reached"),return(done)),
         
        if ant < eps1  then return(done)  ),
        


     asum  )$  
     

/* idek(f,a,kk,fp) computes integral of
    f over [a,inf] using h :1/kk, and with
    fpprec set to fp. This functions calls
    de_sum(f,k) for k = 1,2,...,kk.
    */


     
idek(%fg,a,kk,fp) :=
 block(  [fpprec,asum,h ], local(%ft),   
 
     if ra > fp then
          (print(" ra should be less than fp"), return(done) ),
     fpprec:fp,
     h : 1/2^kk,
     
     /* convert to integration over [0,inf] */
     
     block([fpprec,beta],
       fpprec : fp + 20,         
       beta:bfloat(a),    
       define (%ft(x), %fg(x + beta) )),
       
      /* contribution from u = 0 */
      asum : bfloat(%ft(exp(-1))*2*exp(-1)) ,
      for k:1 thru kk do
        asum : asum + de_sum(%ft,k),
      bfloat(h*asum) )$
      
        
/*  idek_e(f,a,k,fp) uses fdf(f,x,dfp) to 
    return [value-of-integral, estimate-of-error-due-to-floating-point-setting]
    based on Richard Fateman's code
    */

 
idek_e(%fg,a,kk,fp) :=
 block(  [fpprec,asum,h,x0,w0 ], local(%ft),   
 
     if ra > fp then
          (print(" ra should be less than fp"), return(done) ),
     fpprec:fp,
     h : 1/2^kk,
     
     /* convert to integration over [0,inf] */
     
     block([fpprec,beta],
       fpprec : fp + 20,         
       beta:bfloat(a),    
       define (%ft(x), %fg(x + beta) )),
       
      /* contribution from u = 0 */
      
      x0 : bfloat(exp(-1)),
      w0 : 2*x0,
      asum : bfloat( fdf(%ft,x0,10)*w0 ),
      
      /* asum : bfloat(%ft(exp(-1))*2*exp(-1)) ,*/
      
      for k:1 thru kk do
        asum : asum + de_sum(lambda([%z],fdf(%ft,%z,10)),k),
      bfloat(h*asum) )$      
      
/* ide(f,a,ra,fp) integrates f over [a,inf] with
   requested accuracy ra and  with fpprec : fp
   calls de_sum(f,k) for sums due to contributions
   from levels k = 1,2,...,kmax. Search ends when
   the absolute value of the difference of the 
   integral value found is less than
   10^(-ra).
   */      


ide(%fg,a,ra,fp):= 
  block([fpprec,fpprintprec,eps0,asum,h, oldval,newval,kmax,
       vdiff,vdiffold  ], local(%ft),
       
     if ra > fp then
          (print(" ra should be less than fp"), return(done) ),
     fpprec:fp,
     
     /* convert to integration over [0,inf] */
     
     block([fpprec,beta],
       fpprec : fp + 20,         
       beta:bfloat(a),    
       define (%ft(x), %fg(x + beta) )),
     
     eps0:bfloat(10^(-ra)),
     fpprintprec:8,
     kmax : 8,
     
     /* asum from u = 0 */
     
     asum : bfloat(%ft(exp(-1))*2*exp(-1)) ,
     /* print("in ide, asum u = 0  = ",asum ), */
     
     print("  "),
     print("      requested accuracy = ",ra," fpprec = ",fpprec ),
     
     /* start with k = 1, h = 1/2 value */
     h : 1/2,
          
     asum : asum + de_sum(%ft,1),      
         
     
     oldval : bfloat(asum*h),
     
     
     print(" k     value           vdiff  "),
     print("  "),
     printf(true,"~2d    ~14a  ~%",1,string(oldval) ), 
     
     vdiff : 1.0b0,
     vdiffold : 2.0b0,     

     for k:2 thru (kmax+2) do
      ( if k > kmax then 
            (print("k limit reached"),return(done) ),
        h : h/2,
        asum : asum + de_sum(%ft,k),                    
        newval:bfloat(asum*h),       
        vdiff : abs(newval - oldval),        
        printf(true,"~2d    ~14a ~14a  ~%",k,string( expand(newval) ),string(vdiff) ),
        oldval : newval,
        if vdiff > vdiffold then
         ( print("vdiff > vdiffold"), return(done)),
        if vdiff <= eps0 then return(done),
        vdiffold : vdiff ),
        print(" ") )$   
        
        
 /* quad_de based on quad_de1 and ide */
 
  /*    
disp("quad_de(f,a,ra,fp) follows the same route
        as ide(f,a,ra,fp), but instead of printing
        a table, just returns the list [ newval, k-final, vdiff]
         containing the end result. ")$
 */
          
          
quad_de(%%fg,a,ra,fp):= 
  block([fpprec,eps0,oldval,newval,kval,kmax,asum,h,
       vdiff,vdiffold ], local(%ft),
       /* display(a,ra,fp), */
       
     if ra > fp then
          (print(" ra should be less than fp"), return(done) ),
     fpprec:fp,
     block([fpprec,beta],
       fpprec : fp + 20,         
       beta:bfloat(a),    
       define (%ft(x), %%fg(x + beta) )),
       
       /* display(fp,%ft(y)), */
     
     eps0:bfloat(10^(-ra)),
     kmax : 8,
     asum : bfloat(%ft(exp(-1))*2*exp(-1)) ,
     h : 1/2,          
     asum : asum + de_sum(%ft,1),       
     oldval : bfloat(asum*h),
     
     vdiff : 1.0b0,
     vdiffold : 2.0b0,     

     for k:2 thru (kmax+2) do
      ( if k > kmax then 
            (print("k limit reached"),return(done) ),
        kval : k,
        h : h/2,
        asum : asum + de_sum(%ft,k),                    
        newval:bfloat(asum*h),       
        vdiff : abs(newval - oldval),                
        if vdiff > vdiffold then
         ( print(" vdiffnew > vdiffold before vdiff < eps0 reached"),
           print(" return previous level result "),
           newval:oldval,vdiff:vdiffold,kval:kval-1, return(done)),
        if vdiff <= eps0 then return(done),        
        oldval : newval,
        vdiffold : vdiff ),
        
        [ newval, kval,vdiff] )$   
         
      
        
     
/* ide_test(f,a,ra,fp) follows ide(f,a,ra,fp)
    but adds to the table the error compared
    with the true value which must be bound
    to the global symbol tval.
   */    
     
ide_test(%fg,a,ra,fp):= 
  block([fpprec,fpprintprec,eps0,asum,h, oldval,newval,kmax,
       vdiff,vdiffold  ], local(%ft),
       
     if ra > fp then
          (print(" ra should be less than fp"), return(done) ),
     fpprec:fp,
     
     /* convert to integration over [0,inf] */
     
     block([fpprec,beta],
       fpprec : fp + 20,         
       beta:bfloat(a),    
       define (%ft(x), %fg(x + beta) )),
     
     eps0:bfloat(10^(-ra)),
     fpprintprec:8,
     kmax : 8,
     
     /* asum from u = 0 */
     
     asum : bfloat(%ft(exp(-1))*2*exp(-1)) ,
     /* print("in ide, asum u = 0  = ",asum ), */
     
     print("  "),
     print("      requested accuracy = ",ra," fpprec = ",fpprec ),
     
     /* start with k = 1, h = 1/2 value */
     h : 1/2,
          
     asum : asum + de_sum(%ft,1),      
         
     
     oldval : bfloat(asum*h),

     
     verr : abs(oldval - tval),
     print(" k     value           vdiff            verr  "),
     print("  "),
     printf(true,"~2d    ~14a ~14a   ~14a  ~%",1,string(oldval),"              ",string(verr)), 
     
     vdiff : 1.0b0,
     vdiffold : 2.0b0,     

     for k:2 thru (kmax+2) do
      ( if k > kmax then 
            (print("k limit reached"),return(done) ),
        h : h/2,
        asum : asum + de_sum(%ft,k),                    
        newval:bfloat(asum*h),       
        vdiff : abs(newval - oldval),        
        verr : block([fpprec:fp+3],abs(newval - tval)),
           
        printf(true,"~2d    ~14a ~14a   ~14a  ~%",k,string( expand(newval) ),string(vdiff),string(verr)),        
        
        
        oldval : newval,
        if vdiff > vdiffold then
         ( print("vdiff > vdiffold"), return(done)),
        if vdiff <= eps0 then return(done),
        vdiffold : vdiff ),
        print(" ") )$   
        


/***********************************************/     

   
/*
disp(" fdf(f,x,dfp) finds the absolute value
       of the difference of f(x) at the current
     value of fpprec and at the value (fpprec+dfp),
     and returns [f(x), df(x)].
     This code is adapted from Richard
     Fateman's uncert(f,x) function. ")$
 */
 
     
fdf (%ff, %xx, %dfp) :=
   block([fv1,fv2,df],
    fv1:bfloat(%ff(bfloat(%xx))),
    block([fpprec:fpprec + %dfp ],
     fv2:  bfloat(%ff(bfloat(%xx))),
     df: abs(fv2 - fv1) ),
   [bfloat(fv2),bfloat(df)] )$     
   
/******************************************/     
   
 
                
        
bfprint(bf,fpp):=
  block([fpprec, fpprintprec ],
        fpprec : fpp+2,
        fpprintprec:fpp,
        print("   number of digits = ",fpp),
        print("  ",bf) )$