SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
JacobiEllipticFunctions.cpp
Go to the documentation of this file.
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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation