2011年2月12日[应该在2010-2011年间]，我在网上第一次读到这篇文章，这是一篇介绍椭圆积分和椭圆函数极为出色的硕士学位论文。我个人以为，大学数学系本科阶段的分析学学习，应该掌握椭圆函数论的知识，正如中学阶段应该掌握初等函数的知识一样。它是进一步学习其他数学分支的基础。
 
 htt
 p
 ://www.docin.com/p-49451820.html
 
 
 http://d.g.wanfangdata.com.cn/Thesis_Y1519793.aspx
 
 [摘要]：本文以原始文献分析法为研究方法,以研读阿贝尔的法文论文为主,着重考察18-19世纪椭圆积分向椭圆函数转变的过程。论文分为三章。第一章在18世纪微积分进一步发展的背景下,分析欧拉、勒让德、伯努利家族等数学家对椭圆积分进行的研究工作,指出勒让德将椭圆积分发展成三种标准表达形式以及欧拉、伯努利家族得到的椭圆积分加法定理,成为19世纪阿贝尔等人进一步研究椭圆积分和发展椭圆函数论的理论基础。 第二章详细梳理阿贝尔建立和发展椭圆函数论的思想和方法,认为他在试图对椭圆积分进行求解无果的前提下,转变思想,采取与三角函数类比的方式,从反函数的角度入手,由椭圆积分诱导出椭圆函数的概念。另一方面,阿贝尔对R的阶大于4的进行积分,提出了比椭圆积分更广泛的阿贝尔积分,同时进一步推广欧拉的椭圆积分加法定理,得到了现在被称为“阿贝尔大定理”的加法定理,从而创建了定性分析阿贝尔积分的关键性定理. 第三章通过分析雅可比、外尔斯特拉斯、黎曼等人在椭圆函数论方面的成绩,指出由于受到阿贝尔反演思想的影响,上述数学家给出了更一般、简捷的椭圆函数理论形式,并且反演了超椭圆积分,推动了近代复分析函数领域的成型与发展. 文章回溯18-19世纪椭圆函数发展历程,指出阿贝尔在欧拉、勒让德、伯努利家族等数学家工作的基础上,创立椭圆函数理论,由此影响了19世纪特殊函数论的发展,对复变函数论,甚至是近代代数几何学的发展产生深远影响.
 
 [关键词]：椭圆积分 椭圆函数 阿贝尔 函数 积分
 
 中文摘要
 
 绪论
 
 一 椭圆积分的产生及发展
 
 （一） 18世纪的积分技术
 
 （二） 椭圆积分
 
 （三） 勒让德对椭圆积分的贡献
 
 二 阿贝尔在椭圆函数论方面的奠基性工作
 
 （一） 关于阿贝尔
 
 （二） 阿贝尔对椭圆积分的研究
 
 （三） 超椭圆积分的加法定理
 
 （四） 阿贝尔积分
 
 （五） 由反演思想得到的一些性质
 
 （六） 椭圆函数变换的问题
 
 三 阿贝尔思想影响下19世纪椭圆函数论的发展
 
 （一） 雅可比椭圆函数理论
 
 （二） 外尔斯特拉斯的椭圆函数论
 
 （三） 其他数学家的椭圆函数论
 
 结语
 
 部分引文对照
 
 参考文献
 
 致谢 
 
 
 
 
绪论（节选）
 国内单独对阿贝尔的椭圆函数的发展史的研究由于受语言的困难、文化背景的差异、资料的限制等，目前还很少。涉及这方面的研究主要是椭圆函数和椭圆积分解的问题，这是纯粹的数学问题研究。胡作玄的《近代数学史》中对椭圆函数和椭圆积分、阿贝尔函数和阿贝尔积分用归纳总结的方法作了一些介绍，在其著作中涉及了18、19世纪阿贝尔、雅可比、外尔斯特拉斯、黎曼等椭圆函数理念研究者的研究成果，但对其理论的建立方面介绍较少；李文林《数学珍宝》节选翻译了阿贝尔函数和阿贝尔积分论著和雅可比关于雅可比θ函数理论的论文；张致中翻译的《高级超越函数》比较详细的介绍了雅可比椭圆函数和外尔斯特拉斯椭圆函数理论；王竹溪、郭敦仁编写的《特殊函数论》对外氏和雅氏椭圆函数用通俗的数学语言作了介绍。在解恩泽，徐本顺主编的《世界数学家思想方法》对阿贝尔、外尔斯特拉斯等19世纪椭圆函数论研究者建立椭圆函数理论的思想方法作了介绍。梁宗巨、王青建、孙宏安在《世界数学通史》（下）等著作中，对在18、19世纪研究椭圆函数的历史背景、研究者概况作了介绍。吴文俊主编《世界著名科学家传（数学家）》有阿贝尔的专门传记。杜瑞芝主编的《坎坷奇星—阿贝尔》对阿贝尔生平作了详尽介绍。
 国外对阿贝尔的研究很多。可资利用的原始文献有《阿贝尔全集》（法文版）收集了阿贝尔一生中关于椭圆函数论的所有论著（除一篇之外）。关于阿贝尔的传记，B.Holmboe在编辑《阿贝尔全集》时，就在“前言”中做专门介绍，1834年，G.Libri撰写了阿贝尔传，1884年C.A.Bjerknes也撰写了阿贝尔传，最近（2000年），A.Stubhaug出版了一本专门研究阿贝尔的新书，《阿贝尔与他的时代》。吉利斯皮主编的《科学史传记词典》（Dictionary of Scientific Biography）中也有阿贝尔传。一般数学史论著中都对阿贝尔的工作有一定的介绍。其中，[美]M·克莱因的《古今数学思想》比较系统的介绍了从椭圆积分到椭圆函数的发展历程，对阿贝尔和雅克比研究椭圆函数的反演思想作了介绍。
 
 
 
 
 
（三） 勒让德对椭圆积分的贡献
 椭圆积分的历史起点一般公认为由意大利数学家法拉诺（Fagnano,1713-1798）开始研究，1751年发现了双纽线的弧长的积分F(u)=∫[0,u]dw/[sqrt(1-w^4)]的加倍问题而推导出这个公式的。他证明，如果F(w)=2F(u)，那么w^4=4u^2(1-u^4)/[(1+u^4)^2][按：左边应该是w^2吧？]，即w与u之间存在代数关系。欧拉于1761年把倍弧长公式推广成双纽线积分的加法定理，即若F(u)+F(v)=F(w)，则w=[usqrt(1-v^4)+vsqrt(1-u^4)]/(1+u^2v^2)，显然当u=v时，即得出法拉诺关系，欧拉是在1751年12月23日知道法拉诺关系的。后来，雅可比把1751年12月23日这一天称为“椭圆函数的生日。”
 勒让德主要证明了一般椭圆积分∫P(x)/sqrt(R(x))dx（其中P(x)是x的任一有理函数，而R(x)是通常一般的四次多项式），能化成三种类型：
 F=∫dx/[sqrt(1-x^2)sqrt(1-k^2x^2)]
 E=∫x^2dx/[sqrt(1-x^2)sqrt(1-k^2x^2)]
 N=∫dx/[(x-a)sqrt(1-x^2)sqrt(1-k^2x^2)]
 分别称为第一、第二、第三类椭圆积分，显然双纽线积分是第一类椭圆积分。经过变换，这三类积分可化为：
 
 
F(φ,k)=∫[0,φ]dφ/[sqrt(1-k^2sin^2φ)],0<k<1,
 E(φ,k)=∫[0,φ][sqrt(1-k^2sin^2φ)]dφ,0<k<1,
 Π(φ,n,k)=∫[0,φ]dφ/[(1+nsin^2φ)sqrt(1-k^2sin^2φ)],0<k<1,
 这种分法至今仍在使用，其中k称为模。而k'=sqrt(1-k^2)被称为余模。
 当φ=pi/2时，积分F=F(k)=F(pi/2,k),E=E(k)=E(pi/2,k)分别称为第一类和第二类完全椭圆积分。
 通过定义F'=F'(k)=F(pi/2,k'),E'=E'(k)=E(pi/2,k')，勒让德发现它们有如下的关系（1825）：
 FE'+F'E-FF'=pi/2。
 勒让德在他的书中得出了一系列加法公式及变换公式，以及不同的参数n的第三类积分之间的关系。在《椭圆函数研究》第二卷（1826）中，勒让德发表了一个椭圆积分表，它也是今天同类表的基础。
 
 
 
 
 
（三） 超椭圆积分的加法定理
 
 
1826年，阿贝尔撰写了论文《关于一类极为广泛的超越函数的一个一般性质》，在这篇论文中，阿贝尔对R的阶大于4的pdx/sqrt(R)进行积分，提出了比椭圆积分更广泛的阿贝尔积分，并对欧拉的椭圆积分加法定理做了推广，得到了现在称作“阿贝尔大定理”的加法定理。这是阿贝尔积分的一条关键性定理。但是，审查人柯西和勒让德没能使它及时发表，直到1841年重新找到原稿后才面世。阿贝尔在这篇论文中没有对椭圆积分的计算再作处理，而是转向对其性质的研究，他研究椭圆积分的思想发生了转变。
 以往所研究的椭圆积分∫R(x,y)dx中x,y满足y^2=P(x)，且P(x)是多项式。阿贝尔将x,y的关系放宽为x,y由代数方程f(x,y)=0所确定。这样的新积分，显然是椭圆积分的推广（一般化），后人称之为阿贝尔积分。
 阿贝尔对“阿贝尔大定理”的叙述很完整：“我们总能用若干预先确定的其变元是给定函数之变元的代数函数的函数之相思和去表示那特殊个数的函数和，而这特殊中的每个函数被乘上一个有理数并且其变元是任意的。”但证明却不是很严格。
 当f=y^2-P(x)，P(x)是一个6次多项式，而p=(n-2)/2=2时，∫R(x,y)dx是超椭圆积分。
 在后继的论文中，阿贝尔借助于这个定理，发现了椭圆函数的双周期性，从而奠定了椭圆曲线（它们都可以表示成平面中的三次曲线）的理论基础。利用这种性质还可以对椭圆函数做出如下定义：只有极点的双周期解析函数是椭圆函数。
 
 
（五） 由反演思想得到的一些性质
 在《关于椭圆函数的研究》（1827）这篇论文中，他借助于椭圆积分的反函数把椭圆积分的理论归结为椭圆函数的理论。具体地说，阿贝尔所考察的椭圆积分是这样的一些积分，其中被积函数是三次或四次多项式的平方根的有理函数。用现代符号表示，这些积分之中，重要的是函数u(s)=∫[0->s]dx/[sqrt(1-x^2)sqrt(1-k^2x^2)]，其反函数s(u)同样起着重要作用，它恰好是椭圆函数snu。我们使用符号snu是为了表示它是普通的正弦函数的推广。在最基本的情形，即k=0，我们可分别得到u(s)=arcsins=∫[0->s]dx/[sqrt(1-x^2)]和s(u)=sinu。
 u=F(φ,R)=∫[0,φ]dφ/[sqrt(1-R^2sin^2φ)],0<R<1,
 u=∫dx/[sqrt(1-x^2)sqrt(1-R^2x^2)]
 阿贝尔定义了fu=cnu=sqrt(1-sn^2u),Fu=dnu=sqrt(1-R^2sn^2u)，然后发现了snu,cnu,dnu的若干基本性质：
 sn'u=cnudnu,cn'u=-snudnu,dn'u=-k^2snucnu
 还有加法公式：
 sn(u+v)=(snucnvdnv+snvcnudnu)/△
 cn(u+v)=(cnucnv-snusnvdnudnv)/△
 dn(u+v)=(dnudnv-k^2snusnvcnucnv)/△
 这里△=1-R^2sn^2usn^2v。由此，把研究椭圆积分的性质转化为对一类特殊函数椭圆函数的考察，使椭圆函数论成为一种经典的分析工具。
 
 
阿贝尔在下式中定义的量K起着三角函数中pi的作用，
 K=∫[0->1]dx/[sqrt(1-x^2)sqrt(1-k^2x^2)]=∫[0,pi/2]dφ/[sqrt(1-k^2sin^2φ)]=F(k,pi/2),0<k<1,
 与K联系着的是超越量K'，它作为k'的函数相当于K作为k的函数，其中k'由k^2+k'^2=1定义，0<k<1。
 这样，阿贝尔发现，snz的周期是4K及2iK'；cnz的周期是4K及2K+2iK'；dnz的周期是2K及4iK'。这两个周期的比为非实数，因而这些椭圆函数是双周期的，这是伟大的发现之一。
 [
 20110621：利用JT1~JT4计算sn,cn,dn
 K(k=0.985171431009416)=3.16510345444743
 z=0.2
 t=0.5i
 K'=0.5K
 k=kq(q=exp(-0.5*pi)=0.207879576350762)=0.9851714
 k'=ckq(q=exp(-0.5*pi)=0.207879576350762)=0.17157
 错误的：
 u=2Kz/pi=2*3.16510345444743*0.2/pi=0.402993488138034
 正确的：
 u=2Kz=2*3.16510345444743*0.2=1.26604138177897=1.266041
 sn(u=1.266041;k=0.9851714)=0.8564342
 cn(u=1.266041;k=0.9851714)=0.5162562
 dn(u=1.266041;k=0.9851714)=0.5367607
 u=arcsn(sn(u;k);k)=arcsn(0.8564342;0.9851714)=1.266041
 fcomplex retqk=detail::q_from_k(fcomplex(0.985171431,0));
 cout<<retqk<<endl;//(0.20788,0)
 cout<<detail::kq(0.207879576350762)<<endl;//0.985171
 cout<<detail::ckq(0.207879576350762)<<endl;//0.171573
 q=q_from_k(k=0.985171431)=0.20788
 k=kq(q=0.207879576350762)=0.985171431009416
 k'=ckq(q=0.207879576350762)=0.17157287525381
 tau=0.5i=>q=exp(i*pi*tau)=0.207879576350762
    cout<<detail::theta_00(fcomplex(0.2,0),fcomplex(0,0.5))<<endl;//(1.12545,0)
    cout<<detail::theta_00(fcomplex(0.2,0.2),fcomplex(0,0.5))<<endl;//(1.22519,-0.651829)
    cout<<detail::theta_00(fcomplex(0,0),fcomplex(0,0.5))<<endl;//(1.4195,0)=JT3(z=0,tau=0.5i)=1.41949548808377
    cout<<detail::theta_00(fcomplex(3.3416,0),fcomplex(0,0.5))<<endl;//(0.772187,0)
    cout<<detail::theta_00(fcomplex(1.2,0),fcomplex(0,0.5))<<endl;//(1.12545,0)
    cout<<detail::theta_01(fcomplex(0.2,0),fcomplex(0,0.5))<<endl;//(0.868503,0)
    cout<<detail::theta_01(fcomplex(0,0),fcomplex(0,0.5))<<endl;//(0.587974,0)=JT4(z=0,tau=0.5i)=0.587974282891712
    cout<<detail::theta_01(fcomplex(1.5708,0),fcomplex(0,0.5))<<endl;//(1.37765,0)
    cout<<detail::theta_01(fcomplex(1.6708,0),fcomplex(0,0.5))<<endl;//(1.19643,0)
    cout<<detail::theta_01(fcomplex(1.7708,0),fcomplex(0,0.5))<<endl;//(0.942211,0)
    cout<<detail::theta_01(fcomplex(1.8708,0),fcomplex(0,0.5))<<endl;//(0.713678,0)
    cout<<detail::theta_01(fcomplex(1,0),fcomplex(0,0.5))<<endl;//(0.587974,0)
    cout<<detail::theta_01(fcomplex(1.1,0),fcomplex(0,0.5))<<endl;//(0.664798,0)
    cout<<detail::theta_01(fcomplex(1.2,0),fcomplex(0,0.5))<<endl;//(0.868503,0)
    cout<<detail::theta_01(fcomplex(1.3,0),fcomplex(0,0.5))<<endl;//(1.12545,0)
    cout<<detail::theta_01(fcomplex(0.5,0),fcomplex(0,0.5))<<endl;//(1.4195,0)
    cout<<detail::theta_01(fcomplex(0.6,0),fcomplex(0,0.5))<<endl;//(1.33751,0)
    cout<<detail::theta_01(fcomplex(0.7,0),fcomplex(0,0.5))<<endl;//(1.12545,0)
    cout<<detail::theta_01(fcomplex(0.8,0),fcomplex(0,0.5))<<endl;//(0.868503,0)
 JT3(z=0.2,tau=0.5i)=1.12545388495798
 JT3(z=0.2+0.2i,tau=0.5i)=1.22519068532128-0.651828923108493i
 JT3(z=0,tau=0.5i)=1.41949548808377
 验证了公式：JT4(z,tau)=JT3(z+1/2,tau)，JT3(z+1,tau)=JT3(z,tau)
 而不是公式：JT4(z,tau)=JT3(z+pi/2,tau)，JT3(z+pi,tau)=JT3(z,tau)
    cout<<detail::theta_3(fcomplex(0.2,0)*M_PI,fcomplex(0.207879576350762,0))<<endl;//(1.12545,0)=JT3(z=0.2,tau=0.5i)=1.12545388495798
    cout<<detail::theta_3(fcomplex(0.2,0),fcomplex(0.207879576350762,0))<<endl;//(1.38554,0)
    cout<<detail::theta_3(fcomplex(0,0),fcomplex(0.207879576350762,0))<<endl;//(1.4195,0)，JT3(z=0,tau=0.5i)=1.41949548808377
    cout<<detail::theta_4(fcomplex(0.2,0),fcomplex(0.207879576350762,0))<<endl;//(0.619662,0)
    cout<<detail::theta_4(fcomplex(0,0),fcomplex(0.207879576350762,0))<<endl;//(0.587974,0)，JT4(z=0,tau=0.5i)=0.587974282891712
 
 
   cout<<detail::theta_2(fcomplex(0.2,0)*M_PI,fcomplex(0.207879576350762,0))<<endl;//(1.07441,0)
    cout<<detail::theta_10(fcomplex(0.2,0)/M_PI,fcomplex(0,0.5))<<endl;//(1.37177,0)
    cout<<detail::theta_10(fcomplex(0.2,0),fcomplex(0,0.5))<<endl;//(1.07441,-1.11022e-016)
    cout<<detail::theta_10(fcomplex(0,0),fcomplex(0,0.5))<<endl;//(1.40893,0)，某处计算结果：JT2(z=0,tau=0.5i)=0.70446581836561
 
 
   cout<<detail::theta_2(fcomplex(0.2,0),fcomplex(0.207879576350762,0))<<endl;//(1.37177,0)
    cout<<detail::theta_2(fcomplex(0,0),fcomplex(0.207879576350762,0))<<endl;//(1.40893,0)
 
 
   cout<<detail::theta_11(fcomplex(0.2,0),fcomplex(0,0.5))<<endl;//(-0.73828,-5.55112e-017)，某处计算结果：
 JT1(z=0.2,tau=0.5i)=0.369140086660667，JT1无负号，对应的sn的Theta函数表达式中也无负号
    cout<<detail::theta_11(fcomplex(0,0),fcomplex(0,0.5))<<endl;//(7.20051e-017,-4.40904e-033)=JT1(z=0,tau=0.5i)=0
    cout<<detail::theta_1(fcomplex(0.2,0),fcomplex(0.207879576350762,0))<<endl;//(0.235436,0)
    cout<<detail::theta_1(fcomplex(0,0),fcomplex(0.207879576350762,0))<<endl;//(0,0)
 sn(u;k)=JT3(0;tau)*JT1(z;tau)/(JT2(0;tau)*JT4(z;tau))=1.41949548808377*0.369140086660667/(0.70446581836561*0.868502943433161)=0.856434209027962
 cn(u;k)=JT4(0;tau)*JT2(z;tau)/(JT2(0;tau)*JT4(z;tau))=0.587974282891712*0.537202659820004/(0.70446581836561*0.868502943433161)=0.516256182148607
 dn(u;k)=JT4(0;tau)*JT3(z;tau)/(JT3(0;tau)*JT4(z;tau))=0.587974282891712*1.12545388495798/(1.41949548808377*0.868502943433161)=0.516256182148607=0.536760717392963
 
 证明：如果R(x,y)为x,y的有理函数，则∫R(x,sqrt(ax+b))dx和∫R(x,sqrt(ax^2+bx+c))dx为初等函数。
 例如：∫(1/sqrt(1-x^2))dx=arcsinx----multivalued function但对于∫R(x,sqrt(a_3x^3+a_2x^2+a_1x+a_0))dx和∫R(x,sqrt(a_4x^4+a_3x^3+a_2x^2+a_1x+a_0))dx有100年左右找不到什么结果，那时人们把这种积分称作椭圆积分。
 勒让德第一类[不完全]椭圆积分有两种形式：u=F(a,k)和u=F(x=sina;k)，其中，自变量a为幅度的弧度制，自变量x为三角函数值。
 其反函数分别为雅可比幅度函数a=am(u,k)和雅可比椭圆正弦函数x=sn(u;k)=sin(am(u,k))
 u=arcsn(x;k)与x=sn(u;k)的反演关系就是严格的反函数关系。
 u=F(a,k)与x=sn(u;k)存在反演关系：sn(F(a,k);k)=sin(a)。
 u=arcsn(x;k=0)=arcsin(x)<=>x=sn(u;k=0)=sin(u)
 u=arcsn(x=1;k=0)=pi/2<=>x=sn(u=pi/2;k=0)=1
 u=arcsn(x=sqrt(2)/2;k=0)=pi/4<=>x=sn(u=pi/4;k=0)=sqrt(2)/2
 u=arcsn(x=1;k=1)=13.8191004<=>x=sn(u=13.8191004;k=1)=1
 u=arcsn(x=1;k=sqrt(2)/2)=1.8539468=sqrt(2)(ω/2)<=>x=sn(u=sqrt(2)(ω/2);k=sqrt(2)/2)=1
 u=arcsn(x=1;k=0.5)=1.6857504<=>x=sn(u=1.6857504;k=0.5)=1
 u=F(a=0.698132,k=0.5)=0.711647<=>a=am(u=0.711647,k=0.5)=0.698132
 sn(u=0.711647;k=0.5)=sin(0.698132)=0.6427874
 cn(u=0.711647;k=0.5)=cos(0.698132)=0.7660446
 dn(u=0.711647;k=0.5)=0.9469457
 
 
u=F(a=1.396263,k=0.5)=1.48455<=>a=am(u=1.48455,k=0.5)=1.396263
 sn(u=1.48455;k=0.5)=sin(1.396263)=0.9848071
 cn(u=1.48455;k=0.5)=cos(1.396263)=0.1736521
 dn(u=1.48455;k=0.5)=0.870367
 
 
u=F(a=0,k=0.5)=0<=>a=am(u=0,k=0.5)=0
 u=F(a=0,k=1)=0<=>a=am(u=0,k=1)=0
 sn(u=0;k=0.5)=sin(0)=0
 cn(u=0;k=0.5)=cos(0)=1
 dn(u=0;k=0.5)=1
 
 
复数1+2i的正弦为：3.16577851321617+1.95960104142161i
 复椭圆正弦sn(1+2i;k=0)=3.165778+1.95901i
 sin(0.711647)=0.653081908381731
 sn(u=0.711647;k=0)=0.6530819
 cn(u=0.711647;k=0)=0.7572873
 dn(u=0.711647;k=0)=1
 sn(u=0.711647;k=0.5)=0.6427874
 cn(u=0.711647;k=0.5)=0.7660446
 dn(u=0.711647;k=0.5)=0.9469457
 sn(u=0.711647+i;k=0.5)=0.8742257+0.6536986i
 sn(u=0.5+i;k=0.5)=0.6461386+0.7702646i
 sn(u=0.5+0.5i;k=0.5)=0.5193378+0.4214193i
 
 

 2011.6.11
 由虚变换和加法公式计算复椭圆函数对应的调和函数：
 sn(x+yi;k)=?
 sn(x;k)
 sn(yi;k)=isnh(y;k)
 cn(x;k)
 cn(yi;k)=cnh(y;k)
 dn(x;k)
 dn(yi;k)=sqrt(1+k^2snh^2(y;k))
 =>sn(x+yi;k)=(snxcn(yi)dn(yi)+sn(yi)cnxdn(x))/(1-K^2sn^2(x)sn^(yi))
 =(snxcnhysqrt(1+k^2snh^2(y))+isnhycnxdnx)/(1+k^2snh^2(y))
 =>k=0,sin(x+yi)=sinxchy+icosxshy
 
 

 调和函数（二元实函）等式：
 Resin(x+yi)=math.sin(x)*math.cosh(y)
 Imsin(x+yi)=math.cos(x)*math.sinh(y)
 Recos(x+yi)=math.cos(x)*math.cosh(y)
 Imcos(x+yi)=-math.sin(x)*math.sinh(y)
 Resn(x+yi;k)=
 Imsn(x+yi;k)=
 Recn(x+yi;k)=
 Imcn(x+yi;k)=
 Redn(x+yi;k)=
 Imdn(x+yi;k)=
 双曲函数的一种推广：
 sin(yi)=ish(y),cos(yi)=ch(y)=>
 sn(yi;k)=isnh(y;k),cn(yi;k)=cnh(y;k)
 易知snh(y;k)、cnh(y;k)的幂级数为：……，k=0时即分别为sh(y),ch(y)。
 可证明：
 snh(y;k)=sn(y;k')/cn(y;k')
 =>k=0,sh(y)=sn(y;1)/cn(y;1)=[1-2/(e^(2y)+1)]/[2e^y/(e^(2y)+1)]=(e^y-e^(-y))/2=sh(y)
 cnh(y;k)=1/cn(y;k')
 =>k=0,ch(y)=1/cn(y;1)=1/[2e^y/(e^(2y)+1)]=(e^y+e^(-y))/2=ch(y)
 dn(yi;k)=dn(y;k')/cn(y;k')
 =>dn(yi;k=0)=1=dn(y;1)/cn(y;1)=cn(y;1)/cn(y;1)=1
 i/18+pi/9=pi/6
 sn(u=0.711647;e=0.5)=0.6427874
 Resn(u=0.711647;e=0.5)=0.6435
 Imsn(u=0.711647;e=0.5)=0
 sin(pi/18)=0.17364819347858429
 sin(pi/9)=0.34202015399932861
 sin(pi/6)=0.50000000000000000
 cos(pi/18)=0.98480772972106934
 cos(pi/9)=0.93969261646270752
 cos(pi/6)=0.86602538824081421
 Resn(u=pi/18;e=0)=0.173648178645355
 Recn(u=pi/18;e=0)=0.98480775299086
 Redn(u=pi/18;e=0)=1
 Resn(u=pi/9;e=0)=0.342020268405197
 Recn(u=pi/9;e=0)=0.939692615326438
 Redn(u=pi/9;e=0)=1
 Resn(u=pi/6;e=0)=0.500002132588792
 2011.6.9加法公式中的分母应为-，而非+：
 sn(u+v)=(snucnvdnv+snvcnudnu)/(1-(ksnusnv)^2)
 Resn(u=pi/18;e=0.5)=0.173431462562212<sin(pi/18)
 Recn(u=pi/18;e=0.5)=0.984845945129538>cos(pi/18)
 Redn(u=pi/18;e=0.5)=0.996233096336519
 Resn(u=pi/9;e=0.5)=0.340401930678835<sin(pi/9)
 Recn(u=pi/9;e=0.5)=0.940281076753037>cos(pi/9)
 Redn(u=pi/9;e=0.5)=0.985409405212932
 Resn(u=pi/6;e=0.5)=0.495189316792961<sin(pi/6)
 sn(u+v=0.174533+0.349066=0.523599)=(0.173431*0.940281*0.985409+0.340402*0.984846*0.996233)/(1-(0.500000*0.173431*0.340402)^2)=0.495107<sin(pi/6)
 2012.4.7：sc(u,k)=sn(u,k)/cn(u,k)=u+((2-k^2)/3!)u^3+…
 sc'=(sn/dn)'=(cndncn+snsndn)/cn^2=dn/cn^2
 sc(0)=0/1=0
 sc'(0)=1/1^2=1
 sc''=(-k^2sncn^3+2sncndn^2)/cn^4=sncn(-k^2cn^2+2dn^2)/cn^4
 sc''(0)=0
 sc'''=[cndn+cndn^3-k^2cndn-2k^2sn^2cn^4dn+3sn^2cn^2dn(1+dn^2-k^2)]/cn^6
 sc'''(0)=2-k^2
 (1/x)'=-1/x^2
 关键公式1： (f/g)'=(fg^(-1))'=f(g^(-1))'+f'g^(-1)=f(-g'/g^2)+f'(1/g)=(f'g-fg')/g^2，
 关键公式1.1：(1/g)'=-g'/g^2
 (x/x^2)'=(1·x^2-x·2x)/x^4=-1/x^2
 (sinx/cosx)'=(cosxcosx+sinxsinx)/cos^2x=sec^2x
 利用函数的导数、不定积分求反函数的导数、不定积分
 关键公式2.2：∫f^(-1)(x)dx=xf^(-1)(x)-F(f^(-1)(x))+C
 y=f(x)与x=g(y)=f^-1(y)互为反函数（Inverse function），partial inverse部分逆
 关键公式2.1：g'(y)=1/f'(x)
 (f^(-1))'(y)=1/f'(f^(-1)(y))<--x=f^(-1)(y)-->dx/dy=1/(dy/dx)
 例如：(sinx)'=cosx,dx/dsinx=1/cosx=darcsiny/(cosxdx),y'(arcsiny)'=1
 x=siny,y=arcsinx
 y'=1/x'=1/(siny)'=1/sqrt(1-x^2)
 (arcsinx)'=1/sqrt(1-x^2)
 例如：f(x)=(2x+8)^3,f^(-1)(y)=(y^(1/3)-8)/2
 问：f(x)=x-sinx的反函数是初等函数吗？
 微分方程
 z=arcsinw+C
 z'=dz/dw=1/sqrt(1-w^2)
 =>dw/dz=sqrt(1-w^2),(dw/dz)^2=1-w^2
 z=arcsn(w;k)+C
 z'=dz/dw=1/sqrt((1-w^2)(1-k^2w^2))
 dw/dz=sqrt((1-w^2)(1-k^2w^2))
 (dw/dz)^2=(1-w^2)(1-k^2w^2)
 
 
z=P(u)
 u=P^(-1)(z)----？下限为∞的复积分
 du/dz=1/sqrt(4z^3-g_2z-g_3))----不存在g_2=g_3=0的情形
 (dz/du)^2=4z^3-g_2z-g_3
 P'(u)^2=4P(u)^3-g_2P(u)-g_3
 
 2012.4.15
 sin的加法定理<=>arcsin的加法定理
 sin(x+y)=sinxsqrt(1-sin^2y)+sqrt(1-sin^2x)siny
 arcsina+arcsinb=arcsin(asqrt(1-b^2)+bsqrt(1-a^2))
 cosx=sqrt(1-sin^2(x))是单值函数，
 sl'(x)=sqrt(1-sl^4(x))!=cl(x)是单值函数吗？
 sl的加法定理<=>arcsl的加法定理
 sl(x+a)=(sl(x)sl'(a)+sl(a)sl'(x))/(1+sl^2(x)sl^2(a))
 两边求(d/dx)得：
 sl'(x+a)=(f'g-fg')/g^2=[sl'(a)sl'(x)+sl(a)sl''(x)-2sl^3(a)sl(x)sl'^2(x)+(sl^3(a)-sl'(a)sl^2(a))sl^2(x)sl'(x)]/(1+sl^2(x)sl^2(a))^2
 f=sl(x)sl'(a)+sl(a)sl'(x)
 g=1+sl^2(x)sl^2(a)
 f'=sl'(a)sl'(x)+sl(a)sl''(x)
 g'=2sl^2(a)sl(x)sl'(x)
 ]
 
 
 
（一） 雅可比椭圆函数理论
 sn^2u+cn^2u=1,dn^2u+k^2sn^u=1。 
 1835年，雅可比证明任何单变量单值有理型（即亚纯）函数不可能多于两个周期，且周期比必为非实数。
 雅可比在《椭圆函数论新基础》（1835）一书中创立θ函数理论，从而给椭圆函数一个系统的表示，雅可比的θ函数理论将双周期的椭圆函数（也称亚纯函数）定义为θ函数的商，并在此基础上重建椭圆函数理论。
 （二） 外尔斯特拉斯的椭圆函数论
 外尔斯特拉斯把椭圆函数定义为复平面上具有双周期w_1,w_2(Im(w_2/w_1)>0)的亚纯函数，从而具体构造出著名的P函数P(z;w_1,w_2)=1/z^2+∑’[1/(z-mw_1-nw_2)^2-1/(mw_1+nw_2)^2]，其中∑’表示对(m,n)!=(0,0)的所有整数求和。他还定义了δ函数，δ是整函数。他还定义了ζ函数（不要同黎曼的ζ函数混淆），ζ(z)=δ'(z)/δ(z),ζ'(z)=-P(z)。δ函数是表示椭圆函数的有力工具。
 [
 【注解】
 问：实周期w_1=8.626062，虚周期w_2=7.416195i=>g_2=？，g_3=？
 g_2(实周期w_1=8.626062，虚周期w_2=7.416195i)=(0.049881,0)
 g_3(实周期w_1=8.626062，虚周期w_2=7.416195i)=(-0.00112113,0)
 http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp
 P(6.5;g_2=0.049881,g_3=-0.00112113)=0.231868，此为上述网站的jsp计算结果
 以下为C++计算结果：
 相差1个实周期：P(6.5)=P(15.126062)=0.231869
 相差1个虚周期：P(6.5+7.416195i)=0.231869
 相差半个实周期：P(10.813031)=0.220309
 相差半个虚周期：P(6.5+3.708098i)=-0.031417
 http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp
 P(z=0.1,g_2=1,g_3=0)=100.001
 P(z=0.2,g_2=1,g_3=0)=25.0020
 P(z=0.3,g_2=1,g_3=0)=11.1156
 P(z=0.4,g_2=1,g_3=0)=6.25800
 P(z=0.5,g_2=1,g_3=0)=4.01251
 P(z=0.6,g_2=1,g_3=0)=2.79582
 P(z=0.7,g_2=1,g_3=0)=2.06541
 P(z=0.8,g_2=1,g_3=0)=1.59472
 P(z=0.9,g_2=1,g_3=0)=1.27551
 P(z=1,g_2=1,g_3=0)=1.05084
 
 
y''(x)=-2y^3(x),y(0)=0,y'(0)=1<=>y=sl(x)
 P'^2=4P^3-P,g_2=1,g_3=0，双纽线<=>P=？
 sl^(-1)(x) := int {0..x} 1/sqrt(1-s^4)) = x + x^5/10 + x^9/24 + 5*x^13/208 + .... 
 sl(x) = x - 12*x^5/5! + 3024*x^9/9! - 4390848*x^13/13! + .... 
 A104203   Expansion of the sine lemniscate function sl(x).
  1, 0, 0, 0, -12, 0, 0, 0, 3024, 0, 0, 0, -4390848, 0, 0, 0, 21224560896, 0, 0, 0, -257991277243392, 0, 0, 0, 6628234834692624384, 0, 0, 0, -319729080846260095008768, 0, 0, 0, 26571747463798134334265819136, 0, 0, 0 
 ]
 
 （三） 其他数学家的椭圆函数论
 法国数学家刘维尔在1844年最早把双周期性作为刻画椭圆函数的出发点。他把在有限复平面上亚纯的双周期函数定义为椭圆函数。
 椭圆积分及其反演问题到1832年已有一个相当满意的解答，而一般的阿贝尔积分及其反演问题却遇到极大困难。
 关于阿贝尔函数，黎曼发表了两篇文章：一是《阿贝尔函数论》，一是《论θ函数的零点》，这是前一篇的续篇。前一篇由四部分构成，是他生前发表的最深刻且有丰富内容的著作。在这两篇论文中，黎曼具体的反演了阿贝尔积分。在这两篇论文中，包括了以下内容：
 
 
（1）阿贝尔积分的表示及分类，即对由f(z,w)=0定义的黎曼曲面上所有阿贝尔积分进行分类：
 第一类阿贝尔积分，在黎曼曲面上处处有界，线性独立的第一类阿贝尔积分的数目等于曲面的亏格p，如果曲面的连通数N=2p+1，这p个阿贝尔积分被称为基本积分。
 第二类阿贝尔积分，在黎曼曲面上以有限多点为极点。
 第三类阿贝尔积分，在黎曼曲面上具有对数型奇点。
 每一个阿贝尔积分均为上三类积分的和。
 黎曼还引进相伴曲面观念。设黎曼面由多项式F(s,z)=0定义，F对s是n阶，对z是m阶，则相伴曲面由多项式Q(s,z)=0定义，Q对s是n-2阶，对z是m-2阶，这时第一类阿贝尔积分可表为∫Q(s,z)dz/(F_z)。黎曼面上的有理函数也可以借助相伴曲面来表示。在整个黎曼面上的亚纯函数可表示为线性组合S=C_1W_1+…+C_pW_p+C_(p+1)W_(p+1)+…+C_(p+r)W_(p+r)+C_(p+r+1)，其中C_j是常数，W_j(1<=j<=p)是由第一类函数构成的向量空间的基，E_j是任意数目的第二类初等函数，但是C_j(1<=j<=p+r)之间由2p个线性齐次方程相联系。
 /*#include<cmath>*/ 
 #include<iostream>
 #include<complex>
 #include<limits>
 using namespace std; 
 //typedef double real;
 #define real double
 typedef std::complex<real> fcomplex;
 
 
#define M_PI 3.14159265358979323846
 
 
extern "C" _declspec(dllexport)double __stdcall lelp1(double k,double f);
 extern "C" _declspec(dllexport)double __stdcall lelp2(double k,double f);
 
 
//阿贝尔发现：snz的周期是4K及2iK’，cnz的周期是4K及2K+2iK'，dnz的周期是2K+4iK'。
 //约定K叫做1/4实周期，K'叫做1/4虚周期,tau=K'/K是个纯虚数
 /*
 2012.3.11公式P(u)=P(u+2w)=P(u+2w')=P(-u)在代码中是指：
 P(z)=P(z+w_1)=P(z+w_2)=P(-z)!=P(z+w_1/2)!=P(z+w_2/2)，这里w_1,w_2是周期，不是半周期，
 P可以是math::weierstrass_p2(const complex& z,const complex& w_1, const complex& w_2)
 或math::P_theta(const complex& z, const complex& tau, const complex& w1);，
 这两个函数完全一致，它们的输入参数w_1,w_2均是指周期，而不是指半周期。
 20121225注意：P函数、模形式g_2、g_3等数学函数的c++和jsp实现中，参数w_1,w_2是指最小实周期T，最小虚周期T'，即P(z,w_1,w_2)=P(z,T,T')。
 */
 
 
//算术几何平均
 namespace detail{
  fcomplex agm(const fcomplex& a0,const fcomplex& b0,real prec,unsigned m);
 }
 
 
//算术几何平均
 namespace detail{
  fcomplex agm(const fcomplex& a0,const fcomplex& b0,real prec,unsigned m)
  {
   fcomplex a=a0;
   fcomplex b=b0;
   fcomplex t;
   for(unsigned n=0;std::abs(0.5*(a-b))>prec && n<m;++n)
   {
    t=0.5*(a+b);
    b=std::sqrt(a*b);
    a=t;
   }
   return a;
  }
 }
 
 
//勒让德第一类完全椭圆积分
 namespace detail{
  fcomplex K(const fcomplex& k,real e,unsigned m);
  fcomplex K(const fcomplex& z)
  {
   return K(z,std::abs(z)*1e-15,500);
  }
 }
 
 
//勒让德第一类完全椭圆积分
 namespace detail{
  fcomplex K(const fcomplex& k,real e,unsigned m)
  {
   fcomplex a=detail::agm(1.0,std::sqrt(1.0-k*k),e,m);
   return M_PI/2.0/a;
  }
 }
 
 
//雅可比theta函数
 namespace detail{
  //computed to pretty high accuracy(a few epsilon)
  fcomplex theta_q(const fcomplex& z,const fcomplex& q);//JT3
  fcomplex theta_1(const fcomplex& z,const fcomplex& q);//JT1
 }
 
 
//雅可比theta函数
 namespace detail{
  fcomplex theta_q(const fcomplex& z,const fcomplex& q)
  {
   if(std::abs(q)>=1)//XXX invalid q
    return 0;
   fcomplex s=0.0;
   fcomplex t=0.0;
   unsigned n=1;
   do
   {
    t=std::pow(q,n*n)*std::cos(static_cast<real>(2*n)*z);
    s+=t;
    ++n;
   }while(n<100 && std::abs(t)>5*std::numeric_limits<real>::epsilon());
   //compute to high precision because we implement other function using that
      return 1.0+2.0*s;
  }
 
 
 inline fcomplex theta(const fcomplex& z,const fcomplex& tau)
  {
   return theta_q(z*M_PI,std::exp(M_PI*fcomplex(0,1)*tau));
  }
 
 
 fcomplex theta_1(const fcomplex& z,const fcomplex& q)
  {
   if(std::abs(q)>=1)//XXX invalid q
    return 0;
   fcomplex s=0.0;
   fcomplex t=0.0;
   unsigned n=0;
   int f=1;
   do
   {
    t=static_cast<real>(f)*std::pow(q,n*(n+1))*std::sin(static_cast<real>(2*n+1)*z);
    s+=t;
    ++n;
    f*=-1;
   }while(n<100 && std::abs(t)>5*std::numeric_limits<real>::epsilon());
   //compute to high precision because we implement other function using that
   return 2.0*std::pow(q,1.0/4)*s;
  }
 
 
 inline fcomplex theta_2(const fcomplex& z,const fcomplex& q)
  {
   return theta_1(z+M_PI/2,q);
  }
 
 
 inline fcomplex theta_3(const fcomplex& z,const fcomplex& q)
  {
   return theta_q(z,q);
  }
 
 
 inline fcomplex theta_4(const fcomplex& z,const fcomplex& q)
  {
   return theta_3(z,-q);
  }
 
 
 //JT3
  inline fcomplex theta_00(const fcomplex& z,const fcomplex& tau)
  {
   return theta(z,tau);
  }
 
 
 //JT4
  inline fcomplex theta_01(const fcomplex& z,const fcomplex& tau)
  {
   return theta(z+0.5,tau);
  }
 
 
 //JT2
  inline fcomplex theta_10(const fcomplex& z,const fcomplex& tau)
  {
   return std::exp(M_PI*fcomplex(0,1)*(tau/4.0+z))*theta(z+tau/2.0,tau);
  }
 
 
 //JT1
  inline fcomplex theta_11(const fcomplex& z,const fcomplex& tau)
  {
   return std::exp(M_PI*fcomplex(0,1)*(tau/4.0+z+0.5))*theta(z+tau/2.0+0.5,tau);
  }
 }
 
 
namespace detail{
  fcomplex q_from_k(const fcomplex& k);
   fcomplex setup_period(const fcomplex& k,const fcomplex& g);
  //add by Ivan_han
  //计算k(q)
  real kq(real q)
  {
   real s1=0,s2=1;
   for(int i=1;i<=1000;i++)
    s1=s1+2*pow(q,(real)(2*i-1)*(2*i-1)/4.f);
   for(int i=1;i<=1000;i++)
    s2=s2+2*pow(q,(i*i));
   return pow(s1/s2,2);
  }
 
 
 int iseven(int n)
  {
     int ret=1;
     if(n%2==1)
      ret=-1;
     return ret;
  }
 
 
 //计算k'(q)
  real ckq(real q)
  {
   real s1=1,s2=1;
   for(int i=1;i<=1000;i++)
    s1=s1+2*pow(q,(i*i))*iseven(i);
   for(int i=1;i<=1000;i++)
    s2=s2+2*pow(q,(i*i));
   return pow(s1/s2,2);
  }
 }
 
 
namespace detail{
  fcomplex q_from_k(const fcomplex& k)
  {
   fcomplex kprimes=std::pow(1.0-k*k,1.0/4);
   fcomplex l=0.5*(1.0-kprimes)/(1.0+kprimes);
         fcomplex l4=l*l*l*l;
   fcomplex q=l;
   l*=l4;q+=2.0*l;
   l*=l4;q+=15.0*l;
   l*=l4;q+=150.0*l;
   l*=l4;q+=1707.0*l;
   l*=l4;q+=20910.0*l;
   l*=l4;q+=268616.0*l;
   //XXX manually truncated series representation,error should be less than 10^(-4) for |k|<500
   //std::cout<<268616*std::pow(l,25)<<std::endl;
   return q;
  }
 
 
 fcomplex setup_period(const fcomplex& k,const fcomplex& g)
  {
   fcomplex K;//real quarter period
   fcomplex iK_prime;
   fcomplex P;
   K=detail::K(k,1e-25,5000);
   iK_prime=fcomplex(0,1)*detail::K(std::sqrt(1.0-k*k),1e-25,5000);
   if((K/g).imag()>(iK_prime/g).imag())
    P=K;
   else
    P=iK_prime;
   P*=4.0/g;
   return P;
  }
 }
 
 
namespace detail{
  fcomplex sn(const fcomplex& u,const fcomplex& k);
  fcomplex P(const fcomplex& z, const fcomplex& tau);//实周期为1个单位（w1=1）。计算结果同math::detail::weierstrass_e1
  fcomplex P(const fcomplex& z, const fcomplex& tau, const fcomplex& w1);//实周期为w1个单位（w1=w1）。计算结果同math::weierstrass_p2
 
 
 fcomplex e1(const fcomplex& tau);//计算结果同math::detail::weierstrass_e1
  fcomplex e2(const fcomplex& tau);//计算结果同math::detail::weierstrass_e2
  fcomplex e3(const fcomplex& tau);//计算结果约同（≈）math::detail::weierstrass_e3
 
 
 //平稳值e_1,e_2,e_3的计算
  inline fcomplex e1(const fcomplex& tau, const fcomplex& w1)
  {
    return e1(tau)/(w1*w1);
  }
  inline fcomplex e2(const fcomplex& tau, const fcomplex& w1)
  {
    return e2(tau)/(w1*w1);
  }
  inline fcomplex e3(const fcomplex& tau, const fcomplex& w1)
  {
    return e3(tau)/(w1*w1);
  }
 
 
  //不变量g_2,g_3的计算
   inline fcomplex g2(const fcomplex& tau, const fcomplex& w1)
   {
    fcomplex e1=detail::e1(tau,w1);
    fcomplex e2=detail::e2(tau,w1);
    fcomplex e3=detail::e3(tau,w1);
    fcomplex retg2=fcomplex(2.0,0)*(e1*e1+e2*e2+e3*e3);
    return retg2;
   }
   inline fcomplex g3(const fcomplex& tau, const fcomplex& w1)
   {
    fcomplex e1=detail::e1(tau,w1);
    fcomplex e2=detail::e2(tau,w1);
    fcomplex e3=detail::e3(tau,w1);
    fcomplex retg3=fcomplex(4.0,0)*e1*e2*e3;
    return retg3;
   }
 
 
  void Test1()
   {
       std::cout<<"g_2(实周期w_1=8.626062，虚周期w_2=7.416195i)="<<g2(fcomplex(0,7.416195/8.626062),fcomplex(8.626062,0))<<std::endl;
    std::cout<<"g_3(实周期w_1=8.626062，虚周期w_2=7.416195i)="<<g3(fcomplex(0,7.416195/8.626062),fcomplex(8.626062,0))<<std::endl;
   }
 
 
  void Test2()
   {
       std::cout<<"sn(u=1.26604138177897;k=0.985171431009416)="<<sn(fcomplex(1.26604138177897,0),fcomplex(0.985171431009416,0))<<std::endl;
   }
 }
 
 
namespace detail{
  fcomplex sn(const fcomplex& u,const fcomplex& k)
  {
   if(k==0.0)
     return std::sin(u);
   static fcomplex last_k=std::numeric_limits<real>::quiet_NaN();
   static fcomplex f;
   static fcomplex g;
   static fcomplex q;
   static fcomplex P;//period to use for reduction
   if(k!=last_k)
   {
    q=q_from_k(k);
    f=theta_3(0,q)/theta_2(0,q);
    g=theta_3(0,q);
    g*=g;
    P=setup_period(k,g);
    last_k=k;
   }
   fcomplex z=u/g;
 
 
  //reduce imaginary part
   int a=static_cast<int>(z.imag()/P.imag());
   z-=static_cast<real>(a)*P;
 
 
  return f*theta_1(z,q)/theta_4(z,q);
  }
 
 
 fcomplex P(const fcomplex& z, const fcomplex& tau)
  {
     fcomplex z1=fcomplex(M_PI*M_PI,0);
     fcomplex z2=z1*std::pow(theta_00(fcomplex(0,0),tau),2);
     z2=z2*std::pow(theta_10(fcomplex(0,0),tau),2);
     z2=z2*std::pow(theta_01(z,tau),2);
     z2=z2/std::pow(theta_11(z,tau),2);
     fcomplex z3=(z1/fcomplex(3.0,0))*(std::pow(theta_00(fcomplex(0,0),tau),4)+std::pow(theta_10(fcomplex(0,0),tau),4));
     return z2-z3;
  }
 
 
 fcomplex P(const fcomplex& z, const fcomplex& tau, const fcomplex& w1)
  {
     fcomplex retP=P(z/w1,tau)/(w1*w1);
     return retP;
  }
 
 
 fcomplex e1(const fcomplex& tau)
  {
     fcomplex z1=fcomplex(M_PI*M_PI,0);
     fcomplex e1=(z1/fcomplex(3.0,0))*(std::pow(theta_00(fcomplex(0,0),tau),4)+std::pow(theta_01(fcomplex(0,0),tau),4));
     return e1;
  }
 
 
 fcomplex e2(const fcomplex& tau)
  {
     fcomplex z1=fcomplex(M_PI*M_PI,0);
     fcomplex e2=-(z1/fcomplex(3.0,0))*(std::pow(theta_00(fcomplex(0,0),tau),4)+std::pow(theta_10(fcomplex(0,0),tau),4));
     return e2;
  }
 
 
 fcomplex e3(const fcomplex& tau)
  {
     fcomplex z1=fcomplex(M_PI*M_PI,0);
     fcomplex e3=(z1/fcomplex(3.0,0))*(std::pow(theta_10(fcomplex(0,0),tau),4)-std::pow(theta_01(fcomplex(0,0),tau),4));
     return e3;
  }
 
 
}
 
 

 static double __stdcall fk(double k,double f)
 { 
  int m,i,j;
  double s,p,ep,h,aa,bb,w,xx,g,q;
  static double t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  static double c[5]={0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851};
  m=1; g=0.0;
  h=fabs(f); s=fabs(0.0001*h);
  p=1.0e+35; ep=0.000001;
  while ((ep>=0.0000001)&&(fabs(h)>s))
  { 
   g=0.0;
   for(i=1;i<=m;i++)
   { 
    aa=(i-1.0)*h;
    bb=i*h;
    w=0.0;
    for(j=0;j<=4;j++)
    { 
     xx=((bb-aa)*t[j]+(bb+aa))/2.0;
     q=sqrt(1.0-k*k*sin(xx)*sin(xx));
     w=w+c[j]/q;
    }
    g=g+w;
   }
   g=g*h/2.0;
   ep=fabs(g-p)/(1.0+fabs(g));
   p=g; 
   m=m+m; 
   h=0.5*h;
  }
  return(g);
 }
 
 
 
 
 
static double ek(double k,double f)
 { 
  int m,i,j;
  double s,p,ep,h,aa,bb,w,xx,g,q;
  static double t[5]={-0.9061798459,-0.5384693101,0.0,0.5384693101,0.9061798459};
  static double c[5]={0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851};
  m=1;
  g=0.0;
  h=fabs(f);
  s=fabs(0.0001*h);
  p=1.0e+35;
  ep=0.000001;
  while((ep>=0.0000001)&&(fabs(h)>s))
  { 
   g=0.0;
   for(i=1;i<=m;i++)
   { 
    aa=(i-1.0)*h; bb=i*h;
    w=0.0;
    for(j=0;j<=4;j++)
    { 
     xx=((bb-aa)*t[j]+(bb+aa))/2.0;
     q=sqrt(1.0-k*k*sin(xx)*sin(xx));
     w=w+q*c[j];
    }
    g=g+w;
   }
   g=g*h/2.0;
   ep=fabs(g-p)/(1.0+fabs(g));
   p=g;
   m=m+m;
   h=0.5*h;
  }
  return(g);
 }
 
 
double __stdcall lelp1(double k,double f)
 { 
   int n;
   double pi,y,e,ff;
   if(k<0.0)
    k=-k;
   if (k>1.0)
    k=1.0/k;
   pi=3.1415926;
   y=fabs(f);
   n=0;
   while (y>=pi)
   { 
    n=n+1;
    y=y-pi;
   }
   e=1.0;
   if(y>=pi/2.0)
   { 
    n=n+1; 
    e=-e; 
    y=pi-y;
   }
   if(n==0)
    ff=fk(k,y);
   else
   { 
    ff=fk(k,pi/2.0);
    ff=2.0*n*ff+e*fk(k,y);
   }
   if(f<0.0)
    ff=-ff;
   return(ff);
 }
 
 
double __stdcall lelp2(double k,double f)
 { 
   int n;
   double pi,y,e,ff;
   if(k<0.0) k=-k;
   if(k>1.0) k=1.0/k;
   pi=3.1415926; 
   y=fabs(f);
   n=0;
   while(y>=pi)
   { 
    n=n+1;
    y=y-pi;
   }
   e=1.0;
   if(y>=pi/2.0)
   { 
    n=n+1;
    e=-e;
    y=pi-y;
   }
   if(n==0)
    ff=ek(k,y);
   else
   { 
    ff=ek(k,pi/2.0);
    ff=2.0*n*ff+e*ek(k,y);
   }
   if(f<0.0)ff=-ff;
   return(ff);
 }
 namespace Han 
 { 
  //雅可比分号形式的第一类[不完全]椭圆积分Arcsn(e;x)
  double __stdcall lelp1J(double(__stdcall *plelp1)(double e,double x),double e,double x)
  {
    return plelp1(e,asin(x));
  }
 } 
 
 
namespace math 
 { 
  double K(double e)
  {
    return lelp1(e,asin(1.0f));
  }
 
 
 double cK(double e)
  {
    return lelp1(sqrt(1-e*e),asin(1.0f));
  }
 
 
 double P(double x,double e)
  {
    fcomplex ret=detail::P(fcomplex(x,0),fcomplex(0,cK(e)/K(e)),fcomplex(K(e),0));
 #undef real
    return ret.real();
 #define real double
  }
 
 
 //tau=it
  double g2(double t,double w1)
  {
    fcomplex ret=detail::g2(fcomplex(0,t),fcomplex(w1,0));
 #undef real
    return ret.real();
 #define real double
  }
 
 
 //tau=it
  double g3(double t,double w1)
  {
    fcomplex ret=detail::g3(fcomplex(0,t),fcomplex(w1,0));
 #undef real
    return ret.real();
 #define real double
  }
 }
 
 
int main(void) 
 { 
  double ret=Han::lelp1J(lelp1,0.866025,0.466428);//测试命名空间Han，::是作用域解析运算符 
     printf("arcsn(e=0.866025;x=0.466428)=%f\n",ret);//arcsn(e=0.866025;x=0.466428)=0.499999,//sn(u=0.500000+0.000000i;k=0.866025+0.000000i)=0.466428+0.000000i
  
  
  /*
   根据e计算1/4实、虚周期
   real retw1=math::K(complex(0.8660254,0),prec,m).real();
   real retw2=math::K(complex(std::sqrt((real)0.25),0),prec,m).real();
   complex w_1=complex(retw1,0);
   complex w_2=complex(0,retw2);
   tau=w_2/w_1=1.685750i/2.156516=0.781701i
   P(z=0.5,w_1=2.156516,w_2=1.685750i)=4.197128=P(0.5,k=0.8660254)!=4.013195
   P(z=0.5,w_1=1,w_2=tau)=P(z=0.5,tau=0.781701i)=7.750935
   错误的结果：
   g_2(tau=0.781701i)=413.330047
   g_3(tau=0.781701i)=-2086.126608
   正确的结果：
   g2(tau=0.781701i,w1=2.156516)=(17.3333,0)
   g3(tau=0.781701i,w1=2.156516)=(-10.3704,0)
   g2(tau=0.781701i,w1=1)=(374.88,0)
   g3(tau=0.781701i,w1=1)=(-1043.06,0)
   http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp
   jsp计算结果：
   P(0.5;17.3333,-10.3704)=4.19713
 
 
  问：已知w_1,w_2，求e,g_2,g_3？
   P(z=0.5,w_1=2.156516,w_2=1.854049i)=4.151395
   P(z=0.5,w_1=1,tau=0.859743i)=7.294863
   tau=w_2/w_1=0.859743i
   g2(tau=0.859743i,w1=1)=276.176
   g3(tau=0.859743i,w1=1)=-461.88
   g2(tau=0.859743i,w1=2.156516)=12.7695
   g3(tau=0.859743i,w1=2.156516)=-4.59211
   jsp计算结果：
   P(z=0.5,g_2=12.7695,g_3=-4.59211)=4.151395
   P(z=0.231855,g_2=12.7695,g_3=-4.59211)=18.6362
   c++计算结果：
   P(z=0.5,w_1=1,w_2=tau=0.859743i)=7.294862=(7.29486,-7.87764e-031)
   P(z_=z/w_1=0.231855,w_1=1,tau=0.859743i)=19.306386=(19.3064,9.34084e-015)
   P(z=0.5,w_1=2.156516,w_2=w_1*0.859743i)=(4.15139,1.00427e-015)
   P(z_=z/w_1=0.231855,w_1=2.156516,w_2=w_1*0.859743i)=(18.6362,1.02781e-014)
  */
  double retK=math::K(0.8660254);
  printf("K(e=0.8660254)=%f\n",retK);//K(e=0.8660254)=2.156516
  double retcK=math::cK(0.8660254);
  printf("K'(e=0.8660254)=%f\n",retcK);//K'(e=0.8660254)=1.685750
  cout<<"P(z=0.5,e=0.8660254)="<<math::P(0.5,0.8660254)<<endl;//P(z=0.5,e=0.8660254)=4.19713
  cout<<"P(z=0.5,tau=0.781701i)="<<detail::P(fcomplex(0.5,0),fcomplex(0,0.781701))<<endl;//P(z=0.5,tau=0.781701i)=(7.75093,0)
  cout<<"P(z=0.5,tau=0.781701i,w1=2.156516)="<<detail::P(fcomplex(0.5,0),fcomplex(0,0.781701),fcomplex(2.156516,0))<<endl;//P(z=0.5,tau=0.781701i,w1=2.156516)=(4.19713,-1.0154e-015)
  cout<<"g2(tau=0.781701i,w1=2.156516)="<<detail::g2(fcomplex(0,0.781701),fcomplex(2.156516,0))<<endl;//g2(tau=0.781701i,w1=2.156516)=(17.3333,0)
  cout<<"g3(tau=0.781701i,w1=2.156516)="<<detail::g3(fcomplex(0,0.781701),fcomplex(2.156516,0))<<endl;//g3(tau=0.781701i,w1=2.156516)=(-10.3704,0)
  cout<<"g2(tau=0.781701i,w1=1)="<<detail::g2(fcomplex(0,0.781701),fcomplex(1,0))<<endl;//g2(tau=0.781701i,w1=1)=(374.88,0)
  cout<<"g3(tau=0.781701i,w1=1)="<<detail::g3(fcomplex(0,0.781701),fcomplex(1,0))<<endl;//g3(tau=0.781701i,w1=1)=(-1043.06,0)
 
 
    cout<<"g2(tau=0.859743i,w1=1)="<<math::g2(0.859743,1)<<endl;//g2(tau=0.859743i,w1=1)=
     cout<<"g3(tau=0.859743i,w1=1)="<<math::g3(0.859743,1)<<endl;//g3(tau=0.859743i,w1=1)=
     cout<<"g2(tau=0.859743i,w1=2.156516)="<<math::g2(0.859743,2.156516)<<endl;//g2(tau=0.859743i,w1=2.156516)=
     cout<<"g3(tau=0.859743i,w1=2.156516)="<<math::g3(0.859743,2.156516)<<endl;//g3(tau=0.859743i,w1=2.156516)=
 
 
 cout<<detail::P(fcomplex(0.5,0),fcomplex(0,0.859743))<<endl;
  cout<<detail::P(fcomplex(0.231855,0),fcomplex(0,0.859743))<<endl;
 
 
 cout<<detail::P(fcomplex(0.5,0),fcomplex(0,0.859743),fcomplex(2.156516,0))<<endl;
  cout<<detail::P(fcomplex(0.231855,0),fcomplex(0,0.859743),fcomplex(2.156516,0))<<endl;
 
 
 fcomplex x(1,2.1),y;
    std::cout<<"sin(1,2.1)="<<std::sin(x)<<std::endl;
    //复数求模abs(55.6,68.2)=87.992
    std::cout<<"abs(55.6,68.2)="<<std::abs(fcomplex(55.6,68.2))<<std::endl;
    //求实幂指数(16+81i)^3=(-310832,-469233)
    fcomplex cpxZ(16,81);
    fcomplex cpxPow=std::pow(cpxZ,3);
    std::cout<<"(16+81i)^3="<<cpxPow<<std::endl;
    //复数相乘(9+11i)(56+3i)=(471,643)
    fcomplex cpxZ1(9,11),cpxZ2(56,3);
    std::cout<<"(9+11i)(56+3i)="<<cpxZ1*cpxZ2<<std::endl;
    //复数相除(471+643i)/(9+11i)=(56,3)
    cpxZ1=fcomplex(471,643);
    cpxZ2=fcomplex(9,11);
    std::cout<<"(471+643i)/(9+11i)="<<cpxZ1/cpxZ2<<std::endl;
 
 
   detail::Test2();
 
 
 system("pause");
  return 0;
 
 
}