[Maxima-commits] CVS: maxima/share/contrib/levin levin.dem, NONE, 1.1 levin.mac, NONE, 1.1 rtest_levin.mac, NONE, 1.1

[Robert D. <rob...@us...>]

Update of /cvsroot/maxima/maxima/share/contrib/levin
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv13276

Added Files:
	levin.dem levin.mac rtest_levin.mac 
Log Message:
Levin method for acceleration of summations, written by Michel Van den Bergh (mic...@uh...)
and downloaded by me (Robert Dodier) on 2006/12/28 from http://thegate.bertold.org:8081/levin/,
and committed verbatim.

--- NEW FILE: levin.dem ---
fpprec:16$
/* load necessary file */ 
load("levin.mac")$
/* do some random test */
bfloat(bflevin_u_sum(1/n^2,n,1)-%pi^2/6);
/* increase precision */ 
fpprec:50$
bfloat(bflevin_u_sum(1/n^2,n,1)-%pi^2/6);
/* some series with zero terms */
/* this is a bit of a hack */
pi_4(k):=(k-2*fix(k/2))*(-1)^((k-1)/2)/k$
bfloat(bflevin_u_sum(pi_4(k),k,1)-%pi/4);
/* this one's more natural */
pi_4_alt1(k):=if oddp(k) then (-1)^((k-1)/2)/k else 0$
/* but now we need to protect pi_4_alt1(k) against premature evaluation */
bfloat(bflevin_u_sum(funmake(pi_4_alt1,[k]),k,1)-%pi/4);
/* this one works better but is more complicated */
pi_4_alt2(k):=if oddp(k) then (-1)^((k-1)/2)/k else if evenp(k) then 0 else funmake(pi_4_alt2,[k])$
bfloat(bflevin_u_sum(pi_4_alt2(k),k,1)-%pi/4);
/* the following expression shows that we really use three values of a */ 
ratsimp(levin_u_sum(a[n],n,0,3));
/* approximation to the zeta function */
levin_zeta(s):=levin_u_sum(1/n^s,n,1,5)$
plot2d(levin_zeta,[s,1.1,2.0],[y,0.0,11.0])$
/* The sum sum(1/n-log(1+1/n),n,1,inf)  computes the
Euler-Mascheroni constant */
bflevin_u_sum(1/n-log(1+1/n),n,1);
/* oops, we need to use levin_u_sum. */
/* set precision very high */
fpprec:60$
/* compute sum with 30 and 31 terms respectively */
levin_gamma30:bfloat(levin_u_sum(1/n-log(1+1/n),n,1,30));
levin_gamma31:bfloat(levin_u_sum(1/n-log(1+1/n),n,1,31));
/* by the following computation we guess that levin_gamma31 is accurate
up to 27 decimal places */
bfloat(abs(levin_gamma31-levin_gamma30));
/* this turns out to be correct */
bfloat(levin_gamma31-%gamma);


/* reset precision */
fpprec:16$

--- NEW FILE: levin.mac ---
/* levin.mac v0.2.1 (c) Michel Van den Bergh   Licence: GPL.
0.0->0.1 Some suggestions by Richard Fateman incorporated
0.1->0.2 Now much faster. Zero entries are allowed.
0.2->0.2.1 Disable bflevin_u_sum for series with non rational terms
*/

/* see:  "Mathematical Properties of a
New Levin-Type Sequence Transformation"
by Ernst Joachim Weniger
Formula: (2.18)
arXiv:math-ph/0306063 */

levin_debug:false;

/* s     : array with partial sums of a series
   k+1   : number of partial sums to use
   n     : index of first partial sum to be used;
           in other words, we use: s[n],...,s[n+k]
   beta  : translation parameter; this parameter can be used to avoid division
           by zero;
   omega : an array of remainder estimates;
           depending on omega one gets a different type of levin
           transform; omega should not contain zero entries 
           or we get a division by zero */


levin(s,k,n,beta,omega):=block([i,j,numerator,denonimator],
  numerator:sum(((-1)^j)*binomial(k,j)*
             ((beta+n+j)^(k-1)/(beta+n+k)^(k-1))*
             (subvar(s,n+j)/subvar(omega,n+j)),
             j,0,k),
  denominator:sum(((-1)^j)*binomial(k,j)*
             ((beta+n+j)^(k-1)/(beta+n+k)^(k-1))*
             (1/subvar(omega,n+j)),
             j,0,k),
  numerator/denominator
);

/* Here we delete zero entries and compute the partial sums. We also
   set omega such that we get the levin_u transform.  In addition we 
   set beta=0 and n=1. Perhaps we should try several different values
   for beta in case division by zero occurs. */

levin_u(a,terms):=block([i,j,l,n,beta,last,omega_,t_,b_],
  local(omega_,t_,b_),
  n:1,
  beta:0,
  /* remove zeros */
  j:0,
  for i:0 thru terms-1 do (
      if sign(abs(subvar(a,i)))#zero then (
      b_[j]:subvar(a,i),
      j:j+1
    ) else (
      if levin_debug then print("ommitting zero term with index ",i)
    )
  ),
  last:j-1,
  omega_[l]:=l*b_[l],
  t_[l]:=t_[l-1]+b_[l],
  t_[-1]:0,
  levin(t_,last-n,n,beta,omega_)
);

/* User friendly version.

   Example
                               3899836039
   levin_u_sum(1/n^2,n,1,10) = ---------- = 1.644934081345832
                               2370816000
                       1968329
   sum(1/n^2,n,1,10) = ------- = 1.549767731166541
                       1270080

   sum(1/n^2,n,1,inf) = %pi^2/6 = 1.644934066848226  */

levin_u_sum(a_,index_var,start,terms):=block([n,c_,dummy],
  local(c_),
  for n:0 thru terms-1 do (
             dummy:subst(start+n,index_var,a_),
             c_[n]:ev(dummy,eval) ),
   levin_u(c_,terms) 
);

/* Numerical version. This function makes a crude attempt to compute
   an indefinite sum numerically with precision given by fpprec.
   This works (usually) for functions that take rational values only!
   For other functions the intermediate computations much be done in much
   greater precision. In that case use levin_u_sum with an appropriate number of terms.
   If you suspect something is wrong put levin_debug:true.
*/

bflevin_u_sum(a_,index_var,start):=block([n,fpprec_save,old,new,value,first_time,
               count,terms,begin_term,end_term,dummy,non_rational,c_],
  local(c_),
  fpprec_save:fpprec,
  fpprec:fpprec+6,
  first_time:true,
  count:0,
  terms:5,
  non_rational:false,
  value:do (
    if first_time then (
      begin_term:0,
      end_term:terms-1
    ) else (
      begin_term:terms/2,
      end_term:2*terms-1
    ),
    for n:begin_term thru end_term do (
                   dummy:subst(start+n,index_var,a_),
                   c_[n]:ev(dummy,eval),
                   if not(ratnump(c_[n])) and not(non_rational) then (
                     print(
"*** bflevin_u_sum does not work (yet) for series with non rational terms. ***"
                     ),
                     non_rational:true,
                     return(failed)
                   )
               ),
    if non_rational then return(bflevin_u_sum_failed),
    new:levin_u(c_,terms),
    new:bfloat(new),
    if not(first_time) then (
      if levin_debug then print("terms=",terms, " approximation=",new),
      if equal(old,bfloat(0)) and equal(new,bfloat(0)) then return(bfloat(0)),
      if abs(old-new)/max(abs(old),abs(new))<10^(-fpprec_save-3) then return(new),
      if (count:count+1)>10 then return(bflevin_u_sum_failed)
    ),
    terms:2*terms,
    old:new,
    first_time:false
  ),
  fpprec:fpprec_save,
  bfloat(value)
);



--- NEW FILE: rtest_levin.mac ---
kill(all);
done$
block([],load("levin.mac"),done);
done$
fpprec:16;
16$
is(bflevin_u_sum(1/n^2,n,1)-%pi^2/6<10^(-fpprec+1));
true$
fpprec:30;
30$
is(bflevin_u_sum(1/n^2,n,1)-%pi^2/6<10^(-fpprec+1))$
true$
pi_4(k):=(k-2*fix(k/2))*(-1)^((k-1)/2)/k;
pi_4(k):=(k-2*fix(k/2))*(-1)^((k-1)/2)/k$
is(bflevin_u_sum(pi_4(k),k,1)-%pi/4<10^(-fpprec+1));
true$
pi_4_alt1(k):=if oddp(k) then (-1)^((k-1)/2)/k else 0;
pi_4_alt1(k):=if oddp(k) then (-1)^((k-1)/2)/k else 0$
is(bflevin_u_sum(funmake(pi_4_alt1,[k]),k,1)-%pi/4<10^(-fpprec+1));
true$
pi_4_alt2(k):=if oddp(k) then (-1)^((k-1)/2)/k else if evenp(k) then 0 else funmake(pi_4_alt2,[k]);
pi_4_alt2(k):=if oddp(k) then (-1)^((k-1)/2)/k else if evenp(k) then 0 else funmake(pi_4_alt2,[k])$
is(bflevin_u_sum(pi_4_alt2(k),k,1)-%pi/4<10^(-fpprec+1));
true$
sort(listofvars(levin_u_sum(a[n],n,0,10)));
[a[0],a[1],a[2],a[3],a[4],a[5],a[6],a[7],a[8],a[9]]$
levin_zeta(s):=levin_u_sum(1/n^s,n,1,10);
levin_zeta(s):=levin_u_sum(1/n^s,n,1,10)$
block([],load("bffac"),done);
done$
fpprec:16$
16;
is(levin_zeta(1.1)-bfzeta(1.1,10)<10^-5);
true$
is(levin_zeta(1.01)-bfzeta(1.01,10)<10^-4);
true$
bflevin_u_sum(1/n-log(1+1/n),n,1);
bflevin_u_sum_failed$
fpprec:60;
60$
is(levin_u_sum(1/n-log(1+1/n),n,1,31)-%gamma<10^(-29));
true$