share/numeric/rtest_bffac.mac

   1 (kill (all),
   2  if not ?fboundp ('bffac) then load ("bffac"),
   3  float_approx_equal_tolerance : 1e-13,
   4  0);
   5 0;
   6
   7 (float_near_zero (x) := is (abs(x) < float_approx_equal_tolerance),
   8  0);
   9 0$
  10
  11 /* Easy standard values of zeta */
  12 bfzeta (0, 10);
  13 ''(bfloat(-1/2))$
  14
  15 bfzeta (1, 10);
  16 inf$
  17
  18 float_near_zero (bfzeta(-2, 10));
  19 true$
  20
  21 /* zeta(2) = %pi^2/6 (A&S 23.2.24) */
  22 float_near_zero (bfzeta(2, 14) - %pi^2/6);
  23 true$
  24
  25 /* zeta(4) = %pi^4/90 (A&S 23.2.25) */
  26 float_near_zero (bfzeta(4, 15) - %pi^4/90);
  27 true$
  28
  29 /* Don't die horribly with arguments we can't evaluate zeta for */
  30 bfzeta (x, 10);
  31 'bfzeta (x, 10)$
  32
  33 bfzeta (0, fpp);
  34 'bfzeta (0, fpp)$
  35
  36 bfzeta (%i, fpp);
  37 'bfzeta (%i, fpp)$
  38
  39 /* Check the first zero is in the right place */
  40 float_near_zero (bfzeta (1/2+%i*14.13472514173469379045b0, 16));
  41 true$
  42
  43 /* SF bug report #3426: "Algorithm precision problem for small positive values in bfpsi0" */
  44 bfpsi0(1e-10, 1);
  45 -1b10;
  46
  47 bfpsi0(1e-10, 6);
  48 -1b10;
  49
  50 bfpsi0(1e-10, 16);
  51 -10000000000.577215665b0;
  52
  53 bfpsi0(1e-10, 32);
  54 -10000000000.577215665b0;
  55
  56 (reset(float_approx_equal_tolerance),0);
  57 0$