SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JacobiEllipticFunctions.h
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  * NOTE: For higher precision, the methods in this class rely on an external
13  * library, ARPREC (http://crd-legacy.lbl.gov/~dhbailey/mpdist/), in absense of
14  * which they fallback to shogun datatypes. To use it with shogun, configure
15  * ARPREC with `CXX="c++ -fPIC" ./configure' in order to link.
16  */
17 
18 #ifndef JACOBI_ELLIPTIC_FUNCTIONS_H_
19 #define JACOBI_ELLIPTIC_FUNCTIONS_H_
20 
21 #include <shogun/lib/config.h>
22 #include <shogun/base/SGObject.h>
23 #include <limits>
24 #include <math.h>
25 
26 #ifdef HAVE_ARPREC
27 #include <arprec/mp_real.h>
28 #include <arprec/mp_complex.h>
29 #endif //HAVE_ARPREC
30 
31 namespace shogun
32 {
33 
56 {
57 #ifdef HAVE_ARPREC
58  typedef mp_real Real;
59  typedef mp_complex Complex;
60 #else
61  typedef float64_t Real;
62  typedef complex128_t Complex;
63 #endif //HAVE_ARPREC
64 private:
65  static inline Real compute_quarter_period(Real b)
66  {
67 #ifdef HAVE_ARPREC
68  const Real eps=mp_real::_eps;
69  const Real pi=mp_real::_pi;
70 #else
71  const Real eps=std::numeric_limits<Real>::epsilon();
72  const Real pi=M_PI;
73 #endif //HAVE_ARPREC
74  Real a=1.0;
75  Real mm=1.0;
76 
77  int64_t p=2;
78  do
79  {
80  Real a_new=(a+b)*0.5;
81  Real b_new=sqrt(a*b);
82  Real c=(a-b)*0.5;
83  mm=Real(p)*c*c;
84  p<<=1;
85  a=a_new;
86  b=b_new;
87  } while (mm>eps);
88  return pi*0.5/a;
89  }
90 
91  static inline Real poly_six(Real x)
92  {
93  return (132*pow(x,6)+42*pow(x,5)+14*pow(x,4)+5*pow(x,3)+2*pow(x,2)+x);
94  }
95 
96 public:
104  static void ellipKKp(Real L, Real &K, Real &Kp);
105 
114  static void ellipJC(Complex u, Real m, Complex &sn, Complex &cn,
115  Complex &dn);
116 
117 #ifdef HAVE_ARPREC
118 
124  static void ellipKKp(float64_t L, float64_t &K, float64_t &Kp)
125  {
126  mp::mp_init(100, NULL, true);
127  mp_real _K, _Kp;
128  ellipKKp(mp_real(L), _K, _Kp);
129  K=dble(_K);
130  Kp=dble(_Kp);
131  mp::mp_finalize();
132  }
133 
141  static void ellipJC(complex128_t u, float64_t m,
142  complex128_t &sn, complex128_t &cn, complex128_t &dn)
143  {
144  mp::mp_init(100, NULL, true);
145  mp_complex _sn, _cn, _dn;
146  ellipJC(mp_complex(u.real(),u.imag()), mp_real(m), _sn, _cn, _dn);
147  sn=complex128_t(dble(_sn.real),dble(_sn.imag));
148  cn=complex128_t(dble(_cn.real),dble(_cn.imag));
149  dn=complex128_t(dble(_dn.real),dble(_dn.imag));
150  mp::mp_finalize();
151  }
152 #endif //HAVE_ARPREC
153 
155  virtual const char* get_name() const
156  {
157  return "JacobiEllipticFunctions";
158  }
159 };
160 
161 }
162 
163 #endif /* JACOBI_ELLIPTIC_FUNCTIONS_H_ */

SHOGUN Machine Learning Toolbox - Documentation