share/contrib/levin/levin.mac

   1 /*               COPYRIGHT NOTICE
   2
   3 Copyright (C) 2006-2007 Michel Van den Bergh
   4
   5 This program is free software; you can redistribute
   6 it and/or modify it under the terms of the
   7 GNU General Public License as published by
   8 the Free Software Foundation; either version 2
   9 of the License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it
  12 will be useful, but WITHOUT ANY WARRANTY;
  13 without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 GNU General Public License for more details at
  16 http://www.gnu.org/copyleft/gpl.html
  17 */
  18
  19 /* The latest version of this package can be found at
  20           http://www.bertold.org/levin
  21 */
  22
  23
  24
  25 /*
  26 0.0->0.1 Some suggestions by Richard Fateman incorporated
  27 0.1->0.2 Now much faster. Zero entries are allowed.
  28 0.2->0.2.1 Disable bflevin_u_sum for series with non rational terms
  29 0.2->0.3 Argument treatment somewhat more robust. Added GPL header.
  30 0.3->0.4 Extra checks on arguments.
  31          Variable lower bound for levin_u_sum is treated correctly.
  32 0.4->0.5 Got rid of ugly hacks and simply renamed parameters to "non-obvious"
  33          (-)__ names (see discussion on list). This simulates lexical scoping
  34          as long as the user does not use such names.
  35 0.5->0.6 Faster again (for series with rational terms)!
  36          "levin" can now compute progagation of rounding errors.
  37          "bflevin_u_sum" now works for series with non rational terms.
  38 0.6->0.7 Propagation of rounding errors is computed in lower precision.
  39          This gains some time.
  40          Added disclaimer about optimization.
  41 */
  42
  43 /* Semantic note: all functions evaluate all their arguments.
  44    levin_u_sum and bflevin_u_sum perform one round of ev(-) on their
  45    first argument. This is to make (bf)levin_u_sum('(...),...) behave as
  46    expected. Without this, some regression tests fail. Since we simulate
  47    lexical scope using "non-obvious" parameter names, the use of ev seems
  48    ok. I.e. in normal situations there is no risk of ev evaluating things
  49    in the wrong scope.
  50 */
  51
  52 /* Disclaimer: The following code is a proof of concept. It works
  53    reliably but is highly non optimized with regard to speed and memory
  54    usage. The necessary optimizations are quite straight forward however.
  55 */
  56
  57 levin_options["debug"]:false;  /* when true gives some extra output; should
  58                                         really be a loglevel */
  59 levin_options["min_terms"]:5;    /* bflevin_u_sum starts with 5 terms */
  60 levin_options["max_terms"]:640;         /* bflevin_u_sum tries up to 5*2^7 terms
  61                                          before bailing out */
  62 levin_options["min_precision"]:16;  /* bflevin_u_sum will start with this
  63                                            precision */
  64 levin_options["max_precision"]:1000;  /* bflevin_u_sum will not increase the
  65                                            precision up above this value*/
  66
  67 /* see:  "Mathematical Properties of a
  68 New Levin-Type Sequence Transformation"
  69 by Ernst Joachim Weniger
  70 Formula: (2.18)
  71 arXiv:math-ph/0306063 */
  72
  73 /* s     : array with partial sums of a series
  74    k+1   : number of partial sums to use
  75    n     : index of first partial sum to be used;
  76            in other words, we use: s[n],...,s[n+k]
  77    beta  : translation parameter; this parameter can be used to avoid division
  78            by zero; for accurate error estimates it currently has to be
  79            an integer;
  80    omega : an array of remainder estimates;
  81            depending on omega one gets a different type of levin
  82            transform; omega should not contain zero entries
  83            or we get a division by zero
  84    mode  : levin_numeric or levin_algebraic, for levin_numeric the return
  85            value is a list [result,variance]. For levin_algebraic the
  86            return value is just the result.
  87    s__variance     : in numeric mode an array with variances for the elements
  88                      of s
  89    omega__variance : in numeric more an array with variances for the elements
  90                      of omega
  91
  92    Note 1 : All variances are assumed to be in units of 10^(-2*fpprec).
  93
  94    Note 2 : The computation of propagation of variance is not quite
  95    mathematically correct since we are assuming that the errors on the
  96    s[j] and omega[j] are all independent.  However my tests indicate that this
  97    method is reasonably accurate. Note that the correct computation of
  98    propagation of variance is O(N^2).
  99 */
 100
 101
 102 levin(s__,k__,n__,beta__,omega__,mode__,s__variance__,omega__variance__):=
 103        block([i__,j__,l__,numerator__,denominator__,result__,result__variance__,
 104              coeffs__,numerator__variance__,denominator__variance__
 105              ],
 106   local(coeffs__),
 107   coeffs__[l__]:=((-1)^l__)*binomial(k__,l__)*(beta__+n__+l__)^(k__-1),
 108 /* the following sums will be rewritten recursively */
 109   numerator__:sum(coeffs__[j__]*
 110              (
 111                    s__[n__+j__]/omega__[n__+j__]
 112              ),
 113              j__,0,k__),
 114   denominator__:sum(coeffs__[j__]*
 115              (
 116                    1/omega__[n__+j__]
 117              ),
 118              j__,0,k__),
 119   result__:errcatch(numerator__/denominator__),
 120   if result__=[] then (
 121      print (
 122 "*** levin has zero denominator ***"
 123      ),
 124      return(levin_failed)
 125   ),
 126   result__:result__[1],
 127   if mode__=levin_numeric then (
 128         block([fpprec:16], /* start lower precision block */
 129
 130 /* the following formulas are not optimized yet, they are written with an
 131    emphasis on readability, the sums will be rewritten recursively */
 132
 133      numerator__variance__:sum(coeffs__[j__]^2*
 134         (
 135            (1/omega__[n__+j__])^2 * s__variance__[n__+j__] +
 136            (s__[n__+j__]/omega__[n__+j__]^2)^2 * omega__variance__[n__+j__]             ),
 137              j__,0,k__),
 138      denominator__variance__:sum(coeffs__[j__]^2*
 139         (
 140            (1/omega__[n__+j__]^2)^2 * omega__variance__[n__+j__]
 141         ),
 142              j__,0,k__),
 143      result__variance__:
 144               (     (1/denominator__^2) * numerator__variance__   +
 145                    (numerator__/denominator__^2)^2 * denominator__variance__
 146               )
 147      ) /* end lower precision block */
 148   ),
 149   if mode__=levin_algebraic then
 150      return(result__)
 151   else
 152      return([result__,result__variance__])
 153 );
 154
 155
 156
 157
 158
 159 /* Here we delete zero entries and compute the partial sums. We also
 160    set omega such that we get the levin_u transform.  In addition we
 161    set beta=0 and n=1. Perhaps we should try several different values
 162    for beta in case division by zero occurs. */
 163
 164 levin_u(a__,terms__,[mode__]):=block([i__,j__,l__,n__,beta__,last__,
 165                             omega__,t__,b__,t__variance__,omega__variance__
 166                             ],
 167   local(omega__,t__,b__,t__variance__,omega__variance__),
 168   n__:1,
 169   beta__:0,
 170   if mode__#[] then mode__:mode__[1] else mode__:levin_algebraic,
 171   /* remove zeros */
 172   j__:0,
 173   for i__:0 thru terms__-1 do (
 174       if sign(abs(subvar(a__,i__)))#zero then (
 175       b__[j__]:subvar(a__,i__),
 176       j__:j__+1
 177     ) /* else (
 178       if levin_options["debug"] then print("omitting zero term with index ",i)
 179     )*/
 180   ),
 181   last__:j__-1,
 182   if last__<0 then (
 183     if mode__=levin_algebraic then
 184          return(0)            /* best effort */
 185     else
 186          return([bfloat(0),bfloat(0)])
 187   ),
 188   if last__=0 then (
 189      if mode__=levin_algebraic then
 190         return(b__[0])       /* best effort */
 191      else
 192         return([bfloat(b__[0]),bfloat(0)])
 193   ),
 194   if mode__=levin_algebraic then (
 195     omega__[l__]:=l__*b__[l__],
 196     t__[l__]:=t__[l__-1]+b__[l__],
 197     t__[-1]:0,
 198     omega__variance__[l__]:=0,
 199     t__variance__[l__]:=0
 200   ) else (
 201     omega__[l__]:=l__* bfloat(b__[l__]),
 202     t__[l__]:=t__[l__-1]+ bfloat(b__[l__]),
 203     t__[-1]:0,
 204     block([fpprec:16], /* start lower precision block */
 205        omega__variance__[l__]:=l__^2*b__[l__]^2,
 206        t__variance__[l__]:=t__variance__[l__-1]+(abs(bfloat(b__[l__])))^2,
 207        t__variance__[-1]:0
 208     ) /* end lower precision block */
 209   ),
 210   levin(t__,last__-n__,n__,beta__,omega__,mode__,
 211                                t__variance__,omega__variance__)
 212 );
 213
 214 levin_u_complex(a__,terms__,[mode__]):=block([n__,re__,im__,
 215                                         re_result__,im_result__],
 216   local(re__,im__),
 217   if mode__#[] then mode__:mode__[1] else mode__:levin_algebraic,
 218   if mode__=levin_algebraic then (
 219     return(levin_u(a__,terms__,mode__))
 220   ) else (
 221 /*    im__[n__]:=imagpart(arrayapply(a__,[n__])), This give bindstack overflow.
 222                                                   Why???
 223     re__[n__]:=realpart(arrayapply(a__,[n__])), */
 224     for n__:0 thru terms__-1 do (
 225         im__[n__]:imagpart(arrayapply(a__,[n__])),
 226         re__[n__]:realpart(arrayapply(a__,[n__]))
 227     ),
 228     re_result__:levin_u(re__,terms__,mode__),
 229     im_result__:levin_u(im__,terms__,mode__),
 230     return([re_result__[1]+%i*im_result__[1],re_result__[2]+im_result__[2]])
 231   )
 232 );
 233
 234
 235
 236
 237
 238
 239
 240 /* Robert Dodier: */
 241 levin_check_lvalue(e__):=(symbolp(e__) or (subvarp(e__) and levin_check_lvalue(part(e__,0)))) and not stringp(e__);
 242
 243 /* User friendly version.
 244
 245    Example
 246                                3899836039
 247    levin_u_sum(1/n^2,n,1,10) = ---------- = 1.644934081345832
 248                                2370816000
 249                        1968329
 250    sum(1/n^2,n,1,10) = ------- = 1.549767731166541
 251                        1270080
 252
 253    sum(1/n^2,n,1,inf) = %pi^2/6 = 1.644934066848226
 254
 255    If mode=levin_numeric then a floating point error estimate is
 256    returned (the square root of the variance). It is *not* an estimate
 257    for the error of the result.
 258 */
 259
 260
 261 levin_u_sum(a__,index_var__,start__,terms__,[mode__]):=
 262                             block([n__,c__,result__,ratprint:ratprint,
 263                             float2bf:float2bf],
 264   local(c__),
 265   ratprint:false,
 266   float2bf:true,
 267   if mode__#[] then mode__:mode__[1]
 268                               else mode__:levin_algebraic,
 269   if not(integerp(terms__)) or not(constantp(terms__)) or terms__<0 then (
 270      print (
 271 "*** levin_u_sum needs a positive constant integral number of terms ***"
 272      ),
 273      return(levin_u_sum_failed)
 274   ),
 275   if not(levin_check_lvalue(index_var__)) then (
 276      print (
 277 "**** ",index_var," is not an lvalue ***"
 278      ),
 279      return(levin_u_sum_failed)
 280   ),
 281   result__:for n__:0 thru terms__-1 do (
 282              c__[n__]:subst(start__+n__,index_var__,a__),
 283              c__[n__]:ev(c__[n__]),  /* remove one round of quoting */
 284                                     /* this is the plot2d convention */
 285              if mode__=levin_numeric and not(constantp(c__[n__])) then (
 286                      print(
 287 "*** levin_u_sum needs constant terms in numeric mode. ***"
 288                      ),
 289                      return(levin_u_sum_failed)
 290               )
 291    ) ,
 292    if result__=levin_u_sum_failed then return(levin_u_sum_failed),
 293    result__:levin_u_complex(c__,terms__,mode__),
 294    if mode__=levin_algebraic then
 295        return(result__)
 296    else
 297       return([result__[1],bfloat(sqrt(result__[2])*10^(-fpprec))])
 298 );
 299
 300 /* Numerical version. This function makes a crude attempt to compute
 301    an indefinite sum numerically with precision given by fpprec.
 302    If you suspect something is wrong put levin_options["debug"]:true.
 303 */
 304
 305 bflevin_u_sum(a__,index_var__,start__):=block([fpprec:fpprec,fpprec_save__,
 306            precision_overshoot_high__, precision_overshoot_low__,
 307            non_constant__, state__, terms__, new__, error__,
 308            actual_precision__, approx__, approx_old__,
 309            result__,increasing_accuracy__,increasing_terms__,comparing__,
 310            rational__,ratprint:ratprint,float2bf:float2bf,c__],
 311   local(c__,log10__),
 312   log10__(x):=bfloat(log(x)/log(10)),
 313   ratprint:false,
 314   float2bf:true,
 315   if not(levin_check_lvalue(index_var__)) then (
 316      print (
 317 "**** ",index_var__," is not an lvalue ***"
 318      ),
 319      return(bflevin_u_sum_failed)
 320   ),
 321   fpprec_save__:fpprec,
 322   precision_overshoot_high__:fix((3/2)*fpprec_save__), /* I know this is much
 323                                too high in algebraic mode */
 324   precision_overshoot_low__:min(precision_overshoot_high__-1,fpprec_save__+6),
 325   fpprec:min(levin_options["max_precision"],max(levin_options["min_precision"],precision_overshoot_high__)),
 326   terms__:levin_options["min_terms"],
 327   non_constant__:false,
 328   state__:increasing_accuracy__,
 329   result__:do (
 330     rational__:true,
 331     for n__:0 thru terms__-1 do ( /* perhaps we may need to recompute things
 332                               as the user may have used bfloat and the
 333                               precision might have been increased */
 334              c__[n__]:subst(start__+n__,index_var__,a__),
 335              c__[n__]:ev(c__[n__]), /* remove one round of quoting */
 336                                     /* this is the plot2d convention */
 337              if not(constantp(c__[n__])) then (
 338                      print(
 339 "*** bflevin_u_sum does not work  for series with non constant terms. ***"
 340                      ),
 341                    non_constant__:true,
 342                    return(failed)
 343                    ),
 344              if not(ratnump(c__[n__])) then rational__:false
 345     ),
 346     if non_constant__ then return(bflevin_u_sum_failed),
 347     if rational__ then (
 348         approx__:bfloat(levin_u(c__,terms__,levin_algebraic)),
 349         if not constantp(approx__) then return (approx__),
 350         actual_precision__:fpprec
 351     )  else (
 352        new__:levin_u_complex(c__,terms__,levin_numeric),
 353        if not constantp(new__) then return (new__),
 354        approx__:new__[1],
 355        error__:sqrt(new__[2]),
 356        if equal(approx__,bfloat(0)) then actual_precision__:fpprec,
 357                /* assume that a zero is accurate, what else can we do? */
 358        actual_precision__:-fix((log10__(error__/abs(approx__))))+fpprec
 359     ),
 360     if levin_options["debug"] then print (
 361 "==============================","
 362 state=       ",state__,"
 363 precision=   ",fpprec,"
 364 terms=       ",terms__,"
 365 rational=    ",rational__,"
 366 result=      ",approx__,"
 367 fp_accuracy= ",actual_precision__
 368     ),
 369     if state__=increasing_accuracy__ then (
 370          if actual_precision__>=precision_overshoot_low__ then (
 371              state__:comparing__
 372           )
 373            else (
 374              fpprec:fpprec+precision_overshoot_high__-actual_precision__,
 375              if fpprec> levin_options["max_precision"] then (
 376                print (
 377  "*** required precision too high ***"
 378                ),
 379                return(bflevin_u_sum_failed)
 380              )
 381           )
 382     ) else if state__=increasing_terms__ then  (
 383           approx_old__:approx__,
 384           if 2*terms__>levin_options["max_terms"] then (
 385              print (
 386 "*** bflevin_u_sum did not obtain the required precision using ",terms__," terms ***"
 387              ),
 388              return(bflevin_u_sum_failed)
 389           ),
 390           terms__:2*terms__,
 391           state__:increasing_accuracy__
 392     ) else if state__=comparing__ then (
 393        if constantp(approx_old__) then (
 394           if equal(approx_old__,bfloat(0)) and equal(approx__,bfloat(0))
 395                              then return(bfloat(0)),
 396           if levin_options["debug"] then
 397           print("difference=  ",approx_old__-approx__),
 398           if abs(approx_old__-approx__)/
 399                          max(abs(approx_old__),abs(approx__))
 400                   <10^(-fpprec_save__-1) then return(approx__)
 401
 402        ),
 403        state__:increasing_terms__
 404     )
 405   ),
 406   fpprec:fpprec_save__,
 407   bfloat(result__)
 408 );
 409
 410
 411
 412