SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
JacobiEllipticFunctions.h
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  * 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_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation