1#!/usr/bin/env python
2
3import sys
4import polypy
5import polypy_test
6
7polypy_test.init()
8
9def check_roots(p, assignment, roots, expected_roots):
10    ok = len(roots) == len(expected_roots)
11    if ok:
12        for root, expected_root in zip(roots, expected_roots):
13            root_double = root.to_double()
14            ok = abs(root_double - expected_root) < 0.000001
15            if (not ok):
16                print("p = {0}".format(p))
17                print("assignment = {0}".format(assignment))
18                print("root = {0}".format(root))
19                print("expected_root = {0}".format(expected_root))
20                break
21    else:
22        print("p = {0}".format(p))
23        print("assignment = {0}".format(assignment))
24        print("roots = {0}".format(roots))
25        print("expected_roots = {0}".format(expected_roots))
26    polypy_test.check(ok)
27
28# polypy.trace_enable("polynomial")
29# polypy.trace_enable("polynomial::expensive")
30# polypy.trace_enable("factorization")
31# polypy.trace_enable("algebraic_number")
32# polypy.trace_enable("coefficient")
33# polypy.trace_enable("coefficient::sgn")
34# polypy.trace_enable("coefficient::roots")
35# polypy.trace_enable("roots")
36
37
38[x, y, z, w] = [polypy.Variable(name) for name in ['x', 'y', 'z', 'w']]
39[x0, x1, x2, x3, x4, x5, x6, x7] = [polypy.Variable(name) for name in ['x0', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7']]
40
41# polypy_test.start("Polynomial Root Isolation: Speed")
42#
43# polypy.variable_order.set([y, x])
44#
45# y_value_poly = 8192*x**6 + (-14336*x**5) + 75376*x**4 + (-32736*x**3) + 109496*x**2 + 133752*x + (-32441)
46# y_value = polypy.AlgebraicNumber(y_value_poly, 1)
47#
48# assignment = polypy.Assignment()
49# assignment.set_value(y, y_value)
50#
51# p = (1*y**2 + 3)*x**6 + (4*y**2 + 24*y + 20)*x**5 + (36*y**2 + 112*y + 44)*x**4 + (8*y**3 + 114*y**2 + 144*y + 30)*x**3 + (1*y**4 + 24*y**3 + 142*y**2 + 40*y - 15)*x**2 + (2*y**4 + 40*y**3 + 64*y**2 - 32*y - 26)*x + (5*y**4 + 16*y**3 + 7*y**2 - 16*y - 8)
52# roots = p.roots_isolate(assignment)
53
54polypy_test.start("Polynomial Root Isolation: Bugs");
55
56polypy.variable_order.set([x0, x1, x2, x3, x4, x5, x6, x7])
57p = ((16*x0**2 - 32*x0 + 16)*x5**2 + ((32*x0 - 32)*x1)*x5 + (16*x1**2))*x7**2 + ((8*x0**2 - 16*x0 + 8)*x5**2 + ((16*x0**2 - 32*x0 + 16)*x2**2 + (-8*x1**2)))*x7 + ((1*x0**2 - 2*x0 + 1)*x5**2 + ((-2*x0 + 2)*x1)*x5 + (1*x1**2))
58
59p_lc = ((16*x0**2 - 32*x0 + 16)*x5**2 + ((32*x0 - 32)*x1)*x5 + (16*x1**2))
60
61assignment = polypy.Assignment()
62v = polypy.AlgebraicNumber(32*x**2 + (-64*x) + 15, 0)
63assignment.set_value(x0, v)
64assignment.set_value(x1, 28);
65assignment.set_value(x2, -1, 4)
66assignment.set_value(x3, 3)
67assignment.set_value(x4, 277771, 4096)
68assignment.set_value(x5, 786743, 20480)
69assignment.set_value(x6, -6437)
70
71roots = p.roots_isolate(assignment);
72check_roots(p, assignment, roots, [-166318.899794716301953293454440, -8961.48539760581031126063921576])
73
74polypy.variable_order.set([x, w, y, z])
75
76sqrt180 = polypy.AlgebraicNumber(x**2 - 180, 0)
77assignment = polypy.Assignment()
78assignment.set_value(x, 9)
79assignment.set_value(w, sqrt180)
80assignment.set_value(y, 0)
81
82p = ((1*x)*y)*z + (-1*w**2 + (20*x))
83
84roots = p.roots_isolate(assignment)
85check_roots(p, assignment, roots, []);
86
87
88polypy_test.start("Polynomial Root Isolation: Basic")
89
90polypy.variable_order.set([z, y, x])
91
92
93sqrt2 = polypy.AlgebraicNumber(x**2 - 2, 1)
94sqrt3 = polypy.AlgebraicNumber(x**2 - 3, 1)
95
96assignment = polypy.Assignment()
97assignment.set_value(y, sqrt2)
98assignment.set_value(z, sqrt3)
99
100p = ((1*z)*y**2 + (1*z**2)*y)*x**2 + ((1*z**2)*y**2)*x + 1
101roots = p.roots_isolate(assignment)
102check_roots(p, assignment, roots, [-0.536830573034912, -0.241708498946619])
103
104# print assignment
105
106p = (x**2 - y**2)*(x**2 - z**2)
107roots = p.roots_isolate(assignment);
108check_roots(p, assignment, roots, [-1.732050808, -1.414213562, 1.414213562, 1.732050808])
109
110assignment = polypy.Assignment()
111assignment.set_value(y, 0)
112assignment.set_value(z, 0)
113
114# print assignment
115
116p = 10*(x-1)*(x**2 - y - z - 1)
117roots = p.roots_isolate(assignment)
118check_roots(p, assignment, roots, [-1, 1])
119
120p = x**2 - y - z - 1
121roots = p.roots_isolate(assignment)
122check_roots(p, assignment, roots, [-1, 1])
123
124p = x + y + z
125roots = p.roots_isolate(assignment)
126check_roots(p, assignment, roots, [0])
127