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 /* verify that elliptical functions handle non-rectangular complex numbers
 483  * mailing list 2016-07-14 "Jacobi elliptic functions, maxima 5.38.1"
 484  */
 485
 486 /* list of functions here is everything that has a SIMP-[$%](ELLIPTIC|JACOBI)_FOO in sr/ellipt.lisp */
 487 /*
 488 elliptic_fcns :
 489   [ elliptic_e, elliptic_ec, elliptic_eu, elliptic_f, elliptic_kc, elliptic_pi,
 490     inverse_jacobi_cd, inverse_jacobi_cn, inverse_jacobi_cs, inverse_jacobi_dc,
 491     inverse_jacobi_dn, inverse_jacobi_ds, inverse_jacobi_nc, inverse_jacobi_nd,
 492     inverse_jacobi_ns, inverse_jacobi_sc, inverse_jacobi_sd, inverse_jacobi_sn,
 493     jacobi_am,
 494     jacobi_cd, jacobi_cn, jacobi_cs, jacobi_dc, jacobi_dn, jacobi_ds,
 495     jacobi_nc, jacobi_nd, jacobi_ns, jacobi_sc, jacobi_sd, jacobi_sn ];
 496  */
 497
 498 /* code to generate test cases below
 499
 500 elliptic_fcns_1 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a)) # []));
 501 elliptic_fcns_2 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b)) # []));
 502 elliptic_fcns_3 : sublist (elliptic_fcns, lambda ([f], errcatch(f('a, 'b, 'c)) # []));
 503
 504 kill (u1, k1, n1);
 505
 506 with_stdout ("/tmp/foo.mac",
 507
 508   for f in elliptic_fcns_1
 509     do block ([a : f(u1), b : f('rectform(u1))],
 510         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
 511         print (true, ";"),
 512         print ("")),
 513
 514   for f in elliptic_fcns_2
 515     do block ([a : f(u1, k1), b : f('rectform(u1), k1)],
 516         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
 517         print (true, ";"),
 518         print ("")),
 519
 520   for f in elliptic_fcns_3
 521     do block ([a : f(n1, u1, k1), b : f(n1, 'rectform(u1), k1)],
 522         print ('if 'complex_floatp(a) and 'complex_floatp(b) and abs(a - b) < 1e-12 then true else a # b, ";"),
 523         print (true, ";"),
 524         print ("")));
 525  */
 526
 527 (u1 : 0.5*(1.0 - 0.25*%i)^2,
 528  k1 : 0.775,
 529  n1 : -1.375,
 530  complex_floatp(e) := floatnump(realpart(e)) and floatnump(imagpart(e)),
 531  0);
 532 0;
 533
 534 if complex_floatp(elliptic_ec(u1))
 535  and complex_floatp(elliptic_ec(rectform(u1)))
 536  and (abs(elliptic_ec(rectform(u1)) - elliptic_ec(u1)) < 1.0E-12) then true
 537  else elliptic_ec(u1) # elliptic_ec(rectform(u1)) ;
 538 true ;
 539
 540 if complex_floatp(elliptic_kc(u1))
 541  and complex_floatp(elliptic_kc(rectform(u1)))
 542  and (abs(elliptic_kc(rectform(u1)) - elliptic_kc(u1)) < 1.0E-12) then true
 543  else elliptic_kc(u1) # elliptic_kc(rectform(u1)) ;
 544 true ;
 545
 546 if complex_floatp(elliptic_e(u1, k1))
 547  and complex_floatp(elliptic_e(rectform(u1), k1))
 548  and (abs(elliptic_e(rectform(u1), k1) - elliptic_e(u1, k1)) < 1.0E-12)
 549  then true else elliptic_e(u1, k1) # elliptic_e(rectform(u1), k1) ;
 550 true ;
 551
 552 if complex_floatp(elliptic_eu(u1, k1))
 553  and complex_floatp(elliptic_eu(rectform(u1), k1))
 554  and (abs(elliptic_eu(rectform(u1), k1) - elliptic_eu(u1, k1)) < 1.0E-12)
 555  then true else elliptic_eu(u1, k1) # elliptic_eu(rectform(u1), k1) ;
 556 true ;
 557
 558 if complex_floatp(elliptic_f(u1, k1))
 559  and complex_floatp(elliptic_f(rectform(u1), k1))
 560  and (abs(elliptic_f(rectform(u1), k1) - elliptic_f(u1, k1)) < 1.0E-12)
 561  then true else elliptic_f(u1, k1) # elliptic_f(rectform(u1), k1) ;
 562 true ;
 563
 564 if complex_floatp(inverse_jacobi_cd(u1, k1))
 565  and complex_floatp(inverse_jacobi_cd(rectform(u1), k1))
 566  and (abs(inverse_jacobi_cd(rectform(u1), k1) - inverse_jacobi_cd(u1, k1)) < 1.0E-12)
 567  then true else inverse_jacobi_cd(u1, k1) # inverse_jacobi_cd(rectform(u1), k1)
 568  ;
 569 true ;
 570
 571 if complex_floatp(inverse_jacobi_cn(u1, k1))
 572  and complex_floatp(inverse_jacobi_cn(rectform(u1), k1))
 573  and (abs(inverse_jacobi_cn(rectform(u1), k1) - inverse_jacobi_cn(u1, k1)) < 1.0E-12)
 574  then true else inverse_jacobi_cn(u1, k1) # inverse_jacobi_cn(rectform(u1), k1)
 575  ;
 576 true ;
 577
 578 if complex_floatp(inverse_jacobi_cs(u1, k1))
 579  and complex_floatp(inverse_jacobi_cs(rectform(u1), k1))
 580  and (abs(inverse_jacobi_cs(rectform(u1), k1) - inverse_jacobi_cs(u1, k1)) < 1.0E-12)
 581  then true else inverse_jacobi_cs(u1, k1) # inverse_jacobi_cs(rectform(u1), k1)
 582  ;
 583 true ;
 584
 585 if complex_floatp(inverse_jacobi_dc(u1, k1))
 586  and complex_floatp(inverse_jacobi_dc(rectform(u1), k1))
 587  and (abs(inverse_jacobi_dc(rectform(u1), k1) - inverse_jacobi_dc(u1, k1)) < 1.0E-12)
 588  then true else inverse_jacobi_dc(u1, k1) # inverse_jacobi_dc(rectform(u1), k1)
 589  ;
 590 true ;
 591
 592 if complex_floatp(inverse_jacobi_dn(u1, k1))
 593  and complex_floatp(inverse_jacobi_dn(rectform(u1), k1))
 594  and (abs(inverse_jacobi_dn(rectform(u1), k1) - inverse_jacobi_dn(u1, k1)) < 1.0E-12)
 595  then true else inverse_jacobi_dn(u1, k1) # inverse_jacobi_dn(rectform(u1), k1)
 596  ;
 597 true ;
 598
 599 if complex_floatp(inverse_jacobi_ds(u1, k1))
 600  and complex_floatp(inverse_jacobi_ds(rectform(u1), k1))
 601  and (abs(inverse_jacobi_ds(rectform(u1), k1) - inverse_jacobi_ds(u1, k1)) < 1.0E-12)
 602  then true else inverse_jacobi_ds(u1, k1) # inverse_jacobi_ds(rectform(u1), k1)
 603  ;
 604 true ;
 605
 606 if complex_floatp(inverse_jacobi_nc(u1, k1))
 607  and complex_floatp(inverse_jacobi_nc(rectform(u1), k1))
 608  and (abs(inverse_jacobi_nc(rectform(u1), k1) - inverse_jacobi_nc(u1, k1)) < 1.0E-12)
 609  then true else inverse_jacobi_nc(u1, k1) # inverse_jacobi_nc(rectform(u1), k1)
 610  ;
 611 true ;
 612
 613 if complex_floatp(inverse_jacobi_nd(u1, k1))
 614  and complex_floatp(inverse_jacobi_nd(rectform(u1), k1))
 615  and (abs(inverse_jacobi_nd(rectform(u1), k1) - inverse_jacobi_nd(u1, k1)) < 1.0E-12)
 616  then true else inverse_jacobi_nd(u1, k1) # inverse_jacobi_nd(rectform(u1), k1)
 617  ;
 618 true ;
 619
 620 if complex_floatp(inverse_jacobi_ns(u1, k1))
 621  and complex_floatp(inverse_jacobi_ns(rectform(u1), k1))
 622  and (abs(inverse_jacobi_ns(rectform(u1), k1) - inverse_jacobi_ns(u1, k1)) < 1.0E-12)
 623  then true else inverse_jacobi_ns(u1, k1) # inverse_jacobi_ns(rectform(u1), k1)
 624  ;
 625 true ;
 626
 627 if complex_floatp(inverse_jacobi_sc(u1, k1))
 628  and complex_floatp(inverse_jacobi_sc(rectform(u1), k1))
 629  and (abs(inverse_jacobi_sc(rectform(u1), k1) - inverse_jacobi_sc(u1, k1)) < 1.0E-12)
 630  then true else inverse_jacobi_sc(u1, k1) # inverse_jacobi_sc(rectform(u1), k1)
 631  ;
 632 true ;
 633
 634 if complex_floatp(inverse_jacobi_sd(u1, k1))
 635  and complex_floatp(inverse_jacobi_sd(rectform(u1), k1))
 636  and (abs(inverse_jacobi_sd(rectform(u1), k1) - inverse_jacobi_sd(u1, k1)) < 1.0E-12)
 637  then true else inverse_jacobi_sd(u1, k1) # inverse_jacobi_sd(rectform(u1), k1)
 638  ;
 639 true ;
 640
 641 if complex_floatp(inverse_jacobi_sn(u1, k1))
 642  and complex_floatp(inverse_jacobi_sn(rectform(u1), k1))
 643  and (abs(inverse_jacobi_sn(rectform(u1), k1) - inverse_jacobi_sn(u1, k1)) < 1.0E-12)
 644  then true else inverse_jacobi_sn(u1, k1) # inverse_jacobi_sn(rectform(u1), k1)
 645  ;
 646 true ;
 647
 648 if complex_floatp(jacobi_am(u1, k1))
 649  and complex_floatp(jacobi_am(rectform(u1), k1))
 650  and (abs(jacobi_am(rectform(u1), k1) - jacobi_am(u1, k1)) < 1.0E-12) then true
 651  else jacobi_am(u1, k1) # jacobi_am(rectform(u1), k1) ;
 652 true ;
 653
 654 if complex_floatp(jacobi_cd(u1, k1))
 655  and complex_floatp(jacobi_cd(rectform(u1), k1))
 656  and (abs(jacobi_cd(rectform(u1), k1) - jacobi_cd(u1, k1)) < 1.0E-12) then true
 657  else jacobi_cd(u1, k1) # jacobi_cd(rectform(u1), k1) ;
 658 true ;
 659
 660 if complex_floatp(jacobi_cn(u1, k1))
 661  and complex_floatp(jacobi_cn(rectform(u1), k1))
 662  and (abs(jacobi_cn(rectform(u1), k1) - jacobi_cn(u1, k1)) < 1.0E-12) then true
 663  else jacobi_cn(u1, k1) # jacobi_cn(rectform(u1), k1) ;
 664 true ;
 665
 666 if complex_floatp(jacobi_cs(u1, k1))
 667  and complex_floatp(jacobi_cs(rectform(u1), k1))
 668  and (abs(jacobi_cs(rectform(u1), k1) - jacobi_cs(u1, k1)) < 1.0E-12) then true
 669  else jacobi_cs(u1, k1) # jacobi_cs(rectform(u1), k1) ;
 670 true ;
 671
 672 if complex_floatp(jacobi_dc(u1, k1))
 673  and complex_floatp(jacobi_dc(rectform(u1), k1))
 674  and (abs(jacobi_dc(rectform(u1), k1) - jacobi_dc(u1, k1)) < 1.0E-12) then true
 675  else jacobi_dc(u1, k1) # jacobi_dc(rectform(u1), k1) ;
 676 true ;
 677
 678 if complex_floatp(jacobi_dn(u1, k1))
 679  and complex_floatp(jacobi_dn(rectform(u1), k1))
 680  and (abs(jacobi_dn(rectform(u1), k1) - jacobi_dn(u1, k1)) < 1.0E-12) then true
 681  else jacobi_dn(u1, k1) # jacobi_dn(rectform(u1), k1) ;
 682 true ;
 683
 684 if complex_floatp(jacobi_ds(u1, k1))
 685  and complex_floatp(jacobi_ds(rectform(u1), k1))
 686  and (abs(jacobi_ds(rectform(u1), k1) - jacobi_ds(u1, k1)) < 1.0E-12) then true
 687  else jacobi_ds(u1, k1) # jacobi_ds(rectform(u1), k1) ;
 688 true ;
 689
 690 if complex_floatp(jacobi_nc(u1, k1))
 691  and complex_floatp(jacobi_nc(rectform(u1), k1))
 692  and (abs(jacobi_nc(rectform(u1), k1) - jacobi_nc(u1, k1)) < 1.0E-12) then true
 693  else jacobi_nc(u1, k1) # jacobi_nc(rectform(u1), k1) ;
 694 true ;
 695
 696 if complex_floatp(jacobi_nd(u1, k1))
 697  and complex_floatp(jacobi_nd(rectform(u1), k1))
 698  and (abs(jacobi_nd(rectform(u1), k1) - jacobi_nd(u1, k1)) < 1.0E-12) then true
 699  else jacobi_nd(u1, k1) # jacobi_nd(rectform(u1), k1) ;
 700 true ;
 701
 702 if complex_floatp(jacobi_ns(u1, k1))
 703  and complex_floatp(jacobi_ns(rectform(u1), k1))
 704  and (abs(jacobi_ns(rectform(u1), k1) - jacobi_ns(u1, k1)) < 1.0E-12) then true
 705  else jacobi_ns(u1, k1) # jacobi_ns(rectform(u1), k1) ;
 706 true ;
 707
 708 if complex_floatp(jacobi_sc(u1, k1))
 709  and complex_floatp(jacobi_sc(rectform(u1), k1))
 710  and (abs(jacobi_sc(rectform(u1), k1) - jacobi_sc(u1, k1)) < 1.0E-12) then true
 711  else jacobi_sc(u1, k1) # jacobi_sc(rectform(u1), k1) ;
 712 true ;
 713
 714 if complex_floatp(jacobi_sd(u1, k1))
 715  and complex_floatp(jacobi_sd(rectform(u1), k1))
 716  and (abs(jacobi_sd(rectform(u1), k1) - jacobi_sd(u1, k1)) < 1.0E-12) then true
 717  else jacobi_sd(u1, k1) # jacobi_sd(rectform(u1), k1) ;
 718 true ;
 719
 720 if complex_floatp(jacobi_sn(u1, k1))
 721  and complex_floatp(jacobi_sn(rectform(u1), k1))
 722  and (abs(jacobi_sn(rectform(u1), k1) - jacobi_sn(u1, k1)) < 1.0E-12) then true
 723  else jacobi_sn(u1, k1) # jacobi_sn(rectform(u1), k1) ;
 724 true ;
 725
 726 if complex_floatp(elliptic_pi(n1, u1, k1))
 727  and complex_floatp(elliptic_pi(n1, rectform(u1), k1))
 728  and (abs(elliptic_pi(n1, rectform(u1), k1) - elliptic_pi(n1, u1, k1)) < 1.0E-12)
 729  then true else elliptic_pi(n1, u1, k1) # elliptic_pi(n1, rectform(u1), k1) ;
 730 true ;
 731
 732 /*
 733 with_stdout ("/tmp/bar.mac",
 734
 735   for f in elliptic_fcns_1
 736     do block ([a : f(u1), b : f('rectform(u1))],
 737         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
 738         print (true, ";"),
 739         print ("")),
 740
 741   for f in elliptic_fcns_2
 742     do block ([a : f(u1, k1), b : f('rectform(u1), k1)],
 743         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
 744         print (true, ";"),
 745         print ("")),
 746
 747   for f in elliptic_fcns_3
 748     do block ([a : f(n1, u1, k1), b : f(n1, 'rectform(u1), k1)],
 749         print ('if 'complex_bfloatp(a) and 'complex_bfloatp(b) and abs(a - b) < 1b-12 then true else a # b, ";"),
 750         print (true, ";"),
 751         print ("")));
 752  */
 753
 754 (u1 : 0.5b0*(1.0b0 - 0.25b0*%i)^2,
 755  k1 : 0.775b0,
 756  n1 : -1.375b0,
 757  complex_bfloatp(e) := bfloatp(realpart(e)) and bfloatp(imagpart(e)),
 758  0);
 759 0;
 760
 761 if complex_bfloatp(elliptic_ec(u1))
 762  and complex_bfloatp(elliptic_ec(rectform(u1)))
 763  and (abs(elliptic_ec(rectform(u1)) - elliptic_ec(u1)) < 1.0b-12) then true
 764  else elliptic_ec(u1) # elliptic_ec(rectform(u1)) ;
 765 true ;
 766
 767 if complex_bfloatp(elliptic_kc(u1))
 768  and complex_bfloatp(elliptic_kc(rectform(u1)))
 769  and (abs(elliptic_kc(rectform(u1)) - elliptic_kc(u1)) < 1.0b-12) then true
 770  else elliptic_kc(u1) # elliptic_kc(rectform(u1)) ;
 771 true ;
 772
 773 if complex_bfloatp(elliptic_e(u1, k1))
 774  and complex_bfloatp(elliptic_e(rectform(u1), k1))
 775  and (abs(elliptic_e(rectform(u1), k1) - elliptic_e(u1, k1)) < 1.0b-12)
 776  then true else elliptic_e(u1, k1) # elliptic_e(rectform(u1), k1) ;
 777 true ;
 778
 779 if complex_bfloatp(elliptic_eu(u1, k1))
 780  and complex_bfloatp(elliptic_eu(rectform(u1), k1))
 781  and (abs(elliptic_eu(rectform(u1), k1) - elliptic_eu(u1, k1)) < 1.0b-12)
 782  then true else elliptic_eu(u1, k1) # elliptic_eu(rectform(u1), k1) ;
 783 true ;
 784
 785 if complex_bfloatp(elliptic_f(u1, k1))
 786  and complex_bfloatp(elliptic_f(rectform(u1), k1))
 787  and (abs(elliptic_f(rectform(u1), k1) - elliptic_f(u1, k1)) < 1.0b-12)
 788  then true else elliptic_f(u1, k1) # elliptic_f(rectform(u1), k1) ;
 789 true ;
 790
 791 if complex_bfloatp(inverse_jacobi_cd(u1, k1))
 792  and complex_bfloatp(inverse_jacobi_cd(rectform(u1), k1))
 793  and (abs(inverse_jacobi_cd(rectform(u1), k1) - inverse_jacobi_cd(u1, k1)) < 1.0b-12)
 794  then true else inverse_jacobi_cd(u1, k1) # inverse_jacobi_cd(rectform(u1), k1)
 795  ;
 796 true ;
 797
 798 if complex_bfloatp(inverse_jacobi_cn(u1, k1))
 799  and complex_bfloatp(inverse_jacobi_cn(rectform(u1), k1))
 800  and (abs(inverse_jacobi_cn(rectform(u1), k1) - inverse_jacobi_cn(u1, k1)) < 1.0b-12)
 801  then true else inverse_jacobi_cn(u1, k1) # inverse_jacobi_cn(rectform(u1), k1)
 802  ;
 803 true ;
 804
 805 if complex_bfloatp(inverse_jacobi_cs(u1, k1))
 806  and complex_bfloatp(inverse_jacobi_cs(rectform(u1), k1))
 807  and (abs(inverse_jacobi_cs(rectform(u1), k1) - inverse_jacobi_cs(u1, k1)) < 1.0b-12)
 808  then true else inverse_jacobi_cs(u1, k1) # inverse_jacobi_cs(rectform(u1), k1)
 809  ;
 810 true ;
 811
 812 if complex_bfloatp(inverse_jacobi_dc(u1, k1))
 813  and complex_bfloatp(inverse_jacobi_dc(rectform(u1), k1))
 814  and (abs(inverse_jacobi_dc(rectform(u1), k1) - inverse_jacobi_dc(u1, k1)) < 1.0b-12)
 815  then true else inverse_jacobi_dc(u1, k1) # inverse_jacobi_dc(rectform(u1), k1)
 816  ;
 817 true ;
 818
 819 if complex_bfloatp(inverse_jacobi_dn(u1, k1))
 820  and complex_bfloatp(inverse_jacobi_dn(rectform(u1), k1))
 821  and (abs(inverse_jacobi_dn(rectform(u1), k1) - inverse_jacobi_dn(u1, k1)) < 1.0b-12)
 822  then true else inverse_jacobi_dn(u1, k1) # inverse_jacobi_dn(rectform(u1), k1)
 823  ;
 824 true ;
 825
 826 if complex_bfloatp(inverse_jacobi_ds(u1, k1))
 827  and complex_bfloatp(inverse_jacobi_ds(rectform(u1), k1))
 828  and (abs(inverse_jacobi_ds(rectform(u1), k1) - inverse_jacobi_ds(u1, k1)) < 1.0b-12)
 829  then true else inverse_jacobi_ds(u1, k1) # inverse_jacobi_ds(rectform(u1), k1)
 830  ;
 831 true ;
 832
 833 if complex_bfloatp(inverse_jacobi_nc(u1, k1))
 834  and complex_bfloatp(inverse_jacobi_nc(rectform(u1), k1))
 835  and (abs(inverse_jacobi_nc(rectform(u1), k1) - inverse_jacobi_nc(u1, k1)) < 1.0b-12)
 836  then true else inverse_jacobi_nc(u1, k1) # inverse_jacobi_nc(rectform(u1), k1)
 837  ;
 838 true ;
 839
 840 if complex_bfloatp(inverse_jacobi_nd(u1, k1))
 841  and complex_bfloatp(inverse_jacobi_nd(rectform(u1), k1))
 842  and (abs(inverse_jacobi_nd(rectform(u1), k1) - inverse_jacobi_nd(u1, k1)) < 1.0b-12)
 843  then true else inverse_jacobi_nd(u1, k1) # inverse_jacobi_nd(rectform(u1), k1)
 844  ;
 845 true ;
 846
 847 if complex_bfloatp(inverse_jacobi_ns(u1, k1))
 848  and complex_bfloatp(inverse_jacobi_ns(rectform(u1), k1))
 849  and (abs(inverse_jacobi_ns(rectform(u1), k1) - inverse_jacobi_ns(u1, k1)) < 1.0b-12)
 850  then true else inverse_jacobi_ns(u1, k1) # inverse_jacobi_ns(rectform(u1), k1)
 851  ;
 852 true ;
 853
 854 if complex_bfloatp(inverse_jacobi_sc(u1, k1))
 855  and complex_bfloatp(inverse_jacobi_sc(rectform(u1), k1))
 856  and (abs(inverse_jacobi_sc(rectform(u1), k1) - inverse_jacobi_sc(u1, k1)) < 1.0b-12)
 857  then true else inverse_jacobi_sc(u1, k1) # inverse_jacobi_sc(rectform(u1), k1)
 858  ;
 859 true ;
 860
 861 if complex_bfloatp(inverse_jacobi_sd(u1, k1))
 862  and complex_bfloatp(inverse_jacobi_sd(rectform(u1), k1))
 863  and (abs(inverse_jacobi_sd(rectform(u1), k1) - inverse_jacobi_sd(u1, k1)) < 1.0b-12)
 864  then true else inverse_jacobi_sd(u1, k1) # inverse_jacobi_sd(rectform(u1), k1)
 865  ;
 866 true ;
 867
 868 if complex_bfloatp(inverse_jacobi_sn(u1, k1))
 869  and complex_bfloatp(inverse_jacobi_sn(rectform(u1), k1))
 870  and (abs(inverse_jacobi_sn(rectform(u1), k1) - inverse_jacobi_sn(u1, k1)) < 1.0b-12)
 871  then true else inverse_jacobi_sn(u1, k1) # inverse_jacobi_sn(rectform(u1), k1)
 872  ;
 873 true ;
 874
 875 if complex_bfloatp(jacobi_am(u1, k1))
 876  and complex_bfloatp(jacobi_am(rectform(u1), k1))
 877  and (abs(jacobi_am(rectform(u1), k1) - jacobi_am(u1, k1)) < 1.0b-12) then true
 878  else jacobi_am(u1, k1) # jacobi_am(rectform(u1), k1) ;
 879 true ;
 880
 881 if complex_bfloatp(jacobi_cd(u1, k1))
 882  and complex_bfloatp(jacobi_cd(rectform(u1), k1))
 883  and (abs(jacobi_cd(rectform(u1), k1) - jacobi_cd(u1, k1)) < 1.0b-12) then true
 884  else jacobi_cd(u1, k1) # jacobi_cd(rectform(u1), k1) ;
 885 true ;
 886
 887 if complex_bfloatp(jacobi_cn(u1, k1))
 888  and complex_bfloatp(jacobi_cn(rectform(u1), k1))
 889  and (abs(jacobi_cn(rectform(u1), k1) - jacobi_cn(u1, k1)) < 1.0b-12) then true
 890  else jacobi_cn(u1, k1) # jacobi_cn(rectform(u1), k1) ;
 891 true ;
 892
 893 if complex_bfloatp(jacobi_cs(u1, k1))
 894  and complex_bfloatp(jacobi_cs(rectform(u1), k1))
 895  and (abs(jacobi_cs(rectform(u1), k1) - jacobi_cs(u1, k1)) < 1.0b-12) then true
 896  else jacobi_cs(u1, k1) # jacobi_cs(rectform(u1), k1) ;
 897 true ;
 898
 899 if complex_bfloatp(jacobi_dc(u1, k1))
 900  and complex_bfloatp(jacobi_dc(rectform(u1), k1))
 901  and (abs(jacobi_dc(rectform(u1), k1) - jacobi_dc(u1, k1)) < 1.0b-12) then true
 902  else jacobi_dc(u1, k1) # jacobi_dc(rectform(u1), k1) ;
 903 true ;
 904
 905 if complex_bfloatp(jacobi_dn(u1, k1))
 906  and complex_bfloatp(jacobi_dn(rectform(u1), k1))
 907  and (abs(jacobi_dn(rectform(u1), k1) - jacobi_dn(u1, k1)) < 1.0b-12) then true
 908  else jacobi_dn(u1, k1) # jacobi_dn(rectform(u1), k1) ;
 909 true ;
 910
 911 if complex_bfloatp(jacobi_ds(u1, k1))
 912  and complex_bfloatp(jacobi_ds(rectform(u1), k1))
 913  and (abs(jacobi_ds(rectform(u1), k1) - jacobi_ds(u1, k1)) < 1.0b-12) then true
 914  else jacobi_ds(u1, k1) # jacobi_ds(rectform(u1), k1) ;
 915 true ;
 916
 917 if complex_bfloatp(jacobi_nc(u1, k1))
 918  and complex_bfloatp(jacobi_nc(rectform(u1), k1))
 919  and (abs(jacobi_nc(rectform(u1), k1) - jacobi_nc(u1, k1)) < 1.0b-12) then true
 920  else jacobi_nc(u1, k1) # jacobi_nc(rectform(u1), k1) ;
 921 true ;
 922
 923 if complex_bfloatp(jacobi_nd(u1, k1))
 924  and complex_bfloatp(jacobi_nd(rectform(u1), k1))
 925  and (abs(jacobi_nd(rectform(u1), k1) - jacobi_nd(u1, k1)) < 1.0b-12) then true
 926  else jacobi_nd(u1, k1) # jacobi_nd(rectform(u1), k1) ;
 927 true ;
 928
 929 if complex_bfloatp(jacobi_ns(u1, k1))
 930  and complex_bfloatp(jacobi_ns(rectform(u1), k1))
 931  and (abs(jacobi_ns(rectform(u1), k1) - jacobi_ns(u1, k1)) < 1.0b-12) then true
 932  else jacobi_ns(u1, k1) # jacobi_ns(rectform(u1), k1) ;
 933 true ;
 934
 935 if complex_bfloatp(jacobi_sc(u1, k1))
 936  and complex_bfloatp(jacobi_sc(rectform(u1), k1))
 937  and (abs(jacobi_sc(rectform(u1), k1) - jacobi_sc(u1, k1)) < 1.0b-12) then true
 938  else jacobi_sc(u1, k1) # jacobi_sc(rectform(u1), k1) ;
 939 true ;
 940
 941 if complex_bfloatp(jacobi_sd(u1, k1))
 942  and complex_bfloatp(jacobi_sd(rectform(u1), k1))
 943  and (abs(jacobi_sd(rectform(u1), k1) - jacobi_sd(u1, k1)) < 1.0b-12) then true
 944  else jacobi_sd(u1, k1) # jacobi_sd(rectform(u1), k1) ;
 945 true ;
 946
 947 if complex_bfloatp(jacobi_sn(u1, k1))
 948  and complex_bfloatp(jacobi_sn(rectform(u1), k1))
 949  and (abs(jacobi_sn(rectform(u1), k1) - jacobi_sn(u1, k1)) < 1.0b-12) then true
 950  else jacobi_sn(u1, k1) # jacobi_sn(rectform(u1), k1) ;
 951 true ;
 952
 953 if complex_bfloatp(elliptic_pi(n1, u1, k1))
 954  and complex_bfloatp(elliptic_pi(n1, rectform(u1), k1))
 955  and (abs(elliptic_pi(n1, rectform(u1), k1) - elliptic_pi(n1, u1, k1)) < 1.0b-12)
 956  then true else elliptic_pi(n1, u1, k1) # elliptic_pi(n1, rectform(u1), k1) ;
 957 true ;
 958
 959 /* Test derivatives at a few random points using numerical differentiation
 960
 961    func  - a function
 962    df - deriviative of func wrt p-th argument
 963    p - derivative wrt p-th argument
 964    vars - variables in expression deriv
 965    table - table of points to evaluate
 966    eps - report errors > eps
 967    delta - offset for numerical derivative
 968 */
 969 (test_deriv(func,df,vars,p,table,eps,delta) :=
 970 block([badpoints : [],
 971        abserr    : 0,
 972        maxerr    : -1],
 973   for %z in table do
 974     block([%z0, %z1, deriv, nderiv],
 975       %z1:makelist(if i=p then %z[i]+delta else %z[i],i,1,length(%z)),
 976       %z0:makelist(if i=p then %z[i]-delta else %z[i],i,1,length(%z)),
 977       nderiv: float((apply(func,%z1)-apply(func,%z0))/(2*delta)),
 978       deriv : float(subst(maplist("=",vars,%z),df)),
 979       abserr : abs(nderiv-deriv),
 980       maxerr : max(maxerr,abserr),
 981       if abserr > eps then
 982         badpoints : cons ([%z,deriv,nderiv,abserr],badpoints)
 983     ),
 984   if badpoints # [] then
 985     cons(maxerr,badpoints)
 986   else
 987     badpoints
 988     ),done);
 989 done;
 990
 991 /* test points for two-arg functions */
 992 (l2: [[0.2, 0.3], [0.2, 0.5], [1.5, 0.7], [0.5, 0.99], [1.5, 0.8],
 993     [1.57, 0.8], [1.5707, 0.8], [-1.0, 0.8], [2.0, 0.8]],done);
 994 done;
 995
 996 /* test points for three-arg functions */
 997 /* bug #3221 elliptic_pi() wrong for 2nd arg > float(%pi/2) */
 998 (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],
 999     [0.2,1.5,0.8],[0.2,1.57,0.8],[0.2,1.5707,0.8]],done);
1000 done;
1001
1002 test_deriv('elliptic_f,diff(elliptic_f(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
1003 [];
1004
1005 test_deriv('elliptic_f,diff(elliptic_f(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
1006 [];
1007
1008 test_deriv('elliptic_e,diff(elliptic_e(z,m),z),[z,m],1,l2,2.0e-7,1.0e-8);
1009 [];
1010
1011 test_deriv('elliptic_e,diff(elliptic_e(z,m),m),[z,m],2,l2,2.0e-7,1.0e-8);
1012 [];
1013
1014 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),n),[n,z,m],1,l3,2.0e-7,1.0e-8);
1015 [];
1016
1017 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),z),[n,z,m],2,l3,2.0e-7,1.0e-8);
1018 [];
1019
1020 test_deriv('elliptic_pi,diff(elliptic_pi(n,z,m),m),[n,z,m],3,l3,2.0e-7,1.0e-8);
1021 [];
1022
1023 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1024 closeto(e,tol):=block([numer:true,abse],abse:abs(e),if(abse<tol) then true else abse);
1025
1026 /*
1027  * Some additional tests for complex values. Use the fact that
1028  * elliptic_pi(0, phi, m) = elliptic_f(phi,m);
1029  */
1030 closeto(elliptic_pi(0, 0.5, 0.25)- elliptic_f(0.5, 0.25), 1e-15);
1031 true;
1032 closeto(elliptic_pi(0, 0.5*%i, 0.25)- elliptic_f(0.5*%i, 0.25), 1e-15);
1033 true;
1034
1035 closeto(elliptic_pi(0, 0.5b0, 0.25b0)- elliptic_f(0.5b0, 0.25b0), 1e-15);
1036 true;
1037 closeto(elliptic_pi(0, 0.5b0*%i, 0.25b0)- elliptic_f(0.5b0*%i, 0.25b0), 1e-15);
1038 true;
1039
1040 /* Test for Bug 3221 */
1041
1042 closeto(elliptic_pi(0.2,1.57,0.8) - 2.5770799919605668, 1e-14);
1043 true;
1044 closeto(elliptic_pi(0.2,1.58,0.8) - 2.60502920656026151, 1e-14);
1045 true;
1046 closeto(elliptic_pi(0.20,2.00,0.80) - 3.6543835090829025, 7.47e-11);
1047 true;
1048
1049 (oldfpprec:fpprec,fpprec:32);
1050 32;
1051 closeto(elliptic_pi(0.2b0,1.57b0,0.8b0) - 2.5770799919605668058849721013196b0, 3.09b-32);
1052 true;
1053 closeto(elliptic_pi(0.2b0,1.58b0,0.8b0) - 2.6050292065602615145481924917132b0, 1.24b-32);
1054 true;
1055 closeto(elliptic_pi(0.2b0,2.0b0,0.8b0) - 3.6543835090082902501670624830938b0, 1b-32);
1056 true;
1057
1058 (fpprec:oldfpprec,kill(l2,l3,test_deriv,oldfpprec),done);
1059 done;