In [1]:

    

from numpy import roots, random, real, imag
from math import atan, fabs, pi, tan, cos, sin, sqrt
import sympy
from sympy import Symbol

def cubic (Qx, Qy, Ux, Uy):
    
    a = Uy
    b = (Ux + 2*Qy)
    c = (2*Qx - Uy)
    d = -Ux
    
    det = -4*b*b*b*d + b*b*c*c -4*a*c*c*c + 18*a*b*c*d - 27*a*a*d*d
    
    if (det < 0):
        return 'c'
        
    
    if (det > 0):
        
        a = roots([a, b, c, d])
        a = a.real
        a = [atan(a[0]), atan(a[1]), atan(a[2])]
        
        U = [Ux*cos(a[0]) + Uy*sin(a[0]), Ux*cos(a[1]) + Uy*sin(a[1]), Ux*cos(a[2]) + Uy*sin(a[2])]
        rightU = [2*sin(a[0])*cos(a[0]), 2*sin(a[1])*cos(a[1]), 2*sin(a[2])*cos(a[2])]
    
        for i in range(0, 3):
            if (U[i] * rightU[i] < 0):
                a[i] = a[i] + pi
        
        a = sorted(a)
        a = [a[0] - a[0], a[1] - a[0], a[2] - a[0]]
        
        #print a
        
        if (a[2] > pi):
            return 'a'
        else:
            return 'b'


    

In [2]:

    

mu, sigma = 0, 1


    

In [9]:

    

Ux = random.normal(mu, sigma, 1000000)
Uy = random.normal(mu, sigma, 1000000)
Qx = random.normal(mu, sigma, 1000000)
Qy = random.normal(mu, sigma, 1000000)


    

In [4]:

    

a, b, c = 0, 0, 0

for i in range(0, 100000):
    ch = ''
    ch = cubic(Qx[i], Qy[i], Ux[i], Uy[i])
    if (ch == 'a'):
        a = a + 1
    if (ch == 'b'):
        b = b + 1
    if (ch == 'c'):
        c = c + 1
        
print a/100000.0, b/100000.0, c/100000.0


    

    




0.527 0.09225 0.38075

In [5]:

    

def cubic_solver (Qx, Qy, Ux, Uy):
    
    a = (Uy)
    b = (Ux + 2*Qy)
    c = (2*Qx - Uy)
    d = (-Ux)
    
    return (roots([a, b, c, d]))


    

In [10]:

    

# Another test
saddles, beaks, comets = 0, 0, 0

for i in range(0, 1000):
    
    root1, root2, root3 = cubic_solver(Qx[i], Qy[i]+0*i, Ux[i]+0*i, Uy[i]+0*i)
    
    if (Qx[i]*Uy[i] - Qy[i]*Ux[i] > 0):
        saddles = saddles + 1
    elif (imag(root1) == 0 and imag(root2) == 0 and imag(root3) == 0):
        beaks = beaks + 1
        
print saddles/1000.0, beaks/1000.0


    

    




0.5 0.5

In [ ]: