upiterbarg / mpmath

Automatically exported from code.google.com/p/mpmath
Other
0 stars 0 forks source link

Divide by Zero Error in ellippi #228

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago
What steps will reproduce the problem?
1. ellippi(1,1)
2. ellippi(1,0)
3. ellippi(1,>pi/2,0 or 1)

What is the expected output? What do you see instead?
Divide by Zero Error and nan

What version of the product are you using? On what operating system?
Python 3.2.3, mpmath 0.17, win7pro 64bit

Please provide any additional information below.

I have not used either Python or mpmath before, so my code fix may not be 
correct or optimal. The code seems to work and is at the end of the elliptic.py 
module.

        if complete:
            #if m == 0: return ctx.pi/(2*ctx.sqrt(1-n)) #removed by JCB, div/0 error
            if m == 0:                                  #added by JCB
                if n == 1: return ctx.inf               #added by JCB 
                return ctx.pi/(2*ctx.sqrt(1-n))         #added by JCB
            if n == 0: return ctx.ellipk(m)
            if ctx.isinf(n) or ctx.isinf(m): return ctx.zero
        else:
            if z == 0: return z
            if ctx.isinf(n): return ctx.zero
            if ctx.isinf(m): return ctx.zero
        if ctx.isinf(n) or ctx.isinf(z) or ctx.isinf(m):
            raise ValueError
    if complete:
        #if m == 1: return -ctx.inf/ctx.sign(n-1)   #removed by JCB, div/0 error
        if m == 1:                                  #added by JCB
            if n > 1: return -ctx.inf               #added by JCB
            return ctx.inf                          #added by JCB                
        away = False
    else:
        x = z.real
        ctx.prec += max(0, ctx.mag(x))
        pi = +ctx.pi
        away = abs(x) > pi/2
    if away:
        d = ctx.nint(x/pi)
        z = z-pi*d
        P = 2*d*ctx.ellippi(n,m)
        if ctx.isinf(P): return P   #added by JCB, return "(+/-)inf" instead of "nan"
    else:

Original issue reported on code.google.com by JohnBey...@gmail.com on 6 Aug 2012 at 12:43

GoogleCodeExporter commented 9 years ago
Thanks! I have fixed it here: 
https://github.com/fredrik-johansson/mpmath/commit/fc474e10694ea0800ef568afd3e7b
f8f4ab1142c

Original comment by fredrik....@gmail.com on 9 Aug 2012 at 9:56

GoogleCodeExporter commented 9 years ago
Fredrik,
I downloaded and installed your fix. There might still be some issues.
ellippi(1,-2,1) probably should return -inf.
Also, test_linalg.py fails for the print statements and the xrange references.
I have only installed the minimum modules required for mpmath.
Here is an alternate ending for the ellippi function returning either 2 or 3 
terms.
    if away:
        d = ctx.nint(x/pi)
        z = z-pi*d
        #P = 2*d*ctx.ellippi(n,m)
        #if ctx.isinf(P):
            #return ctx.inf
    #else:
        #P = 0
    def terms():
        if complete:
            c, s = ctx.zero, ctx.one
        else:
            c, s = ctx.cos_sin(z)
        x = c**2
        y = 1-m*s**2
        RF = ctx.elliprf(x, y, 1)
        RJ = ctx.elliprj(x, y, 1, 1-n*s**2)
        if away: return s*RF, n*s**3*RJ/3, 2*d*ctx.ellippi(n,m)
        return s*RF, n*s**3*RJ/3
    return ctx.sum_accurately(terms)

Original comment by JohnBey...@gmail.com on 12 Aug 2012 at 3:29