share/contrib/levin/levin.dem

   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 fpprec:16$
  19 /* load necessary file */
  20 load("levin.mac")$
  21 ratprint:false$
  22 /* do a random test (Andrej Vodopivec) */
  23 approx:bfloat(bflevin_u_sum(n/((3*n+1)*(2*n+1)^2),n,1));
  24 exact:-(12*sqrt(3)*log(3)-16*sqrt(3)*log(2)-sqrt(3)*%pi^2+4*%pi)/(8*sqrt(3))$
  25 bfloat(approx-exact);
  26 /* increase precision */
  27 fpprec:50$
  28 approx:bfloat(bflevin_u_sum(n/((3*n+1)*(2*n+1)^2),n,1));
  29 bfloat(approx-exact);
  30 /* The sum sum(1/n-log(1+1/n),n,1,inf)  computes the
  31 Euler-Mascheroni constant */
  32 fpprec:16$
  33 levin_gamma:bflevin_u_sum(1/n-log(1+1/n),n,1);
  34 bfloat(levin_gamma-%gamma);
  35 /* a series with zero terms */
  36 pi4(k):=if oddp(k) then (-1)^((k-1)/2)/k else 0$
  37 /* we need to protect pi4(k) against premature evaluation */
  38 bfloat(bflevin_u_sum('(pi4(k)),k,1)-%pi/4);
  39 /* load package to compute zeta function numerically */
  40 load(bffac)$
  41 float2bf:true$
  42 /* approximation to the zeta function */
  43 levin_zeta(s):=levin_u_sum(1/n^s,n,1,5)$
  44 /* Levin can assign a meaningful value to some divergent series.
  45    Explanation : levin_zeta(s) is an approximate analytic continuation of
  46    zeta(s) for Re s>1. Hence it may be expected to be close to zeta(s) for
  47    (some) s with Re s<1 (by definition zeta(s) is analytic on the complex
  48    plane with 1 excluded).
  49 */
  50 plot2d([levin_zeta(s),zeta(s)],[s,-2.0,2.0],[y,-10.0,10.0])$
  51
  52