SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2013 Soumyajit De 00008 * 00009 * KRYLSTAT Copyright 2011 by Erlend Aune <erlenda@math.ntnu.no> under GPL2+ 00010 * (few parts rewritten and adjusted for shogun) 00011 */ 00012 00013 #include <shogun/mathematics/Math.h> 00014 #include <shogun/mathematics/JacobiEllipticFunctions.h> 00015 00016 using namespace shogun; 00017 00018 void CJacobiEllipticFunctions::ellipKKp(Real L, Real &K, Real &Kp) 00019 { 00020 REQUIRE(L>=0.0, 00021 "CJacobiEllipticFunctions::ellipKKp(): \ 00022 Parameter L should be non-negative\n"); 00023 #ifdef HAVE_ARPREC 00024 const Real eps=Real(std::numeric_limits<float64_t>::epsilon()); 00025 const Real pi=mp_real::_pi; 00026 #else 00027 const Real eps=std::numeric_limits<Real>::epsilon(); 00028 const Real pi=M_PI; 00029 #endif //HAVE_ARPREC 00030 if (L>10.0) 00031 { 00032 K=pi*0.5; 00033 Kp=pi*L+log(4.0); 00034 } 00035 else 00036 { 00037 Real m=exp(-2.0*pi*L); 00038 Real mp=1.0-m; 00039 if (m<eps) 00040 { 00041 K=compute_quarter_period(sqrt(mp)); 00042 Kp=Real(std::numeric_limits<float64_t>::max()); 00043 } 00044 else if (mp<eps) 00045 { 00046 K=Real(std::numeric_limits<float64_t>::max()); 00047 Kp=compute_quarter_period(sqrt(m)); 00048 } 00049 else 00050 { 00051 K=compute_quarter_period(sqrt(mp)); 00052 Kp=compute_quarter_period(sqrt(m)); 00053 } 00054 } 00055 } 00056 00057 void CJacobiEllipticFunctions 00058 ::ellipJC(Complex u, Real m, Complex &sn, Complex &cn, Complex &dn) 00059 { 00060 REQUIRE(m>=0.0 && m<=1.0, 00061 "CJacobiEllipticFunctions::ellipJC(): \ 00062 Parameter m should be >=0 and <=1\n"); 00063 00064 #ifdef HAVE_ARPREC 00065 const Real eps=sqrt(mp_real::_eps); 00066 #else 00067 const Real eps=sqrt(std::numeric_limits<Real>::epsilon()); 00068 #endif //HAVE_ARPREC 00069 if (m>=(1.0-eps)) 00070 { 00071 #ifdef HAVE_ARPREC 00072 complex128_t _u(dble(u.real),dble(u.imag)); 00073 complex128_t t=CMath::tanh(_u); 00074 complex128_t b=CMath::cosh(_u); 00075 complex128_t twon=b*CMath::sinh(_u); 00076 complex128_t ai=0.25*(1.0-dble(m)); 00077 complex128_t _sn=t+ai*(twon-_u)/(b*b); 00078 complex128_t phi=1.0/b; 00079 complex128_t _cn=phi-ai*(twon-_u); 00080 complex128_t _dn=phi+ai*(twon+_u); 00081 sn=mp_complex(_sn.real(),_sn.imag()); 00082 cn=mp_complex(_cn.real(),_cn.imag()); 00083 dn=mp_complex(_dn.real(),_dn.imag()); 00084 #else 00085 Complex t=CMath::tanh(u); 00086 Complex b=CMath::cosh(u); 00087 Complex ai=0.25*(1.0-m); 00088 Complex twon=b*CMath::sinh(u); 00089 sn=t+ai*(twon-u)/(b*b); 00090 Complex phi=Real(1.0)/b; 00091 ai*=t*phi; 00092 cn=phi-ai*(twon-u); 00093 dn=phi+ai*(twon+u); 00094 #endif //HAVE_ARPREC 00095 } 00096 else 00097 { 00098 const Real prec=4.0*eps; 00099 const index_t MAX_ITER=128; 00100 index_t i=0; 00101 Real kappa[MAX_ITER]; 00102 00103 while (i<MAX_ITER && m>prec) 00104 { 00105 Real k; 00106 if (m>0.001) 00107 { 00108 Real mp=sqrt(1.0-m); 00109 k=(1.0-mp)/(1.0+mp); 00110 } 00111 else 00112 k=poly_six(m/4.0); 00113 u/=(1.0+k); 00114 m=k*k; 00115 kappa[i++]=k; 00116 } 00117 Complex sin_u=sin(u); 00118 Complex cos_u=cos(u); 00119 Complex t=Real(0.25*m)*(u-sin_u*cos_u); 00120 sn=sin_u-t*cos_u; 00121 cn=cos_u+t*sin_u; 00122 dn=Real(1.0)+Real(0.5*m)*(cos_u*cos_u); 00123 00124 i--; 00125 while (i>=0) 00126 { 00127 Real k=kappa[i--]; 00128 Complex ksn2=k*(sn*sn); 00129 Complex d=Real(1.0)+ksn2; 00130 sn*=(1.0+k)/d; 00131 cn*=dn/d; 00132 dn=(Real(1.0)-ksn2)/d; 00133 } 00134 } 00135 }