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