1#!/usr/bin/env python 2 3import polypy 4import polypy_test 5 6x = polypy.Variable("x"); 7y = polypy.Variable("y"); 8z = polypy.Variable("z"); 9w = polypy.Variable("w"); 10 11polypy.variable_order.set([w, z, y, x]); 12 13def cmp3(lhs, rhs): 14 return ((lhs > rhs) - (lhs < rhs)) 15 16def check_psc(p, q, expected): 17 psc = p.psc(q) 18 ok = cmp3(psc, expected) == 0 19 if (not ok): 20 print("p = {0}".format(p)) 21 print("q = {0}".format(q)) 22 print("psc = {0}".format(psc)) 23 print("expected = {0}".format(expected)) 24 polypy_test.check(ok) 25 26def check_resultant(p, q, expected): 27 resultant = p.resultant(q) 28 ok = resultant == expected 29 if (not ok): 30 print("p = {0}".format(p)) 31 print("q = {0}".format(q)) 32 print("resultant = {0}".format(resultant)) 33 print("expected = {0}".format(expected)) 34 35polypy_test.init() 36 37# polypy.trace_enable("coefficient::resultant") 38# polypy.trace_enable("polynomial::expensive") 39 40polypy_test.start("Speed") 41 42x5 = polypy.Variable("x5") 43x9 = polypy.Variable("x9") 44x10 = polypy.Variable("x10") 45x11 = polypy.Variable("x11") 46x12 = polypy.Variable("x12") 47 48polypy.variable_order.set([x9, x12, x11, x10, x5]); 49 50A = (((1*x12**6 + (4*x9)*x12**5 + (10*x9**2)*x12**4 + (14*x9**3)*x12**3 + (13*x9**4)*x12**2 + (6*x9**5)*x12 + (1*x9**6))*x11**2)*x10**2 + (((2*x9**2)*x12**6 + (4*x9**3)*x12**5 + (6*x9**4)*x12**4 + (2*x9**5)*x12**3)*x11)*x10 + ((1*x9**4)*x12**6))*x5**2 + ((((-1*x9)*x12**3 + (-1*x9**2)*x12**2 + (-1*x9**3)*x12)*x11**4)*x10**4 + ((-2*x12**6 + (-9*x9)*x12**5 + (-17*x9**2)*x12**4 + (-21*x9**3)*x12**3 + (-12*x9**4)*x12**2 + (-4*x9**5)*x12)*x11**3)*x10**3 + (((-2*x9**2)*x12**6 + (-6*x9**3 - 2*x9)*x12**5 + (-4*x9**4 - 8*x9**2)*x12**4 + (-3*x9**5 - 16*x9**3)*x12**3 + (-1*x9**6 - 18*x9**4)*x12**2 + (-1*x9**7 - 10*x9**5)*x12 + (-2*x9**6))*x11**2)*x10**2 + (((-1*x9**5 - 2*x9**3)*x12**5 + (-1*x9**6 - 4*x9**4)*x12**4 + (-1*x9**7 - 2*x9**5)*x12**3)*x11)*x10)*x5 + ((((1*x9)*x12**3 + (1*x9**2)*x12**2)*x11**5)*x10**5 + ((1*x12**6 + (5*x9)*x12**5 + (7*x9**2)*x12**4 + (6*x9**3)*x12**3 + (3*x9**4)*x12**2 + (1*x9**3)*x12)*x11**4)*x10**4 + (((2*x9**3 + 2*x9)*x12**5 + (2*x9**4 + 7*x9**2)*x12**4 + (1*x9**5 + 13*x9**3)*x12**3 + (1*x9**6 + 8*x9**4)*x12**2 + (4*x9**5)*x12)*x11**3)*x10**3 + (((-1*x9**2)*x12**6 + (1*x9**5)*x12**5 + (1*x9**6 - 2*x9**4 + 1*x9**2)*x12**4 + (2*x9**5 + 4*x9**3)*x12**3 + (6*x9**4)*x12**2 + (1*x9**7 + 4*x9**5)*x12 + (1*x9**6))*x11**2)*x10**2 + (((-2*x9**4)*x12**6 + (-1*x9**6)*x12**4 + (1*x9**7)*x12**3)*x11)*x10 + ((-1*x9**6)*x12**6)) 51B = (((2*x12**6 + (8*x9)*x12**5 + (20*x9**2)*x12**4 + (28*x9**3)*x12**3 + (26*x9**4)*x12**2 + (12*x9**5)*x12 + (2*x9**6))*x11**2)*x10**2 + (((4*x9**2)*x12**6 + (8*x9**3)*x12**5 + (12*x9**4)*x12**4 + (4*x9**5)*x12**3)*x11)*x10 + ((2*x9**4)*x12**6))*x5 + ((((-1*x9)*x12**3 + (-1*x9**2)*x12**2 + (-1*x9**3)*x12)*x11**4)*x10**4 + ((-2*x12**6 + (-9*x9)*x12**5 + (-17*x9**2)*x12**4 + (-21*x9**3)*x12**3 + (-12*x9**4)*x12**2 + (-4*x9**5)*x12)*x11**3)*x10**3 + (((-2*x9**2)*x12**6 + (-6*x9**3 - 2*x9)*x12**5 + (-4*x9**4 - 8*x9**2)*x12**4 + (-3*x9**5 - 16*x9**3)*x12**3 + (-1*x9**6 - 18*x9**4)*x12**2 + (-1*x9**7 - 10*x9**5)*x12 + (-2*x9**6))*x11**2)*x10**2 + (((-1*x9**5 - 2*x9**3)*x12**5 + (-1*x9**6 - 4*x9**4)*x12**4 + (-1*x9**7 - 2*x9**5)*x12**3)*x11)*x10) 52A.psc(B) 53 54polypy.variable_order.set([w, z, y, x]); 55 56A = w - 2*x**3*y**2*z - x**2*y*z**3 - x*y**3*z**2 - z - 4 57B = 2304*x**12 - 5184*x**10 + 2304*x**9 - 3996*x**8 - 2592*x**7 + 5376*x**6 - 3456*x**5 + 3132*x**4 + 2832*x**3 - 2736*x**2 + 169 58resultant = A.resultant(B) 59 60polypy_test.start("Bugs") 61 62polypy.variable_order.set([w, z, y, x]); 63 64p = 7*x**6 + 14*x**5 + (7*y - 1) 65q = 42*x**5 + 70*x**4 66expected = [92254156521408*y**5-461361174686720*y**4+244807578081920*y**3-51113953954560*y**2+4803956911040*y-170197631744, 67 174327582240*y**3-74711820960*y**2+10673117280*y-508243680, 68 0, 69 0, 70 -6860, 71 42] 72check_psc(p, q, expected) 73 74polypy_test.start("Resultants") 75 76p_sqrt2 = x**2 - 2 77p_sqrt3 = y**2 - 3 78 79add = z - (x + y) 80res1 = add.resultant(p_sqrt2) 81check_resultant(res1, p_sqrt3, 1*z**4 - 10*z**2 + 1) 82 83sub = z - (x - y) 84res1 = sub.resultant(p_sqrt2) 85check_resultant(res1, p_sqrt3, 1*z**4 - 10*z**2 + 1) 86 87mul = z - x*y 88res1 = mul.resultant(p_sqrt2) 89check_resultant(res1, p_sqrt3, 1*z**4 - 12*z**2 + 36) 90 91polypy_test.start("Principal Sub-resultant Coefficients") 92 93p = (y-1)*z*x**2 + y*(z-1)*x + y*z 94q = z*(y-1)*x + y*z 95expected = [-y**2*z**2 + y**3*z**2 + y*z**3 - 2*y**2*z**3 + y**3*z**3, -z + y*z] 96check_psc(p, q, expected) 97 98p = (y-1)*z*x**2 + y*(z-1)*x + y*z 99q = (z-1)*y*x**2 + z*(y-1)*x + y*z 100expected = [-y**4*z + y**3*z**2 + 3*y**4*z**2 + y**2*z**3 - 6*y**3*z**3 - y*z**4 + 3*y**2*z**4, 101 -y**2 + 2*y**2*z + z**2 - 2*y*z**2, 102 1] 103check_psc(p, q, expected) 104 105p = (y-3)*x**3 + (y-2)*x**2 + (y-1)*x + y 106q = (z-3)*x**3 + (z-2)*x**2 + (z-1)*x + z 107expected = [16*y**3 - 48*y**2*z + 48*y*z**2 - 16*z**3, 108 0, 109 y - z, 110 1] 111check_psc(p, q, expected) 112expected = [-16*y**3 + 48*y**2*z - 48*y*z**2 + 16*z**3, 113 0, 114 z - y, 115 1] 116check_psc(q, p, expected) 117 118p = x**2 + y**2 - 1 119q = p.derivative() 120expected = [-4 + 4*y**2, 2] 121check_psc(p, q, expected) 122check_psc(q, p, expected); 123 124p = (y-1)*x**2 125q = (y-2)*x**2 126expected = [0, 0, 1] 127check_psc(p, q, expected) 128 129p = x**2 130q = x**2 131expected = [0, 0, 1] 132check_psc(p, q, expected) 133 134p = x**2 135q = x 136expected = [0, 1] 137check_psc(p, q, expected) 138 139p = x**2 + y**2 + z**2 - 1 140q = p.derivative() 141expected = [-4 + 4*y**2 + 4*z**2, 2] 142check_psc(p, q, expected) 143 144p = (y-3)*x**3 + (y-2)*x**2 + (y-1)*x + y 145q = (z-2)*x**2 + (z-1)*x + z 146expected = [-3*y - 5*y**2 + 3*z - y*z + 7*y**2*z + 6*z**2 - y*z**2 - 3*y**2*z**2 + 3*z**3 - 3*y*z**3 + y**2*z**3, 147 -3 + 3*y - 2*z - y*z + z**2, 148 -2 + z] 149check_psc(p, q, expected) 150 151p = z*x**3 + (y - 1)*x**2 + z*y 152q = y*x**3 + (z - 1)*x**2 + z*y 153expected = [y**5*z**2 - y**6*z**2 - 3*y**4*z**3 + 2*y**5*z**3 - y**6*z**3 + 3*y**3*z**4 + 3*y**5*z**4 - y**2*z**5 - 2*y**3*z**5 - 3*y**4*z**5 + y**2*z**6 + y**3*z**6, 154 y**3*z - y**4*z - 2*y**2*z**2 + y**3*z**2 + y*z**3 + y**2*z**3 - y*z**4, 155 y - y**2 - z + z**2, 156 1] 157check_psc(p, q, expected) 158 159p = z*x**4 + (y - 1)*x**3 + z*y 160q = y*x**4 + (z - 1)*x**3 + z*y 161expected = [-y**7*z**3 + y**8*z**3 + 4*y**6*z**4 - 3*y**7*z**4 + y**8*z**4 - 6*y**5*z**5 + 2*y**6*z**5 - 4*y**7*z**5 + 4*y**4*z**6 + 2*y**5*z**6 + 6*y**6*z**6 - y**3*z**7 - 3*y**4*z**7 - 4*y**5*z**7 + y**3*z**8 + y**4*z**8, 162 y**5*z**2 - y**6*z**2 - 3*y**4*z**3 + 2*y**5*z**3 + 3*y**3*z**4 - y**2*z**5 - 2*y**3*z**5 + y**2*z**6, 163 0, 164 y - y**2 - z + z**2, 165 1] 166check_psc(p, q, expected) 167