SHOGUN  4.1.0
JacobiEllipticFunctions.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
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
16 using namespace shogun;
17
18 void CJacobiEllipticFunctions::ellipKKp(Real L, Real &K, Real &Kp)
19 {
20  REQUIRE(L>=0.0,
21  "CJacobiEllipticFunctions::ellipKKp(): \
22  Parameter L should be non-negative\n");
23 #ifdef HAVE_ARPREC
24  const Real eps=Real(std::numeric_limits<float64_t>::epsilon());
25  const Real pi=mp_real::_pi;
26 #else
27  const Real eps=std::numeric_limits<Real>::epsilon();
28  const Real pi=M_PI;
29 #endif //HAVE_ARPREC
30  if (L>10.0)
31  {
32  K=pi*0.5;
33  Kp=pi*L+log(4.0);
34  }
35  else
36  {
37  Real m=exp(-2.0*pi*L);
38  Real mp=1.0-m;
39  if (m<eps)
40  {
41  K=compute_quarter_period(sqrt(mp));
43  }
44  else if (mp<eps)
45  {
47  Kp=compute_quarter_period(sqrt(m));
48  }
49  else
50  {
51  K=compute_quarter_period(sqrt(mp));
52  Kp=compute_quarter_period(sqrt(m));
53  }
54  }
55 }
56
58  ::ellipJC(Complex u, Real m, Complex &sn, Complex &cn, Complex &dn)
59 {
60  REQUIRE(m>=0.0 && m<=1.0,
61  "CJacobiEllipticFunctions::ellipJC(): \
62  Parameter m should be >=0 and <=1\n");
63
64 #ifdef HAVE_ARPREC
65  const Real eps=sqrt(mp_real::_eps);
66 #else
67  const Real eps=sqrt(std::numeric_limits<Real>::epsilon());
68 #endif //HAVE_ARPREC
69  if (m>=(1.0-eps))
70  {
71 #ifdef HAVE_ARPREC
72  complex128_t _u(dble(u.real),dble(u.imag));
75  complex128_t twon=b*CMath::sinh(_u);
76  complex128_t ai=0.25*(1.0-dble(m));
77  complex128_t _sn=t+ai*(twon-_u)/(b*b);
78  complex128_t phi=1.0/b;
79  complex128_t _cn=phi-ai*(twon-_u);
80  complex128_t _dn=phi+ai*(twon+_u);
81  sn=mp_complex(_sn.real(),_sn.imag());
82  cn=mp_complex(_cn.real(),_cn.imag());
83  dn=mp_complex(_dn.real(),_dn.imag());
84 #else
85  Complex t=CMath::tanh(u);
86  Complex b=CMath::cosh(u);
87  Complex ai=0.25*(1.0-m);
88  Complex twon=b*CMath::sinh(u);
89  sn=t+ai*(twon-u)/(b*b);
90  Complex phi=Real(1.0)/b;
91  ai*=t*phi;
92  cn=phi-ai*(twon-u);
93  dn=phi+ai*(twon+u);
94 #endif //HAVE_ARPREC
95  }
96  else
97  {
98  const Real prec=4.0*eps;
99  const index_t MAX_ITER=128;
100  index_t i=0;
101  Real kappa[MAX_ITER];
102
103  while (i<MAX_ITER && m>prec)
104  {
105  Real k;
106  if (m>0.001)
107  {
108  Real mp=sqrt(1.0-m);
109  k=(1.0-mp)/(1.0+mp);
110  }
111  else
112  k=poly_six(m/4.0);
113  u/=(1.0+k);
114  m=k*k;
115  kappa[i++]=k;
116  }
117  Complex sin_u=sin(u);
118  Complex cos_u=cos(u);
119  Complex t=Real(0.25*m)*(u-sin_u*cos_u);
120  sn=sin_u-t*cos_u;
121  cn=cos_u+t*sin_u;
122  dn=Real(1.0)+Real(0.5*m)*(cos_u*cos_u);
123
124  i--;
125  while (i>=0)
126  {
127  Real k=kappa[i--];
128  Complex ksn2=k*(sn*sn);
129  Complex d=Real(1.0)+ksn2;
130  sn*=(1.0+k)/d;
131  cn*=dn/d;
132  dn=(Real(1.0)-ksn2)/d;
133  }
134  }
135 }
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)
static const float64_t epsilon
Definition: libbmrm.cpp:25
#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:66

SHOGUN Machine Learning Toolbox - Documentation