q-special-functions/Sokal.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
// Author: Daisuke Kanaizumi
// Affiliation: Department of Applied Mathematics, Waseda University

// Verified computation of infinite q-Pochhammer symbol by the Sokal's algorithm


// Reference
// Sokal, Alan. D. (2002). arXiv preprint math/0212035.

#ifndef SOKAL_HPP
#define SOKAL_HPP

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/complex.hpp>
#include <kv/constants.hpp>
#include <kv/Heine.hpp>
#include <kv/Pochhammer.hpp>
#include <cmath>

namespace kv{
  template <class T> complex<interval<T> >Sokal_qp(const complex<interval<T> >& z,const complex<interval<T> >& q, int N=50){
    if (abs(q)>=1){
      throw std::domain_error("value of q must be under 1");
    }
    complex<interval<T> >res,sum,zz;
    sum=0.;
    zz=1.;
    interval<T> gamma,pi;
    T rad;
    pi=kv::constants<interval<T> >::pi();
    gamma=1.;
    while(abs(q)<=exp(-gamma)){
      gamma=gamma+1;
    }
    while(abs(z)>=exp((N+1)*gamma)){
      N=N+10;
    }
    while(abs(z*pow(q,N+1))>=1){
      N=N+10;
    }

    for(int n=0;n<=N-1;n++){
      sum=sum+zz*pow(q,n*(n-1)*0.5)/qPochhammer(q,q,n);

      zz*=-z;
    }
    rad=(pow(abs(z),N)*exp(pi*pi/6./gamma-N*(N-1)*gamma*0.5)/(1-abs(z)*exp(-(N+1))*gamma)).upper();
    res=complex_nbd(sum,rad);
    return res;
  }

  template <class T> complex<interval<T> >Sokal_qp(const complex<interval<T> >& z,const interval<T>& q, int N=100){
    if (q>=1){
      throw std::domain_error("value of q must be under 1");
    }
    if (q<=0){
      throw std::domain_error("value of q must be positive");
    }
    complex<interval<T> >res,sum,zz;
    sum=0.;
    zz=1.;
    interval<T> gamma, pi;
    pi=kv::constants<interval<T> >::pi();
    T rad;
    gamma=1.;
    while(abs(q)<=exp(-gamma)){
      gamma=gamma+1;
    }
    while(abs(z)>=exp((N+1)*gamma)){
      N=N+10;
    }
    while(abs(z)*pow(q,N+1)>=1){
      N=N+10;
    }

    for(int n=0;n<=N-1;n++){
      sum=sum+zz*pow(q,n*(n-1)*0.5)/qPochhammer(q,q,n);

      zz=-z*zz;
    }
    rad=(pow(abs(z),N)*exp(pi*pi/6./gamma-N*(N-1)*gamma*0.5)/(1-abs(z)*exp(-(N+1))*gamma)).upper();
    res=complex_nbd(sum,rad);
    return res;
  }

}
#endif