SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
JacobiEllipticFunctions.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Soumyajit De
8  *
9  * KRYLSTAT Copyright 2011 by Erlend Aune <erlenda@math.ntnu.no> under GPL2+
10  * (few parts rewritten and adjusted for shogun)
11  */
12 
15 #include <shogun/lib/config.h>
16 
17 using namespace shogun;
18 
19 void CJacobiEllipticFunctions::ellipKKp(Real L, Real &K, Real &Kp)
20 {
21  REQUIRE(L>=0.0,
22  "CJacobiEllipticFunctions::ellipKKp(): \
23  Parameter L should be non-negative\n");
24 #if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
25  const Real eps=Real(std::numeric_limits<float64_t>::epsilon());
26  const Real pi=mp_real::_pi;
27 #else
28  const Real eps=std::numeric_limits<Real>::epsilon();
29  const Real pi=M_PI;
30 #endif //(HAVE_ARPREC && USE_GPL_SHOGUN)
31  if (L>10.0)
32  {
33  K=pi*0.5;
34  Kp=pi*L+log(4.0);
35  }
36  else
37  {
38  Real m=exp(-2.0*pi*L);
39  Real mp=1.0-m;
40  if (m<eps)
41  {
42  K=compute_quarter_period(sqrt(mp));
44  }
45  else if (mp<eps)
46  {
48  Kp=compute_quarter_period(sqrt(m));
49  }
50  else
51  {
52  K=compute_quarter_period(sqrt(mp));
53  Kp=compute_quarter_period(sqrt(m));
54  }
55  }
56 }
57 
59  ::ellipJC(Complex u, Real m, Complex &sn, Complex &cn, Complex &dn)
60 {
61  REQUIRE(m>=0.0 && m<=1.0,
62  "CJacobiEllipticFunctions::ellipJC(): \
63  Parameter m should be >=0 and <=1\n");
64 
65 #if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
66  const Real eps=sqrt(mp_real::_eps);
67 #else
68  const Real eps=sqrt(std::numeric_limits<Real>::epsilon());
69 #endif //(HAVE_ARPREC && USE_GPL_SHOGUN)
70  if (m>=(1.0-eps))
71  {
72 #if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
73  complex128_t _u(dble(u.real),dble(u.imag));
76  complex128_t twon=b*CMath::sinh(_u);
77  complex128_t ai=0.25*(1.0-dble(m));
78  complex128_t _sn=t+ai*(twon-_u)/(b*b);
79  complex128_t phi=1.0/b;
80  complex128_t _cn=phi-ai*(twon-_u);
81  complex128_t _dn=phi+ai*(twon+_u);
82  sn=mp_complex(_sn.real(),_sn.imag());
83  cn=mp_complex(_cn.real(),_cn.imag());
84  dn=mp_complex(_dn.real(),_dn.imag());
85 #else
86  Complex t=CMath::tanh(u);
87  Complex b=CMath::cosh(u);
88  Complex ai=0.25*(1.0-m);
89  Complex twon=b*CMath::sinh(u);
90  sn=t+ai*(twon-u)/(b*b);
91  Complex phi=Real(1.0)/b;
92  ai*=t*phi;
93  cn=phi-ai*(twon-u);
94  dn=phi+ai*(twon+u);
95 #endif //(HAVE_ARPREC && USE_GPL_SHOGUN)
96  }
97  else
98  {
99  const Real prec=4.0*eps;
100  const index_t MAX_ITER=128;
101  index_t i=0;
102  Real kappa[MAX_ITER];
103 
104  while (i<MAX_ITER && m>prec)
105  {
106  Real k;
107  if (m>0.001)
108  {
109  Real mp=sqrt(1.0-m);
110  k=(1.0-mp)/(1.0+mp);
111  }
112  else
113  k=poly_six(m/4.0);
114  u/=(1.0+k);
115  m=k*k;
116  kappa[i++]=k;
117  }
118  Complex sin_u=sin(u);
119  Complex cos_u=cos(u);
120  Complex t=Real(0.25*m)*(u-sin_u*cos_u);
121  sn=sin_u-t*cos_u;
122  cn=cos_u+t*sin_u;
123  dn=Real(1.0)+Real(0.5*m)*(cos_u*cos_u);
124 
125  i--;
126  while (i>=0)
127  {
128  Real k=kappa[i--];
129  Complex ksn2=k*(sn*sn);
130  Complex d=Real(1.0)+ksn2;
131  sn*=(1.0+k)/d;
132  cn*=dn/d;
133  dn=(Real(1.0)-ksn2)/d;
134  }
135  }
136 }
std::complex< float64_t > complex128_t
Definition: common.h:67
int32_t index_t
Definition: common.h:62
#define REQUIRE(x,...)
Definition: SGIO.h:206
static float64_t cosh(float64_t x)
acos(x), x being a complex128_t not implemented
Definition: Math.h:880
static void ellipKKp(Real L, Real &K, Real &Kp)
static void ellipJC(Complex u, Real m, Complex &sn, Complex &cn, Complex &dn)
#define M_PI
workaround for log2 being a define on cygwin
Definition: Math.h:59
static float64_t sinh(float64_t x)
asin(x), x being a complex128_t not implemented
Definition: Math.h:844
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
static float64_t tanh(float64_t x)
atan2(x), x being a complex128_t not implemented
Definition: Math.h:808
Matrix::Scalar max(Matrix m)
Definition: Redux.h:68

SHOGUN Machine Learning Toolbox - Documentation