home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / ELLIPTIC.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  48.4 KB  |  1,295 lines
  1.  
  2. (* :Name: StartUp`Elliptic` *)
  3.  
  4. (* :Title: Elliptic Functions *)
  5.  
  6. (* :Author: Daniel R. Grayson, Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Extends the coverage of elliptic functions from that in the kernel
  10.     to include the Jacobi elliptic functions and their inverses, the
  11.     Neville theta functions, the amplitude function am, the nome
  12.     function q, and the Weierstrass P function and its derivative.
  13. *)
  14.  
  15. (* :Context: System` *)
  16.  
  17. (* :Package Version: 2.0 *)
  18.  
  19. (* :Copyright: Copyright 1990  Wolfram Research, Inc.
  20.  
  21.     Permission is hereby granted to modify and/or make copies of
  22.     this file for any purpose other than direct profit, or as part
  23.     of a commercial product, provided this copyright notice is left
  24.     intact.  Sale, other than for the cost of media, is prohibited.
  25.  
  26.     Permission is hereby granted to reproduce part or all of
  27.     this file, provided that the source is acknowledged.
  28. *)
  29.  
  30. (* :History:
  31.     Extensively revised by Jerry B. Keiper, November 1990.
  32.     Version 1.0 by Daniel R. Grayson, 1988.
  33. *)
  34.  
  35. (* :Keywords:
  36.     elliptic functions, theta functions, Weierstrass P function
  37. *)
  38.  
  39. (* :Source:
  40.     Abramowitz and Stegun: "Handbook of Mathematical Functions",
  41.         Dover Publications, New York. chap. 16, 17.
  42.     Whittaker and Watson: "A Course of Modern Analysis",
  43.         Cambridge University Press.  chap. 20, 21, 22.
  44. *)
  45.  
  46. (* :Warning:
  47.     Adds rules to EllipticExp[ ], EllipticLog[ ], and InverseFunction[ ]
  48. *)
  49.  
  50. (* :Mathematica Version: 2.0 *)
  51.  
  52. (* :Limitation:
  53.     Affected by inexactness in input numbers especially near
  54.     branch cuts in the inverse functions.
  55.  
  56.     Ranges of certain inverse Jacobi elliptic functions differ from
  57.     the ranges of the corresponding limiting case inverse trigonometric
  58.     (both circular and hyperbolic) functions.
  59.  
  60.     This package is not complete; there are many more features which
  61.     could be added.  The theory of elliptic functions is too rich to
  62.     be able to cover everything without a great deal of work.    The
  63.     Weierstrass functions are much less complete than are the Jacobi
  64.     functions.
  65.  
  66.     Weierstrass functions do not work with complex invariants g2 and g3.
  67. *)
  68.  
  69. (* :Discussion:
  70.     This package defines the mathematical functions JacobiSN, JacobiCN,
  71.     JacobiDN, JacobiCD, JacobiSD, JacobiND, JacobiDC, JacobiNC, JacobiSC,
  72.     JacobiNS, JacobiDS, JacobiCS, InverseJacobiSN, InverseJacobiCN,
  73.     InverseJacobiDN, InverseJacobiCD, InverseJacobiSD, InverseJacobiND,
  74.     InverseJacobiDC, InverseJacobiNC, InverseJacobiSC, InverseJacobiNS,
  75.     InverseJacobiDS, InverseJacobiCS, JacobiAmplitude, EllipticNomeQ, 
  76.     EllipticThetaS, EllipticThetaC, EllipticThetaD, EllipticThetaN,
  77.     WeierstrassP, WeierstrassPPrime, and InverseWeierstrassP as well
  78.     as the derivatives of most of these functions.
  79.  
  80.     The Jacobi elliptic functions are evaluated in terms of the Neville
  81.     theta functions after argument reductions on the first argument.
  82.     The inverse Jacobi elliptic functions are evaluated in terms
  83.     EllipticLog and then the result is shifted to the appropriate
  84.     parallelogram.  JacobiAmplitude is evaluated as the ArcSin of
  85.     JacobiSN and then the result is shifted to the appropriate value.
  86.     EllipticNomeQ is evaluated in terms of EllipticK[m] and EllipticK[1-m].
  87.     The Neville theta functions are evaluated in terms of the Jacobi
  88.     theta functions.  WeierstrassP and WeierstrassPPrime are evaluated
  89.     in terms of EllipticExp and InverseWeierstrassP is evaluated in
  90.     terms of EllipticLog.  The result of InverseWeierstrassP is not 
  91.     shifted to any particular period parallelogram since the basic
  92.     periods are not calulated.
  93. *)
  94.  
  95. Needs["System`", "StartUp`Attributes.m"]
  96.  
  97. Begin["System`"]
  98.  
  99. Unprotect[ JacobiSN, JacobiSD, JacobiSC,
  100.            JacobiCS, JacobiCN, JacobiCD,
  101.            JacobiNC, JacobiND, JacobiNS,
  102.            JacobiDN, JacobiDS, JacobiDC,
  103.            JacobiAmplitude,
  104.            InverseJacobiSN, InverseJacobiSD, InverseJacobiSC,
  105.            InverseJacobiCS, InverseJacobiCN, InverseJacobiCD,
  106.            InverseJacobiNS, InverseJacobiNC, InverseJacobiND,
  107.            InverseJacobiDS, InverseJacobiDC, InverseJacobiDN,
  108.            EllipticNomeQ, EllipticThetaS, EllipticThetaC,
  109.            EllipticThetaN, EllipticThetaD,
  110.            WeierstrassP, WeierstrassPPrime,
  111.            InverseWeierstrassP,
  112.            EllipticExp, EllipticLog]
  113.  
  114. Map[ Clear,
  115.           {JacobiSN, JacobiSD, JacobiSC,
  116.            JacobiCS, JacobiCN, JacobiCD,
  117.            JacobiNC, JacobiNS, JacobiND,
  118.            JacobiDN, JacobiDS, JacobiDC,
  119.            JacobiAmplitude,
  120.            InverseJacobiSN, InverseJacobiSD, InverseJacobiSC,
  121.            InverseJacobiCS, InverseJacobiCN, InverseJacobiCD,
  122.            InverseJacobiNS, InverseJacobiNC, InverseJacobiND,
  123.            InverseJacobiDS, InverseJacobiDC, InverseJacobiDN,
  124.            EllipticNomeQ, EllipticThetaS, EllipticThetaC,
  125.            EllipticThetaN, EllipticThetaD,
  126.            WeierstrassP, WeierstrassPPrime,
  127.            InverseWeierstrassP,
  128.            EllipticExp, EllipticLog}]
  129.  
  130. (* ---------------------- Jacobi elliptic functions ---------------------- *)
  131.  
  132. JacobiSN::usage =
  133. "JacobiSN[u, m] gives the Jacobi elliptic function sn at u for the parameter m."
  134. JacobiSD::usage =
  135. "JacobiSD[u, m] gives the Jacobi elliptic function sd at u for the parameter m."
  136. JacobiSC::usage =
  137. "JacobiSC[u, m] gives the Jacobi elliptic function sc at u for the parameter m."
  138. JacobiCS::usage =
  139. "JacobiCS[u, m] gives the Jacobi elliptic function cs at u for the parameter m."
  140. JacobiCN::usage =
  141. "JacobiCN[u, m] gives the Jacobi elliptic function cn at u for the parameter m."
  142. JacobiCD::usage =
  143. "JacobiCD[u, m] gives the Jacobi elliptic function cd at u for the parameter m."
  144. JacobiNS::usage =
  145. "JacobiNS[u, m] gives the Jacobi elliptic function ns at u for the parameter m."
  146. JacobiNC::usage =
  147. "JacobiNC[u, m] gives the Jacobi elliptic function nc at u for the parameter m."
  148. JacobiND::usage =
  149. "JacobiND[u, m] gives the Jacobi elliptic function nd at u for the parameter m."
  150. JacobiDN::usage =
  151. "JacobiDN[u, m] gives the Jacobi elliptic function dn at u for the parameter m."
  152. JacobiDS::usage =
  153. "JacobiDS[u, m] gives the Jacobi elliptic function ds at u for the parameter m."
  154. JacobiDC::usage =
  155. "JacobiDC[u, m] gives the Jacobi elliptic function dc at u for the parameter m."
  156.  
  157. JacobiAmplitude::usage =
  158. "JacobiAmplitude[u, m] gives the amplitude for Jacobi elliptic functions."
  159.  
  160. InverseJacobiSN::usage =
  161. "InverseJacobiSN[v, m] gives the inverse of the Jacobi elliptic function sn at
  162. v for the parameter m.  The result will be in the parallelogram with sides 2K
  163. and 2I K' centered at 0."
  164. InverseJacobiSD::usage =
  165. "InverseJacobiSD[v, m] gives the inverse of the Jacobi elliptic function sd at
  166. v for the parameter m.  The result will be in the parallelogram with sides 2K
  167. and 2I K' centered at 0."
  168. InverseJacobiSC::usage =
  169. "InverseJacobiSC[v, m] gives the inverse of the Jacobi elliptic function sc at
  170. v for the parameter m.  The result will be in the parallelogram with sides 2K
  171. and 2I K' centered at 0."
  172. InverseJacobiCS::usage =
  173. "InverseJacobiCS[v, m] gives the inverse of the Jacobi elliptic function cs at
  174. v for the parameter m.  The result will be in the parallelogram with sides 2K
  175. and 2I K' centered at K."
  176. InverseJacobiCN::usage =
  177. "InverseJacobiCN[v, m] gives the inverse of the Jacobi elliptic function cn at
  178. v for the parameter m.  The result will be in the parallelogram with sides 2K
  179. and 2I K' centered at K."
  180. InverseJacobiCD::usage =
  181. "InverseJacobiCD[v, m] gives the inverse of the Jacobi elliptic function cd at
  182. v for the parameter m.  The result will be in the parallelogram with sides 2K
  183. and 2I K' centered at K."
  184. InverseJacobiNS::usage =
  185. "InverseJacobiNS[v, m] gives the inverse of the Jacobi elliptic function ns at
  186. v for the parameter m.  The result will be in the parallelogram with sides 2K
  187. and 2I K' centered at I K'."
  188. InverseJacobiNC::usage =
  189. "InverseJacobiNC[v, m] gives the inverse of the Jacobi elliptic function nc at
  190. v for the parameter m.  The result will be in the parallelogram with sides 2K
  191. and 2I K' centered at I K'."
  192. InverseJacobiDS::usage =
  193. "InverseJacobiDS[v, m] gives the inverse of the Jacobi elliptic function ds at
  194. v for the parameter m.  The result will be in the parallelogram with sides 2K
  195. and 2I K' centered at K + I K'."
  196. InverseJacobiDC::usage =
  197. "InverseJacobiDC[v, m] gives the inverse of the Jacobi elliptic function dc at
  198. v for the parameter m.  The result will be in the parallelogram with sides 2K
  199. and 2I K' centered at K + I K'."
  200. InverseJacobiND::usage =
  201. "InverseJacobiND[v, m] gives the inverse of the Jacobi elliptic function nd at
  202. v for the parameter m.  The result will be in the parallelogram with sides 2K
  203. and 2I K' centered at I K'."
  204. InverseJacobiDN::usage =
  205. "InverseJacobiDN[v, m] gives the inverse of the Jacobi elliptic function dn at
  206. v for the parameter m.  The result will be in the parallelogram with sides 2K
  207. and 2I K' centered at K + I K'."
  208.  
  209. EllipticNomeQ::usage =
  210. "EllipticNomeQ[m] gives the nome q == Exp[-Pi EllipticK[1-m]/EllipticK[m]]."
  211.  
  212. EllipticThetaS::usage =
  213. "EllipticThetaS[u, m] gives Neville's elliptic theta function theta_s(u, m)."
  214. EllipticThetaC::usage =
  215. "EllipticThetaC[u, m] gives Neville's elliptic theta function theta_c(u, m)."
  216. EllipticThetaN::usage =
  217. "EllipticThetaN[u, m] gives Neville's elliptic theta function theta_n(u, m)."
  218. EllipticThetaD::usage =
  219. "EllipticThetaD[u, m] gives Neville's elliptic theta function theta_d(u, m)."
  220.  
  221.  (* ---------------------- Weierstrass functions ---------------------- *)
  222.  
  223. WeierstrassP::usage =
  224. "WeierstrassP[u, g2, g3] gives the Weierstrass elliptic function P."
  225.  
  226. WeierstrassPPrime::usage =
  227. "WeierstrassPPrime[u, g2, g3] gives the derivative with respect to u of the
  228. Weierstrass elliptic function P'."
  229.  
  230. InverseWeierstrassP::usage =
  231. "InverseWeierstrassP[{P, P'}, g2, g3] gives the value of u such that
  232. P = WeierstrassP[u, g2, g3] and P' = WeierstrassPPrime[u, g2, g3].  Note
  233. that P and P' are not independent.  No attempt is made to ensure that the
  234. result is in any particular period parallelogram."
  235.  
  236.  (* ---------------------- end of introduction ---------------------- *)
  237.  
  238. Begin["Elliptic`Private`"]
  239.  
  240.  (* ------------------ derivative of EllipticLog -----------------------*)
  241.  
  242. Unprotect[EllipticExp, EllipticLog, InverseFunction];
  243.  
  244. Off[EllipticLog::ellnp];
  245. EllipticLog/: D[EllipticLog[{x_,y_}, _List], t_] := D[x, t]/(2y);
  246. On[EllipticLog::ellnp];
  247.  
  248.  (* ----------------------- private functions -------------------------*)
  249.  
  250. (* :Note:
  251.     The following variables are used as storage for previous values
  252.     that are likely to be needed again.  We thus avoid the need to
  253.     actually re-evaluate them.  They are paired as input and output
  254.     values for various functions.
  255. *)
  256.  
  257. $LastEllipticmK;    $LastEllipticK;
  258. $LastEllipticmKp;    $LastEllipticKp;
  259. $LastEllipticmT;    $LastEllipticT;
  260. $LastEllipticmQ;    $LastEllipticQ;
  261. $LastOneIn;        $LastOneOut;
  262.  
  263. $EllipK[m_?NumberQ] :=
  264.     Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmK]},
  265.     If[m === $LastEllipticmK && prec <= oldprec,
  266.         N[$LastEllipticK, prec],
  267.       (* else *)
  268.         prec = EllipticK[m];
  269.         AbortProtect[If[NumberQ[prec],
  270.             $LastEllipticmK = m; $LastEllipticK = prec]];
  271.         prec]
  272.     ]
  273.  
  274. $EllipKp[m_?NumberQ] :=
  275.     Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmKp]},
  276.     If[m === $LastEllipticmKp && prec <= oldprec,
  277.         N[$LastEllipticKp, prec],
  278.       (* else *)
  279.         prec = EllipticK[1-m];
  280.         AbortProtect[If[NumberQ[prec],
  281.             $LastEllipticmKp = m; $LastEllipticKp = prec]];
  282.         prec]
  283.     ]
  284.  
  285. EllipticNomeQ[m_?NumberQ] :=
  286.     Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmQ]},
  287.     If[TrueQ[m == 0.], Return[m]];
  288.     If[m === $LastEllipticmQ && prec <= oldprec,
  289.         N[$LastEllipticQ, prec],
  290.       (* else *)
  291.         prec = Exp[-N[Pi, prec] $EllipKp[m]/$EllipK[m]];
  292.         AbortProtect[If[NumberQ[prec],
  293.             $LastEllipticmQ = m; $LastEllipticQ = prec]];
  294.         prec]
  295.     ]
  296.  
  297. $EllipticTheta0[n_, m_] :=
  298.     Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmT]},
  299.     If[m === $LastEllipticmT && prec <= oldprec,
  300.         N[$LastEllipticT[[n]], prec],
  301.       (* else *)
  302.         prec = {EllipticThetaPrime[1, 0, m], EllipticTheta[2, 0, m],
  303.         EllipticTheta[3, 0, m], EllipticTheta[4, 0, m]};
  304.         AbortProtect[If[Apply[And, Map[NumberQ, prec]],
  305.         $LastEllipticmT = m; $LastEllipticT = prec]];
  306.         prec[[n]]]
  307.     ]
  308.  
  309. $EllipticTheta[p_, u_, m_] :=
  310.     Module[{prec = Precision[{u, m}], K2, mm, q},
  311.     mm = N[m, prec];
  312.     K2 = 2 $EllipK[mm];
  313.     q = EllipticNomeQ[mm];
  314.     If[p == 0,
  315.         K2 EllipticTheta[1, N[Pi u/K2, prec], q]/
  316.         (N[Pi, prec] $EllipticTheta0[1, q]),
  317.       (* else *)
  318.         EllipticTheta[p+1, N[Pi u/K2, prec], q]/$EllipticTheta0[p+1, q]
  319.         ]
  320.     ]
  321.  
  322. EllipticThetaS[u_?NumberQ, m_?NumberQ] :=
  323.     $EllipticTheta[0, u, m] /; Precision[{u, m}] < Infinity
  324.  
  325. EllipticThetaC[u_?NumberQ, m_?NumberQ] :=
  326.     $EllipticTheta[1, u, m] /; Precision[{u, m}] < Infinity
  327.  
  328. EllipticThetaD[u_?NumberQ, m_?NumberQ] :=
  329.     $EllipticTheta[2, u, m] /; Precision[{u, m}] < Infinity
  330.  
  331. EllipticThetaN[u_?NumberQ, m_?NumberQ] :=
  332.     $EllipticTheta[3, u, m] /; Precision[{u, m}] < Infinity
  333.  
  334. EllipticThetaS[0, m_] := 0
  335. EllipticThetaC[0, m_] := 1
  336. EllipticThetaD[0, m_] := 1
  337. EllipticThetaN[0, m_] := 1
  338.  
  339. Fix[eqn_Equal,x_,y_,form_] :=
  340.     Module[{newx=x,newy=y,oldx=x,oldy=y,c,f,sub,newform = form,newout,
  341.         null = {x->0,y->0},a,b,disc,e1,precision,c4,c6,rtdisc,s},
  342.     If[ $LastOneIn === {eqn,x,y,form}, Return[$LastOneOut]];
  343.  
  344.     f = Expand[eqn[[1]]-eqn[[2]]];
  345.     precision = Precision[f];
  346.     If [ precision == Infinity , precision = 16 ];
  347.  
  348.     (* remove coeff of y^2 *)
  349.     f = Expand[f / (D[f,{y,2}]/2)];
  350.  
  351.     (* remove coeff of x^3 *)
  352.     c = -(D[f,{x,3}]/6);
  353.     sub = {x -> x/c, y -> y/c};
  354.     oldx = Expand[oldx c];
  355.     oldy = Expand[oldy c];
  356.     f = Expand[(f /. sub) c^2];
  357.     newform = Expand[newform /. sub];
  358.     newx = Expand[newx /. sub];
  359.     newy = Expand[newy /. sub];
  360.  
  361.     (* remove linear y terms *)
  362.     c = ((D[f,y] /. null) + (D[f,x,y] /. null) x) / 2;
  363.     sub = {y -> y-c};
  364.     oldy = Expand[oldy + c];
  365.     f = Expand[(f /. sub)];
  366.     newform = Expand[newform /. sub];
  367.     newx = Expand[newx /. sub];
  368.     newy = Expand[newy /. sub];
  369.  
  370.     (* remove x^2 term *)
  371.     c = (D[f,{x,2}]/2 /. null)/ 3;
  372.     sub = {x -> x+c};
  373.     oldx = Expand[oldx - c];
  374.     f = Expand[(f /. sub)];
  375.     newform = Expand[newform /. sub];
  376.     newx = Expand[newx /. sub];
  377.     newy = Expand[newy /. sub];
  378.  
  379.     (* convert to inexact *)
  380.     f = N[f,precision];
  381.  
  382.     (* compute the largest real root *)
  383.     disc = 4 (D[f,x] /. null)^3 - 27 (f /. null)^2 ;
  384.     c4 = (D[f,x] /. null) 48;
  385.     c6 = (f /. null) 864;
  386.     If[Head[disc] =!= Real && disc =!=0, Return[$Failed]];
  387.     If[ disc < 0,
  388.         rtdisc = 96 Sqrt [ - 3 disc ];
  389.         e1 = ((c6+rtdisc)^(1/3) + (c6-rtdisc)^(1/3))/12,
  390.         (* else *)
  391.         s = (c6/1728 + I Sqrt[ disc/108 ])^(1/3);
  392.         e1 = 2 Re[s],
  393.         (* undecided *)
  394.         Return[$Failed]
  395.     ];
  396.     sub = {x->x+e1};
  397.     f = Expand[(f /. sub)];
  398.     newform = Expand[newform /. sub];
  399.     newx = Expand[newx /. sub];
  400.     newy = Expand[newy /. sub];
  401.     oldx = Expand[oldx - e1];
  402.     a = - (D[f,{x,2}] /. null)/2;
  403.     b = - D[f,x] /. null;
  404.     c = Expand[newform / (Dt[x] / (2y))];
  405.     (* sigh, what if c hasn t simplified down to a number here ? *)
  406.     If[c==0,c=1];
  407.     newout = {{a,b}, Apply[Function,{{newx, newy} /. {x->#1, y->#2}}], 
  408.             Apply[Function,{{oldx, oldy} /. {x->#1, y->#2}}], c};
  409.  
  410.     (* don't let interrupts break the memory function *)
  411.     AbortProtect[{$LastOneIn, $LastOneOut} = {in, newout}];
  412.     newout
  413.     ]
  414.  
  415. EqnEllipticExp[u_, {eqn_Equal, x_, y_, form_}] :=
  416.     Module[{c = Fix[eqn,x,y,form]},
  417.     If[Head[c] === List && Length[c] == 4,
  418.         c = Apply[c[[2]], EllipticExp[u/c[[4]], c[[1]]]]];
  419.     c
  420.     ]
  421.  
  422. EqnEllipticLog[xy_List, {eqn_Equal, x_, y_, form_}] :=
  423.     Module[{c = Fix[eqn,x,y,form]},
  424.     If[Head[c] === List && Length[c] == 4,
  425.         c = c[[4]] EllipticLog[Apply[c[[3]], xy], c[[1]]]];
  426.     c
  427.     ]
  428.  
  429. EqnEllipticExp[u_, {eqn_Equal, x_, y_}] := EqnEllipticExp[u, {eqn, x, y, 0}] 
  430. EqnEllipticLog[u_, {eqn_Equal, x_, y_}] := EqnEllipticLog[u, {eqn, x, y, 0}] 
  431.  
  432. (* this one fills in the square root *)
  433. $EllipticLog[xi_,a_,b_] :=
  434.     EllipticLog[{xi, Sqrt[xi(xi(xi + a) + b)]}, {a,b}]
  435.  
  436. (* this one is useful for integrals with two arbitrary limits *)
  437. $EllipticLog[x0_,x1_,a_,b_] := $EllipticLog[x1,a,b] - $EllipticLog[x0,a,b]
  438.  
  439. (* this converts from the parameter m to the pair a,b *)
  440. $EllipticLogm[x_, m_] := $EllipticLog[x, m+1, m]
  441.  
  442.   (* ---------------- Weierstrass elliptic functions ---------------- *)
  443.  
  444. Weierstrass[u_?NumberQ, g2_?NumberQ, g3_?NumberQ] := 
  445.     EqnEllipticExp[u, {Y^2 == 4X^3 - g2 X - g3 , X, Y, Dt[X]/Y}]
  446.  
  447. WeierstrassP[u_, g2_, g3_] :=
  448.     Module[{p}, p[[1]] /; (Head[p = Weierstrass[u, g2, g3]] === List)];
  449.  
  450. WeierstrassPPrime[u_, g2_, g3_] :=
  451.     Module[{p}, p[[2]] /; (Head[p = Weierstrass[u, g2, g3]] === List)];
  452.  
  453. InverseWeierstrassP[{p_, pp_}, g2_, g3_] := 
  454.     Module[{u, Y, X},
  455.     u /; (u = EqnEllipticLog[{p, pp}, {Y^2 == 4X^3-g2 X-g3, X, Y, Dt[X]/Y}];
  456.         NumberQ[u])
  457.     ];
  458.  
  459.   (* ---------------- Weierstrass  derivatives ----------------------- *)
  460.  
  461. Derivative[n_Integer, 0, 0][WeierstrassP] ^:=
  462.     Derivative[n-1, 0, 0][WeierstrassPPrime] /; n > 0
  463.  
  464. Derivative[1, 0, 0][WeierstrassPPrime] ^=
  465.     (6 WeierstrassP[#1, #2, #3]^2 - #2/2)&;        (* AS 18.6.4, p640 *)
  466.  
  467. Derivative[n_, 0, 0][WeierstrassPPrime] ^:=
  468.     Module[{i, r},
  469.         r = Evaluate[12 Sum[Binomial[n-2, i]
  470.             Derivative[i, 0, 0][WeierstrassPPrime][#1, #2, #3] 
  471.             Derivative[n-2-i, 0, 0][WeierstrassP][#1, #2, #3],
  472.             {i, 0, n-2}]] &;
  473.         r = Share[r];
  474.         AbortProtect[
  475.         Unprotect[WeierstrassPPrime];
  476.         Literal[Derivative[n, 0, 0][WeierstrassPPrime]] ^= r;
  477.         Protect[WeierstrassPPrime]
  478.         ];
  479.         r
  480.         ] /; n > 1
  481.  
  482.   (* -----------elementary special cases, AS, p571, 16.6 ---------------- *)
  483.  
  484. JacobiSN[u_, 0|0.] := Sin[u];    JacobiSN[u_, 1|1.] := Tanh[u];
  485. JacobiCN[u_, 0|0.] := Cos[u];    JacobiCN[u_, 1|1.] := Sech[u];
  486. JacobiDN[u_, 0|0.] := 1;    JacobiDN[u_, 1|1.] := Sech[u];
  487.  
  488. JacobiCD[u_, 0|0.] := Cos[u];    JacobiCD[u_, 1|1.] := 1;
  489. JacobiSD[u_, 0|0.] := Sin[u];    JacobiSD[u_, 1|1.] := Sinh[u];
  490. JacobiND[u_, 0|0.] := 1;    JacobiND[u_, 1|1.] := Cosh[u];
  491.  
  492. JacobiDC[u_, 0|0.] := Sec[u];    JacobiDC[u_, 1|1.] := 1;
  493. JacobiNC[u_, 0|0.] := Sec[u];    JacobiNC[u_, 1|1.] := Cosh[u];
  494. JacobiSC[u_, 0|0.] := Tan[u];    JacobiSC[u_, 1|1.] := Sinh[u];
  495.  
  496. JacobiNS[u_, 0|0.] := Csc[u];    JacobiNS[u_, 1|1.] := Coth[u];
  497. JacobiDS[u_, 0|0.] := Csc[u];    JacobiDS[u_, 1|1.] := Csch[u];
  498. JacobiCS[u_, 0|0.] := Cot[u];    JacobiCS[u_, 1|1.] := Csch[u];
  499.  
  500.  (* -------- special symbolic and exact cases, AS, p571, 16.5 ------------*)
  501.  
  502.     (* note: row 7 has been changed to deal with the singularity *)
  503.  
  504. $JacobiSpecialValues =
  505.     {{0&, 1&, 1&},
  506.     {(1 + (1-#)^(1/2))^(-1/2)&, (1-#)^(1/4)(1 + (1-#)^(1/2))^(-1/2)&,
  507.         (1-#)^(1/4)&},
  508.     {1&, 0&, (1-#)^(1/2)&},
  509.     {I #^(-1/4)&, (1 + #^(1/2))^(1/2)(#^(-1/4))&, (1 + #^(1/2))^(1/2)&},
  510.     {2^(-1/2) #^(-1/4) ((1 + #^(1/2))^(1/2) + I (1 - #^(1/2))^(1/2))&,
  511.         2^(-1/2) (1-I) (1-#)^(1/4) (#^(-1/4))&,
  512.         2^(-1/2) (1-#)^(1/4) ((1 + (1-#)^(1/2))^(1/2) -
  513.             I (1 - (1-#)^(1/2))^(1/2))&},
  514.     {#^(-1/4)&, -I (1 - #^(1/2))^(1/2) (#^(-1/4))&, (1 - #^(1/2))^(1/2)&},
  515.     {1&, -I&, -I #^(1/2)&}  (* this row is ss, cs, ds at I K' *),
  516.     {(1 - (1-#)^(1/2))^(-1/2)&, -I (1-(1-#)^(1/2))^(-1/2) ((1-#)^(1/4))&,
  517.         -I (1-#)^(1/4)&},
  518.     {#^(-1/2)&, -I (1-#)^(1/2) (#^(-1/2))&, 0&},
  519.     (* row 10 is sn, cn, dn at 3K/2 + I K'/2.  It is not in AS *)
  520.     {((1-I)((1+(1-#)^(1/2))^(1/2)+I(1-(1-#)^(1/2))^(1/2))/(2 #^(1/4)))&,
  521.        (-((1-#)/#)^(1/4)((1+#^(1/2))^(1/2)+I(1-#^(1/2))^(1/2))/
  522.         ((1+(1-#)^(1/2))^(1/2)-I(1-(1-#)^(1/2))^(1/2)))&,
  523.        ((1-#)^(1/4)2^(-1/2)((1+(1-#)^(1/2))^(1/2)+
  524.         I(1-(1-#)^(1/2))^(1/2)))&}};
  525.  
  526. SymbJacobi[p_, q_, u_, m_] :=
  527.     Module[{K, K1, r, s, sign = 1},
  528.     K = EllipticK[m];
  529.     K1 = EllipticK[1 - m];
  530.     If[Head[K] =!= EllipticK || Head[K1] =!= EllipticK, Return[$Failed]];
  531.     r = u /. {K -> 1, K1 -> 0};
  532.     s = u /. {K -> 0, K1 -> 1};
  533.     If[r K + s K1 =!= u, Return[$Failed]];
  534.     r *= 2; s *= -2 I;
  535.     If[!IntegerQ[r] || !IntegerQ[s], Return[$Failed]];
  536.     r = Mod[r, 8]; s = Mod[s, 8];    (* 1/8 period phase shifts *)
  537.  
  538.     (* argument reductions: AS, p572, 16.8 *)
  539.     (* shifts by periods and half periods *)
  540.     If[r > 3, r -= 4; If[Mod[p+q, 4] != 1, sign = -sign]];     (* u - 2K *)
  541.     If[s > 3, s -= 4; If[p+q != 3, sign = -sign]];         (* u - 2iK' *)
  542.     If[s == 0 && r == 3,                     (* 2K - u *)
  543.         r = 1; If[p==1 || q==1, sign = -sign]];
  544.     If[r == 0 && s == 3,                     (* 2iK' - u *)
  545.         s = 1; If[p==3 || q==3, sign = -sign]];
  546.     If[r + 2s > 6,                      (* u = 2(K+iK`)-u *)
  547.         s = 4 - s; r = 4 - r; If[p==2 || q==2, sign = -sign]];
  548.  
  549.     (* row in table 16.5, AS, p571.  row 10 is a special added case *)
  550.     If[r == 3 && s == 1, r = 10, r += 3s + 1];
  551.  
  552.     (* dispose of singularities in row 7 first *)
  553.     If[r == 7, If[p==3, Return[0]]; If[q==3, Return[ComplexInfinity]]];
  554.  
  555.     (* dispose of pn and nq *)
  556.     If[p == 3, Return[sign/$JacobiSpecialValues[[r, q+1]][m]]];
  557.     If[q == 3, Return[sign $JacobiSpecialValues[[r, p+1]][m]]];
  558.     sign $JacobiSpecialValues[[r,p+1]][m]/$JacobiSpecialValues[[r,q+1]][m]
  559.     ];
  560.  
  561. JacobiSN[u_, m_] :=
  562.     Module[{jac = SymbJacobi[0,3,u,m]},
  563.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  564. JacobiCN[u_, m_] :=
  565.     Module[{jac = SymbJacobi[1,3,u,m]},
  566.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  567. JacobiDN[u_, m_] :=
  568.     Module[{jac = SymbJacobi[2,3,u,m]},
  569.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  570.  
  571. JacobiCD[u_, m_] :=
  572.     Module[{jac = SymbJacobi[1,2,u,m]},
  573.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  574. JacobiSD[u_, m_] :=
  575.     Module[{jac = SymbJacobi[0,2,u,m]},
  576.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  577. JacobiND[u_, m_] :=
  578.     Module[{jac = SymbJacobi[3,2,u,m]},
  579.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  580.  
  581. JacobiDC[u_, m_] :=
  582.     Module[{jac = SymbJacobi[2,1,u,m]},
  583.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  584. JacobiNC[u_, m_] :=
  585.     Module[{jac = SymbJacobi[3,1,u,m]},
  586.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  587. JacobiSC[u_, m_] :=
  588.     Module[{jac = SymbJacobi[0,1,u,m]},
  589.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  590.  
  591. JacobiNS[u_, m_] :=
  592.     Module[{jac = SymbJacobi[3,0,u,m]},
  593.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  594. JacobiDS[u_, m_] :=
  595.     Module[{jac = SymbJacobi[2,0,u,m]},
  596.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  597. JacobiCS[u_, m_] :=
  598.     Module[{jac = SymbJacobi[1,0,u,m]},
  599.     jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
  600.  
  601.  (*-------------------- general numerical cases ---------------------------*)
  602.  
  603. $NJacobi[pp_, qq_, u_, m_] :=
  604.     Module[{p = pp, q = qq, K, K1, r, s, den, sign = 1, ru, mm},
  605.     mm = N[m, Precision[u]];
  606.     K = $EllipK[mm];
  607.     K1 = $EllipKp[mm];
  608.     den = Re[K] Re[K1] + Im[K] Im[K1];
  609.     If[den == 0,
  610.         r = Floor[2(Re[K] Re[u] + Im[K] Im[u])/Abs[K]^2];
  611.         s = 0,
  612.       (* else *)
  613.         r = Floor[2(Re[K1] Re[u] + Im[K1] Im[u])/den];
  614.         s = Floor[2(Re[K] Im[u] - Im[K] Re[u])/den]
  615.         ];
  616.     ru = u - (2 Floor[r/4] K + 2 I Floor[s/4] K1);
  617.     r = Mod[r, 8]; s = Mod[s, 8];
  618.     If[r > 3, If[Mod[p+q, 4] != 1, sign = -sign]];
  619.     If[s > 3, s -= 4; If[p+q != 3, sign = -sign]];
  620.     If[s > 1, ru = 2(K+I K1) - ru; If[p==2 || q==2, sign = -sign]];
  621.     sign $EllipticTheta[p, ru, mm]/$EllipticTheta[q, ru, mm]
  622.     ]
  623.  
  624. (* here are the twelve functions : one for each ordered pair of letters
  625.      from the set S,D,N,C *)
  626.  
  627. JacobiSN[u_?NumberQ,m_?NumberQ] :=
  628.     $NJacobi[0, 3, u, m] /; Precision[{u, m}] < Infinity
  629. JacobiCN[u_?NumberQ,m_?NumberQ] :=
  630.     $NJacobi[1, 3, u, m] /; Precision[{u, m}] < Infinity
  631. JacobiDN[u_?NumberQ,m_?NumberQ] :=
  632.     $NJacobi[2, 3, u, m] /; Precision[{u, m}] < Infinity
  633.  
  634. JacobiCD[u_?NumberQ,m_?NumberQ] :=
  635.     $NJacobi[1, 2, u, m] /; Precision[{u, m}] < Infinity
  636. JacobiSD[u_?NumberQ,m_?NumberQ] :=
  637.     $NJacobi[0, 2, u, m] /; Precision[{u, m}] < Infinity
  638. JacobiND[u_?NumberQ,m_?NumberQ] :=
  639.     $NJacobi[3, 2, u, m] /; Precision[{u, m}] < Infinity
  640.  
  641. JacobiDC[u_?NumberQ,m_?NumberQ] :=
  642.     $NJacobi[2, 1, u, m] /; Precision[{u, m}] < Infinity
  643. JacobiNC[u_?NumberQ,m_?NumberQ] :=
  644.     $NJacobi[3, 1, u, m] /; Precision[{u, m}] < Infinity
  645. JacobiSC[u_?NumberQ,m_?NumberQ] :=
  646.     $NJacobi[0, 1, u, m] /; Precision[{u, m}] < Infinity
  647.  
  648. JacobiNS[u_?NumberQ,m_?NumberQ] :=
  649.     $NJacobi[3, 0, u, m] /; Precision[{u, m}] < Infinity
  650. JacobiDS[u_?NumberQ,m_?NumberQ] :=
  651.     $NJacobi[2, 0, u, m] /; Precision[{u, m}] < Infinity
  652. JacobiCS[u_?NumberQ,m_?NumberQ] :=
  653.     $NJacobi[1, 0, u, m] /; Precision[{u, m}] < Infinity
  654.  
  655.   (* ------------------ Derivatives: AS, p574, 16.16 ---------------------- *)
  656.  
  657. JacobiSN/: Derivative[n_Integer,0][JacobiSN] :=
  658.     Module[{j, r},
  659.         r = Evaluate[
  660.         Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCN][#1, #2] *
  661.             Derivative[n-1-j, 0][JacobiDN][#1, #2], {j, 0, n-1}]]&;
  662.         r = Share[r];
  663.         AbortProtect[
  664.         Unprotect[JacobiSN];
  665.         Literal[Derivative[n,0][JacobiSN]] ^= r;
  666.         Protect[JacobiSN]
  667.         ];
  668.         r] /; n > 0
  669. JacobiCN/: Derivative[n_Integer,0][JacobiCN] :=
  670.     Module[{j, r},
  671.         r = Evaluate[
  672.         -Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSN][#1, #2] *
  673.             Derivative[n-1-j, 0][JacobiDN][#1, #2], {j, 0, n-1}]]&;
  674.         r = Share[r];
  675.         AbortProtect[
  676.         Unprotect[JacobiCN];
  677.         Literal[Derivative[n,0][JacobiCN]] ^= r;
  678.         Protect[JacobiCN]
  679.         ];
  680.         r] /; n > 0
  681. JacobiDN/: Derivative[n_Integer,0][JacobiDN] :=
  682.     Module[{j, r},
  683.         r = Evaluate[
  684.         -#2 Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSN][#1, #2] *
  685.             Derivative[n-1-j, 0][JacobiCN][#1, #2], {j, 0, n-1}]]&;
  686.         r = Share[r];
  687.         AbortProtect[
  688.         Unprotect[JacobiDN];
  689.         Literal[Derivative[n,0][JacobiDN]] ^= r;
  690.         Protect[JacobiDN]
  691.         ];
  692.         r] /; n > 0
  693.  
  694. JacobiCD/: Derivative[n_Integer,0][JacobiCD] :=
  695.     Module[{j, r},
  696.         r = Evaluate[
  697.         (#2-1) Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSD][#1, #2]*
  698.             Derivative[n-1-j, 0][JacobiND][#1, #2], {j, 0, n-1}]]&;
  699.         r = Share[r];
  700.         AbortProtect[
  701.         Unprotect[JacobiCD];
  702.         Literal[Derivative[n,0][JacobiCD]] ^= r;
  703.         Protect[JacobiCD]
  704.         ];
  705.         r] /; n > 0
  706. JacobiSD/: Derivative[n_Integer,0][JacobiSD] :=
  707.     Module[{j, r},
  708.         r = Evaluate[
  709.         Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCD][#1, #2] *
  710.             Derivative[n-1-j, 0][JacobiND][#1, #2], {j, 0, n-1}]]&;
  711.         r = Share[r];
  712.         AbortProtect[
  713.         Unprotect[JacobiSD];
  714.         Literal[Derivative[n,0][JacobiSD]] ^= r;
  715.         Protect[JacobiSD]
  716.         ];
  717.         r] /; n > 0
  718. JacobiND/: Derivative[n_Integer,0][JacobiND] :=
  719.     Module[{j, r},
  720.         r = Evaluate[
  721.         #2 Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSD][#1, #2] *
  722.             Derivative[n-1-j, 0][JacobiCD][#1, #2], {j, 0, n-1}]]&;
  723.         r = Share[r];
  724.         AbortProtect[
  725.         Unprotect[JacobiND];
  726.         Literal[Derivative[n,0][JacobiND]] ^= r;
  727.         Protect[JacobiND]
  728.         ];
  729.         r] /; n > 0
  730.  
  731. JacobiDC/: Derivative[n_Integer,0][JacobiDC] :=
  732.     Module[{j, r},
  733.         r = Evaluate[
  734.         (1-#2) Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSC][#1, #2] *
  735.             Derivative[n-1-j, 0][JacobiNC][#1, #2], {j, 0, n-1}]]&;
  736.         r = Share[r];
  737.         AbortProtect[
  738.         Unprotect[JacobiDC];
  739.         Literal[Derivative[n,0][JacobiDC]] ^= r;
  740.         Protect[JacobiDC]
  741.         ];
  742.         r] /; n > 0
  743. JacobiNC/: Derivative[n_Integer,0][JacobiNC] :=
  744.     Module[{j, r},
  745.         r = Evaluate[
  746.         Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSC][#1, #2] *
  747.             Derivative[n-1-j, 0][JacobiDC][#1, #2], {j, 0, n-1}]]&;
  748.         r = Share[r];
  749.         AbortProtect[
  750.         Unprotect[JacobiNC];
  751.         Literal[Derivative[n,0][JacobiNC]] ^= r;
  752.         Protect[JacobiNC]
  753.         ];
  754.         r] /; n > 0
  755. JacobiSC/: Derivative[n_Integer,0][JacobiSC] :=
  756.     Module[{j, r},
  757.         r = Evaluate[
  758.         Sum[Binomial[n-1, j] Derivative[j, 0][JacobiDC][#1, #2] *
  759.             Derivative[n-1-j, 0][JacobiNC][#1, #2], {j, 0, n-1}]]&;
  760.         r = Share[r];
  761.         AbortProtect[
  762.         Unprotect[JacobiSC];
  763.         Literal[Derivative[n,0][JacobiSC]] ^= r;
  764.         Protect[JacobiSC]
  765.         ];
  766.         r] /; n > 0
  767.  
  768. JacobiNS/: Derivative[n_Integer,0][JacobiNS] :=
  769.     Module[{j, r},
  770.         r = Evaluate[
  771.         -Sum[Binomial[n-1, j] Derivative[j, 0][JacobiDS][#1, #2] *
  772.             Derivative[n-1-j, 0][JacobiCS][#1, #2], {j, 0, n-1}]]&;
  773.         r = Share[r];
  774.         AbortProtect[
  775.         Unprotect[JacobiNS];
  776.         Literal[Derivative[n,0][JacobiNS]] ^= r;
  777.         Protect[JacobiNS]
  778.         ];
  779.         r] /; n > 0
  780. JacobiDS/: Derivative[n_Integer,0][JacobiDS] :=
  781.     Module[{j, r},
  782.         r = Evaluate[
  783.         -Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCS][#1, #2] *
  784.             Derivative[n-1-j, 0][JacobiNS][#1, #2], {j, 0, n-1}]]&;
  785.         r = Share[r];
  786.         AbortProtect[
  787.         Unprotect[JacobiDS];
  788.         Literal[Derivative[n,0][JacobiDS]] ^= r;
  789.         Protect[JacobiDS]
  790.         ];
  791.         r] /; n > 0
  792. JacobiCS/: Derivative[n_Integer,0][JacobiCS] :=
  793.     Module[{j, r},
  794.         r = Evaluate[
  795.         -Sum[Binomial[n-1, j] Derivative[j, 0][JacobiNS][#1, #2] *
  796.             Derivative[n-1-j, 0][JacobiDS][#1, #2], {j, 0, n-1}]]&;
  797.         r = Share[r];
  798.         AbortProtect[
  799.         Unprotect[JacobiCS];
  800.         Literal[Derivative[n,0][JacobiCS]] ^= r;
  801.         Protect[JacobiCS]
  802.         ];
  803.         r] /; n > 0
  804.  
  805.   (* --------------------Inverse Jacobi functions ----------------------- *)
  806.  
  807. InverseFunction[JacobiSN, 1, 2] = InverseJacobiSN
  808. InverseFunction[JacobiCN, 1, 2] = InverseJacobiCN
  809. InverseFunction[JacobiDN, 1, 2] = InverseJacobiDN
  810.  
  811. InverseFunction[JacobiCD, 1, 2] = InverseJacobiCD
  812. InverseFunction[JacobiSD, 1, 2] = InverseJacobiSD
  813. InverseFunction[JacobiND, 1, 2] = InverseJacobiND
  814.  
  815. InverseFunction[JacobiDC, 1, 2] = InverseJacobiDC
  816. InverseFunction[JacobiNC, 1, 2] = InverseJacobiNC
  817. InverseFunction[JacobiSC, 1, 2] = InverseJacobiSC
  818.  
  819. InverseFunction[JacobiNS, 1, 2] = InverseJacobiNS
  820. InverseFunction[JacobiDS, 1, 2] = InverseJacobiDS
  821. InverseFunction[JacobiCS, 1, 2] = InverseJacobiCS
  822.  
  823. InverseFunction[InverseJacobiSN, 1, 2] = JacobiSN
  824. InverseFunction[InverseJacobiCN, 1, 2] = JacobiCN
  825. InverseFunction[InverseJacobiDN, 1, 2] = JacobiDN
  826.  
  827. InverseFunction[InverseJacobiCD, 1, 2] = JacobiCD
  828. InverseFunction[InverseJacobiSD, 1, 2] = JacobiSD
  829. InverseFunction[InverseJacobiND, 1, 2] = JacobiND
  830.  
  831. InverseFunction[InverseJacobiDC, 1, 2] = JacobiDC
  832. InverseFunction[InverseJacobiNC, 1, 2] = JacobiNC
  833. InverseFunction[InverseJacobiSC, 1, 2] = JacobiSC
  834.  
  835. InverseFunction[InverseJacobiNS, 1, 2] = JacobiNS
  836. InverseFunction[InverseJacobiDS, 1, 2] = JacobiDS
  837. InverseFunction[InverseJacobiCS, 1, 2] = JacobiCS
  838.  
  839.  (* ------------ elementary special cases, AS, p571, 16.6 -------------- *)
  840.  
  841. InverseJacobiSN[u_, 0|0.] := ArcSin[u];    InverseJacobiSN[u_, 1|1.] := ArcTanh[u];
  842. InverseJacobiCN[u_, 0|0.] := ArcCos[u];    InverseJacobiCN[u_, 1|1.] := ArcSech[u];
  843. (* InverseJacobiDN[_, 0] has no inv *)    InverseJacobiDN[u_, 1|1.] := ArcSech[u];
  844.  
  845. InverseJacobiCD[u_, 0|0.] := ArcCos[u];    (* InverseJacobiCD[_, 1] has no inv *)
  846. InverseJacobiSD[u_, 0|0.] := ArcSin[u];    InverseJacobiSD[u_, 1|1.] := ArcSinh[u];
  847. (* InverseJacobiND[_, 0] has no inv *)    InverseJacobiND[u_, 1|1.] := ArcCosh[u];
  848.  
  849. InverseJacobiDC[u_, 0|0.] := ArcSec[u];    (* InverseJacobiDC[_, 1] has no inv *)
  850. InverseJacobiNC[u_, 0|0.] := ArcSec[u];    InverseJacobiNC[u_, 1|1.] := ArcCosh[u];
  851. InverseJacobiSC[u_, 0|0.] := ArcTan[u];    InverseJacobiSC[u_, 1|1.] := ArcSinh[u];
  852.  
  853. InverseJacobiNS[u_, 0|0.] := ArcCsc[u];    InverseJacobiNS[u_, 1|1.] := ArcCoth[u];
  854. InverseJacobiDS[u_, 0|0.] := ArcCsc[u];    InverseJacobiDS[u_, 1|1.] := ArcCsch[u];
  855. InverseJacobiCS[u_, 0|0.] := ArcCot[u];    InverseJacobiCS[u_, 1|1.] := ArcCsch[u];
  856.  
  857. InverseFunction::noinv = "The function `1` is not invertible."
  858.  
  859. $NoInverse[f_, m_] := (Message[InverseFunction::noinv, f[#, m]&]; False)
  860.  
  861. InverseJacobiDN[_, m_] := 0 /; (TrueQ[m == 0.0] && $NoInverse[JacobiDN, 0])
  862. InverseJacobiND[_, m_] := 0 /; (TrueQ[m == 0.0] && $NoInverse[JacobiND, 0])
  863.  
  864. InverseJacobiCD[_, m_] := 0 /; (TrueQ[m-1 == 0.0] && $NoInverse[JacobiCD, 1])
  865. InverseJacobiDC[_, m_] := 0 /; (TrueQ[m-1 == 0.0] && $NoInverse[JacobiDC, 1])
  866.  
  867.  (* ---------------- general numerical inverses --------------------------*)
  868.  
  869. $InverseJacobi[jac_,a_,b_] :=
  870.     Module[{j = -jac/2, s, t, A, B, CC},
  871.     If[j === 0,
  872.         0,
  873.     (* else *)
  874.         A = j^3;
  875.         B = a j^2 - 1;
  876.         CC = b j;
  877.         t = Sqrt[B^2-4 A CC];
  878.         If[(B == 0) || Re[t/B] <= 0,
  879.         t = (t-B)/(2 A),
  880.           (* else *)
  881.         t = (-2 CC)/(t+B)];
  882.         s = j t;
  883.         EllipticLog[{s,t},{a,b}]
  884.         ]]
  885.  
  886. (* the Jacobi functions don't REALLY have inverses *)
  887.  
  888. (*                    .        .
  889.                     .         .
  890.                     .       .
  891.                     N D N D
  892.                     S C S C
  893.    Here is the diagram of zeros:    N D N D
  894.                     S C S C . . .
  895.  *)
  896.  
  897. $InvJacArgRed[u_, m_, n_, o_, p_] :=
  898.     Module[{ru, r, s, d, K, IK1},
  899.     K = $EllipK[N[m, Precision[u]]];
  900.     IK1 = I $EllipKp[N[m, Precision[u]]];
  901.     ru = u - (n K + o IK1);
  902.     d = 2(Re[K] Im[IK1] - Im[K] Re[IK1]);
  903.     (* what do we do if K and IK1 are not independent? *)
  904.     r = Round[(Im[IK1] Re[ru] - Re[IK1] Im[ru])/d];
  905.     s = Round[(Re[K] Im[ru] - Im[K] Re[ru])/d];
  906.     ru = u - 2(r K + s IK1);
  907.     If[Switch[p, 1, OddQ[r], 2, OddQ[s], 3, Xor[OddQ[r], OddQ[s]]],
  908.         ru = 2(n K + o IK1) - ru];
  909.     ru
  910.     ]
  911.  
  912. InverseJacobiSN[v_?NumberQ,m_?NumberQ] :=
  913.     Module[{u},
  914.     u = 2 $InverseJacobi[v, 2 + 2m, (1 - m)^2];
  915.     $InvJacArgRed[u, m, 0, 0, 1]
  916.     ] /; Precision[{v, m}] < Infinity
  917. InverseJacobiCN[v_?NumberQ,m_?NumberQ] :=
  918.     Module[{u},
  919.     u = $EllipK[m] - 2 $InverseJacobi[v/Sqrt[1-m], 2 - 4m, 1 ];
  920.     $InvJacArgRed[u, m, 1, 0, 3]
  921.     ] /; Precision[{v, m}] < Infinity
  922. InverseJacobiDN[v_?NumberQ,m_?NumberQ] :=
  923.     Module[{u},
  924.     u = 2 $InverseJacobi[I/v, -4 + 2m, m^2] - I $EllipKp[m];
  925.     $InvJacArgRed[u, m, 1, 1, 2]
  926.     ] /; (Precision[{v, m}] < Infinity && m != 0.0)
  927.  
  928. InverseJacobiCD[v_?NumberQ,m_?NumberQ] :=
  929.     Module[{u},
  930.     u = $EllipK[m] - 2 $InverseJacobi[v, 2 + 2m, (1 - m)^2];
  931.     $InvJacArgRed[u, m, 1, 0, 1]
  932.     ] /; (Precision[{v, m}] < Infinity && (m-1) != 0.0)
  933. InverseJacobiSD[v_?NumberQ,m_?NumberQ] :=
  934.     Module[{u},
  935.     u = 2 $InverseJacobi[v, 2 - 4m, 1 ];
  936.     $InvJacArgRed[u, m, 0, 0, 3]
  937.     ] /; Precision[{v, m}] < Infinity
  938. InverseJacobiND[v_?NumberQ,m_?NumberQ] :=
  939.     Module[{u},
  940.     u = 2 $InverseJacobi[I v, -4 + 2m, m^2] - I $EllipKp[m];
  941.     $InvJacArgRed[u, m, 0, 1, 2]
  942.     ] /; (Precision[{v, m}] < Infinity && m != 0.0)
  943.  
  944. InverseJacobiDC[v_?NumberQ,m_?NumberQ] :=
  945.     Module[{u},
  946.     u = $EllipK[m] - 2 $InverseJacobi[1/v, 2 + 2m, (1 - m)^2];
  947.     $InvJacArgRed[u, m, 1, 1, 1]
  948.     ] /; (Precision[{v, m}] < Infinity && (m-1) != 0.0)
  949. InverseJacobiNC[v_?NumberQ,m_?NumberQ] :=
  950.     Module[{u},
  951.     u = $EllipK[m] - 2 $InverseJacobi[1/(v Sqrt[1-m]), 2 - 4m, 1 ];
  952.     $InvJacArgRed[u, m, 0, 1, 3]
  953.     ] /; Precision[{v, m}] < Infinity
  954. InverseJacobiSC[v_?NumberQ,m_?NumberQ] :=
  955.     Module[{u},
  956.     u = 2 $InverseJacobi[v, -4 + 2m, m^2];
  957.     $InvJacArgRed[u, m, 0, 0, 2]
  958.     ] /; Precision[{v, m}] < Infinity
  959.  
  960. InverseJacobiNS[v_?NumberQ,m_?NumberQ] :=
  961.     Module[{u},
  962.     u = 2 $InverseJacobi[1/v, 2 + 2m, (1 - m)^2];;
  963.     $InvJacArgRed[u, m, 0, 1, 1]
  964.     ] /; Precision[{v, m}] < Infinity
  965. InverseJacobiDS[v_?NumberQ,m_?NumberQ] :=
  966.     Module[{u},
  967.     u = 2 $InverseJacobi[1/v, 2 - 4m, 1 ];
  968.     $InvJacArgRed[u, m, 1, 1, 3]
  969.     ] /; Precision[{v, m}] < Infinity
  970. InverseJacobiCS[v_?NumberQ,m_?NumberQ] :=
  971.     Module[{u},
  972.     u = 2 $InverseJacobi[1/v, -4 + 2m, m^2];
  973.     $InvJacArgRed[u, m, 1, 0, 2]
  974.     ] /; Precision[{v, m}] < Infinity
  975.  
  976.  (* ------------ derivatives of the inverse functions ----------------- *)
  977.  
  978. InverseJacobiSN/: Derivative[1,0][InverseJacobiSN] :=
  979.         (1/Derivative[1,0][JacobiSN][InverseJacobiSN[#1, #2], #2])&;
  980. InverseJacobiCN/: Derivative[1,0][InverseJacobiCN] :=
  981.         (1/Derivative[1,0][JacobiCN][InverseJacobiCN[#1, #2], #2])&;
  982. InverseJacobiDN/: Derivative[1,0][InverseJacobiDN] :=
  983.         (1/Derivative[1,0][JacobiDN][InverseJacobiDN[#1, #2], #2])&;
  984.  
  985. InverseJacobiCD/: Derivative[1,0][InverseJacobiCD] :=
  986.         (1/Derivative[1,0][JacobiCD][InverseJacobiCD[#1, #2], #2])&;
  987. InverseJacobiSD/: Derivative[1,0][InverseJacobiSD] :=
  988.         (1/Derivative[1,0][JacobiSD][InverseJacobiSD[#1, #2], #2])&;
  989. InverseJacobiND/: Derivative[1,0][InverseJacobiND] :=
  990.         (1/Derivative[1,0][JacobiND][InverseJacobiND[#1, #2], #2])&;
  991.  
  992. InverseJacobiDC/: Derivative[1,0][InverseJacobiDC] :=
  993.         (1/Derivative[1,0][JacobiDC][InverseJacobiDC[#1, #2], #2])&;
  994. InverseJacobiNC/: Derivative[1,0][InverseJacobiNC] :=
  995.         (1/Derivative[1,0][JacobiNC][InverseJacobiNC[#1, #2], #2])&;
  996. InverseJacobiSC/: Derivative[1,0][InverseJacobiSC] :=
  997.         (1/Derivative[1,0][JacobiSC][InverseJacobiSC[#1, #2], #2])&;
  998.  
  999. InverseJacobiNS/: Derivative[1,0][InverseJacobiNS] :=
  1000.         (1/Derivative[1,0][JacobiNS][InverseJacobiNS[#1, #2], #2])&;
  1001. InverseJacobiDS/: Derivative[1,0][InverseJacobiDS] :=
  1002.         (1/Derivative[1,0][JacobiDS][InverseJacobiDS[#1, #2], #2])&;
  1003. InverseJacobiCS/: Derivative[1,0][InverseJacobiCS] :=
  1004.         (1/Derivative[1,0][JacobiCS][InverseJacobiCS[#1, #2], #2])&;
  1005.  
  1006.  (* ----------------------- JacobiAmplitude ------------------------ *)
  1007.  
  1008. JacobiAmplitude[u_, 0] := u;
  1009. JacobiAmplitude[u_, 1] := 2 ArcTan[Exp[u]] - Pi/2;
  1010. JacobiAmplitude[u_?NumberQ, m_?NumberQ] :=
  1011.     Module[{am = ArcSin[JacobiSN[u,m]], pi, s, k2, kp, km},
  1012.         k2 = 2 $EllipK[m];
  1013.         s = EllipticF[am, m];
  1014.         pi = N[Pi, Accuracy[am]+2];
  1015.         km = Re[(u - s)/k2];
  1016.         kp = Re[(u + s)/k2];
  1017.         If[Abs[km - Round[km]] > Abs[kp - Round[kp]],
  1018.             Round[kp] pi - am,
  1019.            (* else *)
  1020.             Round[km] pi + am]
  1021.         ];
  1022.  
  1023. JacobiAmplitude/: Derivative[1,0][JacobiAmplitude] = JacobiDN
  1024.  
  1025. InverseFunction[JacobiAmplitude, 1, 2] = EllipticF
  1026. InverseFunction[EllipticF, 1, 2] = JacobiAmplitude 
  1027.  
  1028. End[] (* Private *)
  1029.  
  1030. SetAttributes[JacobiSN, ReadProtected];
  1031. SetAttributes[JacobiSD, ReadProtected];
  1032. SetAttributes[JacobiSC, ReadProtected];
  1033. SetAttributes[JacobiCS, ReadProtected];
  1034. SetAttributes[JacobiCN, ReadProtected];
  1035. SetAttributes[JacobiCD, ReadProtected];
  1036. SetAttributes[JacobiNS, ReadProtected];
  1037. SetAttributes[JacobiNC, ReadProtected];
  1038. SetAttributes[JacobiND, ReadProtected];
  1039. SetAttributes[JacobiDN, ReadProtected];
  1040. SetAttributes[JacobiDS, ReadProtected];
  1041. SetAttributes[JacobiDC, ReadProtected];
  1042. SetAttributes[InverseJacobiSN, ReadProtected];
  1043. SetAttributes[InverseJacobiSD, ReadProtected];
  1044. SetAttributes[InverseJacobiSC, ReadProtected];
  1045. SetAttributes[InverseJacobiCS, ReadProtected];
  1046. SetAttributes[InverseJacobiCN, ReadProtected];
  1047. SetAttributes[InverseJacobiCD, ReadProtected];
  1048. SetAttributes[InverseJacobiNS, ReadProtected];
  1049. SetAttributes[InverseJacobiNC, ReadProtected];
  1050. SetAttributes[InverseJacobiND, ReadProtected];
  1051. SetAttributes[InverseJacobiDN, ReadProtected];
  1052. SetAttributes[InverseJacobiDS, ReadProtected];
  1053. SetAttributes[InverseJacobiDC, ReadProtected];
  1054. SetAttributes[JacobiAmplitude, ReadProtected];
  1055. SetAttributes[EllipticNomeQ, ReadProtected];
  1056. SetAttributes[EllipticThetaS, ReadProtected];
  1057. SetAttributes[EllipticThetaC, ReadProtected];
  1058. SetAttributes[EllipticThetaN, ReadProtected];
  1059. SetAttributes[EllipticThetaD, ReadProtected];
  1060. SetAttributes[WeierstrassP, ReadProtected];
  1061. SetAttributes[WeierstrassPPrime, ReadProtected];
  1062. SetAttributes[InverseWeierstrassP, ReadProtected];
  1063. SetAttributes[EllipticExp, ReadProtected];
  1064. SetAttributes[EllipticLog, ReadProtected];
  1065.  
  1066. Protect[JacobiSN, JacobiSD, JacobiSC, JacobiCS, JacobiCN, JacobiCD,
  1067.     JacobiNS, JacobiNC, JacobiND, JacobiDN, JacobiDS, JacobiDC]
  1068. Protect[InverseJacobiSN, InverseJacobiSD, InverseJacobiSC, InverseJacobiCS,
  1069.     InverseJacobiCN, InverseJacobiCD, InverseJacobiNS, InverseJacobiNC,
  1070.     InverseJacobiDS, InverseJacobiDC, InverseJacobiND, InverseJacobiDN]
  1071. Protect[JacobiAmplitude, EllipticNomeQ, EllipticThetaS, EllipticThetaC,
  1072.     EllipticThetaN, EllipticThetaD, WeierstrassP, WeierstrassPPrime,
  1073.     InverseWeierstrassP]
  1074. Protect[EllipticExp, EllipticLog, InverseFunction]
  1075.  
  1076. End[] (* Package *)
  1077.  
  1078. (* :Tests: *)
  1079.  
  1080. (* :Note:
  1081.     To run the tests just type `check', `complexcheck', `complexcheck2',
  1082.     or `symbcheckall'.  The test symbcheck[n] checks the symbolic values
  1083.     at (n K + I K' Range[0,15])/2.
  1084. *)
  1085.  
  1086. (*
  1087. see[result_,parms___] := If[Chop[result] === 0, Okay, Error[result,parms]]
  1088.  
  1089. see2[result_,parms___] :=
  1090.     If[Chop[result, .1^5] === 0, Okay, Error[result,parms]]
  1091.  
  1092. check1[ x_, m_ ] := see2 [ x - JacobiAmplitude[EllipticF[x, m], m]]
  1093. check2[u_,m_] := see[ - JacobiDN[u,m]^2 + 1 - m + m JacobiCN[u,m]^2 ]
  1094. check3[u_,m_] := see[ - m JacobiCN[u,m]^2 - m JacobiSN[u,m]^2 + m ]
  1095. check4[u_,m_] := see[ - (1-m) JacobiND[u,m]^2 + (1-m) + m(1-m)JacobiSD[u,m]^2 ]
  1096. check5[u_,m_] := see[ -m(1-m)JacobiSD[u,m]^2 - m JacobiCD[u,m]^2 + m ]
  1097. checkP[u_, g2_, g3_] :=
  1098.     Module[{p, pp, uu, err},
  1099.     p = WeierstrassP[u, g2, g3];
  1100.     pp = WeierstrassPPrime[u, g2, g3];
  1101.     uu = InverseWeierstrassP[{p, pp}, g2, g3];
  1102.     err = Abs[p-WeierstrassP[uu, g2, g3]];
  1103.     err += Abs[pp-WeierstrassPPrime[uu, g2, g3]];
  1104.     see[err, P, u, g2, g3]
  1105.     ];
  1106.  
  1107. checkinv[m_] :=
  1108.   Module[{K, K1, u, r},
  1109.      K = Elliptic`Private`$EllipK[m];
  1110.      K1 = Elliptic`Private`$EllipKp[m];
  1111.      r = Table[, {12}];
  1112.      u = Random[Real, {-1,1}] K + Random[Real, {-1,1}] I K1;
  1113.      r[[1]] = see[ u - InverseJacobiSN[JacobiSN[u,m],m] , u , m, SN ];
  1114.      r[[2]] = see[ u - InverseJacobiSD[JacobiSD[u,m],m] , u , m, SD ];
  1115.      r[[3]] = see[ u - InverseJacobiSC[JacobiSC[u,m],m] , u , m, SC ];
  1116.      u += K;
  1117.      r[[4]] = see[ u - InverseJacobiCD[JacobiCD[u,m],m] , u , m, CD ];
  1118.      r[[5]] = see[ u - InverseJacobiCN[JacobiCN[u,m],m] , u , m, CN ];
  1119.      r[[6]] = see[ u - InverseJacobiCS[JacobiCS[u,m],m] , u , m, CS ];
  1120.      u += I K1;
  1121.      r[[7]] = see[ u - InverseJacobiDC[JacobiDC[u,m],m] , u , m, DC ];
  1122.      r[[8]] = see[ u - InverseJacobiDS[JacobiDS[u,m],m] , u , m, DS ];
  1123.      r[[9]] = see[ u - InverseJacobiDN[JacobiDN[u,m],m] , u , m, DN ];
  1124.      u -= K;
  1125.      r[[10]] = see[ u - InverseJacobiNS[JacobiNS[u,m],m] , u , m, NS ];
  1126.      r[[11]] = see[ u - InverseJacobiNC[JacobiNC[u,m],m] , u , m, NC ];
  1127.      r[[12]] = see[ u - InverseJacobiND[JacobiND[u,m],m] , u , m, ND ];
  1128.      r
  1129.      ]
  1130.  
  1131. checkinv2[u_,m_] := {
  1132.      see[ u - JacobiSD[InverseJacobiSD[u,m],m] , u , m, SD ],
  1133.      see[ u - JacobiSN[InverseJacobiSN[u,m],m] , u , m, SN ],
  1134.      see[ u - JacobiSC[InverseJacobiSC[u,m],m] , u , m, SC ],
  1135.      see[ u - JacobiCS[InverseJacobiCS[u,m],m] , u , m, CS ],
  1136.      see[ u - JacobiCD[InverseJacobiCD[u,m],m] , u , m, CD ],
  1137.      see[ u - JacobiCN[InverseJacobiCN[u,m],m] , u , m, CN ],
  1138.      see[ u - JacobiDS[InverseJacobiDS[u,m],m] , u , m, DS ],
  1139.      see[ u - JacobiDC[InverseJacobiDC[u,m],m] , u , m, DC ],
  1140.      see[ u - JacobiNS[InverseJacobiNS[u,m],m] , u , m, NS ],
  1141.      see[ u - JacobiND[InverseJacobiND[u,m],m] , u , m, ND ],
  1142.      see[ u - JacobiDN[InverseJacobiDN[u,m],m] , u , m, DN ],
  1143.      see[ u - JacobiNC[InverseJacobiNC[u,m],m] , u , m, NC ]
  1144.      }
  1145.  
  1146. checkder[u_, m_] := {
  1147.      see2[(JacobiSD[u + .00001, m] - JacobiSD[u - .00001, m])/
  1148.         (.00002 Derivative[1, 0][JacobiSD][u, m]) - 1, u, m, SD],
  1149.      see2[(JacobiSN[u + .00001, m] - JacobiSN[u - .00001, m])/
  1150.         (.00002 Derivative[1, 0][JacobiSN][u, m]) - 1, u, m, SN],
  1151.      see2[(JacobiSC[u + .00001, m] - JacobiSC[u - .00001, m])/
  1152.         (.00002 Derivative[1, 0][JacobiSC][u, m]) - 1, u, m, SC],
  1153.      see2[(JacobiCS[u + .00001, m] - JacobiCS[u - .00001, m])/
  1154.         (.00002 Derivative[1, 0][JacobiCS][u, m]) - 1, u, m, CS],
  1155.      see2[(JacobiCD[u + .00001, m] - JacobiCD[u - .00001, m])/
  1156.         (.00002 Derivative[1, 0][JacobiCD][u, m]) - 1, u, m, CD],
  1157.      see2[(JacobiCN[u + .00001, m] - JacobiCN[u - .00001, m])/
  1158.         (.00002 Derivative[1, 0][JacobiCN][u, m]) - 1, u, m, CN],
  1159.      see2[(JacobiDS[u + .00001, m] - JacobiDS[u - .00001, m])/
  1160.         (.00002 Derivative[1, 0][JacobiDS][u, m]) - 1, u, m, DS],
  1161.      see2[(JacobiDC[u + .00001, m] - JacobiDC[u - .00001, m])/
  1162.         (.00002 Derivative[1, 0][JacobiDC][u, m]) - 1, u, m, DC],
  1163.      see2[(JacobiNS[u + .00001, m] - JacobiNS[u - .00001, m])/
  1164.         (.00002 Derivative[1, 0][JacobiNS][u, m]) - 1, u, m, NS],
  1165.      see2[(JacobiND[u + .00001, m] - JacobiND[u - .00001, m])/
  1166.         (.00002 Derivative[1, 0][JacobiND][u, m]) - 1, u, m, ND],
  1167.      see2[(JacobiDN[u + .00001, m] - JacobiDN[u - .00001, m])/
  1168.         (.00002 Derivative[1, 0][JacobiDN][u, m]) - 1, u, m, DN],
  1169.      see2[(JacobiNC[u + .00001, m] - JacobiNC[u - .00001, m])/
  1170.         (.00002 Derivative[1, 0][JacobiNC][u, m]) - 1, u, m, NC]}
  1171.  
  1172. checkPder[u_, g2_, g3_] := {
  1173.      see2[(WeierstrassP[u+.00001,g2,g3]-WeierstrassP[u-.00001,g2,g3])/
  1174.     (.00002 Derivative[1, 0, 0][WeierstrassP][u,g2,g3]) - 1,u,g2,g3,P],
  1175.      see2[(WeierstrassPPrime[u+.00001,g2,g3]-WeierstrassPPrime[u-.00001,g2,g3])/
  1176.     (.00002 Derivative[1,0,0][WeierstrassPPrime][u,g2,g3])-1,u,g2,g3,PP]}
  1177.  
  1178. checkderinv[u_, m_] := {
  1179.      see2[(InverseJacobiSD[u + .00001, m] - InverseJacobiSD[u - .00001, m])/
  1180.     (.00002 Derivative[1, 0][InverseJacobiSD][u, m]) - 1, u, m, SD],
  1181.      see2[(InverseJacobiSN[u + .00001, m] - InverseJacobiSN[u - .00001, m])/
  1182.     (.00002 Derivative[1, 0][InverseJacobiSN][u, m]) - 1, u, m, SN],
  1183.      see2[(InverseJacobiSC[u + .00001, m] - InverseJacobiSC[u - .00001, m])/
  1184.     (.00002 Derivative[1, 0][InverseJacobiSC][u, m]) - 1, u, m, SC],
  1185.      see2[(InverseJacobiCS[u + .00001, m] - InverseJacobiCS[u - .00001, m])/
  1186.     (.00002 Derivative[1, 0][InverseJacobiCS][u, m]) - 1, u, m, CS],
  1187.      see2[(InverseJacobiCD[u + .00001, m] - InverseJacobiCD[u - .00001, m])/
  1188.     (.00002 Derivative[1, 0][InverseJacobiCD][u, m]) - 1, u, m, CD],
  1189.      see2[(InverseJacobiCN[u + .00001, m] - InverseJacobiCN[u - .00001, m])/
  1190.     (.00002 Derivative[1, 0][InverseJacobiCN][u, m]) - 1, u, m, CN],
  1191.      see2[(InverseJacobiDS[u + .00001, m] - InverseJacobiDS[u - .00001, m])/
  1192.     (.00002 Derivative[1, 0][InverseJacobiDS][u, m]) - 1, u, m, DS],
  1193.      see2[(InverseJacobiDC[u + .00001, m] - InverseJacobiDC[u - .00001, m])/
  1194.     (.00002 Derivative[1, 0][InverseJacobiDC][u, m]) - 1, u, m, DC],
  1195.      see2[(InverseJacobiNS[u + .00001, m] - InverseJacobiNS[u - .00001, m])/
  1196.     (.00002 Derivative[1, 0][InverseJacobiNS][u, m]) - 1, u, m, NS],
  1197.      see2[(InverseJacobiND[u + .00001, m] - InverseJacobiND[u - .00001, m])/
  1198.     (.00002 Derivative[1, 0][InverseJacobiND][u, m]) - 1, u, m, ND],
  1199.      see2[(InverseJacobiDN[u + .00001, m] - InverseJacobiDN[u - .00001, m])/
  1200.     (.00002 Derivative[1, 0][InverseJacobiDN][u, m]) - 1, u, m, DN],
  1201.      see2[(InverseJacobiNC[u + .00001, m] - InverseJacobiNC[u - .00001, m])/
  1202.     (.00002 Derivative[1, 0][InverseJacobiNC][u, m]) - 1, u, m, NC]}
  1203.  
  1204. check := {
  1205.     check1[Random[Real, {-25, 25}],Random[]] ,
  1206.     check2[Random[Real, {-25, 25}],Random[]] ,
  1207.     check3[Random[Real, {-25, 25}],Random[]] ,
  1208.     check4[Random[Real, {-25, 25}],Random[]] ,
  1209.     check5[Random[Real, {-25, 25}],Random[]] ,
  1210.     checkP[Random[Real, {-5, 5}],
  1211.         Random[Real, {-10, 10}], Random[Real, {-10, 10}]] ,
  1212.     checkinv[Random[]],
  1213.     checkinv2[Random[],Random[]],
  1214.     checkder[Random[Real, {-25, 25}],Random[]],
  1215.     checkPder[Random[Real, {-5, 5}],
  1216.         Random[Real, {-5, 5}],Random[Real, {-5, 5}]]
  1217.     }
  1218.  
  1219. complexcheck := {
  1220.     check1[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
  1221.     check2[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
  1222.     check3[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
  1223.     check4[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
  1224.     check5[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
  1225.     checkP[Random[Complex, {-5-5I, 5+5I}],
  1226.         Random[Real, {-10, 10}], Random[Real, {-10, 10}]] ,
  1227.     checkinv[Random[]],
  1228.     checkinv2[Random[Complex, {-25-25I, 25+25I}],Random[]],
  1229.     checkder[Random[Complex, {-25-25I, 25+25I}],Random[]],
  1230.     checkderinv[Random[Complex, {-25-25I, 25+25I}],Random[]],
  1231.     checkPder[Random[Complex, {-5-5I,5+5I}],
  1232.         Random[Real, {-5, 5}], Random[Real, {-5, 5}]]
  1233.     }
  1234.  
  1235. complexcheck2 := {
  1236.     check1[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
  1237.     check2[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
  1238.     check3[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
  1239.     check4[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
  1240.     check5[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
  1241.     (*  Weierstrass doesn't support complex invariants.
  1242.     checkP[Random[Complex, {-5-5I, 5+5I}],
  1243.         Random[Complex, {-10-10I, 10+10I}],
  1244.         Random[Complex, {-10-10I, 10+10I}]] ,
  1245.     *)
  1246.     checkinv[Random[Complex]],
  1247.     checkinv2[Random[Complex, {-25-25I, 25+25I}],Random[Complex]],
  1248.     checkder[Random[Complex, {-25-25I, 25+25I}],Random[Complex]],
  1249.     checkderinv[Random[Complex, {-25-25I, 25+25I}],Random[Complex]]
  1250.     }
  1251.  
  1252. argtab[nK_, K1_] := Table[nK + I K1 j/2, {j, 0, 15}];
  1253. compare[{jac_, njac_}] :=
  1254.     Module[{},
  1255.     If[Head[jac] === DirectedInfinity && Abs[njac] > 10.^5, Return[0]];
  1256.     If[jac ==  0,
  1257.         If[Abs[njac] < 0.1^5, Return[0], Return[{jac, njac}]]];
  1258.     If[NumberQ[jac] && NumberQ[njac], Return[Chop[njac/jac - 1, .1^5]]];
  1259.     {jac, njac}];
  1260.  
  1261. scheck[f_, n_] :=
  1262.     Module[{m, nm = Random[ ], jac, njac, err},
  1263.     jac = Map[f[#, m]&, argtab[EllipticK[m] n/2, EllipticK[1-m]]];
  1264.     njac = Map[f[#, nm]&,
  1265.         .1^13 + argtab[EllipticK[nm] n/2, EllipticK[1-nm]]];
  1266.     jac = jac /. EllipticK[m] -> EllipticK[nm];
  1267.     jac = jac /. EllipticK[1-m] -> EllipticK[1-nm];
  1268.     jac = N[(jac /. m -> nm)];
  1269.     jac = Transpose[{jac, njac}];
  1270.     err = Map[compare, jac];
  1271.     Print[err];
  1272.     err]
  1273.  
  1274. symbcheck[n_] := {
  1275.     (* check the symbolic values at (n K + I K' Range[0,15])/2 *)
  1276.     {scheck[JacobiSN, n], SN},
  1277.     {scheck[JacobiCN, n], CN},
  1278.     {scheck[JacobiDN, n], DN},
  1279.     {scheck[JacobiCD, n], CD},
  1280.     {scheck[JacobiSD, n], SD},
  1281.     {scheck[JacobiND, n], ND},
  1282.     {scheck[JacobiDC, n], DC},
  1283.     {scheck[JacobiNC, n], NC},
  1284.     {scheck[JacobiSC, n], SC},
  1285.     {scheck[JacobiNS, n], NS},
  1286.     {scheck[JacobiDS, n], DS},
  1287.     {scheck[JacobiCS, n], CS}
  1288.     }
  1289.  
  1290. symbcheckall :=
  1291.     (* check the symbolic values at (n K + m I K')/2, 0 <= m, n <= 15 *)
  1292.     Table[Print[n," * K / 2"]; symbcheck[n], {n, 0, 15}]
  1293.  
  1294. *)
  1295.