tests/rtest_elliptic.mac

   1 /* Tests of Jacobi elliptic functions and elliptic integrals */
   2
   3 kill(all);
   4 done$
   5
   6 /* derivatives */
   7 diff(jacobi_sn(u,m),u);
   8 jacobi_cn(u,m)*jacobi_dn(u,m);
   9
  10 diff(jacobi_sn(u,m),m);
  11 jacobi_cn(u,m)*jacobi_dn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))
  12   /(2*m)
  13  +jacobi_cn(u,m)^2*jacobi_sn(u,m)/(2*(1-m));
  14
  15 diff(jacobi_cn(u,m),u);
  16 -jacobi_dn(u,m)*jacobi_sn(u,m);
  17
  18 diff(jacobi_cn(u,m),m);
  19 -(jacobi_dn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/(2*m))
  20   -(jacobi_cn(u,m)*jacobi_sn(u,m)^2/(2*(1-m)));
  21
  22 diff(jacobi_dn(u,m),u);
  23 -m*jacobi_cn(u,m)*jacobi_sn(u,m);
  24
  25 diff(jacobi_dn(u,m),m);
  26 -(jacobi_cn(u,m)*jacobi_sn(u,m)*(u-elliptic_e(asin(jacobi_sn(u,m)),m)/(1-m))/2)
  27   -(jacobi_dn(u,m)*'jacobi_sn(u,m)^2/(2*(1-m)));
  28
  29 diff(inverse_jacobi_sn(u,m),u);
  30 1/(sqrt(1-u^2)*sqrt(1-m*u^2));
  31
  32 diff(inverse_jacobi_sn(u,m),m);
  33 ((elliptic_e(asin(u),m)-(1-m)*elliptic_f(asin(u),m))/m-(u*sqrt(1-u^2)/sqrt(1-m*u^2)))/(1-m);
  34
  35 diff(inverse_jacobi_cn(u,m),u);
  36 -(1/(sqrt(1-u^2)*sqrt(m*u^2-m+1)));
  37
  38 diff(inverse_jacobi_cn(u,m),m);
  39 ((elliptic_e(asin(sqrt(1-u^2)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)),m))/m
  40    -(sqrt(1-u^2)*abs(u)/sqrt(1-m*(1-u^2))))/(1-m);
  41
  42 diff(inverse_jacobi_dn(u,m),u);
  43 1/(sqrt(1-u^2)*sqrt(u^2+m-1));
  44
  45 diff(inverse_jacobi_dn(u,m),m);
  46 ((elliptic_e(asin(sqrt(1-u^2)/sqrt(m)),m)-(1-m)*elliptic_f(asin(sqrt(1-u^2)/sqrt(m)),m))/m
  47   -sqrt(1-(1-u^2)/m)*sqrt(1-u^2)/(sqrt(m)*abs(u)))/(1-m)
  48 -(sqrt(1-u^2)/(2*m^(3/2)*sqrt(1-(1-u^2)/m)*abs(u)));
  49
  50 diff(elliptic_e(phi,m),phi);
  51 sqrt(1-m*sin(phi)^2);
  52
  53 diff(elliptic_e(phi,m),m);
  54 (elliptic_e(phi,m)-elliptic_f(phi,m))/(2*m);
  55
  56 diff(elliptic_f(phi,m),phi);
  57 1/sqrt(1-m*sin(phi)^2);
  58
  59 diff(elliptic_f(phi,m),m);
  60 ((elliptic_e(phi,m)-(1-m)*elliptic_f(phi,m))/m
  61    -(cos(phi)*sin(phi)/sqrt(1-m*sin(phi)^2)))
  62  /(2*(1-m));
  63
  64 diff(elliptic_pi(n,phi,m),n);
  65 ((-(n*sqrt(1-m*sin(phi)^2)*sin(2*phi))/(2*(1-n*sin(phi)^2)))
  66  +((m-n)*elliptic_f(phi,m))/n+elliptic_e(phi,m)
  67  +((n^2-m)*elliptic_pi(n,phi,m))/n)
  68  /(2*(m-n)*(n-1));
  69
  70 diff(elliptic_pi(n,phi,m),phi);
  71 1/(sqrt(1-m*sin(phi)^2)*(1-n*sin(phi)^2));
  72
  73 diff(elliptic_pi(n,phi,m),m);
  74 ((-(m*sin(2*phi))/(2*(m-1)*sqrt(1-m*sin(phi)^2)))
  75  +elliptic_e(phi,m)/(m-1)+elliptic_pi(n,phi,m))
  76  /(2*(n-m));
  77
  78 diff(elliptic_kc(m),m);
  79 (elliptic_ec(m)-(1-m)*elliptic_kc(m))/(2*(1-m)*m);
  80
  81 diff(elliptic_ec(m),m);
  82 (elliptic_ec(m)-elliptic_kc(m))/(2*m);
  83
  84 /* Integrals */
  85
  86 integrate(jacobi_sn(u,m),u); /* A&S 16.24.1 */
  87 log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m);
  88
  89 integrate(jacobi_cn(u,m),u); /* A&S 16.24.2 */
  90 acos(jacobi_dn(u,m))/sqrt(m);
  91
  92 integrate(jacobi_dn(u,m),u); /* A&S 16.24.3 */
  93 asin(jacobi_sn(u,m));
  94
  95 integrate(jacobi_cd(u,m),u); /* A&S 16.24.4 */
  96 log(sqrt(m)*jacobi_sd(u,m)+jacobi_nd(u,m))/sqrt(m);
  97
  98 /* Use functions.wolfram.com 09.35.21.001.01, not A&S 16.24.5 */
  99 integrate(jacobi_sd(u,m),u);
 100 -(sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)*asin(sqrt(m)*jacobi_cd(u,m))
 101   /((1-m)*sqrt(m)));
 102
 103 /* Use functions.wolfram.com 09.32.21.0001.01, not A&S 16.24.6 */
 104 integrate(jacobi_nd(u,m),u);
 105 sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m));
 106
 107 integrate(jacobi_dc(u,m),u); /* A&S 16.24.7 */
 108 log(jacobi_sc(u,m)+jacobi_nc(u,m));
 109
 110 integrate(jacobi_nc(u,m),u); /* A&S 16.24.8 */
 111 log(sqrt(1-m)*jacobi_sc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
 112
 113 integrate(jacobi_sc(u,m),u); /* A&S 16.24.9 */
 114 log(sqrt(1-m)*jacobi_nc(u,m)+jacobi_dc(u,m))/sqrt(1-m);
 115
 116 integrate(jacobi_ns(u,m),u); /* A&S 16.24.10 */
 117 log(jacobi_ds(u,m)-jacobi_cs(u,m));
 118
 119 /* Use functions.wolfram.com 09.30.21.0001.01, not A&S 16.24.11 */
 120 integrate(jacobi_ds(u,m),u);
 121 log((1-jacobi_cn(u,m))/jacobi_sn(u,m));
 122
 123 integrate(jacobi_cs(u,m),u); /* A&S 16.24.12 */
 124 log(jacobi_ns(u,m)-jacobi_ds(u,m));
 125
 126 /* Check the integrals and derivatives by confirming
 127
 128        f(x,m)-diff(integral(x,m),x),x) = constant
 129
 130   Look at Taylor expansion about zero, rather than messing about with
 131   elliptic function.  This was sufficient to find several errors  */
 132 (te(f,n):=ratsimp(
 133            taylor(f(x,m)
 134            -diff(integrate(f(x,m),x),x),x,0,n)),
 135 ti(f,n):=ratsimp(
 136            taylor(integrate(f(x,m),x),x,0,n)
 137             -integrate(taylor(f(x,m),x,0,n-1),x)),
 138 td(f,n):=ratsimp(
 139            taylor(diff(f(x,m),x),x,0,n)
 140             -diff(taylor(f(x,m),x,0,n+1),x)),
 141 /* Compare analytic and numerical integral */
 142 ni(f,x1,x2,m):= block(
 143   [x,n,I,In,Ia,key:1,eps:1.0e-14],
 144   I:integrate(f(x,n),x),
 145   In:quad_qag(f(x,m),x,x1,x2,key),
 146   Ia:subst([x=float(x2),n=float(m)],I)-subst([x=float(x1),n=float(m)],I),
 147   if ( abs(Ia-In[1]) < eps ) then
 148     true
 149   else
 150     [Ia,In[1]]
 151 ),
 152 done);
 153 done;
 154
 155 te(jacobi_sn,4);
 156 0;
 157
 158 te(jacobi_cn,4);
 159 0;
 160
 161 te(jacobi_dn,4);
 162 0;
 163
 164 te(jacobi_cd,4);
 165 0;
 166
 167 assume(m>0,m<1); /* also ok for m<0 and m>0 */
 168 [m > 0, m < 1];
 169 te(jacobi_sd,4);
 170 0;
 171 forget(m>0,m<1);
 172 [m > 0, m < 1];
 173
 174 te(jacobi_nd,4);
 175 0;
 176
 177 te(jacobi_dc,4);
 178 0;
 179
 180 te(jacobi_nc,4);
 181 0;
 182
 183 te(jacobi_sc,4);
 184 0;
 185
 186 /* jacobi_ns, jacobi_ds and jacobi_cs are singular at x=0
 187    Compare numerical and analytic integrals for a single case.
 188 */
 189
 190 ni(jacobi_ns,1,2,1/2);
 191 true;
 192
 193 ni(jacobi_ds,1,2,1/2);
 194 true;
 195
 196 ni(jacobi_cs,1,2,1/2);
 197 true;
 198
 199 kill(te,ti,td,ni);
 200 done;
 201
 202 /* Slightly modified version of test_table taken from rtest_expintegral.mac */
 203
 204 (test_table(func,table,eps) :=
 205 block([badpoints : [],
 206        abserr    : 0,
 207        maxerr    : -1],
 208   for entry in table do
 209     block([z : entry[1],
 210            result, answer],
 211       z : expand(rectform(bfloat(entry[1]))),
 212       result : rectform(apply(func, z)),
 213       answer : expand(rectform(bfloat(entry[2]))),
 214       abserr : abs(result-answer),
 215       maxerr : max(maxerr,abserr),
 216       if abserr > eps then
 217         badpoints : cons ([z,result,answer,abserr],badpoints)
 218     ),
 219   if badpoints # [] then
 220     cons(maxerr,badpoints)
 221   else
 222     badpoints
 223 ),done);
 224 done;
 225
 226 /* These test values come from http://getnet.net/~cherry/testrf.mac */
 227
 228 (mrf:[[[1,2,0],1.3110287771461b0],
 229      [[%i,-%i,0],1.8540746773014b0],
 230      [[0.5,1,0],1.8540746773014b0],
 231      [[%i-1,%i,0],0.79612586584234b0-%i*(1.2138566698365b0)],
 232      [[2,3,4],0.58408284167715b0],
 233      [[%i,-%i,2],1.0441445654064b0],
 234      [[%i-1,%i,1-%i],0.93912050218619b0-%i*(0.53296252018635b0)]],
 235 done);
 236 done;
 237
 238 test_table('carlson_rf, mrf, 1.5b-13);
 239 [];
 240
 241 (mrc:[[[0,1/4],bfloat(%pi)],
 242       [[9/4,2],log(2b0)],
 243       [[0,%i],(1-%i)*(1.1107207345396b0)],
 244       [[-%i,%i],1.2260849569072b0-%i*(0.34471136988768b0)],
 245       [[1/4,-2],log(2b0)/3],
 246       [[%i,-1],0.77778596920447b0+%i*(0.19832484993429b0)],
 247       [[0,1/4],%pi],
 248       [[9/4,2],log(2)],
 249       [[2,1],-log(sqrt(2)-1)],
 250       [[-%i,%i],-log(sqrt(2)-1)/2+ %pi/4-%i*(log(sqrt(2)-1)/2+%pi/4)],
 251       [[1/4,-2],log(2)/3],
 252       [[%i,-1],sqrt(sqrt(2)/4-1/4)*atan(sqrt(sqrt(2)-1))-
 253                sqrt(sqrt(2)/16+1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1)+
 254                %i*(sqrt(sqrt(2)/4+1/4)*atan(sqrt(sqrt(2)-1))+sqrt(sqrt(2)/16-
 255                1/16)*log(-sqrt(2*sqrt(2)+2)+sqrt(2)+1))],
 256       [[0,1],%pi/2],
 257       [[%i,%i+1],%pi/4+%i*log(sqrt(2)-1)/2]],
 258 done);
 259 done;
 260
 261 test_table('carlson_rc, mrc, 2.5b-14);
 262 [];
 263
 264 (mrj:[[[0,1,2,3],0.77688623778582b0],
 265       [[2,3,4,5],0.14297579667157b0],
 266       [[2,3,4,-1+%i],0.13613945827771b0-%i*(0.38207561624427b0)],
 267       [[%i,-%i,0,2],1.6490011662711b0],
 268       [[-1+%i,-1-%i,1,2],0.9414835884122b0],
 269       [[%i,-%i,0,1-%i],1.8260115229009b0+%i*(1.2290661908643b0)],
 270       [[-1+%i,-1-%i,1,-3+%i],-0.61127970812028b0-%i*(1.0684038390007b0)],
 271       [[-1+%i,-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)],
 272       [[2,3,4,-0.5],0.24723819703052b0],
 273       [[2,3,4,-5],-0.12711230042964b0]],
 274 done);
 275 done;
 276
 277 test_table('carlson_rj, mrj, 1b-13);
 278 [];
 279
 280 (mrd:[[[0,2,1],1.7972103521034b0],
 281       [[2,3,4],0.16510527294261b0],
 282       [[%i,-%i,2],0.6593385415422b0],
 283       [[0,%i,-%i],1.270819627191b0+%i*(2.7811120159521b0)],
 284       [[0,%i-1,%i],-1.8577235439239b0-%i*(0.96193450888839b0)],
 285       [[-2-%i,-%i,-1+%i],1.8249027393704b0-%i*(1.2218475784827b0)]],
 286 done);
 287 done;
 288
 289 test_table('carlson_rd, mrd, 6e-14);
 290 [];
 291
 292 /* Some tests of the Jacobian elliptic functions.
 293  * Just some tests at random points, to verify that we are accurate.
 294  * The reference values were obtained from Mathematica, but we could
 295  * also just compute the values using a much larger fpprec.
 296  */
 297
 298 test_table('jacobi_sn,
 299            [[[1b0+%i*1b0, .7b0], 1.134045971912365274954b0 + 0.352252346922494477621b0*%i],
 300             [[1b0+%i*1b0, 2b0], 0.98613318109123804740b0 + 0.09521910819727230780b0*%i],
 301             [[1b0+%i*1b0, 2b0+3b0*%i], 0.94467077879445294981b0 - 0.19586410083100945528b0* %i],
 302             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 303               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 304              1.345233534539675700312281b0 - 7.599023926615176284214056b-2 * %i]],
 305            1b-15);
 306 [];
 307
 308 test_table('jacobi_cn,
 309            [[[1b0+%i*1b0, .7b0], 0.571496591371764254029b0 - 0.698989917271916772991b0*%i ],
 310             [[1b0+%i*1b0, 2b0], 0.33759463268642431412b0 - 0.27814044708010806377b0*%i],
 311             [[1b0+%i*1b0, 2b0+3b0*%i], -0.52142079485827170824b0 - 0.35485177134179601850b0*%i],
 312             [[100b0, .7b0], 0.93004753815774770476196b0]],
 313            6b-14);
 314 [];
 315
 316 test_table('jacobi_dn,
 317            [[[1b0+%i*1b0, .7b0], 0.62297154372331777630564880787568b0 - 0.448863598031509643267241389621738b0 *%i ],
 318             [[1b0+%i*1b0, 2b0], 0.1913322443206041462495606602242b0 - 0.9815253294150083432282549919753b0 * %i],
 319             [[1b0+%i*1b0, 2b0+3b0*%i], 0.6147387452173944656984656771134b0 - 1.4819401302071697495918834416787b0 * %i],
 320             [[100b0, .7b0], 0.95157337933724324055428565654872978b0],
 321             [[1.785063352082689082581887b0 *%i + 8.890858759528509578661528b-1,
 322               9.434463451695984398149033b-1 - 1.476052973708684178844821b-1 * %i],
 323              -8.617730683333292717095686b-1 *%i - 2.663978258141280808361839b-1]],
 324            2b-15);
 325 [];
 326
 327 /* These routines for cn and dn work well for small (<= 1?) values of
 328  * u and m.  They have known issues for large real values of u.
 329  */
 330 (ascending_transform(u,m) :=
 331   block([root_m : expand(rectform(sqrt(m))), mu, root_mu1, v],
 332     mu : expand(rectform(4*root_m/(1+root_m)^2)),
 333     root_mu1 : expand(rectform((1-root_m)/(1+root_m))),
 334     v : expand(rectform(u/(1+root_mu1))),
 335     [v, mu, root_mu1]),
 336  elliptic_dn_ascending(u,m) :=
 337   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 338   else
 339     block([v, mu, root_mu1, new],
 340       [v, mu, root_mu1] : ascending_transform(u,m),
 341       new : elliptic_dn_ascending(v, mu),
 342       expand(rectform((1-root_mu1)/mu*(new^2 + root_mu1)/new))),
 343  elliptic_cn_ascending(u,m) :=
 344   if is(abs(m-1) < 4*10^(-fpprec)) then sech(u)
 345   else
 346     block([v, mu, root_mu1, new],
 347       [v, mu, root_mu1] : ascending_transform(u,m),
 348       new : elliptic_dn_ascending(v, mu),
 349       expand(rectform((1+root_mu1)/mu*(new^2-root_mu1)/new))),
 350  done);
 351 done;
 352
 353 /* Test with random values for the argument and parameter. */
 354 (test_random(n, eps, testf, truef) :=
 355   block([badpoints : [], maxerr : -1],
 356     for k : 1 thru n do
 357       block([z : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 358              m : bfloat(1-2*random(1d0) + %i * (1-2*random(1d0))),
 359              ans, expected, abserr],
 360         ans : testf(z, m),
 361         expected : truef(z, m),
 362         abserr : abs(ans-expected),
 363         maxerr : max(maxerr, abserr),
 364         if abserr > eps then
 365           badpoints : cons([[z,m], ans, expected, abserr], badpoints)),
 366     if badpoints # [] then
 367       cons(maxerr, badpoints)
 368     else
 369       badpoints),
 370  done);
 371 done;
 372
 373 test_random(50, 2b-15, 'jacobi_dn, 'elliptic_dn_ascending);
 374 [];
 375
 376 test_random(50, 2b-15, 'jacobi_cn, 'elliptic_cn_ascending);
 377 [];
 378
 379 /* Test elliptic_f by using the fact that
 380
 381 inverse_jacobi_sn(x,m) = elliptic_f(asin(x), m)
 382
 383 */
 384 (test_ef(x,m) := jacobi_sn(elliptic_f(asin(x),m), m), id(x,m):=x, done);
 385 done;
 386
 387 test_random(100, 6b-13, 'test_ef, 'id);
 388 [];
 389
 390 /* Test elliptic_kc These values are from
 391  *
 392  * http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
 393  */
 394
 395 block([oldfpprec : fpprec, fpprec:100],
 396   test_table('elliptic_kc,
 397              [[[1/2], 8*%pi^(3/2)/gamma(-1/4)^2],
 398               [[17-12*sqrt(2)], 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2],
 399               [[-1], gamma(1/4)^2/4/sqrt(2*%pi)]],
 400              2b-100));
 401 [];
 402
 403 /* Some tests for specific values */
 404 inverse_jacobi_sn(1,m);
 405 elliptic_kc(m);
 406
 407 inverse_jacobi_sn(x,0);
 408 asin(x);
 409
 410 inverse_jacobi_sn(x,1);
 411 log(tan(asin(x)/2 + %pi/4));
 412
 413 inverse_jacobi_sd(1/sqrt(1-m),m);
 414 elliptic_kc(m);
 415
 416 inverse_jacobi_ds(sqrt(1-m),m);
 417 elliptic_kc(m);
 418
 419 /* elliptic_kc(1) is undefined */
 420 errcatch(elliptic_kc(1));
 421 [];
 422 errcatch(elliptic_kc(1.0));
 423 [];
 424 errcatch(elliptic_kc(1.0b0));
 425 [];
 426
 427 /* Test noun/verb results from elliptic functions */
 428 diff(inverse_jacobi_sn(x,0),x);
 429 1/sqrt(1-x^2);
 430
 431 diff(elliptic_pi(4/3,phi,0),phi);
 432 sec(phi)^2/(1-tan(phi)^2/3);
 433
 434 diff(elliptic_pi(3/4,phi,0),phi);
 435 sec(phi)^2/(1+tan(phi)^2/4);
 436
 437 diff(elliptic_pi(1,phi,0),phi);
 438 sec(phi)^2;
 439
 440 /* signbfloat:false provokes "Is 0 positive, negative, or zero?" */
 441 (signbfloat:false,
 442  diff(elliptic_pi(1,phi,0),phi));
 443 sec(phi)^2;
 444
 445 (reset (signbfloat), 0);
 446 0;
 447
 448 /* Special values */
 449 elliptic_ec(1);
 450 1;
 451
 452 elliptic_ec(0);
 453 %pi/2;
 454
 455 elliptic_ec(1/2);
 456 gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2);
 457
 458 elliptic_ec(-1);
 459 sqrt(2)*elliptic_ec(1/2);
 460
 461 /* Test periodicity of elliptic_e */
 462 elliptic_e(x, 1);
 463 2*round(x/%pi)+sin(x);
 464
 465 elliptic_e(3,1/3);
 466 elliptic_e(3-%pi,1/3)+2*elliptic_ec(1/3);
 467
 468 /* Bug #2629: elliptic_kc(3.0) not accurate */
 469 test_table('elliptic_kc, [[[3.0], elliptic_kc(3b0)]], 1b-15);
 470 [];
 471
 472 /* Bug #2630: inverse_jacobi_cn(-2.0, 3.0) generates an error */
 473 test_table('inverse_jacobi_cn, [[[-2.0, 3.0], 2.002154760912212-3.202503914656527*%i]],
 474   1b-15);
 475 [];
 476
 477 /* Bug #2615: Numeric evaluation of inverse Jacobi elliptic functions is wrong for some inputs
 478 */
 479 is(abs(jacobi_dn(inverse_jacobi_dn(-2.0,3.0), 3.0) + 2) < 1d-14);
 480 true;
 481
 482 /* elliptical functions handling of non-rectangular complex numbers
 483  * mailing list 2016-07-14 "Jacobi elliptic functions, maxima 5.38.1"
 484  * mailing list 2020-05-28 "Unexpected result from ev"
 485  */
 486
 487 /* list of functions here is everything that has a SIMP-[$%](ELLIPTIC|JACOBI)_FOO in sr/ellipt.lisp */
 488 /*
 489 elliptic_fcns :
 490   [ elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi,
 491     inverse_jacobi_cd, inverse_jacobi_cn, inverse_jacobi_cs, inverse_jacobi_dc,
 492     inverse_jacobi_dn, inverse_jacobi_ds, inverse_jacobi_nc, inverse_jacobi_nd,
 493     inverse_jacobi_ns, inverse_jacobi_sc, inverse_jacobi_sd, inverse_jacobi_sn,
 494     jacobi_am,
 495     jacobi_cd, jacobi_cn, jacobi_cs, jacobi_dc, jacobi_dn, jacobi_ds,
 496     jacobi_nc, jacobi_nd, jacobi_ns, jacobi_sc, jacobi_sd, jacobi_sn ];
 497  */
 498
 499 /* code to generate test cases below
 500
 501 elliptic_fcns_1 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a)) # []));
 502 elliptic_fcns_2 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b)) # []));
 503 elliptic_fcns_3 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b, 'c)) # []));
 504
 505 kill (u1, k1, n1);
 506
 507 with_stdout ("/tmp/foo.mac",
 508
 509   for f in elliptic_fcns_1
 510     do block ([a : f('rectform(u1))],
 511         print ('complex_floatp(a), ";"),
 512         print (true, ";"),
 513         print ("")),
 514
 515   for f in elliptic_fcns_2
 516     do block ([a : f('rectform(u1), k1)],
 517         print ('complex_floatp(a), ";"),
 518         print (true, ";"),
 519         print ("")),
 520
 521   for f in elliptic_fcns_3
 522     do block ([a : f(n1, 'rectform(u1), k1)],
 523         print ('complex_floatp(a), ";"),
 524         print (true, ";"),
 525         print ("")));
 526  */
 527
 528 (u1 : 0.5*(1.0 - 0.25*%i)^2,
 529  k1 : 0.775,
 530  n1 : -1.375,
 531  complex_floatp(e) := floatnump(realpart(e)) and floatnump(imagpart(e)),
 532  0);
 533 0;
 534
 535 complex_floatp(elliptic_ec(rectform(u1))) ;
 536 true ;
 537
 538 complex_floatp(elliptic_kc(rectform(u1))) ;
 539 true ;
 540
 541 complex_floatp(elliptic_e(rectform(u1), k1)) ;
 542 true ;
 543
 544 complex_floatp(elliptic_eu(rectform(u1), k1)) ;
 545 true ;
 546
 547 complex_floatp(elliptic_f(rectform(u1), k1)) ;
 548 true ;
 549
 550 complex_floatp(inverse_jacobi_cd(rectform(u1), k1)) ;
 551 true ;
 552
 553 complex_floatp(inverse_jacobi_cn(rectform(u1), k1)) ;
 554 true ;
 555
 556 complex_floatp(inverse_jacobi_cs(rectform(u1), k1)) ;
 557 true ;
 558
 559 complex_floatp(inverse_jacobi_dc(rectform(u1), k1)) ;
 560 true ;
 561
 562 complex_floatp(inverse_jacobi_dn(rectform(u1), k1)) ;
 563 true ;
 564
 565 complex_floatp(inverse_jacobi_ds(rectform(u1), k1)) ;
 566 true ;
 567
 568 complex_floatp(inverse_jacobi_nc(rectform(u1), k1)) ;
 569 true ;
 570
 571 complex_floatp(inverse_jacobi_nd(rectform(u1), k1)) ;
 572 true ;
 573
 574 complex_floatp(inverse_jacobi_ns(rectform(u1), k1)) ;
 575 true ;
 576
 577 complex_floatp(inverse_jacobi_sc(rectform(u1), k1)) ;
 578 true ;
 579
 580 complex_floatp(inverse_jacobi_sd(rectform(u1), k1)) ;
 581 true ;
 582
 583 complex_floatp(inverse_jacobi_sn(rectform(u1), k1)) ;
 584 true ;
 585
 586 complex_floatp(jacobi_am(rectform(u1), k1)) ;
 587 true ;
 588
 589 complex_floatp(jacobi_cd(rectform(u1), k1)) ;
 590 true ;
 591
 592 complex_floatp(jacobi_cn(rectform(u1), k1)) ;
 593 true ;
 594
 595 complex_floatp(jacobi_cs(rectform(u1), k1)) ;
 596 true ;
 597
 598 complex_floatp(jacobi_dc(rectform(u1), k1)) ;
 599 true ;
 600
 601 complex_floatp(jacobi_dn(rectform(u1), k1)) ;
 602 true ;
 603
 604 complex_floatp(jacobi_ds(rectform(u1), k1)) ;
 605 true ;
 606
 607 complex_floatp(jacobi_nc(rectform(u1), k1)) ;
 608 true ;
 609
 610 complex_floatp(jacobi_nd(rectform(u1), k1)) ;
 611 true ;
 612
 613 complex_floatp(jacobi_ns(rectform(u1), k1)) ;
 614 true ;
 615
 616 complex_floatp(jacobi_sc(rectform(u1), k1)) ;
 617 true ;
 618
 619 complex_floatp(jacobi_sd(rectform(u1), k1)) ;
 620 true ;
 621
 622 complex_floatp(jacobi_sn(rectform(u1), k1)) ;
 623 true ;
 624
 625 complex_floatp(elliptic_pi(n1, rectform(u1), k1)) ;
 626 true ;
 627
 628 /*
 629 with_stdout ("/tmp/bar.mac",
 630
 631   for f in elliptic_fcns_1
 632     do block ([a : f('rectform(u1))],
 633         print ('complex_bfloatp(a), ";"),
 634         print (true, ";"),
 635         print ("")),
 636
 637   for f in elliptic_fcns_2
 638     do block ([a : f('rectform(u1), k1)],
 639         print ('complex_bfloatp(a), ";"),
 640         print (true, ";"),
 641         print ("")),
 642
 643   for f in elliptic_fcns_3
 644     do block ([a : f(n1, 'rectform(u1), k1)],
 645         print ('complex_bfloatp(a), ";"),
 646         print (true, ";"),
 647         print ("")));
 648  */
 649
 650 (u1 : 0.5b0*(1.0b0 - 0.25b0*%i)^2,
 651  k1 : 0.775b0,
 652  n1 : -1.375b0,
 653  complex_bfloatp(e) := bfloatp(realpart(e)) and bfloatp(imagpart(e)),
 654  0);
 655 0;
 656
 657 complex_bfloatp(elliptic_ec(rectform(u1))) ;
 658 true ;
 659
 660 complex_bfloatp(elliptic_kc(rectform(u1))) ;
 661 true ;
 662
 663 complex_bfloatp(elliptic_e(rectform(u1), k1)) ;
 664 true ;
 665
 666 complex_bfloatp(elliptic_eu(rectform(u1), k1)) ;
 667 true ;
 668
 669 complex_bfloatp(elliptic_f(rectform(u1), k1)) ;
 670 true ;
 671
 672 complex_bfloatp(inverse_jacobi_cd(rectform(u1), k1)) ;
 673 true ;
 674
 675 complex_bfloatp(inverse_jacobi_cn(rectform(u1), k1)) ;
 676 true ;
 677
 678 complex_bfloatp(inverse_jacobi_cs(rectform(u1), k1)) ;
 679 true ;
 680
 681 complex_bfloatp(inverse_jacobi_dc(rectform(u1), k1)) ;
 682 true ;
 683
 684 complex_bfloatp(inverse_jacobi_dn(rectform(u1), k1)) ;
 685 true ;
 686
 687 complex_bfloatp(inverse_jacobi_ds(rectform(u1), k1)) ;
 688 true ;
 689
 690 complex_bfloatp(inverse_jacobi_nc(rectform(u1), k1)) ;
 691 true ;
 692
 693 complex_bfloatp(inverse_jacobi_nd(rectform(u1), k1)) ;
 694 true ;
 695
 696 complex_bfloatp(inverse_jacobi_ns(rectform(u1), k1)) ;
 697 true ;
 698
 699 complex_bfloatp(inverse_jacobi_sc(rectform(u1), k1)) ;
 700 true ;
 701
 702 complex_bfloatp(inverse_jacobi_sd(rectform(u1), k1)) ;
 703 true ;
 704
 705 complex_bfloatp(inverse_jacobi_sn(rectform(u1), k1)) ;
 706 true ;
 707
 708 complex_bfloatp(jacobi_am(rectform(u1), k1)) ;
 709 true ;
 710
 711 complex_bfloatp(jacobi_cd(rectform(u1), k1)) ;
 712 true ;
 713
 714 complex_bfloatp(jacobi_cn(rectform(u1), k1)) ;
 715 true ;
 716
 717 complex_bfloatp(jacobi_cs(rectform(u1), k1)) ;
 718 true ;
 719
 720 complex_bfloatp(jacobi_dc(rectform(u1), k1)) ;
 721 true ;
 722
 723 complex_bfloatp(jacobi_dn(rectform(u1), k1)) ;
 724 true ;
 725
 726 complex_bfloatp(jacobi_ds(rectform(u1), k1)) ;
 727 true ;
 728
 729 complex_bfloatp(jacobi_nc(rectform(u1), k1)) ;
 730 true ;
 731
 732 complex_bfloatp(jacobi_nd(rectform(u1), k1)) ;
 733 true ;
 734
 735 complex_bfloatp(jacobi_ns(rectform(u1), k1)) ;
 736 true ;
 737
 738 complex_bfloatp(jacobi_sc(rectform(u1), k1)) ;
 739 true ;
 740
 741 complex_bfloatp(jacobi_sd(rectform(u1), k1)) ;
 742 true ;
 743
 744 complex_bfloatp(jacobi_sn(rectform(u1), k1)) ;
 745 true ;
 746
 747 complex_bfloatp(elliptic_pi(n1, rectform(u1), k1)) ;
 748 true ;
 749
 750 /* Test derivatives at a few random points using numerical differentiation
 751
 752    func  - a function
 753    df - deriviative of func wrt p-th argument
 754    p - derivative wrt p-th argument
 755    vars - variables in expression deriv
 756    table - table of points to evaluate
 757    eps - report errors > eps
 758    delta - offset for numerical derivative
 759 */
 760 (test_deriv(func,df,vars,p,table,eps,delta) :=
 761 block([badpoints : [],
 762        abserr    : 0,
 763        maxerr    : -1],
 764   for %z in table do
 765     block([%z0, %z1, deriv, nderiv],
 766       %z1:makelist(if i=p then %z[i]+delta else %z[i],i,1,length(%z)),
 767       %z0:makelist(if i=p then %z[i]-delta else %z[i],i,1,length(%z)),
 768       nderiv: float((apply(func,%z1)-apply(func,%z0))/(2*delta)),
 769       deriv : float(subst(maplist("=",vars,%z),df)),
 770       abserr : abs(nderiv-deriv),
 771       maxerr : max(maxerr,abserr),
 772       if abserr > eps then
 773         badpoints : cons ([%z,deriv,nderiv,abserr],badpoints)
 774     ),
 775   if badpoints # [] then
 776     cons(maxerr,badpoints)
 777   else
 778     badpoints
 779     ),done);
 780 done;
 781
 782 /* test points for two-arg functions */
 783 (l2: [[0.2, 0.3], [0.2, 0.5], [1.5, 0.7], [0.5, 0.99], [1.5, 0.8],
 784     [1.57, 0.8], [1.5707, 0.8], [-1.0, 0.8], [2.0, 0.8]],done);
 785 done;
 786
 787 /* test points for three-arg functions */
 788 /* bug #3221 elliptic_pi() wrong for 2nd arg > float(%pi/2) */
 789 (l3: [[0.1,0.2,0.3],[0.1,0.2,0.5],[0.2,1.5,0.7],[0.001,0.5,0.99],
 790     [0.2,1.5,0.8],[0.2,1.57,0.8],[0.2,1.5707,0.8]],done);
 791 done;
 792
 793 test_deriv('elliptic_f,diff(elliptic_f(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
 794 [];
 795
 796 test_deriv('elliptic_f,diff(elliptic_f(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
 797 [];
 798
 799 test_deriv('elliptic_e,diff(elliptic_e(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
 800 [];
 801
 802 test_deriv('elliptic_e,diff(elliptic_e(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
 803 [];
 804
 805 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),n),[n,z,m],1,l3,2.0e-7,1.0e-8);
 806 [];
 807
 808 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),z),[n,z,m],2,l3,2.0e-7,1.0e-8);
 809 [];
 810
 811 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),m),[n,z,m],3,l3,2.0e-7,1.0e-8);
 812 [];
 813
 814 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 815 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
 816
 817 /*
 818  * Some additional tests for complex values. Use the fact that
 819  * elliptic_pi(0, phi, m) = elliptic_f(phi,m);
 820  */
 821 closeto(elliptic_pi(0, 0.5, 0.25)- elliptic_f(0.5, 0.25), 1e-15);
 822 true;
 823 closeto(elliptic_pi(0, 0.5*%i, 0.25)- elliptic_f(0.5*%i, 0.25), 1e-15);
 824 true;
 825
 826 closeto(elliptic_pi(0, 0.5b0, 0.25b0)- elliptic_f(0.5b0, 0.25b0), 1e-15);
 827 true;
 828 closeto(elliptic_pi(0, 0.5b0*%i, 0.25b0)- elliptic_f(0.5b0*%i, 0.25b0), 1e-15);
 829 true;
 830
 831 /* Test for Bug 3221 */
 832
 833 closeto(elliptic_pi(0.2,1.57,0.8) - 2.5770799919605668, 1e-14);
 834 true;
 835 closeto(elliptic_pi(0.2,1.58,0.8) - 2.60502920656026151, 1e-14);
 836 true;
 837 closeto(elliptic_pi(0.20,2.00,0.80) - 3.6543835090829025, 7.47e-11);
 838 true;
 839
 840 (oldfpprec:fpprec,fpprec:32);
 841 32;
 842 closeto(elliptic_pi(0.2b0,1.57b0,0.8b0) - 2.5770799919605668058849721013196b0, 3.09b-32);
 843 true;
 844 closeto(elliptic_pi(0.2b0,1.58b0,0.8b0) - 2.6050292065602615145481924917132b0, 1.24b-32);
 845 true;
 846 closeto(elliptic_pi(0.2b0,2.0b0,0.8b0) - 3.6543835090082902501670624830938b0, 1b-32);
 847 true;
 848
 849 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
 850 done;