q-special-functions/Heine.hpp at master · Daisuke-Kanaizumi/q-special-functions · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
// Author: Daisuke Kanaizumi
// Affiliation: Department of Applied Mathematics, Waseda University
// Date: March 3rd, 2017

// This program was made in order to calculate the Heine hypergeometric function
// 2\phi1(a,b,c,q,z), this is a q extension of Gaussian hypergeometric function

// reference
// Fredrik Johansson, Computing hypergeometric functions rigorously, arXiv, 2016

#ifndef HEINE_HPP
#define HEINE_HPP

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/geoseries.hpp>
#include <algorithm>
#include <limits>
#include <cmath>
#include <kv/convert.hpp> // this was included to use complex numbers
#include <kv/complex.hpp>
namespace kv {

  //0<q<1, |z|<1
  template <class T> interval<T> Heine(const interval<T>& a, const interval<T>& b, const interval<T>& c,const interval<T>& q, const interval<T> & z) {
    int N;
    N=1000;
    interval<T> mid,rad,res,ratio,first;
    mid=1.;
    while(abs(c)>pow(1/q,N)){
      N=N+500;
      // throw std::domain_error("value of N not large enough");
    }
   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
    if (abs(z)>=1){
      throw std::domain_error("absolute value of z must be under 1");
    }
    for(int n=1;n<=N-1;n++){
      mid=mid+qPochhammer(interval<T>(a),interval<T>(q),int (n))*qPochhammer(interval<T>(b),interval<T>(q),int (n))*pow(z,n)
	/qPochhammer(interval<T>(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
    }
    first=abs(qPochhammer(interval<T>(a),interval<T>(q),int (N))*qPochhammer(interval<T>(b),interval<T>(q),int (N))*pow(z,N)
	      /qPochhammer(interval<T>(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*(1+abs(pow(q,N)*(c-a)/(1-c*pow(q,N))))*(1+abs(pow(q,N)*(q-b)/(1-pow(q,N+1))));
    if(ratio<1){
      rad=(first/(1-ratio)).upper();
      res=mid+rad*interval<T>(-1.,1.);
      return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }
 }
  template <class T> interval<T> q_log(const interval<T>& (q),const interval<T> &(z)){
    // verification program for q-logarithm l_q(1-z)
    // reference
    // Adrienne W. Kemp, C. David Kemp, The q-cluster distribution
    // Journal of Statistical Plannning and Inference, 2009

    // 0<q<1, 0<z<1
   if (q>=1){
     throw std::domain_error("value of q must be under 1");
   }
   if (q<=0){
     throw std::domain_error("q must be positive");
   }
   if (z>=1){
     throw std::domain_error("value of z must be under 1");
   }
   if (z<=0){
     throw std::domain_error("z must be positive");
   }
    interval<T>res;
    res=-z*Heine(interval<T>(q),interval<T>(q),interval<T>(q*q),interval<T>(q),interval<T>(z));
    return res;
  }
template <class T> complex<interval<T> >little_q_Jacobi(const complex<interval<T> >& lambda,const complex<interval<T> >& x,const complex<interval<T> >& a, const complex<interval<T> >& b,const interval<T>& q) {
    // verification program for little q-Jacobi function
    // reference
    // Ahmed Fitouhi, Fethi Bouzeffour, Wafa Binous
    // Expansion and asymptotic in terms of basic Bessel functions
    // Applied Mathematics and Computation 188 (2007) 2034–2044
    complex<interval<T> >res;
    res=Heine(complex<interval<T> >(a*lambda),complex<interval<T> >(a/lambda),complex<interval<T> >(a*b),interval<T>(q),complex<interval<T> >(-b*x));
    return res;
 }
 template <class T> complex<interval<T> >Heine(const complex<interval<T> >& a, const complex<interval<T> >& b, const complex<interval<T> >& c,const interval<T>& q, const complex<interval<T> >& z) {
    int N;
    N=1000;
    complex<interval<T> > mid,res;
    interval<T>ratio,first;
    T rad;
    mid=1.;

    while(abs(c)>pow(1/q,N)){
      N=N+500;
      // throw std::domain_error("value of N not large enough");
    }
    if (q>=1){
      throw std::domain_error("value of q must be under 1");
    }
    if (q<=0){
      throw std::domain_error("q must be positive");
    }
    if (abs(z)>=1){
      throw std::domain_error("absolute value of z must be under 1");
    }
    for(int n=1;n<=N-1;n++){
      mid=mid+qPochhammer(complex<interval<T> >(a),interval<T>(q) ,int (n))*qPochhammer(complex<interval<T> >(b),interval<T>(q),int (n))*pow(z,n)
	/qPochhammer(complex<interval<T> >(c),interval<T>(q),int (n))/qPochhammer(interval<T>(q),interval<T>(q),int (n));
    }
    first=abs(qPochhammer(complex<interval<T> >(a),interval<T>(q),int (N))*qPochhammer(complex<interval<T> >(b),interval<T>(q),int (N))*pow(z,N)
	      /qPochhammer(complex<interval<T> >(c),interval<T>(q),int (N))/qPochhammer(interval<T>(q),interval<T>(q),int (N)));
    ratio=abs(z)*(1+abs(pow(q,N)*(c-a)/(1-c*pow(q,N))))*(1+abs(pow(q,N)*(q-b)/(1-pow(q,N+1))));
    if(abs(ratio)<1){
      rad=(first/(1-ratio)).upper();
      res=complex_nbd(mid,rad);
      return res;
    }
    else{
      throw std::domain_error("ratio is more than 1");
    }
}
  template <class T> complex<interval<T> >complex_nbd(const complex<interval<T> >& a,const T & rad){
    interval<T> rc,ic,rcc,icc;
    complex<interval<T> >res;
    rc=a.real();
    ic=a.imag();
    rcc=rc+rad*interval<T>(-1.,1.);
    icc=ic+rad*interval<T>(-1.,1.);
    res=complex<interval<T> >(rcc,icc);
    return res;
  }
}
#endif