Pari-Python update

I have uploaded my current state of the project to sourceforge:

http://sourceforge.net/projects/pari-python/

Sorry, it is only available through CVS. I don’t have much time to make a binary distribution. I started to use it myself, here I post some example module with computations of some mock theta functions:

from pari import *
ser_prec(300)
q = var('q')
# q is in fact q^(1/8)
theta1 = eta(-q^4,0)^2/eta(q^8,0)
theta2 = eta(q^4,0)^2/eta(q^8,0)
theta3 = 2*q*eta(q^16,0)^2/eta(q^8,0)
eta3 = q*eta(q^8,0)^3
 def th1(x):
    return exp(-Pi*I/12)*eta((x+1)/2,1)^2/eta(x,1)
def th2(x):
    return eta(x/2,1)^2/eta(x,1)
def th3(x):
    return 2*eta(2*x,1)^2/eta(x,1)
def modprecision():
    return ser_prec()
def modG(k):
    res = -bernfrac(k)/2/k
    for i in range(1, ser_prec()-1):
        res += sigma(i,k-1)*q^i
    res += O(q^ser_prec())
    return res
def modE(k):
    return modG(k)*2*k/-bernfrac(k)

S = var('x')
# exp(Pi*I*s) in terms of 2*Pi*I*s
Es=exp(S/2)+O(S^5)
# qz = exp(Pi*I*z)

def th(qz):
    bound = int(float(sqrt(ser_prec())))
    res = O(q^ser_prec())
    for n in range(-bound, bound):
        res += cast(-1)^n*q^(4*(n^2+n))*qz^(2*n)
    res *= q*I*qz
    return res

def mu(qu, qv):
    bound = int(float(sqrt(ser_prec())))
    res = O(q^ser_prec())
    for n in range(-bound, bound):
        res += cast(-1)^n*q^(4*(n^2+n))*qv^(2*n)/(1-q^(8*n)*qu^2)
    res*=qu/th(qv)
    return res

mu1 = polcoeff(mu(I,q^2*Es),1)
mu2 = polcoeff(mu(I,I*q^2*Es),1)
mu3 = polcoeff(mu(q^2,I*q^2*Es),1)

mu12 = O(q^ser_prec())+sum(-10,10,lambda n: cast(-1)^n*(n^2/2-n/2+1/8)*q^(4*(n^2+n))/(1+q^(8*n-4)))
mu22 = O(q^ser_prec())+sum(-10,10,lambda n: cast(-1)^n*(n^2/2-n/2+1/8)*q^(4*(n^2+n))/(1-q^(8*n-4)))
mu32 = O(q^ser_prec())+sum(-10,10,lambda n: cast(-1)^n*(n^2)/2*q^(4*(n^2+n))/(1+q^(8*n)))

E2 = modE(2)

tau0 = I*1.231241+2.234235 # random point
print 'Verifying transformation rules of theta functions...'
print abs(th1(tau0+1) - th2(tau0))
print abs(th2(tau0+1) - th1(tau0))
print abs(th3(tau0+1) - exp(I*Pi/4)*th3(tau0))
print abs(th1(-1/tau0) - exp(-I*Pi/4)*sqrt(tau0)*th1(tau0))
print abs(th2(-1/tau0) - exp(-I*Pi/4)*sqrt(tau0)*th3(tau0))
print abs(th3(-1/tau0) - exp(-I*Pi/4)*sqrt(tau0)*th2(tau0))

def modeval(f,x,ord):
    return subst(Pol(f,q),q,exp(2*Pi*I*x/ord))

h1=mu12/q^3/eta3
h2=-mu22/q^3/eta3
h3=-mu32/eta3

def H1(x):
    return modeval(h1,x,8)
def H2(x):
    return modeval(h2,x,8)
def H3(x):
    return modeval(h3,x,8)

def EE(x):
    return 2*intnum(0,x,lambda t: exp(-Pi*t^2))

def beta(x):
    return erfc(sqrt(Pi*x))

def RR(u,x):
    y = imag(x)
    a = imag(u)/imag(x)
    res = 0
    for n in range(-20,20):
        nu=n+1/2
        s1=sign(nu)
        s2=sign(nu+a)
        k = s1-EE((nu+a)*sqrt(2*y))
        if s1==s2:
            k = s1*beta((nu+a)^2*2*y)
        res += cast(-1)^n*exp(-Pi*I*nu^2*x-2*Pi*I*nu*u)*k
    return res

def RRd(x):
    u=1/2-x/2
    y = imag(x)
    a = -1/2
    res = 0
    for n in range(-20,20):
        nu=n+1/2
        s1=sign(n)
        k = s1*beta(n^2*2*y)
        res += I*exp(-Pi*I*n^2*x)*(-k*n-sqrt(2*y)/2/I/y*2*exp(-Pi*2*y*(nu+a)^2)/2/Pi/I)
    return res

def R1(x):
    y = imag(x)
    res = -1/2/Pi/sqrt(y*2)*conj(th1(x))
    for n in range(-20,20):
        res += abs(n)*beta(2*y*n^2)*exp(-Pi*I*n^2*x)/2
    return res

def R2(x):
    y = imag(x)
    res = -1/2/Pi/sqrt(y*2)*conj(th2(x))
    for n in range(-20,20):
        res += cast(-1)^n*abs(n)*beta(2*y*n^2)*exp(-Pi*I*n^2*x)/2
    return res

def R3(x):
    y = imag(x)
    res = 1/2/Pi/sqrt(y*2)*conj(th3(x))
    for n in range(-20,20):
        res -= abs(n+1/2)*beta(2*y*(n+1/2)^2)*exp(-Pi*I*(n+1/2)^2*x)/2
    return res

def M1(x):
    return H1(x)+R1(x)
def M2(x):
    return H2(x)+R2(x)
def M3(x):
    return H3(x)+R3(x)
def M(x):
    return M1(x)*th1(x)+M2(x)*th2(x)-M3(x)*th3(x)

print 'Verifying transformation rules of M functions...'
print abs(M1(tau0+1) - M2(tau0))
print abs(M2(tau0+1) - M1(tau0))
print abs(M3(tau0+1) - exp(-I*Pi/4)*M3(tau0))
print abs(M1(-1/tau0) - exp(I*Pi/4)*sqrt(tau0)^3*M1(tau0))
print abs(M2(-1/tau0) + exp(I*Pi/4)*sqrt(tau0)^3*M3(tau0))
print abs(M3(-1/tau0) + exp(I*Pi/4)*sqrt(tau0)^3*M2(tau0))

 

About these ads