hlliu20 / matrix Goto Github PK

def cp_x4_1(a, b, c, d, e): if a == 0: re = cp_x4(b, c, d, e) else : b = b/a c = c/a d = d/a e = e/a y_li = cp_x3(-1, c, 4*e-b*d, b*b*e-4*c*e+d*d) print(y_li) y = y_li[2] x12 = cp_x2(1, b/2-(b*b/4+c+y)**(1/2), y/2-(y*y/4-e)**(1/2)) x34 = cp_x2(1, b/2+(b*b/4+c+y)**(1/2), y/2+(y*y/4-e)**(1/2)) re = [x12[0]+x34[0], x12[1]+x34[1], x12[2], x12[3], x34[2], x34[3]] return re

def cp_x4(a, b, c, d, e): if a == 0: re = cp_x3(b, c, d, e) else : P=(c*c+12*a*e-3*b*d)/9 Q=(27*a*d*d+2*c*c*c+27*b*b*e-72*a*c*e-9*b*c*d)/54 D=(Q*Q-P*P*P)**(1/2) u=(Q+D)**(1/3) if abs((Q+D)**(1/3)) > abs((Q-D)**(1/3)) else (Q-D)**(1/3) if u==0: v=0 else : v=P/u w = -0.5+3**(1/2)/2j mst={} for k in [1, 2, 3]: m = (b*b-8*a*c/3+4*a*(w**(k-1)*u+w**(4-k)*v))**(1/2) if m==0: s, t = b*b-8*a*c/3, 0 else : s = 2*b*b-16*a*c/3-4*a*(w**(k-1)*u+w**(4-k)*v) t = (8*a*b*c-16*a*a*d-2*b*b*b)/m mst[m]=(s, t) m_li = mst.keys() m=0 for it in m_li: if abs(it)>abs(m): m= it s, t = mst[m][0], mst[m][1] x=[] for i in [1, 2, 3, 4]: x.append((-1*b+(-1)**(math.ceil(i/2))*m+(-1)**(i+1)*(s+(-1)**(math.ceil(i/2))*t)**(1/2))/(4*a)) re = [0, 0] for j in x: if type(j)==complex: if abs(j.imag) > 1e-6: re[1]+=1 re.append(j) else : re[0]+=1 re.append(j.real) else : re[0]+=1 re.append(j) return re