Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

convert function to mpmath

I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.

Here is the function:

from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
    num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
    denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
    return num/denom

def myfisher_mp(a1,a2,b1,b2):
    p = hyp_mp(a1,a2,b1,b2)
    outp = p
    for i in range(0,a1):
        temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
        if temp > p:
            break
        else:
            outp += temp
    for i in range(0,a2):
        temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
        if temp > p:
            break
        else:
            outp += temp
    return outp

A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.

click to hide/show revision 2
No.2 Revision

updated 2011-03-04 12:16:14 -0600

DSM gravatar image

convert mpmath function to mpmathcython

I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.

Here is the function:

from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
    num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
    denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
    return num/denom

def myfisher_mp(a1,a2,b1,b2):
    p = hyp_mp(a1,a2,b1,b2)
    outp = p
    for i in range(0,a1):
        temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
        if temp > p:
            break
        else:
            outp += temp
    for i in range(0,a2):
        temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
        if temp > p:
            break
        else:
            outp += temp
    return outp

A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.