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