
(* :Name: StartUp`Elliptic` *)

(* :Title: Elliptic Functions *)

(* :Author: Daniel R. Grayson, Jerry B. Keiper *)

(* :Summary:
	Extends the coverage of elliptic functions from that in the kernel
	to include the Jacobi elliptic functions and their inverses, the
	Neville theta functions, the amplitude function am, the nome
	function q, and the Weierstrass P function and its derivative.
*)

(* :Context: System` *)

(* :Package Version: 2.0 *)

(* :Copyright: Copyright 1990  Wolfram Research, Inc.

	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Extensively revised by Jerry B. Keiper, November 1990.
	Version 1.0 by Daniel R. Grayson, 1988.
*)

(* :Keywords:
	elliptic functions, theta functions, Weierstrass P function
*)

(* :Source:
	Abramowitz and Stegun: "Handbook of Mathematical Functions",
		Dover Publications, New York. chap. 16, 17.
	Whittaker and Watson: "A Course of Modern Analysis",
		Cambridge University Press.  chap. 20, 21, 22.
*)

(* :Warning:
	Adds rules to EllipticExp[ ], EllipticLog[ ], and InverseFunction[ ]
*)

(* :Mathematica Version: 2.0 *)

(* :Limitation:
	Affected by inexactness in input numbers especially near
	branch cuts in the inverse functions.

	Ranges of certain inverse Jacobi elliptic functions differ from
	the ranges of the corresponding limiting case inverse trigonometric
	(both circular and hyperbolic) functions.

	This package is not complete; there are many more features which
	could be added.  The theory of elliptic functions is too rich to
	be able to cover everything without a great deal of work.    The
	Weierstrass functions are much less complete than are the Jacobi
	functions.

	Weierstrass functions do not work with complex invariants g2 and g3.
*)

(* :Discussion:
	This package defines the mathematical functions JacobiSN, JacobiCN,
	JacobiDN, JacobiCD, JacobiSD, JacobiND, JacobiDC, JacobiNC, JacobiSC,
	JacobiNS, JacobiDS, JacobiCS, InverseJacobiSN, InverseJacobiCN,
	InverseJacobiDN, InverseJacobiCD, InverseJacobiSD, InverseJacobiND,
	InverseJacobiDC, InverseJacobiNC, InverseJacobiSC, InverseJacobiNS,
	InverseJacobiDS, InverseJacobiCS, JacobiAmplitude, EllipticNomeQ, 
	EllipticThetaS, EllipticThetaC, EllipticThetaD, EllipticThetaN,
	WeierstrassP, WeierstrassPPrime, and InverseWeierstrassP as well
	as the derivatives of most of these functions.

	The Jacobi elliptic functions are evaluated in terms of the Neville
	theta functions after argument reductions on the first argument.
	The inverse Jacobi elliptic functions are evaluated in terms
	EllipticLog and then the result is shifted to the appropriate
	parallelogram.  JacobiAmplitude is evaluated as the ArcSin of
	JacobiSN and then the result is shifted to the appropriate value.
	EllipticNomeQ is evaluated in terms of EllipticK[m] and EllipticK[1-m].
	The Neville theta functions are evaluated in terms of the Jacobi
	theta functions.  WeierstrassP and WeierstrassPPrime are evaluated
	in terms of EllipticExp and InverseWeierstrassP is evaluated in
	terms of EllipticLog.  The result of InverseWeierstrassP is not 
	shifted to any particular period parallelogram since the basic
	periods are not calulated.
*)

Needs["System`", "StartUp`Attributes.m"]

Begin["System`"]

Unprotect[ JacobiSN, JacobiSD, JacobiSC,
           JacobiCS, JacobiCN, JacobiCD,
           JacobiNC, JacobiND, JacobiNS,
           JacobiDN, JacobiDS, JacobiDC,
           JacobiAmplitude,
           InverseJacobiSN, InverseJacobiSD, InverseJacobiSC,
           InverseJacobiCS, InverseJacobiCN, InverseJacobiCD,
           InverseJacobiNS, InverseJacobiNC, InverseJacobiND,
           InverseJacobiDS, InverseJacobiDC, InverseJacobiDN,
           EllipticNomeQ, EllipticThetaS, EllipticThetaC,
           EllipticThetaN, EllipticThetaD,
           WeierstrassP, WeierstrassPPrime,
           InverseWeierstrassP,
           EllipticExp, EllipticLog]

Map[ Clear,
          {JacobiSN, JacobiSD, JacobiSC,
           JacobiCS, JacobiCN, JacobiCD,
           JacobiNC, JacobiNS, JacobiND,
           JacobiDN, JacobiDS, JacobiDC,
           JacobiAmplitude,
           InverseJacobiSN, InverseJacobiSD, InverseJacobiSC,
           InverseJacobiCS, InverseJacobiCN, InverseJacobiCD,
           InverseJacobiNS, InverseJacobiNC, InverseJacobiND,
           InverseJacobiDS, InverseJacobiDC, InverseJacobiDN,
           EllipticNomeQ, EllipticThetaS, EllipticThetaC,
           EllipticThetaN, EllipticThetaD,
           WeierstrassP, WeierstrassPPrime,
           InverseWeierstrassP,
           EllipticExp, EllipticLog}]

(* ---------------------- Jacobi elliptic functions ---------------------- *)

JacobiSN::usage =
"JacobiSN[u, m] gives the Jacobi elliptic function sn at u for the parameter m."
JacobiSD::usage =
"JacobiSD[u, m] gives the Jacobi elliptic function sd at u for the parameter m."
JacobiSC::usage =
"JacobiSC[u, m] gives the Jacobi elliptic function sc at u for the parameter m."
JacobiCS::usage =
"JacobiCS[u, m] gives the Jacobi elliptic function cs at u for the parameter m."
JacobiCN::usage =
"JacobiCN[u, m] gives the Jacobi elliptic function cn at u for the parameter m."
JacobiCD::usage =
"JacobiCD[u, m] gives the Jacobi elliptic function cd at u for the parameter m."
JacobiNS::usage =
"JacobiNS[u, m] gives the Jacobi elliptic function ns at u for the parameter m."
JacobiNC::usage =
"JacobiNC[u, m] gives the Jacobi elliptic function nc at u for the parameter m."
JacobiND::usage =
"JacobiND[u, m] gives the Jacobi elliptic function nd at u for the parameter m."
JacobiDN::usage =
"JacobiDN[u, m] gives the Jacobi elliptic function dn at u for the parameter m."
JacobiDS::usage =
"JacobiDS[u, m] gives the Jacobi elliptic function ds at u for the parameter m."
JacobiDC::usage =
"JacobiDC[u, m] gives the Jacobi elliptic function dc at u for the parameter m."

JacobiAmplitude::usage =
"JacobiAmplitude[u, m] gives the amplitude for Jacobi elliptic functions."

InverseJacobiSN::usage =
"InverseJacobiSN[v, m] gives the inverse of the Jacobi elliptic function sn at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at 0."
InverseJacobiSD::usage =
"InverseJacobiSD[v, m] gives the inverse of the Jacobi elliptic function sd at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at 0."
InverseJacobiSC::usage =
"InverseJacobiSC[v, m] gives the inverse of the Jacobi elliptic function sc at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at 0."
InverseJacobiCS::usage =
"InverseJacobiCS[v, m] gives the inverse of the Jacobi elliptic function cs at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K."
InverseJacobiCN::usage =
"InverseJacobiCN[v, m] gives the inverse of the Jacobi elliptic function cn at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K."
InverseJacobiCD::usage =
"InverseJacobiCD[v, m] gives the inverse of the Jacobi elliptic function cd at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K."
InverseJacobiNS::usage =
"InverseJacobiNS[v, m] gives the inverse of the Jacobi elliptic function ns at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at I K'."
InverseJacobiNC::usage =
"InverseJacobiNC[v, m] gives the inverse of the Jacobi elliptic function nc at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at I K'."
InverseJacobiDS::usage =
"InverseJacobiDS[v, m] gives the inverse of the Jacobi elliptic function ds at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K + I K'."
InverseJacobiDC::usage =
"InverseJacobiDC[v, m] gives the inverse of the Jacobi elliptic function dc at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K + I K'."
InverseJacobiND::usage =
"InverseJacobiND[v, m] gives the inverse of the Jacobi elliptic function nd at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at I K'."
InverseJacobiDN::usage =
"InverseJacobiDN[v, m] gives the inverse of the Jacobi elliptic function dn at
v for the parameter m.  The result will be in the parallelogram with sides 2K
and 2I K' centered at K + I K'."

EllipticNomeQ::usage =
"EllipticNomeQ[m] gives the nome q == Exp[-Pi EllipticK[1-m]/EllipticK[m]]."

EllipticThetaS::usage =
"EllipticThetaS[u, m] gives Neville's elliptic theta function theta_s(u, m)."
EllipticThetaC::usage =
"EllipticThetaC[u, m] gives Neville's elliptic theta function theta_c(u, m)."
EllipticThetaN::usage =
"EllipticThetaN[u, m] gives Neville's elliptic theta function theta_n(u, m)."
EllipticThetaD::usage =
"EllipticThetaD[u, m] gives Neville's elliptic theta function theta_d(u, m)."

 (* ---------------------- Weierstrass functions ---------------------- *)

WeierstrassP::usage =
"WeierstrassP[u, g2, g3] gives the Weierstrass elliptic function P."

WeierstrassPPrime::usage =
"WeierstrassPPrime[u, g2, g3] gives the derivative with respect to u of the
Weierstrass elliptic function P'."

InverseWeierstrassP::usage =
"InverseWeierstrassP[{P, P'}, g2, g3] gives the value of u such that
P = WeierstrassP[u, g2, g3] and P' = WeierstrassPPrime[u, g2, g3].  Note
that P and P' are not independent.  No attempt is made to ensure that the
result is in any particular period parallelogram."

 (* ---------------------- end of introduction ---------------------- *)

Begin["Elliptic`Private`"]

 (* ------------------ derivative of EllipticLog -----------------------*)

Unprotect[EllipticExp, EllipticLog, InverseFunction];

Off[EllipticLog::ellnp];
EllipticLog/: D[EllipticLog[{x_,y_}, _List], t_] := D[x, t]/(2y);
On[EllipticLog::ellnp];

 (* ----------------------- private functions -------------------------*)

(* :Note:
	The following variables are used as storage for previous values
	that are likely to be needed again.  We thus avoid the need to
	actually re-evaluate them.  They are paired as input and output
	values for various functions.
*)

$LastEllipticmK;	$LastEllipticK;
$LastEllipticmKp;	$LastEllipticKp;
$LastEllipticmT;	$LastEllipticT;
$LastEllipticmQ;	$LastEllipticQ;
$LastOneIn;		$LastOneOut;

$EllipK[m_?NumberQ] :=
    Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmK]},
	If[m === $LastEllipticmK && prec <= oldprec,
	    N[$LastEllipticK, prec],
	  (* else *)
	    prec = EllipticK[m];
	    AbortProtect[If[NumberQ[prec],
			$LastEllipticmK = m; $LastEllipticK = prec]];
	    prec]
	]

$EllipKp[m_?NumberQ] :=
    Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmKp]},
	If[m === $LastEllipticmKp && prec <= oldprec,
	    N[$LastEllipticKp, prec],
	  (* else *)
	    prec = EllipticK[1-m];
	    AbortProtect[If[NumberQ[prec],
			$LastEllipticmKp = m; $LastEllipticKp = prec]];
	    prec]
	]

EllipticNomeQ[m_?NumberQ] :=
    Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmQ]},
	If[TrueQ[m == 0.], Return[m]];
	If[m === $LastEllipticmQ && prec <= oldprec,
	    N[$LastEllipticQ, prec],
	  (* else *)
	    prec = Exp[-N[Pi, prec] $EllipKp[m]/$EllipK[m]];
	    AbortProtect[If[NumberQ[prec],
			$LastEllipticmQ = m; $LastEllipticQ = prec]];
	    prec]
	]

$EllipticTheta0[n_, m_] :=
    Module[{prec = Precision[m], oldprec = Precision[$LastEllipticmT]},
	If[m === $LastEllipticmT && prec <= oldprec,
	    N[$LastEllipticT[[n]], prec],
	  (* else *)
	    prec = {EllipticThetaPrime[1, 0, m], EllipticTheta[2, 0, m],
		EllipticTheta[3, 0, m], EllipticTheta[4, 0, m]};
	    AbortProtect[If[Apply[And, Map[NumberQ, prec]],
		$LastEllipticmT = m; $LastEllipticT = prec]];
	    prec[[n]]]
	]

$EllipticTheta[p_, u_, m_] :=
    Module[{prec = Precision[{u, m}], K2, mm, q},
	mm = N[m, prec];
	K2 = 2 $EllipK[mm];
	q = EllipticNomeQ[mm];
	If[p == 0,
	    K2 EllipticTheta[1, N[Pi u/K2, prec], q]/
		(N[Pi, prec] $EllipticTheta0[1, q]),
	  (* else *)
	    EllipticTheta[p+1, N[Pi u/K2, prec], q]/$EllipticTheta0[p+1, q]
	    ]
	]

EllipticThetaS[u_?NumberQ, m_?NumberQ] :=
	$EllipticTheta[0, u, m] /; Precision[{u, m}] < Infinity

EllipticThetaC[u_?NumberQ, m_?NumberQ] :=
	$EllipticTheta[1, u, m] /; Precision[{u, m}] < Infinity

EllipticThetaD[u_?NumberQ, m_?NumberQ] :=
	$EllipticTheta[2, u, m] /; Precision[{u, m}] < Infinity

EllipticThetaN[u_?NumberQ, m_?NumberQ] :=
	$EllipticTheta[3, u, m] /; Precision[{u, m}] < Infinity

EllipticThetaS[0, m_] := 0
EllipticThetaC[0, m_] := 1
EllipticThetaD[0, m_] := 1
EllipticThetaN[0, m_] := 1

Fix[eqn_Equal,x_,y_,form_] :=
    Module[{newx=x,newy=y,oldx=x,oldy=y,c,f,sub,newform = form,newout,
		null = {x->0,y->0},a,b,disc,e1,precision,c4,c6,rtdisc,s},
	If[ $LastOneIn === {eqn,x,y,form}, Return[$LastOneOut]];

	f = Expand[eqn[[1]]-eqn[[2]]];
	precision = Precision[f];
	If [ precision == Infinity , precision = 16 ];

	(* remove coeff of y^2 *)
	f = Expand[f / (D[f,{y,2}]/2)];

	(* remove coeff of x^3 *)
	c = -(D[f,{x,3}]/6);
	sub = {x -> x/c, y -> y/c};
	oldx = Expand[oldx c];
	oldy = Expand[oldy c];
	f = Expand[(f /. sub) c^2];
	newform = Expand[newform /. sub];
	newx = Expand[newx /. sub];
	newy = Expand[newy /. sub];

	(* remove linear y terms *)
	c = ((D[f,y] /. null) + (D[f,x,y] /. null) x) / 2;
	sub = {y -> y-c};
	oldy = Expand[oldy + c];
	f = Expand[(f /. sub)];
	newform = Expand[newform /. sub];
	newx = Expand[newx /. sub];
	newy = Expand[newy /. sub];

	(* remove x^2 term *)
	c = (D[f,{x,2}]/2 /. null)/ 3;
	sub = {x -> x+c};
	oldx = Expand[oldx - c];
	f = Expand[(f /. sub)];
	newform = Expand[newform /. sub];
	newx = Expand[newx /. sub];
	newy = Expand[newy /. sub];

	(* convert to inexact *)
	f = N[f,precision];

	(* compute the largest real root *)
	disc = 4 (D[f,x] /. null)^3 - 27 (f /. null)^2 ;
	c4 = (D[f,x] /. null) 48;
	c6 = (f /. null) 864;
	If[Head[disc] =!= Real && disc =!=0, Return[$Failed]];
	If[ disc < 0,
		rtdisc = 96 Sqrt [ - 3 disc ];
		e1 = ((c6+rtdisc)^(1/3) + (c6-rtdisc)^(1/3))/12,
	    (* else *)
		s = (c6/1728 + I Sqrt[ disc/108 ])^(1/3);
		e1 = 2 Re[s],
	    (* undecided *)
		Return[$Failed]
	];
	sub = {x->x+e1};
	f = Expand[(f /. sub)];
	newform = Expand[newform /. sub];
	newx = Expand[newx /. sub];
	newy = Expand[newy /. sub];
	oldx = Expand[oldx - e1];
	a = - (D[f,{x,2}] /. null)/2;
	b = - D[f,x] /. null;
	c = Expand[newform / (Dt[x] / (2y))];
	(* sigh, what if c hasn t simplified down to a number here ? *)
	If[c==0,c=1];
	newout = {{a,b}, Apply[Function,{{newx, newy} /. {x->#1, y->#2}}], 
			Apply[Function,{{oldx, oldy} /. {x->#1, y->#2}}], c};

	(* don't let interrupts break the memory function *)
	AbortProtect[{$LastOneIn, $LastOneOut} = {in, newout}];
	newout
	]

EqnEllipticExp[u_, {eqn_Equal, x_, y_, form_}] :=
    Module[{c = Fix[eqn,x,y,form]},
	If[Head[c] === List && Length[c] == 4,
		c = Apply[c[[2]], EllipticExp[u/c[[4]], c[[1]]]]];
	c
	]

EqnEllipticLog[xy_List, {eqn_Equal, x_, y_, form_}] :=
    Module[{c = Fix[eqn,x,y,form]},
	If[Head[c] === List && Length[c] == 4,
		c = c[[4]] EllipticLog[Apply[c[[3]], xy], c[[1]]]];
	c
	]

EqnEllipticExp[u_, {eqn_Equal, x_, y_}] := EqnEllipticExp[u, {eqn, x, y, 0}] 
EqnEllipticLog[u_, {eqn_Equal, x_, y_}] := EqnEllipticLog[u, {eqn, x, y, 0}] 

(* this one fills in the square root *)
$EllipticLog[xi_,a_,b_] :=
    EllipticLog[{xi, Sqrt[xi(xi(xi + a) + b)]}, {a,b}]

(* this one is useful for integrals with two arbitrary limits *)
$EllipticLog[x0_,x1_,a_,b_] := $EllipticLog[x1,a,b] - $EllipticLog[x0,a,b]

(* this converts from the parameter m to the pair a,b *)
$EllipticLogm[x_, m_] := $EllipticLog[x, m+1, m]

  (* ---------------- Weierstrass elliptic functions ---------------- *)

Weierstrass[u_?NumberQ, g2_?NumberQ, g3_?NumberQ] := 
    EqnEllipticExp[u, {Y^2 == 4X^3 - g2 X - g3 , X, Y, Dt[X]/Y}]

WeierstrassP[u_, g2_, g3_] :=
    Module[{p}, p[[1]] /; (Head[p = Weierstrass[u, g2, g3]] === List)];

WeierstrassPPrime[u_, g2_, g3_] :=
    Module[{p}, p[[2]] /; (Head[p = Weierstrass[u, g2, g3]] === List)];

InverseWeierstrassP[{p_, pp_}, g2_, g3_] := 
    Module[{u, Y, X},
	u /; (u = EqnEllipticLog[{p, pp}, {Y^2 == 4X^3-g2 X-g3, X, Y, Dt[X]/Y}];
		NumberQ[u])
	];

  (* ---------------- Weierstrass  derivatives ----------------------- *)

Derivative[n_Integer, 0, 0][WeierstrassP] ^:=
	Derivative[n-1, 0, 0][WeierstrassPPrime] /; n > 0

Derivative[1, 0, 0][WeierstrassPPrime] ^=
	(6 WeierstrassP[#1, #2, #3]^2 - #2/2)&;		(* AS 18.6.4, p640 *)

Derivative[n_, 0, 0][WeierstrassPPrime] ^:=
	Module[{i, r},
	    r = Evaluate[12 Sum[Binomial[n-2, i]
			Derivative[i, 0, 0][WeierstrassPPrime][#1, #2, #3] 
			Derivative[n-2-i, 0, 0][WeierstrassP][#1, #2, #3],
			{i, 0, n-2}]] &;
	    r = Share[r];
	    AbortProtect[
		Unprotect[WeierstrassPPrime];
		Literal[Derivative[n, 0, 0][WeierstrassPPrime]] ^= r;
		Protect[WeierstrassPPrime]
		];
	    r
	    ] /; n > 1

  (* -----------elementary special cases, AS, p571, 16.6 ---------------- *)

JacobiSN[u_, 0|0.] := Sin[u];	JacobiSN[u_, 1|1.] := Tanh[u];
JacobiCN[u_, 0|0.] := Cos[u];	JacobiCN[u_, 1|1.] := Sech[u];
JacobiDN[u_, 0|0.] := 1;	JacobiDN[u_, 1|1.] := Sech[u];

JacobiCD[u_, 0|0.] := Cos[u];	JacobiCD[u_, 1|1.] := 1;
JacobiSD[u_, 0|0.] := Sin[u];	JacobiSD[u_, 1|1.] := Sinh[u];
JacobiND[u_, 0|0.] := 1;	JacobiND[u_, 1|1.] := Cosh[u];

JacobiDC[u_, 0|0.] := Sec[u];	JacobiDC[u_, 1|1.] := 1;
JacobiNC[u_, 0|0.] := Sec[u];	JacobiNC[u_, 1|1.] := Cosh[u];
JacobiSC[u_, 0|0.] := Tan[u];	JacobiSC[u_, 1|1.] := Sinh[u];

JacobiNS[u_, 0|0.] := Csc[u];	JacobiNS[u_, 1|1.] := Coth[u];
JacobiDS[u_, 0|0.] := Csc[u];	JacobiDS[u_, 1|1.] := Csch[u];
JacobiCS[u_, 0|0.] := Cot[u];	JacobiCS[u_, 1|1.] := Csch[u];

 (* -------- special symbolic and exact cases, AS, p571, 16.5 ------------*)

	(* note: row 7 has been changed to deal with the singularity *)

$JacobiSpecialValues =
	{{0&, 1&, 1&},
	{(1 + (1-#)^(1/2))^(-1/2)&, (1-#)^(1/4)(1 + (1-#)^(1/2))^(-1/2)&,
		(1-#)^(1/4)&},
	{1&, 0&, (1-#)^(1/2)&},
	{I #^(-1/4)&, (1 + #^(1/2))^(1/2)(#^(-1/4))&, (1 + #^(1/2))^(1/2)&},
	{2^(-1/2) #^(-1/4) ((1 + #^(1/2))^(1/2) + I (1 - #^(1/2))^(1/2))&,
		2^(-1/2) (1-I) (1-#)^(1/4) (#^(-1/4))&,
		2^(-1/2) (1-#)^(1/4) ((1 + (1-#)^(1/2))^(1/2) -
			I (1 - (1-#)^(1/2))^(1/2))&},
	{#^(-1/4)&, -I (1 - #^(1/2))^(1/2) (#^(-1/4))&, (1 - #^(1/2))^(1/2)&},
	{1&, -I&, -I #^(1/2)&}  (* this row is ss, cs, ds at I K' *),
	{(1 - (1-#)^(1/2))^(-1/2)&, -I (1-(1-#)^(1/2))^(-1/2) ((1-#)^(1/4))&,
		-I (1-#)^(1/4)&},
	{#^(-1/2)&, -I (1-#)^(1/2) (#^(-1/2))&, 0&},
	(* row 10 is sn, cn, dn at 3K/2 + I K'/2.  It is not in AS *)
	{((1-I)((1+(1-#)^(1/2))^(1/2)+I(1-(1-#)^(1/2))^(1/2))/(2 #^(1/4)))&,
	   (-((1-#)/#)^(1/4)((1+#^(1/2))^(1/2)+I(1-#^(1/2))^(1/2))/
		((1+(1-#)^(1/2))^(1/2)-I(1-(1-#)^(1/2))^(1/2)))&,
	   ((1-#)^(1/4)2^(-1/2)((1+(1-#)^(1/2))^(1/2)+
		I(1-(1-#)^(1/2))^(1/2)))&}};

SymbJacobi[p_, q_, u_, m_] :=
    Module[{K, K1, r, s, sign = 1},
	K = EllipticK[m];
	K1 = EllipticK[1 - m];
	If[Head[K] =!= EllipticK || Head[K1] =!= EllipticK, Return[$Failed]];
	r = u /. {K -> 1, K1 -> 0};
	s = u /. {K -> 0, K1 -> 1};
	If[r K + s K1 =!= u, Return[$Failed]];
	r *= 2; s *= -2 I;
	If[!IntegerQ[r] || !IntegerQ[s], Return[$Failed]];
	r = Mod[r, 8]; s = Mod[s, 8];	(* 1/8 period phase shifts *)

	(* argument reductions: AS, p572, 16.8 *)
	(* shifts by periods and half periods *)
	If[r > 3, r -= 4; If[Mod[p+q, 4] != 1, sign = -sign]];	 (* u - 2K *)
	If[s > 3, s -= 4; If[p+q != 3, sign = -sign]];		 (* u - 2iK' *)
	If[s == 0 && r == 3,					 (* 2K - u *)
		r = 1; If[p==1 || q==1, sign = -sign]];
	If[r == 0 && s == 3,					 (* 2iK' - u *)
		s = 1; If[p==3 || q==3, sign = -sign]];
	If[r + 2s > 6,					  (* u = 2(K+iK`)-u *)
		s = 4 - s; r = 4 - r; If[p==2 || q==2, sign = -sign]];

	(* row in table 16.5, AS, p571.  row 10 is a special added case *)
	If[r == 3 && s == 1, r = 10, r += 3s + 1];

	(* dispose of singularities in row 7 first *)
	If[r == 7, If[p==3, Return[0]]; If[q==3, Return[ComplexInfinity]]];

	(* dispose of pn and nq *)
	If[p == 3, Return[sign/$JacobiSpecialValues[[r, q+1]][m]]];
	If[q == 3, Return[sign $JacobiSpecialValues[[r, p+1]][m]]];
	sign $JacobiSpecialValues[[r,p+1]][m]/$JacobiSpecialValues[[r,q+1]][m]
	];

JacobiSN[u_, m_] :=
    Module[{jac = SymbJacobi[0,3,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiCN[u_, m_] :=
    Module[{jac = SymbJacobi[1,3,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiDN[u_, m_] :=
    Module[{jac = SymbJacobi[2,3,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]

JacobiCD[u_, m_] :=
    Module[{jac = SymbJacobi[1,2,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiSD[u_, m_] :=
    Module[{jac = SymbJacobi[0,2,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiND[u_, m_] :=
    Module[{jac = SymbJacobi[3,2,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]

JacobiDC[u_, m_] :=
    Module[{jac = SymbJacobi[2,1,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiNC[u_, m_] :=
    Module[{jac = SymbJacobi[3,1,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiSC[u_, m_] :=
    Module[{jac = SymbJacobi[0,1,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]

JacobiNS[u_, m_] :=
    Module[{jac = SymbJacobi[3,0,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiDS[u_, m_] :=
    Module[{jac = SymbJacobi[2,0,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]
JacobiCS[u_, m_] :=
    Module[{jac = SymbJacobi[1,0,u,m]},
	jac /; (jac =!= $Failed)] /; u == 0 || !FreeQ[u, EllipticK]

 (*-------------------- general numerical cases ---------------------------*)

$NJacobi[pp_, qq_, u_, m_] :=
    Module[{p = pp, q = qq, K, K1, r, s, den, sign = 1, ru, mm},
	mm = N[m, Precision[u]];
	K = $EllipK[mm];
	K1 = $EllipKp[mm];
	den = Re[K] Re[K1] + Im[K] Im[K1];
	If[den == 0,
	    r = Floor[2(Re[K] Re[u] + Im[K] Im[u])/Abs[K]^2];
	    s = 0,
	  (* else *)
	    r = Floor[2(Re[K1] Re[u] + Im[K1] Im[u])/den];
	    s = Floor[2(Re[K] Im[u] - Im[K] Re[u])/den]
	    ];
	ru = u - (2 Floor[r/4] K + 2 I Floor[s/4] K1);
	r = Mod[r, 8]; s = Mod[s, 8];
	If[r > 3, If[Mod[p+q, 4] != 1, sign = -sign]];
	If[s > 3, s -= 4; If[p+q != 3, sign = -sign]];
	If[s > 1, ru = 2(K+I K1) - ru; If[p==2 || q==2, sign = -sign]];
	sign $EllipticTheta[p, ru, mm]/$EllipticTheta[q, ru, mm]
	]

(* here are the twelve functions : one for each ordered pair of letters
     from the set S,D,N,C *)

JacobiSN[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[0, 3, u, m] /; Precision[{u, m}] < Infinity
JacobiCN[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[1, 3, u, m] /; Precision[{u, m}] < Infinity
JacobiDN[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[2, 3, u, m] /; Precision[{u, m}] < Infinity

JacobiCD[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[1, 2, u, m] /; Precision[{u, m}] < Infinity
JacobiSD[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[0, 2, u, m] /; Precision[{u, m}] < Infinity
JacobiND[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[3, 2, u, m] /; Precision[{u, m}] < Infinity

JacobiDC[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[2, 1, u, m] /; Precision[{u, m}] < Infinity
JacobiNC[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[3, 1, u, m] /; Precision[{u, m}] < Infinity
JacobiSC[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[0, 1, u, m] /; Precision[{u, m}] < Infinity

JacobiNS[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[3, 0, u, m] /; Precision[{u, m}] < Infinity
JacobiDS[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[2, 0, u, m] /; Precision[{u, m}] < Infinity
JacobiCS[u_?NumberQ,m_?NumberQ] :=
	$NJacobi[1, 0, u, m] /; Precision[{u, m}] < Infinity

  (* ------------------ Derivatives: AS, p574, 16.16 ---------------------- *)

JacobiSN/: Derivative[n_Integer,0][JacobiSN] :=
	Module[{j, r},
	    r = Evaluate[
		Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCN][#1, #2] *
		    Derivative[n-1-j, 0][JacobiDN][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiSN];
		Literal[Derivative[n,0][JacobiSN]] ^= r;
		Protect[JacobiSN]
		];
	    r] /; n > 0
JacobiCN/: Derivative[n_Integer,0][JacobiCN] :=
	Module[{j, r},
	    r = Evaluate[
		-Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSN][#1, #2] *
		    Derivative[n-1-j, 0][JacobiDN][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiCN];
		Literal[Derivative[n,0][JacobiCN]] ^= r;
		Protect[JacobiCN]
		];
	    r] /; n > 0
JacobiDN/: Derivative[n_Integer,0][JacobiDN] :=
	Module[{j, r},
	    r = Evaluate[
		-#2 Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSN][#1, #2] *
		    Derivative[n-1-j, 0][JacobiCN][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiDN];
		Literal[Derivative[n,0][JacobiDN]] ^= r;
		Protect[JacobiDN]
		];
	    r] /; n > 0

JacobiCD/: Derivative[n_Integer,0][JacobiCD] :=
	Module[{j, r},
	    r = Evaluate[
		(#2-1) Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSD][#1, #2]*
		    Derivative[n-1-j, 0][JacobiND][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiCD];
		Literal[Derivative[n,0][JacobiCD]] ^= r;
		Protect[JacobiCD]
		];
	    r] /; n > 0
JacobiSD/: Derivative[n_Integer,0][JacobiSD] :=
	Module[{j, r},
	    r = Evaluate[
		Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCD][#1, #2] *
		    Derivative[n-1-j, 0][JacobiND][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiSD];
		Literal[Derivative[n,0][JacobiSD]] ^= r;
		Protect[JacobiSD]
		];
	    r] /; n > 0
JacobiND/: Derivative[n_Integer,0][JacobiND] :=
	Module[{j, r},
	    r = Evaluate[
		#2 Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSD][#1, #2] *
		    Derivative[n-1-j, 0][JacobiCD][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiND];
		Literal[Derivative[n,0][JacobiND]] ^= r;
		Protect[JacobiND]
		];
	    r] /; n > 0

JacobiDC/: Derivative[n_Integer,0][JacobiDC] :=
	Module[{j, r},
	    r = Evaluate[
		(1-#2) Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSC][#1, #2] *
		    Derivative[n-1-j, 0][JacobiNC][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiDC];
		Literal[Derivative[n,0][JacobiDC]] ^= r;
		Protect[JacobiDC]
		];
	    r] /; n > 0
JacobiNC/: Derivative[n_Integer,0][JacobiNC] :=
	Module[{j, r},
	    r = Evaluate[
		Sum[Binomial[n-1, j] Derivative[j, 0][JacobiSC][#1, #2] *
		    Derivative[n-1-j, 0][JacobiDC][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiNC];
		Literal[Derivative[n,0][JacobiNC]] ^= r;
		Protect[JacobiNC]
		];
	    r] /; n > 0
JacobiSC/: Derivative[n_Integer,0][JacobiSC] :=
	Module[{j, r},
	    r = Evaluate[
		Sum[Binomial[n-1, j] Derivative[j, 0][JacobiDC][#1, #2] *
		    Derivative[n-1-j, 0][JacobiNC][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiSC];
		Literal[Derivative[n,0][JacobiSC]] ^= r;
		Protect[JacobiSC]
		];
	    r] /; n > 0

JacobiNS/: Derivative[n_Integer,0][JacobiNS] :=
	Module[{j, r},
	    r = Evaluate[
		-Sum[Binomial[n-1, j] Derivative[j, 0][JacobiDS][#1, #2] *
		    Derivative[n-1-j, 0][JacobiCS][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiNS];
		Literal[Derivative[n,0][JacobiNS]] ^= r;
		Protect[JacobiNS]
		];
	    r] /; n > 0
JacobiDS/: Derivative[n_Integer,0][JacobiDS] :=
	Module[{j, r},
	    r = Evaluate[
		-Sum[Binomial[n-1, j] Derivative[j, 0][JacobiCS][#1, #2] *
		    Derivative[n-1-j, 0][JacobiNS][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiDS];
		Literal[Derivative[n,0][JacobiDS]] ^= r;
		Protect[JacobiDS]
		];
	    r] /; n > 0
JacobiCS/: Derivative[n_Integer,0][JacobiCS] :=
	Module[{j, r},
	    r = Evaluate[
		-Sum[Binomial[n-1, j] Derivative[j, 0][JacobiNS][#1, #2] *
		    Derivative[n-1-j, 0][JacobiDS][#1, #2], {j, 0, n-1}]]&;
	    r = Share[r];
	    AbortProtect[
		Unprotect[JacobiCS];
		Literal[Derivative[n,0][JacobiCS]] ^= r;
		Protect[JacobiCS]
		];
	    r] /; n > 0

  (* --------------------Inverse Jacobi functions ----------------------- *)

InverseFunction[JacobiSN, 1, 2] = InverseJacobiSN
InverseFunction[JacobiCN, 1, 2] = InverseJacobiCN
InverseFunction[JacobiDN, 1, 2] = InverseJacobiDN

InverseFunction[JacobiCD, 1, 2] = InverseJacobiCD
InverseFunction[JacobiSD, 1, 2] = InverseJacobiSD
InverseFunction[JacobiND, 1, 2] = InverseJacobiND

InverseFunction[JacobiDC, 1, 2] = InverseJacobiDC
InverseFunction[JacobiNC, 1, 2] = InverseJacobiNC
InverseFunction[JacobiSC, 1, 2] = InverseJacobiSC

InverseFunction[JacobiNS, 1, 2] = InverseJacobiNS
InverseFunction[JacobiDS, 1, 2] = InverseJacobiDS
InverseFunction[JacobiCS, 1, 2] = InverseJacobiCS

InverseFunction[InverseJacobiSN, 1, 2] = JacobiSN
InverseFunction[InverseJacobiCN, 1, 2] = JacobiCN
InverseFunction[InverseJacobiDN, 1, 2] = JacobiDN

InverseFunction[InverseJacobiCD, 1, 2] = JacobiCD
InverseFunction[InverseJacobiSD, 1, 2] = JacobiSD
InverseFunction[InverseJacobiND, 1, 2] = JacobiND

InverseFunction[InverseJacobiDC, 1, 2] = JacobiDC
InverseFunction[InverseJacobiNC, 1, 2] = JacobiNC
InverseFunction[InverseJacobiSC, 1, 2] = JacobiSC

InverseFunction[InverseJacobiNS, 1, 2] = JacobiNS
InverseFunction[InverseJacobiDS, 1, 2] = JacobiDS
InverseFunction[InverseJacobiCS, 1, 2] = JacobiCS

 (* ------------ elementary special cases, AS, p571, 16.6 -------------- *)

InverseJacobiSN[u_, 0|0.] := ArcSin[u];	InverseJacobiSN[u_, 1|1.] := ArcTanh[u];
InverseJacobiCN[u_, 0|0.] := ArcCos[u];	InverseJacobiCN[u_, 1|1.] := ArcSech[u];
(* InverseJacobiDN[_, 0] has no inv *)	InverseJacobiDN[u_, 1|1.] := ArcSech[u];

InverseJacobiCD[u_, 0|0.] := ArcCos[u];	(* InverseJacobiCD[_, 1] has no inv *)
InverseJacobiSD[u_, 0|0.] := ArcSin[u];	InverseJacobiSD[u_, 1|1.] := ArcSinh[u];
(* InverseJacobiND[_, 0] has no inv *)	InverseJacobiND[u_, 1|1.] := ArcCosh[u];

InverseJacobiDC[u_, 0|0.] := ArcSec[u];	(* InverseJacobiDC[_, 1] has no inv *)
InverseJacobiNC[u_, 0|0.] := ArcSec[u];	InverseJacobiNC[u_, 1|1.] := ArcCosh[u];
InverseJacobiSC[u_, 0|0.] := ArcTan[u];	InverseJacobiSC[u_, 1|1.] := ArcSinh[u];

InverseJacobiNS[u_, 0|0.] := ArcCsc[u];	InverseJacobiNS[u_, 1|1.] := ArcCoth[u];
InverseJacobiDS[u_, 0|0.] := ArcCsc[u];	InverseJacobiDS[u_, 1|1.] := ArcCsch[u];
InverseJacobiCS[u_, 0|0.] := ArcCot[u];	InverseJacobiCS[u_, 1|1.] := ArcCsch[u];

InverseFunction::noinv = "The function `1` is not invertible."

$NoInverse[f_, m_] := (Message[InverseFunction::noinv, f[#, m]&]; False)

InverseJacobiDN[_, m_] := 0 /; (TrueQ[m == 0.0] && $NoInverse[JacobiDN, 0])
InverseJacobiND[_, m_] := 0 /; (TrueQ[m == 0.0] && $NoInverse[JacobiND, 0])

InverseJacobiCD[_, m_] := 0 /; (TrueQ[m-1 == 0.0] && $NoInverse[JacobiCD, 1])
InverseJacobiDC[_, m_] := 0 /; (TrueQ[m-1 == 0.0] && $NoInverse[JacobiDC, 1])

 (* ---------------- general numerical inverses --------------------------*)

$InverseJacobi[jac_,a_,b_] :=
    Module[{j = -jac/2, s, t, A, B, CC},
	If[j === 0,
	    0,
	(* else *)
	    A = j^3;
	    B = a j^2 - 1;
	    CC = b j;
	    t = Sqrt[B^2-4 A CC];
	    If[(B == 0) || Re[t/B] <= 0,
		t = (t-B)/(2 A),
	      (* else *)
		t = (-2 CC)/(t+B)];
	    s = j t;
	    EllipticLog[{s,t},{a,b}]
	    ]]

(* the Jacobi functions don't REALLY have inverses *)

(*					.	    .
					.         .
					.       .
					N D N D
					S C S C
   Here is the diagram of zeros:	N D N D
					S C S C . . .
 *)

$InvJacArgRed[u_, m_, n_, o_, p_] :=
    Module[{ru, r, s, d, K, IK1},
	K = $EllipK[N[m, Precision[u]]];
	IK1 = I $EllipKp[N[m, Precision[u]]];
	ru = u - (n K + o IK1);
	d = 2(Re[K] Im[IK1] - Im[K] Re[IK1]);
	(* what do we do if K and IK1 are not independent? *)
	r = Round[(Im[IK1] Re[ru] - Re[IK1] Im[ru])/d];
	s = Round[(Re[K] Im[ru] - Im[K] Re[ru])/d];
	ru = u - 2(r K + s IK1);
	If[Switch[p, 1, OddQ[r], 2, OddQ[s], 3, Xor[OddQ[r], OddQ[s]]],
	    ru = 2(n K + o IK1) - ru];
	ru
	]

InverseJacobiSN[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[v, 2 + 2m, (1 - m)^2];
	$InvJacArgRed[u, m, 0, 0, 1]
	] /; Precision[{v, m}] < Infinity
InverseJacobiCN[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = $EllipK[m] - 2 $InverseJacobi[v/Sqrt[1-m], 2 - 4m, 1 ];
	$InvJacArgRed[u, m, 1, 0, 3]
	] /; Precision[{v, m}] < Infinity
InverseJacobiDN[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[I/v, -4 + 2m, m^2] - I $EllipKp[m];
	$InvJacArgRed[u, m, 1, 1, 2]
	] /; (Precision[{v, m}] < Infinity && m != 0.0)

InverseJacobiCD[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = $EllipK[m] - 2 $InverseJacobi[v, 2 + 2m, (1 - m)^2];
	$InvJacArgRed[u, m, 1, 0, 1]
	] /; (Precision[{v, m}] < Infinity && (m-1) != 0.0)
InverseJacobiSD[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[v, 2 - 4m, 1 ];
	$InvJacArgRed[u, m, 0, 0, 3]
	] /; Precision[{v, m}] < Infinity
InverseJacobiND[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[I v, -4 + 2m, m^2] - I $EllipKp[m];
	$InvJacArgRed[u, m, 0, 1, 2]
	] /; (Precision[{v, m}] < Infinity && m != 0.0)

InverseJacobiDC[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = $EllipK[m] - 2 $InverseJacobi[1/v, 2 + 2m, (1 - m)^2];
	$InvJacArgRed[u, m, 1, 1, 1]
	] /; (Precision[{v, m}] < Infinity && (m-1) != 0.0)
InverseJacobiNC[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = $EllipK[m] - 2 $InverseJacobi[1/(v Sqrt[1-m]), 2 - 4m, 1 ];
	$InvJacArgRed[u, m, 0, 1, 3]
	] /; Precision[{v, m}] < Infinity
InverseJacobiSC[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[v, -4 + 2m, m^2];
	$InvJacArgRed[u, m, 0, 0, 2]
	] /; Precision[{v, m}] < Infinity

InverseJacobiNS[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[1/v, 2 + 2m, (1 - m)^2];;
	$InvJacArgRed[u, m, 0, 1, 1]
	] /; Precision[{v, m}] < Infinity
InverseJacobiDS[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[1/v, 2 - 4m, 1 ];
	$InvJacArgRed[u, m, 1, 1, 3]
	] /; Precision[{v, m}] < Infinity
InverseJacobiCS[v_?NumberQ,m_?NumberQ] :=
    Module[{u},
	u = 2 $InverseJacobi[1/v, -4 + 2m, m^2];
	$InvJacArgRed[u, m, 1, 0, 2]
	] /; Precision[{v, m}] < Infinity

 (* ------------ derivatives of the inverse functions ----------------- *)

InverseJacobiSN/: Derivative[1,0][InverseJacobiSN] :=
		(1/Derivative[1,0][JacobiSN][InverseJacobiSN[#1, #2], #2])&;
InverseJacobiCN/: Derivative[1,0][InverseJacobiCN] :=
		(1/Derivative[1,0][JacobiCN][InverseJacobiCN[#1, #2], #2])&;
InverseJacobiDN/: Derivative[1,0][InverseJacobiDN] :=
		(1/Derivative[1,0][JacobiDN][InverseJacobiDN[#1, #2], #2])&;

InverseJacobiCD/: Derivative[1,0][InverseJacobiCD] :=
		(1/Derivative[1,0][JacobiCD][InverseJacobiCD[#1, #2], #2])&;
InverseJacobiSD/: Derivative[1,0][InverseJacobiSD] :=
		(1/Derivative[1,0][JacobiSD][InverseJacobiSD[#1, #2], #2])&;
InverseJacobiND/: Derivative[1,0][InverseJacobiND] :=
		(1/Derivative[1,0][JacobiND][InverseJacobiND[#1, #2], #2])&;

InverseJacobiDC/: Derivative[1,0][InverseJacobiDC] :=
		(1/Derivative[1,0][JacobiDC][InverseJacobiDC[#1, #2], #2])&;
InverseJacobiNC/: Derivative[1,0][InverseJacobiNC] :=
		(1/Derivative[1,0][JacobiNC][InverseJacobiNC[#1, #2], #2])&;
InverseJacobiSC/: Derivative[1,0][InverseJacobiSC] :=
		(1/Derivative[1,0][JacobiSC][InverseJacobiSC[#1, #2], #2])&;

InverseJacobiNS/: Derivative[1,0][InverseJacobiNS] :=
		(1/Derivative[1,0][JacobiNS][InverseJacobiNS[#1, #2], #2])&;
InverseJacobiDS/: Derivative[1,0][InverseJacobiDS] :=
		(1/Derivative[1,0][JacobiDS][InverseJacobiDS[#1, #2], #2])&;
InverseJacobiCS/: Derivative[1,0][InverseJacobiCS] :=
		(1/Derivative[1,0][JacobiCS][InverseJacobiCS[#1, #2], #2])&;

 (* ----------------------- JacobiAmplitude ------------------------ *)

JacobiAmplitude[u_, 0] := u;
JacobiAmplitude[u_, 1] := 2 ArcTan[Exp[u]] - Pi/2;
JacobiAmplitude[u_?NumberQ, m_?NumberQ] :=
	Module[{am = ArcSin[JacobiSN[u,m]], pi, s, k2, kp, km},
		k2 = 2 $EllipK[m];
		s = EllipticF[am, m];
		pi = N[Pi, Accuracy[am]+2];
		km = Re[(u - s)/k2];
		kp = Re[(u + s)/k2];
		If[Abs[km - Round[km]] > Abs[kp - Round[kp]],
			Round[kp] pi - am,
		   (* else *)
			Round[km] pi + am]
		];

JacobiAmplitude/: Derivative[1,0][JacobiAmplitude] = JacobiDN

InverseFunction[JacobiAmplitude, 1, 2] = EllipticF
InverseFunction[EllipticF, 1, 2] = JacobiAmplitude 

End[] (* Private *)

SetAttributes[JacobiSN, ReadProtected];
SetAttributes[JacobiSD, ReadProtected];
SetAttributes[JacobiSC, ReadProtected];
SetAttributes[JacobiCS, ReadProtected];
SetAttributes[JacobiCN, ReadProtected];
SetAttributes[JacobiCD, ReadProtected];
SetAttributes[JacobiNS, ReadProtected];
SetAttributes[JacobiNC, ReadProtected];
SetAttributes[JacobiND, ReadProtected];
SetAttributes[JacobiDN, ReadProtected];
SetAttributes[JacobiDS, ReadProtected];
SetAttributes[JacobiDC, ReadProtected];
SetAttributes[InverseJacobiSN, ReadProtected];
SetAttributes[InverseJacobiSD, ReadProtected];
SetAttributes[InverseJacobiSC, ReadProtected];
SetAttributes[InverseJacobiCS, ReadProtected];
SetAttributes[InverseJacobiCN, ReadProtected];
SetAttributes[InverseJacobiCD, ReadProtected];
SetAttributes[InverseJacobiNS, ReadProtected];
SetAttributes[InverseJacobiNC, ReadProtected];
SetAttributes[InverseJacobiND, ReadProtected];
SetAttributes[InverseJacobiDN, ReadProtected];
SetAttributes[InverseJacobiDS, ReadProtected];
SetAttributes[InverseJacobiDC, ReadProtected];
SetAttributes[JacobiAmplitude, ReadProtected];
SetAttributes[EllipticNomeQ, ReadProtected];
SetAttributes[EllipticThetaS, ReadProtected];
SetAttributes[EllipticThetaC, ReadProtected];
SetAttributes[EllipticThetaN, ReadProtected];
SetAttributes[EllipticThetaD, ReadProtected];
SetAttributes[WeierstrassP, ReadProtected];
SetAttributes[WeierstrassPPrime, ReadProtected];
SetAttributes[InverseWeierstrassP, ReadProtected];
SetAttributes[EllipticExp, ReadProtected];
SetAttributes[EllipticLog, ReadProtected];

Protect[JacobiSN, JacobiSD, JacobiSC, JacobiCS, JacobiCN, JacobiCD,
	JacobiNS, JacobiNC, JacobiND, JacobiDN, JacobiDS, JacobiDC]
Protect[InverseJacobiSN, InverseJacobiSD, InverseJacobiSC, InverseJacobiCS,
	InverseJacobiCN, InverseJacobiCD, InverseJacobiNS, InverseJacobiNC,
	InverseJacobiDS, InverseJacobiDC, InverseJacobiND, InverseJacobiDN]
Protect[JacobiAmplitude, EllipticNomeQ, EllipticThetaS, EllipticThetaC,
	EllipticThetaN, EllipticThetaD, WeierstrassP, WeierstrassPPrime,
	InverseWeierstrassP]
Protect[EllipticExp, EllipticLog, InverseFunction]

End[] (* Package *)

(* :Tests: *)

(* :Note:
	To run the tests just type `check', `complexcheck', `complexcheck2',
	or `symbcheckall'.  The test symbcheck[n] checks the symbolic values
	at (n K + I K' Range[0,15])/2.
*)

(*
see[result_,parms___] := If[Chop[result] === 0, Okay, Error[result,parms]]

see2[result_,parms___] :=
	If[Chop[result, .1^5] === 0, Okay, Error[result,parms]]

check1[ x_, m_ ] := see2 [ x - JacobiAmplitude[EllipticF[x, m], m]]
check2[u_,m_] := see[ - JacobiDN[u,m]^2 + 1 - m + m JacobiCN[u,m]^2 ]
check3[u_,m_] := see[ - m JacobiCN[u,m]^2 - m JacobiSN[u,m]^2 + m ]
check4[u_,m_] := see[ - (1-m) JacobiND[u,m]^2 + (1-m) + m(1-m)JacobiSD[u,m]^2 ]
check5[u_,m_] := see[ -m(1-m)JacobiSD[u,m]^2 - m JacobiCD[u,m]^2 + m ]
checkP[u_, g2_, g3_] :=
    Module[{p, pp, uu, err},
	p = WeierstrassP[u, g2, g3];
	pp = WeierstrassPPrime[u, g2, g3];
	uu = InverseWeierstrassP[{p, pp}, g2, g3];
	err = Abs[p-WeierstrassP[uu, g2, g3]];
	err += Abs[pp-WeierstrassPPrime[uu, g2, g3]];
	see[err, P, u, g2, g3]
	];

checkinv[m_] :=
  Module[{K, K1, u, r},
     K = Elliptic`Private`$EllipK[m];
     K1 = Elliptic`Private`$EllipKp[m];
     r = Table[, {12}];
     u = Random[Real, {-1,1}] K + Random[Real, {-1,1}] I K1;
     r[[1]] = see[ u - InverseJacobiSN[JacobiSN[u,m],m] , u , m, SN ];
     r[[2]] = see[ u - InverseJacobiSD[JacobiSD[u,m],m] , u , m, SD ];
     r[[3]] = see[ u - InverseJacobiSC[JacobiSC[u,m],m] , u , m, SC ];
     u += K;
     r[[4]] = see[ u - InverseJacobiCD[JacobiCD[u,m],m] , u , m, CD ];
     r[[5]] = see[ u - InverseJacobiCN[JacobiCN[u,m],m] , u , m, CN ];
     r[[6]] = see[ u - InverseJacobiCS[JacobiCS[u,m],m] , u , m, CS ];
     u += I K1;
     r[[7]] = see[ u - InverseJacobiDC[JacobiDC[u,m],m] , u , m, DC ];
     r[[8]] = see[ u - InverseJacobiDS[JacobiDS[u,m],m] , u , m, DS ];
     r[[9]] = see[ u - InverseJacobiDN[JacobiDN[u,m],m] , u , m, DN ];
     u -= K;
     r[[10]] = see[ u - InverseJacobiNS[JacobiNS[u,m],m] , u , m, NS ];
     r[[11]] = see[ u - InverseJacobiNC[JacobiNC[u,m],m] , u , m, NC ];
     r[[12]] = see[ u - InverseJacobiND[JacobiND[u,m],m] , u , m, ND ];
     r
     ]

checkinv2[u_,m_] := {
     see[ u - JacobiSD[InverseJacobiSD[u,m],m] , u , m, SD ],
     see[ u - JacobiSN[InverseJacobiSN[u,m],m] , u , m, SN ],
     see[ u - JacobiSC[InverseJacobiSC[u,m],m] , u , m, SC ],
     see[ u - JacobiCS[InverseJacobiCS[u,m],m] , u , m, CS ],
     see[ u - JacobiCD[InverseJacobiCD[u,m],m] , u , m, CD ],
     see[ u - JacobiCN[InverseJacobiCN[u,m],m] , u , m, CN ],
     see[ u - JacobiDS[InverseJacobiDS[u,m],m] , u , m, DS ],
     see[ u - JacobiDC[InverseJacobiDC[u,m],m] , u , m, DC ],
     see[ u - JacobiNS[InverseJacobiNS[u,m],m] , u , m, NS ],
     see[ u - JacobiND[InverseJacobiND[u,m],m] , u , m, ND ],
     see[ u - JacobiDN[InverseJacobiDN[u,m],m] , u , m, DN ],
     see[ u - JacobiNC[InverseJacobiNC[u,m],m] , u , m, NC ]
     }

checkder[u_, m_] := {
     see2[(JacobiSD[u + .00001, m] - JacobiSD[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiSD][u, m]) - 1, u, m, SD],
     see2[(JacobiSN[u + .00001, m] - JacobiSN[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiSN][u, m]) - 1, u, m, SN],
     see2[(JacobiSC[u + .00001, m] - JacobiSC[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiSC][u, m]) - 1, u, m, SC],
     see2[(JacobiCS[u + .00001, m] - JacobiCS[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiCS][u, m]) - 1, u, m, CS],
     see2[(JacobiCD[u + .00001, m] - JacobiCD[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiCD][u, m]) - 1, u, m, CD],
     see2[(JacobiCN[u + .00001, m] - JacobiCN[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiCN][u, m]) - 1, u, m, CN],
     see2[(JacobiDS[u + .00001, m] - JacobiDS[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiDS][u, m]) - 1, u, m, DS],
     see2[(JacobiDC[u + .00001, m] - JacobiDC[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiDC][u, m]) - 1, u, m, DC],
     see2[(JacobiNS[u + .00001, m] - JacobiNS[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiNS][u, m]) - 1, u, m, NS],
     see2[(JacobiND[u + .00001, m] - JacobiND[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiND][u, m]) - 1, u, m, ND],
     see2[(JacobiDN[u + .00001, m] - JacobiDN[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiDN][u, m]) - 1, u, m, DN],
     see2[(JacobiNC[u + .00001, m] - JacobiNC[u - .00001, m])/
		(.00002 Derivative[1, 0][JacobiNC][u, m]) - 1, u, m, NC]}

checkPder[u_, g2_, g3_] := {
     see2[(WeierstrassP[u+.00001,g2,g3]-WeierstrassP[u-.00001,g2,g3])/
	(.00002 Derivative[1, 0, 0][WeierstrassP][u,g2,g3]) - 1,u,g2,g3,P],
     see2[(WeierstrassPPrime[u+.00001,g2,g3]-WeierstrassPPrime[u-.00001,g2,g3])/
	(.00002 Derivative[1,0,0][WeierstrassPPrime][u,g2,g3])-1,u,g2,g3,PP]}

checkderinv[u_, m_] := {
     see2[(InverseJacobiSD[u + .00001, m] - InverseJacobiSD[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiSD][u, m]) - 1, u, m, SD],
     see2[(InverseJacobiSN[u + .00001, m] - InverseJacobiSN[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiSN][u, m]) - 1, u, m, SN],
     see2[(InverseJacobiSC[u + .00001, m] - InverseJacobiSC[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiSC][u, m]) - 1, u, m, SC],
     see2[(InverseJacobiCS[u + .00001, m] - InverseJacobiCS[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiCS][u, m]) - 1, u, m, CS],
     see2[(InverseJacobiCD[u + .00001, m] - InverseJacobiCD[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiCD][u, m]) - 1, u, m, CD],
     see2[(InverseJacobiCN[u + .00001, m] - InverseJacobiCN[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiCN][u, m]) - 1, u, m, CN],
     see2[(InverseJacobiDS[u + .00001, m] - InverseJacobiDS[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiDS][u, m]) - 1, u, m, DS],
     see2[(InverseJacobiDC[u + .00001, m] - InverseJacobiDC[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiDC][u, m]) - 1, u, m, DC],
     see2[(InverseJacobiNS[u + .00001, m] - InverseJacobiNS[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiNS][u, m]) - 1, u, m, NS],
     see2[(InverseJacobiND[u + .00001, m] - InverseJacobiND[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiND][u, m]) - 1, u, m, ND],
     see2[(InverseJacobiDN[u + .00001, m] - InverseJacobiDN[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiDN][u, m]) - 1, u, m, DN],
     see2[(InverseJacobiNC[u + .00001, m] - InverseJacobiNC[u - .00001, m])/
	(.00002 Derivative[1, 0][InverseJacobiNC][u, m]) - 1, u, m, NC]}

check := {
	check1[Random[Real, {-25, 25}],Random[]] ,
	check2[Random[Real, {-25, 25}],Random[]] ,
	check3[Random[Real, {-25, 25}],Random[]] ,
	check4[Random[Real, {-25, 25}],Random[]] ,
	check5[Random[Real, {-25, 25}],Random[]] ,
	checkP[Random[Real, {-5, 5}],
		Random[Real, {-10, 10}], Random[Real, {-10, 10}]] ,
	checkinv[Random[]],
	checkinv2[Random[],Random[]],
	checkder[Random[Real, {-25, 25}],Random[]],
	checkPder[Random[Real, {-5, 5}],
		Random[Real, {-5, 5}],Random[Real, {-5, 5}]]
	}

complexcheck := {
	check1[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
	check2[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
	check3[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
	check4[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
	check5[Random[Complex, {-25-25I, 25+25I}],Random[]] ,
	checkP[Random[Complex, {-5-5I, 5+5I}],
		Random[Real, {-10, 10}], Random[Real, {-10, 10}]] ,
	checkinv[Random[]],
	checkinv2[Random[Complex, {-25-25I, 25+25I}],Random[]],
	checkder[Random[Complex, {-25-25I, 25+25I}],Random[]],
	checkderinv[Random[Complex, {-25-25I, 25+25I}],Random[]],
	checkPder[Random[Complex, {-5-5I,5+5I}],
		Random[Real, {-5, 5}], Random[Real, {-5, 5}]]
	}

complexcheck2 := {
	check1[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
	check2[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
	check3[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
	check4[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
	check5[Random[Complex, {-25-25I, 25+25I}],Random[Complex]] ,
	(*  Weierstrass doesn't support complex invariants.
	checkP[Random[Complex, {-5-5I, 5+5I}],
		Random[Complex, {-10-10I, 10+10I}],
		Random[Complex, {-10-10I, 10+10I}]] ,
	*)
	checkinv[Random[Complex]],
	checkinv2[Random[Complex, {-25-25I, 25+25I}],Random[Complex]],
	checkder[Random[Complex, {-25-25I, 25+25I}],Random[Complex]],
	checkderinv[Random[Complex, {-25-25I, 25+25I}],Random[Complex]]
	}

argtab[nK_, K1_] := Table[nK + I K1 j/2, {j, 0, 15}];
compare[{jac_, njac_}] :=
    Module[{},
	If[Head[jac] === DirectedInfinity && Abs[njac] > 10.^5, Return[0]];
	If[jac ==  0,
	    If[Abs[njac] < 0.1^5, Return[0], Return[{jac, njac}]]];
	If[NumberQ[jac] && NumberQ[njac], Return[Chop[njac/jac - 1, .1^5]]];
	{jac, njac}];

scheck[f_, n_] :=
    Module[{m, nm = Random[ ], jac, njac, err},
	jac = Map[f[#, m]&, argtab[EllipticK[m] n/2, EllipticK[1-m]]];
	njac = Map[f[#, nm]&,
		.1^13 + argtab[EllipticK[nm] n/2, EllipticK[1-nm]]];
	jac = jac /. EllipticK[m] -> EllipticK[nm];
	jac = jac /. EllipticK[1-m] -> EllipticK[1-nm];
	jac = N[(jac /. m -> nm)];
	jac = Transpose[{jac, njac}];
	err = Map[compare, jac];
	Print[err];
	err]

symbcheck[n_] := {
	(* check the symbolic values at (n K + I K' Range[0,15])/2 *)
	{scheck[JacobiSN, n], SN},
	{scheck[JacobiCN, n], CN},
	{scheck[JacobiDN, n], DN},
	{scheck[JacobiCD, n], CD},
	{scheck[JacobiSD, n], SD},
	{scheck[JacobiND, n], ND},
	{scheck[JacobiDC, n], DC},
	{scheck[JacobiNC, n], NC},
	{scheck[JacobiSC, n], SC},
	{scheck[JacobiNS, n], NS},
	{scheck[JacobiDS, n], DS},
	{scheck[JacobiCS, n], CS}
	}

symbcheckall :=
	(* check the symbolic values at (n K + m I K')/2, 0 <= m, n <= 15 *)
	Table[Print[n," * K / 2"]; symbcheck[n], {n, 0, 15}]

*)
