Weierstrass elliptic function

"P-function" redirects here. For the phase-space function representing a quantum state, see Glauber–Sudarshan P representation.

In mathematics, the Weierstrass elliptic functions are elliptic functions that take a particularly simple form. They are named for Karl Weierstrass. This class of functions is also referred to as ℘-functions and they are usually denoted by the symbol ℘, a uniquely fancy script p. They play an important role in the theory of elliptic functions, i.e., meromorphic functions that are doubly periodic. A ℘-function together with its derivative can be used to parameterize elliptic curves and they generate the field of elliptic functions with respect to a given period lattice.

Symbol for Weierstrass $\mathrm{\wp }$-function

Model of Weierstrass $\mathrm{\wp }$-function

Motivation

A cubic of the form $C_{g_2,g_3}^{ℂ}=\{(x,y)\in {ℂ}^2:y^2=4x^3-g_{2}x-g_3\}$, where $g_2,g_3\in ℂ$ are complex numbers with $g_2^3-27g_3^2\ne 0$, cannot be rationally parameterized. Yet one still wants to find a way to parameterize it.

For the quadric $K=\{(x,y)\in {ℝ}^2:x^2+y^2=1\}$; the unit circle, there exists a (non-rational) parameterization using the sine function and its derivative the cosine function: $$\psi :ℝ/2\pi \mathbb{Z}\to K,t\mapsto (\sin t,\cos t).$$ Because of the periodicity of the sine and cosine $ℝ/2\pi \mathbb{Z}$ is chosen to be the domain, so the function is bijective.

In a similar way one can get a parameterization of $C_{g_2,g_3}^{ℂ}$ by means of the doubly periodic $\mathrm{\wp }$-function and its derivative, namely via $(x,y)=(\mathrm{\wp }(z),{\mathrm{\wp }}^{\prime }(z))$. This parameterization has the domain $ℂ/\Lambda$, which is topologically equivalent to a torus.

There is another analogy to the trigonometric functions. Consider the integral function $$a(x)=\underset{0}{\overset{x}{\int }}\frac{dy}{\sqrt{1-y^2}}.$$ It can be simplified by substituting $y=\sin t$ and $s=\arcsin x$: $$a(x)=\underset{0}{\overset{s}{\int }}dt=s=\arcsin x.$$ That means $a^{-1}(x)=\sin x$. So the sine function is an inverse function of an integral function.

Elliptic functions are the inverse functions of elliptic integrals. In particular, let: $$u(z)=\underset{z}{\overset{\mathrm{\infty }}{\int }}\frac{ds}{\sqrt{4s^3-g_{2}s-g_3}}.$$ Then the extension of $u^{-1}$ to the complex plane equals the $\mathrm{\wp }$-function. This invertibility is used in complex analysis to provide a solution to certain nonlinear differential equations satisfying the Painlevé property, i.e., those equations that admit poles as their only movable singularities.

Definition

Visualization of the $\mathrm{\wp }$-function with invariants $g_2=1+i$ and $g_3=2-3i$ in which white corresponds to a pole, black to a zero.

Let ${\omega }_1,{\omega }_2\in ℂ$ be two complex numbers that are linearly independent over $ℝ$ and let $\Lambda :=\mathbb{Z}{\omega }_1+\mathbb{Z}{\omega }_2:=\{m{\omega }_1+n{\omega }_2:m,n\in \mathbb{Z}\}$ be the period lattice generated by those numbers. Then the $\mathrm{\wp }$-function is defined as follows:

$\mathrm{\wp }(z,{\omega }_1,{\omega }_2):=\mathrm{\wp }(z)=\frac{1}{z^2}+\sum _{\lambda \in \Lambda \setminus \{0\}}(\frac{1}{(z-\lambda {)}^2}-\frac{1}{{\lambda }^2}).$

This series converges locally uniformly absolutely in the complex torus $ℂ/\Lambda$.

It is common to use $1$ and $\tau$ in the upper half-plane $ℍ:=\{z\in ℂ:\mathrm{Im}(z)>0\}$ as generators of the lattice. Dividing by ${\omega }_1$ maps the lattice $\mathbb{Z}{\omega }_1+\mathbb{Z}{\omega }_2$ isomorphically onto the lattice $\mathbb{Z}+\mathbb{Z}\tau$ with $\tau =\frac{{\omega }_2}{{\omega }_1}$. Because $-\tau$ can be substituted for $\tau$, without loss of generality we can assume $\tau \in ℍ$, and then define $\mathrm{\wp }(z,\tau ):=\mathrm{\wp }(z,1,\tau )$. With that definition, we have $\mathrm{\wp }(z,{\omega }_1,{\omega }_2)={\omega }_1^{-2}\mathrm{\wp }(z/{\omega }_1,{\omega }_2/{\omega }_1)$.

Properties

• $\mathrm{\wp }$ is a meromorphic function with a pole of order 2 at each period $\lambda$ in $\Lambda$.
• $\mathrm{\wp }$ is a homogeneous function in that:

$\mathrm{\wp }(\lambda z,\lambda {\omega }_1,\lambda {\omega }_2)={\lambda }^{-2}\mathrm{\wp }(z,{\omega }_1,{\omega }_2).$

• $\mathrm{\wp }$ is an even function. That means $\mathrm{\wp }(z)=\mathrm{\wp }(-z)$ for all $z\in ℂ\setminus \Lambda$, which can be seen in the following way:

$\begin{array}{cc}\mathrm{\wp }(-z)& =\frac{1}{(-z{)}^2}+\sum _{\lambda \in \Lambda \setminus \{0\}}(\frac{1}{(-z-\lambda {)}^2}-\frac{1}{{\lambda }^2})\\ & =\frac{1}{z^2}+\sum _{\lambda \in \Lambda \setminus \{0\}}(\frac{1}{(z+\lambda {)}^2}-\frac{1}{{\lambda }^2})\\ & =\frac{1}{z^2}+\sum _{\lambda \in \Lambda \setminus \{0\}}(\frac{1}{(z-\lambda {)}^2}-\frac{1}{{\lambda }^2})=\mathrm{\wp }(z).\end{array}$

The second last equality holds because $\{-\lambda :\lambda \in \Lambda \}=\Lambda$. Since the sum converges absolutely this rearrangement does not change the limit.

• The derivative of $\mathrm{\wp }$ is given by: $${\mathrm{\wp }}^{\prime }(z)=-2\sum _{\lambda \in \Lambda }\frac{1}{(z-\lambda {)}^3}.$$
• $\mathrm{\wp }$ and ${\mathrm{\wp }}^{\prime }$ are doubly periodic with the periods ${\omega }_1$ and ${\omega }_2$. This means: $$\begin{array}{cc}\mathrm{\wp }(z+{\omega }_1)& =\mathrm{\wp }(z)=\mathrm{\wp }(z+{\omega }_2),\text{and}\\ {\mathrm{\wp }}^{\prime }(z+{\omega }_1)& ={\mathrm{\wp }}^{\prime }(z)={\mathrm{\wp }}^{\prime }(z+{\omega }_2).\end{array}$$ It follows that $\mathrm{\wp }(z+\lambda )=\mathrm{\wp }(z)$ and ${\mathrm{\wp }}^{\prime }(z+\lambda )={\mathrm{\wp }}^{\prime }(z)$ for all $\lambda \in \Lambda$.

Laurent expansion

Let $r:=\min \{|\lambda |:0\ne \lambda \in \Lambda \}$. Then for $0<|z|<r$ the $\mathrm{\wp }$-function has the following Laurent expansion $$\mathrm{\wp }(z)=\frac{1}{z^2}+\sum _{n=1}^{\mathrm{\infty }}(2n+1)G_{2n+2}{z}^{2n}$$ where $$G_n=\sum _{0\ne \lambda \in \Lambda }{\lambda }^{-n}$$ for $n\ge 3$ are so called Eisenstein series.

Differential equation

Set $g_2=60G_4$ and $g_3=140G_6$. Then the $\mathrm{\wp }$-function satisfies the differential equation $$\mathrm{\wp }{\text{'}}^2(z)=4{\mathrm{\wp }}^3(z)-g_2\mathrm{\wp }(z)-g_3.$$ This relation can be verified by forming a linear combination of powers of $\mathrm{\wp }$ and ${\mathrm{\wp }}^{\prime }$ to eliminate the pole at $z=0$. This yields an entire elliptic function that has to be constant by Liouville's theorem.

Invariants

The real part of the invariant g₃ as a function of the square of the nome q on the unit disk.

The imaginary part of the invariant g₃ as a function of the square of the nome q on the unit disk.

The coefficients of the above differential equation $g_2$ and $g_3$ are known as the invariants. Because they depend on the lattice $\Lambda$ they can be viewed as functions in ${\omega }_1$ and ${\omega }_2$.

The series expansion suggests that $g_2$ and $g_3$ are homogeneous functions of degree $-4$ and $-6$. That is $$g_2(\lambda {\omega }_1,\lambda {\omega }_2)={\lambda }^{-4}{g}_2({\omega }_1,{\omega }_2)$$ $$g_3(\lambda {\omega }_1,\lambda {\omega }_2)={\lambda }^{-6}{g}_3({\omega }_1,{\omega }_2)$$ for $\lambda \ne 0$.

If ${\omega }_1$ and ${\omega }_2$ are chosen in such a way that $\mathrm{Im}\left(\frac{{\omega }_2}{{\omega }_1}\right)>0$, $g_2$ and $g_3$ can be interpreted as functions on the upper half-plane $ℍ:=\{z\in ℂ:\mathrm{Im}(z)>0\}$.

Let $\tau =\frac{{\omega }_2}{{\omega }_1}$. One has: $$g_2(1,\tau )={\omega }_1^4{g}_2({\omega }_1,{\omega }_2),$$ $$g_3(1,\tau )={\omega }_1^6{g}_3({\omega }_1,{\omega }_2).$$ That means g₂ and g₃ are only scaled by doing this. Set $$g_2(\tau ):=g_2(1,\tau )$$ and $$g_3(\tau ):=g_3(1,\tau ).$$ As functions of $\tau \in ℍ$, $g_2$ and $g_3$ are so called modular forms.

The Fourier series for $g_2$ and $g_3$ are given as follows: $$g_2(\tau )=\frac{4}{3}{\pi }^4[1+240\sum _{k=1}^{\mathrm{\infty }}{\sigma }_3(k)q^{2k}]$$ $$g_3(\tau )=\frac{8}{27}{\pi }^6[1-504\sum _{k=1}^{\mathrm{\infty }}{\sigma }_5(k)q^{2k}]$$ where $${\sigma }_m(k):=\sum _{d\mid k}{d}^m$$ is the divisor function and $q=e^{\pi i\tau }$ is the nome.

Modular discriminant

The real part of the discriminant as a function of the square of the nome q on the unit disk.

The modular discriminant $\Delta$ is defined as the discriminant of the characteristic polynomial of the differential equation $$\mathrm{\wp }{\text{'}}^2(z)=4{\mathrm{\wp }}^3(z)-g_2\mathrm{\wp }(z)-g_3$$ as follows: $$\Delta =g_2^3-27g_3^2.$$ The discriminant is a modular form of weight $12$. That is, under the action of the modular group, it transforms as $$\Delta \left(\frac{a\tau +b}{c\tau +d}\right)={(c\tau +d)}^{12}\Delta (\tau )$$ where $a,b,d,c\in \mathbb{Z}$ with $ad-bc=1$.

Note that $\Delta =(2\pi {)}^{12}{\eta }^{24}$ where $\eta$ is the Dedekind eta function.

For the Fourier coefficients of $\Delta$, see Ramanujan tau function.

The constants e₁, e₂ and e₃

$e_1$, $e_2$ and $e_3$ are usually used to denote the values of the $\mathrm{\wp }$-function at the half-periods. $$e_1\equiv \mathrm{\wp }\left(\frac{{\omega }_1}{2}\right)$$ $$e_2\equiv \mathrm{\wp }\left(\frac{{\omega }_2}{2}\right)$$ $$e_3\equiv \mathrm{\wp }\left(\frac{{\omega }_1+{\omega }_2}{2}\right)$$ They are pairwise distinct and only depend on the lattice $\Lambda$ and not on its generators.

$e_1$, $e_2$ and $e_3$ are the roots of the cubic polynomial $4\mathrm{\wp }(z{)}^3-g_2\mathrm{\wp }(z)-g_3$ and are related by the equation: $$e_1+e_2+e_3=0.$$ Because those roots are distinct the discriminant $\Delta$ does not vanish on the upper half plane. Now we can rewrite the differential equation: $$\mathrm{\wp }{\text{'}}^2(z)=4(\mathrm{\wp }(z)-e_1)(\mathrm{\wp }(z)-e_2)(\mathrm{\wp }(z)-e_3).$$ That means the half-periods are zeros of ${\mathrm{\wp }}^{\prime }$.

The invariants $g_2$ and $g_3$ can be expressed in terms of these constants in the following way: $$g_2=-4(e_1{e}_2+e_1{e}_3+e_2{e}_3)$$ $$g_3=4e_1{e}_2{e}_3$$ $e_1$, $e_2$ and $e_3$ are related to the modular lambda function: $$\lambda (\tau )=\frac{e_3-e_2}{e_1-e_2},\tau =\frac{{\omega }_2}{{\omega }_1}.$$

Relation to Jacobi's elliptic functions

For numerical work, it is often convenient to calculate the Weierstrass elliptic function in terms of Jacobi's elliptic functions.

The basic relations are: $$\mathrm{\wp }(z)=e_3+\frac{e_1-e_3}{{\mathrm{sn}}^{2}w}=e_2+(e_1-e_3)\frac{{\mathrm{dn}}^{2}w}{{\mathrm{sn}}^{2}w}=e_1+(e_1-e_3)\frac{{\mathrm{cn}}^{2}w}{{\mathrm{sn}}^{2}w}$$ where $e_1,e_2$ and $e_3$ are the three roots described above and where the modulus k of the Jacobi functions equals $$k=\sqrt{\frac{e_2-e_3}{e_1-e_3}}$$ and their argument w equals $$w=z\sqrt{e_1-e_3}.$$

Relation to Jacobi's theta functions

The function $\mathrm{\wp }(z,\tau )=\mathrm{\wp }(z,1,{\omega }_2/{\omega }_1)$ can be represented by Jacobi's theta functions: $$\mathrm{\wp }(z,\tau )={(\pi {\theta }_2(0,q){\theta }_3(0,q)\frac{{\theta }_4(\pi z,q)}{{\theta }_1(\pi z,q)})}^2-\frac{{\pi }^2}{3}({\theta }_2^4(0,q)+{\theta }_3^4(0,q))$$ where $q=e^{\pi i\tau }$ is the nome and $\tau$ is the period ratio $(\tau \in ℍ)$. This also provides a very rapid algorithm for computing $\mathrm{\wp }(z,\tau )$.

Relation to elliptic curves

See also: Elliptic curve § Elliptic curves over the complex numbers

Consider the embedding of the cubic curve in the complex projective plane

${\overline{C}}_{g_2,g_3}^{ℂ}=\{(x,y)\in {ℂ}^2:y^2=4x^3-g_{2}x-g_3\}\cup \{O\}\subset {ℂ}^2\cup {\mathbb{P}}_1(ℂ)={\mathbb{P}}_2(ℂ).$

where $O$ is a point lying on the line at infinity ${\mathbb{P}}_1(ℂ)$. For this cubic there exists no rational parameterization, if $\Delta \ne 0$. In this case it is also called an elliptic curve. Nevertheless there is a parameterization in homogeneous coordinates that uses the $\mathrm{\wp }$-function and its derivative ${\mathrm{\wp }}^{\prime }$:

$\phi (\mathrm{\wp },{\mathrm{\wp }}^{\prime }):ℂ/\Lambda \to {\overline{C}}_{g_2,g_3}^{ℂ},z\mapsto \{\begin{array}{cc}[\mathrm{\wp }(z):{\mathrm{\wp }}^{\prime }(z):1]& z\notin \Lambda \\ [0:1:0]& z\in \Lambda \end{array}$

Now the map $\phi$ is bijective and parameterizes the elliptic curve ${\overline{C}}_{g_2,g_3}^{ℂ}$.

$ℂ/\Lambda$ is an abelian group and a topological space, equipped with the quotient topology.

It can be shown that every Weierstrass cubic is given in such a way. That is to say that for every pair $g_2,g_3\in ℂ$ with $\Delta =g_2^3-27g_3^2\ne 0$ there exists a lattice $\mathbb{Z}{\omega }_1+\mathbb{Z}{\omega }_2$, such that

$g_2=g_2({\omega }_1,{\omega }_2)$ and $g_3=g_3({\omega }_1,{\omega }_2)$.

The statement that elliptic curves over $\mathbb{Q}$ can be parameterized over $\mathbb{Q}$, is known as the modularity theorem. This is an important theorem in number theory. It was part of Andrew Wiles' proof (1995) of Fermat's Last Theorem.

Addition theorem

The addition theorem states that if $z,w,$ and $z+w$ do not belong to $\Lambda$, then $$\det \left[\begin{array}{ccc}1& \mathrm{\wp }(z)& {\mathrm{\wp }}^{\prime }(z)\\ 1& \mathrm{\wp }(w)& {\mathrm{\wp }}^{\prime }(w)\\ 1& \mathrm{\wp }(z+w)& -{\mathrm{\wp }}^{\prime }(z+w)\end{array}\right]=0.$$ This states that the points $P=(\mathrm{\wp }(z),{\mathrm{\wp }}^{\prime }(z)),$ $Q=(\mathrm{\wp }(w),{\mathrm{\wp }}^{\prime }(w)),$ and $R=(\mathrm{\wp }(z+w),-{\mathrm{\wp }}^{\prime }(z+w))$ are collinear, the geometric form of the group law of an elliptic curve.

This can be proven by considering constants $A,B$ such that $${\mathrm{\wp }}^{\prime }(z)=A\mathrm{\wp }(z)+B,{\mathrm{\wp }}^{\prime }(w)=A\mathrm{\wp }(w)+B.$$ Then the elliptic function $${\mathrm{\wp }}^{\prime }(\zeta )-A\mathrm{\wp }(\zeta )-B$$ has a pole of order three at zero, and therefore three zeros whose sum belongs to $\Lambda$. Two of the zeros are $z$ and $w$, and thus the third is congruent to $-z-w$.

Alternative form

The addition theorem can be put into the alternative form, for $z,w,z-w,z+w\in ̸\Lambda$: $$\mathrm{\wp }(z+w)=\frac{1}{4}{\left[\frac{{\mathrm{\wp }}^{\prime }(z)-{\mathrm{\wp }}^{\prime }(w)}{\mathrm{\wp }(z)-\mathrm{\wp }(w)}\right]}^2-\mathrm{\wp }(z)-\mathrm{\wp }(w).$$

As well as the duplication formula: $$\mathrm{\wp }(2z)=\frac{1}{4}{\left[\frac{{\mathrm{\wp }}^{″}(z)}{{\mathrm{\wp }}^{\prime }(z)}\right]}^2-2\mathrm{\wp }(z).$$

Proofs

This can be proven from the addition theorem shown above. The points $P=(\mathrm{\wp }(u),{\mathrm{\wp }}^{\prime }(u)),Q=(\mathrm{\wp }(v),{\mathrm{\wp }}^{\prime }(v)),$ and $R=(\mathrm{\wp }(u+v),-{\mathrm{\wp }}^{\prime }(u+v))$ are collinear and lie on the curve $y^2=4x^3-g_{2}x-g_3$. The slope of that line is $$m=\frac{y_P-y_Q}{x_P-x_Q}=\frac{{\mathrm{\wp }}^{\prime }(u)-{\mathrm{\wp }}^{\prime }(v)}{\mathrm{\wp }(u)-\mathrm{\wp }(v)}.$$ So $x=x_P=\mathrm{\wp }(u)$, $x=x_Q=\mathrm{\wp }(v)$, and $x=x_R=\mathrm{\wp }(u+v)$ all satisfy a cubic $$(mx+q{)}^2=4x^3-g_{2}x-g_3,$$ where $q$ is a constant. This becomes $$4x^3-m^2{x}^2-(2mq+g_2)x-g_3-q^2=0.$$ Thus $x_P+x_Q+x_R=\frac{m^2}{4}$ which provides the wanted formula $\mathrm{\wp }(u+v)+\mathrm{\wp }(u)+\mathrm{\wp }(v)=\frac{1}{4}{\left[\frac{{\mathrm{\wp }}^{\prime }(u)-{\mathrm{\wp }}^{\prime }(v)}{\mathrm{\wp }(u)-\mathrm{\wp }(v)}\right]}^2.$

A direct proof is as follows. Any elliptic function $f$ can be expressed as: $$f(u)=c\prod _{i=1}^n\frac{\sigma (u-a_i)}{\sigma (u-b_i)}c\in ℂ$$ where $\sigma$ is the Weierstrass sigma function and $a_i,b_i$ are the respective zeros and poles in the period parallelogram. Considering the function $\mathrm{\wp }(u)-\mathrm{\wp }(v)$ as a function of $u$, we have $$\mathrm{\wp }(u)-\mathrm{\wp }(v)=c\frac{\sigma (u+v)\sigma (u-v)}{\sigma (u{)}^2}.$$ Multiplying both sides by $u^2$ and letting $u\to 0$, we have $1=-c\sigma (v{)}^2$, so $c=-\frac{1}{\sigma (v{)}^2}⟹\mathrm{\wp }(u)-\mathrm{\wp }(v)=-\frac{\sigma (u+v)\sigma (u-v)}{\sigma (u{)}^2\sigma (v{)}^2}.$

By definition the Weierstrass zeta function: $\frac{d}{dz}\ln \sigma (z)=\zeta (z)$ therefore we logarithmically differentiate both sides with respect to $u$ obtaining: $$\frac{{\mathrm{\wp }}^{\prime }(u)}{\mathrm{\wp }(u)-\mathrm{\wp }(v)}=\zeta (u+v)-2\zeta (u)-\zeta (u-v)$$ Once again by definition ${\zeta }^{\prime }(z)=-\mathrm{\wp }(z)$ thus by differentiating once more on both sides and rearranging the terms we obtain $$-\mathrm{\wp }(u+v)=-\mathrm{\wp }(u)+\frac{1}{2}\frac{{\mathrm{\wp }}^{″}(v)[\mathrm{\wp }(u)-\mathrm{\wp }(v)]-{\mathrm{\wp }}^{\prime }(u)[{\mathrm{\wp }}^{\prime }(u)-{\mathrm{\wp }}^{\prime }(v)]}{[\mathrm{\wp }(u)-\mathrm{\wp }(v){]}^2}$$ Knowing that ${\mathrm{\wp }}^{″}$ has the following differential equation $2{\mathrm{\wp }}^{″}=12{\mathrm{\wp }}^2-g_2$ and rearranging the terms one gets the wanted formula $$\mathrm{\wp }(u+v)=\frac{1}{4}{\left[\frac{{\mathrm{\wp }}^{\prime }(u)-{\mathrm{\wp }}^{\prime }(v)}{\mathrm{\wp }(u)-\mathrm{\wp }(v)}\right]}^2-\mathrm{\wp }(u)-\mathrm{\wp }(v).$$

Typography

The Weierstrass's elliptic function is usually written with a rather special, lower case script letter ℘, which was Weierstrass's own notation introduced in his lectures of 1862–1863. It should not be confused with the normal mathematical script letters P: 𝒫 and 𝓅.

In computing, the letter ℘ is available as \wp in TeX. In Unicode the code point is U+2118 ℘ <reserved-2118>, with the more correct alias weierstrass elliptic function. In HTML, it can be escaped as ℘ or ℘.

Character information

Preview	℘
Unicode name	SCRIPT CAPITAL P / WEIERSTRASS ELLIPTIC FUNCTION
Encodings	decimal	hex
Unicode	8472	U+2118
UTF-8	226 132 152	E2 84 98
Numeric character reference	℘	℘
Named character reference	℘, ℘

See also

• Weierstrass functions
• Jacobi elliptic functions
• Lemniscate elliptic functions

References

• N. I. Akhiezer, Elements of the Theory of Elliptic Functions, (1970) Moscow, translated into English as AMS Translations of Mathematical Monographs Volume 79 (1990) AMS, Rhode Island
• Tom M. Apostol, Modular Functions and Dirichlet Series in Number Theory, Second Edition (1990), Springer, New York (See chapter 1.)
• K. Chandrasekharan, Elliptic functions (1980), Springer-Verlag
• Konrad Knopp, Funktionentheorie II (1947), Dover Publications; Republished in English translation as Theory of Functions (1996), Dover Publications
• Serge Lang, Elliptic Functions (1973), Addison-Wesley,
• E. T. Whittaker and G. N. Watson, A Course of Modern Analysis, Cambridge University Press, 1952, chapters 20 and 21

External links

• Weierstrass's elliptic functions on Mathworld.
• Chapter 23, Weierstrass Elliptic and Modular Functions in DLMF (Digital Library of Mathematical Functions) by W. P. Reinhardt and P. L. Walker.
• Weierstrass P function and its derivative implemented in C by David Dumas