JacobiEllipticFunctions.h

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 2013 Soumyajit De
 *
 * KRYLSTAT Copyright 2011 by Erlend Aune <erlenda@math.ntnu.no> under GPL2+
 * (few parts rewritten and adjusted for shogun)
 *
 * NOTE: For higher precision, the methods in this class rely on an external
 * library, ARPREC (http://crd-legacy.lbl.gov/~dhbailey/mpdist/), in absense of
 * which they fallback to shogun datatypes. To use it with shogun, configure
 * ARPREC with `CXX="c++ -fPIC" ./configure' in order to link.
 */

#ifndef JACOBI_ELLIPTIC_FUNCTIONS_H_
#define JACOBI_ELLIPTIC_FUNCTIONS_H_

#include <shogun/lib/config.h>
#include <shogun/base/SGObject.h>
#include <shogun/mathematics/Math.h>
#include <limits>

#if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
#include <arprec/mp_real.h>
#include <arprec/mp_complex.h>
#endif //(HAVE_ARPREC && USE_GPL_SHOGUN)

namespace shogun
{

class CJacobiEllipticFunctions: public CSGObject
{
#if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
    typedef mp_real Real;
    typedef mp_complex Complex;
#else
    typedef float64_t Real;
    typedef complex128_t Complex;
#endif //(HAVE_ARPREC && USE_GPL_SHOGUN)
private:
    static inline Real compute_quarter_period(Real b)
    {
#if defined HAVE_ARPREC && defined USE_GPL_SHOGUN
        const Real eps=mp_real::_eps;
        const Real pi=mp_real::_pi;
#else
        const Real eps=std::numeric_limits<Real>::epsilon();
        const Real pi=M_PI;
#endif //(HAVE_ARPREC && USE_GPL_SHOGUN)
        Real a=1.0;
        Real mm=1.0;

        int64_t p=2;
        do
        {
            Real a_new=(a+b)*0.5;
            Real b_new=sqrt(a*b);
            Real c=(a-b)*0.5;
            mm=Real(p)*c*c;
            p<<=1;
            a=a_new;
            b=b_new;
        } while (mm>eps);
        return pi*0.5/a;
    }

    static inline Real poly_six(Real x)
    {
        return (132*pow(x,6)+42*pow(x,5)+14*pow(x,4)+5*pow(x,3)+2*pow(x,2)+x);
    }

public:
    static void ellipKKp(Real L, Real &K, Real &Kp);

    static void ellipJC(Complex u, Real m, Complex &sn, Complex &cn,
        Complex &dn);

#if defined HAVE_ARPREC && defined USE_GPL_SHOGUN

    static void ellipKKp(float64_t L, float64_t &K, float64_t &Kp)
    {
        mp::mp_init(100, NULL, true);
        mp_real _K, _Kp;
        ellipKKp(mp_real(L), _K, _Kp);
        K=dble(_K);
        Kp=dble(_Kp);
        mp::mp_finalize();
    }

    static void ellipJC(complex128_t u, float64_t m,
        complex128_t &sn, complex128_t &cn, complex128_t &dn)
    {
        mp::mp_init(100, NULL, true);
        mp_complex _sn, _cn, _dn;
        ellipJC(mp_complex(u.real(),u.imag()), mp_real(m), _sn, _cn, _dn);
        sn=complex128_t(dble(_sn.real),dble(_sn.imag));
        cn=complex128_t(dble(_cn.real),dble(_cn.imag));
        dn=complex128_t(dble(_dn.real),dble(_dn.imag));
        mp::mp_finalize();
    }
#endif //(HAVE_ARPREC && USE_GPL_SHOGUN)

    virtual const char* get_name() const
    {
        return "JacobiEllipticFunctions";
    }
};

}

#endif /* JACOBI_ELLIPTIC_FUNCTIONS_H_ */