test_sf.c

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

{
  gsl_sf_result r;
  double ga[100];
  int s = 0;
  int sa;

  TEST_SF(s,  gsl_sf_gegenpoly_1_e, (-0.2,   1.0, &r), -0.4, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_1_e, ( 0.0,   1.0, &r), 2.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_1_e, ( 1.0,   1.0, &r), 2.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_1_e, ( 1.0,   0.5, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_1_e, ( 5.0,   1.0, &r), 10.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_1_e, ( 100.0, 0.5, &r), 100.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gegenpoly_2_e, (-0.2,   0.5, &r), 0.12, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_2_e, ( 0.0,   1.0, &r), 1.00, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_2_e, ( 1.0,   1.0, &r), 3.00, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_2_e, ( 1.0,   0.1, &r), -0.96, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_2_e, ( 5.0,   1.0, &r), 55.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_2_e, ( 100.0, 0.5, &r), 4950.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gegenpoly_3_e, (-0.2,   0.5, &r), 0.112, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_3_e, ( 0.0,   1.0, &r), -2.0/3.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_3_e, ( 1.0,   1.0, &r), 4.000, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_3_e, ( 1.0,   0.1, &r), -0.392, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_3_e, ( 5.0,   1.0, &r), 220.000, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_3_e, ( 100.0, 0.5, &r), 161600.000, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (1,       1.0, 1.0, &r), 2.000              , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (10,      1.0, 1.0, &r), 11.000             , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (10,      1.0, 0.1, &r), -0.4542309376      , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (10,      5.0, 1.0, &r), 9.23780e+4         , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (10,    100.0, 0.5, &r), 1.5729338392690000e+13, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (1000,  100.0, 1.0, &r), 3.3353666135627322e+232, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (100,  2000.0, 1.0, &r), 5.8753432034937579e+202, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (103,   207.0, 2.0, &r), 1.4210272202235983e+145, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gegenpoly_n_e, (103,    -0.4, 0.3, &r), -1.64527498094522e-04, TEST_TOL1, GSL_SUCCESS);

  sa = 0;
  gsl_sf_gegenpoly_array(99, 5.0, 1.0, ga);
  sa += ( test_sf_frac_diff( ga[1],     10.0    ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff( ga[10], 9.23780e+4 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_gegenpoly_array");
  s += sa;

  return s;
}


int test_jac(void)
{
  double u, m;
  double sn, cn, dn;
  int stat_ej;
  int s = 0;
  int sa;

  u = 0.5;
  m = 0.5;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn, 0.4707504736556572833, TEST_TOL0, "gsl_sf_elljac_e(0.5|0.5) sn");
  sa += test_sf_val(cn, 0.8822663948904402865, TEST_TOL0, "gsl_sf_elljac_e(0.5|0.5) cn");
  sa += test_sf_val(dn, 0.9429724257773856873, TEST_TOL0, "gsl_sf_elljac_e(0.5|0.5) dn");
  gsl_test(s, "  gsl_sf_elljac_e(0.5|0.5)");
  s += sa;

  u = 1.0;
  m = 0.3;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn, 0.8187707145344889190, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.3) sn");
  sa += test_sf_val(cn, 0.5741206467465548795, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.3) cn");
  sa += test_sf_val(dn, 0.8938033089590823040, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.3) dn");
  gsl_test(sa, "  gsl_sf_elljac_e(1.0|0.3)");
  s += sa;

  u = 1.0;
  m = 0.6;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn, 0.7949388393365780943, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.6) sn");
  sa += test_sf_val(cn, 0.6066895760718277578, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.6) cn");
  sa += test_sf_val(dn, 0.7879361300438814425, TEST_TOL0, "gsl_sf_elljac_e(1.0|0.6) dn");
  gsl_test(sa, "  gsl_sf_elljac_e(1.0|0.6)");
  s += sa;

  u = 3.0;
  m = 0.6;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn,  0.7432676860864044186, TEST_TOL0, " gsl_sf_elljac_e(3.0|0.6) sn");
  sa += test_sf_val(cn, -0.6689941306317733154, TEST_TOL0, " gsl_sf_elljac_e(3.0|0.6) cn");
  sa += test_sf_val(dn,  0.8176379933025723259, TEST_TOL0, " gsl_sf_elljac_e(3.0|0.6) dn");
  gsl_test(sa, "  gsl_sf_elljac_e(3.0|0.6)");
  s += sa;

  u = 2.0;
  m = 0.999999;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn, 0.96402778575700186570, TEST_TOL1, "gsl_sf_elljac_e(2.0|0.999999) sn");
  sa += test_sf_val(cn, 0.26580148285600686381, TEST_TOL1, "gsl_sf_elljac_e(2.0|0.999999) cn");
  sa += test_sf_val(dn, 0.26580323105264131136, TEST_TOL1, "gsl_sf_elljac_e(2.0|0.999999) dn");
  gsl_test(sa, "  gsl_sf_elljac_e(2.0|0.999999)");
  s += sa;

  /* test supplied by Ivan Panchenko */
  u = 1.69695970624443;
  m = 0.270378013104138;
  sa = 0;
  stat_ej = gsl_sf_elljac_e(u, m, &sn, &cn, &dn);
  sa += test_sf_val(sn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(1.69..|0.27..) sn");
  sa += test_sf_val(cn, 0.0, TEST_TOL1, "gsl_sf_elljac_e(1.69..|0.27..) cn");
  sa += test_sf_val(dn, 0.8541791304497336, TEST_TOL1, "gsl_sf_elljac_e(1.69..|0.27..) dn");
  gsl_test(sa, "  gsl_sf_elljac_e(1.69695970624443|0.270378013104138)");
  s += sa;



  /* Check known values from Abramowitz & Stegun, Table 16.5 */
  u = 0;
  m = 0.1;

  {
    double mc = 1 - m;
    /* quarter period K is (pi/2)/agm(1,mc) */
    double K = (M_PI_2)/ 0.9741726903999478375938128316;

    double A = 1.0 / sqrt(1+sqrt(mc));
    double B = pow(mc, 0.25) / sqrt(1+sqrt(mc));
    double C = pow(mc, 0.25);
    double C2 = sqrt(mc);

    double eps = 1e-10;
    
    sa = 0;
    stat_ej = gsl_sf_elljac_e(0.0, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 0.0, TEST_TOL0, "gsl_sf_elljac_e(0|0.1) sn");
    sa += test_sf_val(cn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(0|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(0|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(0|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(-eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, -eps, TEST_TOL0, "gsl_sf_elljac_e(-1e-10|0.1) sn");
    sa += test_sf_val(cn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(-1e-10|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(-1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(-1e-10|0.1)");
    s += sa;
    
    sa = 0;
    stat_ej = gsl_sf_elljac_e(eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, eps, TEST_TOL0, "gsl_sf_elljac_e(1e-10|0.1) sn");
    sa += test_sf_val(cn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(1e-10|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(1e-30, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 1e-30, TEST_TOL0, "gsl_sf_elljac_e(1e-30|0.1) sn");
    sa += test_sf_val(cn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(1e-30|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL0, "gsl_sf_elljac_e(1e-30|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(1e-30|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K / 2.0 - eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, A - eps*B*C, TEST_TOL2, "gsl_sf_elljac_e(K/2-1e-10|0.1) sn");
    sa += test_sf_val(cn, B + eps*A*C, TEST_TOL2, "gsl_sf_elljac_e(K/2-1e-10|0.1) cn");
    sa += test_sf_val(dn, C + m*eps*A*B, TEST_TOL2, "gsl_sf_elljac_e(K/2-1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K/2-1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K / 2.0, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, A, TEST_TOL2, "gsl_sf_elljac_e(K/2|0.1) sn");
    sa += test_sf_val(cn, B, TEST_TOL2, "gsl_sf_elljac_e(K/2|0.1) cn");
    sa += test_sf_val(dn, C, TEST_TOL2, "gsl_sf_elljac_e(K/2|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K/2|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K / 2.0 + eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, A + eps*B*C, TEST_TOL2, "gsl_sf_elljac_e(K/2+1e-10|0.1) sn");
    sa += test_sf_val(cn, B - eps*A*C, TEST_TOL2, "gsl_sf_elljac_e(K/2+1e-10|0.1) cn");
    sa += test_sf_val(dn, C - m*eps*A*B, TEST_TOL2, "gsl_sf_elljac_e(K/2+1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K/2+1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K - eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 1.0, TEST_TOL1, "gsl_sf_elljac_e(K-1e-10|0.1) sn");
    sa += test_sf_val(cn, eps*C2, 10*TEST_SNGL, "gsl_sf_elljac_e(K-1e-10|0.1) cn");
    sa += test_sf_val(dn, C2, TEST_TOL2, "gsl_sf_elljac_e(K-1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K-1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 1.0, TEST_TOL1, "gsl_sf_elljac_e(K|0.1) sn");
    sa += test_sf_val(cn, 0.0, TEST_TOL1, "gsl_sf_elljac_e(K|0.1) cn");
    sa += test_sf_val(dn, C2, TEST_TOL2, "gsl_sf_elljac_e(K|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(K + eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 1.0, TEST_TOL1, "gsl_sf_elljac_e(K+1e-10|0.1) sn");
    sa += test_sf_val(cn, -eps*C2, 10*TEST_SNGL, "gsl_sf_elljac_e(K+1e-10|0.1) cn");
    sa += test_sf_val(dn, C2, TEST_TOL2, "gsl_sf_elljac_e(K+1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(K+1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(3.0*K / 2.0, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, A, TEST_TOL2, "gsl_sf_elljac_e(3K/2|0.1) sn");
    sa += test_sf_val(cn, -B, TEST_TOL2, "gsl_sf_elljac_e(3K/2|0.1) cn");
    sa += test_sf_val(dn, C, TEST_TOL2, "gsl_sf_elljac_e(3K/2|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(3K/2|0.1)");
    s += sa;


    sa = 0;
    stat_ej = gsl_sf_elljac_e(2.0*K - eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, eps, 10*TEST_SNGL, "gsl_sf_elljac_e(2K-1e-10|0.1) sn");
    sa += test_sf_val(cn, -1.0, TEST_TOL1, "gsl_sf_elljac_e(2K-1e-10|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL1, "gsl_sf_elljac_e(2K-1e-10|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(2K-1e-10|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(2.0*K, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, 0.0, TEST_TOL1, "gsl_sf_elljac_e(2K|0.1) sn");
    sa += test_sf_val(cn, -1.0, TEST_TOL1, "gsl_sf_elljac_e(2K|0.1) cn");
    sa += test_sf_val(dn, 1.0, TEST_TOL1, "gsl_sf_elljac_e(2K|0.1) dn");
    gsl_test(sa, "  gsl_sf_elljac_e(2K|0.1)");
    s += sa;

    sa = 0;
    stat_ej = gsl_sf_elljac_e(2.0*K + eps, m, &sn, &cn, &dn);
    sa += test_sf_val(sn, -eps, 10*TEST_SNGL, "gsl_sf_elljac_e(2K+1e-10|0.1) sn");
    sa += test_sf_val(cn, -1.0, TEST_TOL1, "gsl_sf_elljac_e(