Hyp.cpp

#include "Headers.h"
#include "Hyp.h"
#include "Gamma.h"

#define PRECISION   1E-16
#define ROUND(x)    ((x >= 0)?floor(x+0.5):ceil(x-0.5))

double Hypergeometric2F1Series(double a, double b, double c, double z)
{
    unsigned int i = 1;
    double term = a*b*z/c;
    double sum = 1+term;

    do
    {
        term *= ((a+i)*(b+i)*z)/((c+i)*(i+1));
        sum += term;
        i++;
    }
    while(fabs(term) > PRECISION);
    return sum;
}

double Hypergeometric2F1(double a, double b, double c, double z)
{
//  double delta = 0.051776695296636881100211090526;

    if(fabs(z) < PRECISION) return 1;           // z = 0;
    else if(fabs(z-1) < PRECISION)              // z = 1;
    {
        if(fabs(c-a-b) < PRECISION) return 1E300;
        else return Gamma(c)*Gamma(c-a-b)/(Gamma(c-a)*Gamma(c-b));
    }

    if(z < -1)
    {
        double w = 1/(1-z);
        return pow(w,a) * Hypergeometric2F1Series(a,c-b,a-b+1,w) * exp(LnGamma(c) + LnGamma(b-a) - LnGamma(b) - LnGamma(c-a)) +
               pow(w,b) * Hypergeometric2F1Series(b,c-a,b-a+1,w) * exp(LnGamma(c) + LnGamma(a-b) - LnGamma(a) - LnGamma(c-b));
    }
    else if(z < 0)
    {
        double w = z/(z-1);
        return pow(1-w,a) * Hypergeometric2F1Series(a,c-b,c,w);
    }
    else if(z <= 0.5)
    {
        return Hypergeometric2F1Series(a,b,c,z);
    }
    else if(z <= 1)
    {
        double w = 1-z;
        return Hypergeometric2F1Series(a,b,a+b-c+1,w) * exp(LnGamma(c) + LnGamma(c-a-b) - LnGamma(c-a) - LnGamma(c-b)) +
               pow(w,c-a-b) * Hypergeometric2F1Series(c-a,c-b,c-a-b+1,w) * exp(LnGamma(c) + LnGamma(a+b-c) - LnGamma(a) - LnGamma(b));
    }
    else if(z <= 2)
    {
        double w = 1 - (1/z);
        double re = pow(1-w,a) * Hypergeometric2F1Series(a,a-c+1,a+b-c+1,w) * exp(LnGamma(c) + LnGamma(c-a-b) - LnGamma(c-a) - LnGamma(c-b)) +
                    pow(1-w,a) * pow(w,c-a-b) * cos(M_PI*(c-a-b)) * Hypergeometric2F1Series(1-b,c-b,c-a-b+1,w) * exp(LnGamma(c) + LnGamma(a+b-c) - LnGamma(a) - LnGamma(b));
        //double im = -pow(1-w,a) * pow(w,c-a-b) * sin(M_PI*(c-a-b)) * Hypergeometric2F1Series(c-a,1-a,c-a-b+1,w) * exp(LnGamma(c) + LnGamma(a+b-c) - LnGamma(a) - LnGamma(b));
        return re;
    }
    else
    {
        double w = 1/z;
        double re = pow(w,a) * cos(M_PI*a) * Hypergeometric2F1Series(a,a-c+1,a-b+1,w) * exp(LnGamma(c) + LnGamma(b-a) - LnGamma(b) - LnGamma(c-a)) +
                    pow(w,b) * cos(M_PI*b) * Hypergeometric2F1Series(b-c+1,b,b-a+1,w) * exp(LnGamma(c) + LnGamma(a-b) - LnGamma(a) - LnGamma(c-b));
        //double im = pow(w,a) * sin(M_PI*a) * Hypergeometric2F1Series(a,a-c+1,a-b+1,w) * exp(LnGamma(c) + LnGamma(b-a) - LnGamma(b) - LnGamma(c-a)) +
        //            pow(w,b) * sin(M_PI*b) * Hypergeometric2F1Series(b-c+1,b,b-a+1,w) * exp(LnGamma(c) + LnGamma(a-b) - LnGamma(a) - LnGamma(c-b));
        return re;
    }
}

Hyp2F1::Hyp2F1(double a, double b, double c)
{
    this->a = a;
    this->b = b;
    this->c = c;

    this->c1 = fabs(c-a-b);
    this->c2 = Gamma(c)*Gamma(c-a-b)/(Gamma(c-a)*Gamma(c-b));
    this->c3 = exp(LnGamma(c) + LnGamma(b-a) - LnGamma(b) - LnGamma(c-a));
    this->c4 = exp(LnGamma(c) + LnGamma(a-b) - LnGamma(a) - LnGamma(c-b));
    this->c5 = exp(LnGamma(c) + LnGamma(c-a-b) - LnGamma(c-a) - LnGamma(c-b));
    this->c6 = exp(LnGamma(c) + LnGamma(a+b-c) - LnGamma(a) - LnGamma(b));
    this->c7 = cos(M_PI*(c-a-b));
    this->c8 = sin(M_PI*(c-a-b));
    this->c9 = cos(M_PI*a);
    this->c10 = cos(M_PI*b);
    this->c11 = sin(M_PI*a);
    this->c12 = sin(M_PI*b);
}

double Hyp2F1::operator()(double z)
{
    if(fabs(z) < PRECISION) return 1;           // z = 0;
    else if(fabs(z-1) < PRECISION)              // z = 1;
    {
        if(this->c1 < PRECISION) return 1E300;
        else return this->c2;
    }

    if(z < -1)
    {
        double w = 1/(1-z);
        return pow(w,a) * Hypergeometric2F1Series(a,c-b,a-b+1,w) * this->c3 +
               pow(w,b) * Hypergeometric2F1Series(b,c-a,b-a+1,w) * this->c4;
    }
    else if(z < 0)
    {
        double w = z/(z-1);
        return pow(1-w,a) * Hypergeometric2F1Series(a,c-b,c,w);
    }
    else if(z <= 0.5)
    {
        return Hypergeometric2F1Series(a,b,c,z);
    }
    else if(z <= 1)
    {
        double w = 1-z;
        return Hypergeometric2F1Series(a,b,a+b-c+1,w) * this->c5 +
               pow(w,c-a-b) * Hypergeometric2F1Series(c-a,c-b,c-a-b+1,w) * this->c6;
    }
    else if(z <= 2)
    {
        double w = 1 - (1/z);
        double re = pow(1-w,a) * Hypergeometric2F1Series(a,a-c+1,a+b-c+1,w) * this->c5 +
                    pow(1-w,a) * pow(w,c-a-b) * this->c7 * Hypergeometric2F1Series(1-b,c-b,c-a-b+1,w) * this->c6;
        //double im = -pow(1-w,a) * pow(w,c-a-b) * this->c8 * Hypergeometric2F1Series(c-a,1-a,c-a-b+1,w) * this->c6;
        return re;
    }
    else
    {
        double w = 1/z;

        double h1 = pow(w,a) * Hypergeometric2F1Series(a,a-c+1,a-b+1,w) * this->c3;
        double h2 = pow(w,b) * Hypergeometric2F1Series(b-c+1,b,b-a+1,w) * this->c4;

        double re = this->c9 * h1 + this->c10 * h2;
        //double im = this->c11 * h1 + this->c12 * h2;
        
        return re;
    }
}