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 * NOTE: For higher precision, the methods in this class rely on an external 00013 * library, ARPREC (http://crd-legacy.lbl.gov/~dhbailey/mpdist/), in absense of 00014 * which they fallback to shogun datatypes. To use it with shogun, configure 00015 * ARPREC with `CXX="c++ -fPIC" ./configure' in order to link. 00016 */ 00017 00018 #ifndef JACOBI_ELLIPTIC_FUNCTIONS_H_ 00019 #define JACOBI_ELLIPTIC_FUNCTIONS_H_ 00020 00021 #include <shogun/lib/config.h> 00022 #include <shogun/base/SGObject.h> 00023 #include <limits> 00024 #include <math.h> 00025 00026 #ifdef HAVE_ARPREC 00027 #include <arprec/mp_real.h> 00028 #include <arprec/mp_complex.h> 00029 #endif //HAVE_ARPREC 00030 00031 namespace shogun 00032 { 00033 00055 class CJacobiEllipticFunctions: public CSGObject 00056 { 00057 #ifdef HAVE_ARPREC 00058 typedef mp_real Real; 00059 typedef mp_complex Complex; 00060 #else 00061 typedef float64_t Real; 00062 typedef complex128_t Complex; 00063 #endif //HAVE_ARPREC 00064 private: 00065 static inline Real compute_quarter_period(Real b) 00066 { 00067 #ifdef HAVE_ARPREC 00068 const Real eps=mp_real::_eps; 00069 const Real pi=mp_real::_pi; 00070 #else 00071 const Real eps=std::numeric_limits<Real>::epsilon(); 00072 const Real pi=M_PI; 00073 #endif //HAVE_ARPREC 00074 Real a=1.0; 00075 Real mm=1.0; 00076 00077 int64_t p=2; 00078 do 00079 { 00080 Real a_new=(a+b)*0.5; 00081 Real b_new=sqrt(a*b); 00082 Real c=(a-b)*0.5; 00083 mm=Real(p)*c*c; 00084 p<<=1; 00085 a=a_new; 00086 b=b_new; 00087 } while (mm>eps); 00088 return pi*0.5/a; 00089 } 00090 00091 static inline Real poly_six(Real x) 00092 { 00093 return (132*pow(x,6)+42*pow(x,5)+14*pow(x,4)+5*pow(x,3)+2*pow(x,2)+x); 00094 } 00095 00096 public: 00104 static void ellipKKp(Real L, Real &K, Real &Kp); 00105 00114 static void ellipJC(Complex u, Real m, Complex &sn, Complex &cn, 00115 Complex &dn); 00116 00117 #ifdef HAVE_ARPREC 00118 00124 static void ellipKKp(float64_t L, float64_t &K, float64_t &Kp) 00125 { 00126 mp::mp_init(100, NULL, true); 00127 mp_real _K, _Kp; 00128 ellipKKp(mp_real(L), _K, _Kp); 00129 K=dble(_K); 00130 Kp=dble(_Kp); 00131 mp::mp_finalize(); 00132 } 00133 00141 static void ellipJC(complex128_t u, float64_t m, 00142 complex128_t &sn, complex128_t &cn, complex128_t &dn) 00143 { 00144 mp::mp_init(100, NULL, true); 00145 mp_complex _sn, _cn, _dn; 00146 ellipJC(mp_complex(u.real(),u.imag()), mp_real(m), _sn, _cn, _dn); 00147 sn=complex128_t(dble(_sn.real),dble(_sn.imag)); 00148 cn=complex128_t(dble(_cn.real),dble(_cn.imag)); 00149 dn=complex128_t(dble(_dn.real),dble(_dn.imag)); 00150 mp::mp_finalize(); 00151 } 00152 #endif //HAVE_ARPREC 00153 00155 virtual const char* get_name() const 00156 { 00157 return "JacobiEllipticFunctions"; 00158 } 00159 }; 00160 00161 } 00162 00163 #endif /* JACOBI_ELLIPTIC_FUNCTIONS_H_ */