1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
#!
# -*- coding: utf-8 -*-
#
# Polynomial generator using sage
import subprocess, random, sys
# We need imaginary unit over Q
I = QuadraticField (-1, 'I').gen ()
random.seed()
class PolynomialGenerator():
def __init__(self):
self._degree = 0
self._x = var('x')
self._pol = 1
self._eps = Rational("1 / 100000000")
self._comment = "! Polygen generated polynomial\n"
self._comment += "! This polynomial has the following roots: \n"
def add_complex_root (self, complex_root, mult = 1):
self._pol *= (self._x - complex_root)**mult
self._comment += "!# (%s)\n" % str (N(complex_root, digits = 200)).replace (" +", "e0, ").replace (" -", "e0, - ").replace("*I", "e0")
def add_complex_conjugate_roots (self, complex_root, mult = 1):
self._pol *= (self._x - complex_root)**mult
self._pol *= (self._x - complex_root.conjugate ())**mult
def add_cluster (self, cluster_center, n):
for i in range(n):
gamma = Rational(RR(cos(2*i*pi/n))) + Rational(RR(sin(2*i*pi/n))) * I
gamma = gamma * self._eps
self.add_complex_root (cluster_center + gamma)
def get_pol_file (self):
self._pol = expand (self._pol)
coeffs = list (map (lambda x : x[0], self._pol.coefficients()))
imag_coeffs = list (map (lambda x : x.imag (), coeffs))
real = reduce (lambda x, y : x and y, map (lambda x : x == 0, imag_coeffs))
pol_file = self._comment + "\n"
pol_file += "Degree = %d;\n" % (len (coeffs) - 1)
# Check if the coefficients are real
if real:
pol_file += "Real;\n"
else:
pol_file += "Complex;\n"
# Determine poly structure
pol_file += "Rational;\n"
pol_file += "Monomial;\n"
pol_file += "\n"
for c in coeffs:
r = c.real ()
i = c.imag ()
pol_file += str (r) + " "
if not real:
pol_file += str (i) + " "
pol_file += "\n"
return pol_file
def usage ():
print "%s degree [type]" % sys.argv[0]
print ""
print " type can be one of: "
print " - separate: Well separate integer roots (Wilkinson-style), default choice"
print " - clustered: Polynomial with clusters of roots"
print " - mixed: Polynomial with both clusters and isolated roots"
print " - multiple: Polynomial with multiple roots"
print " - all: Polynomial with all the above"
print ""
sys.exit (1)
if __name__ == "__main__":
if (len (sys.argv) < 2 or len (sys.argv) > 3):
usage ()
polygen = PolynomialGenerator ()
n = int(sys.argv[1])
t = "separate"
if (len (sys.argv) == 3):
t = sys.argv[2]
integer_grid = range (1, 100)
if t == "separate":
for i in range(n):
polygen.add_complex_root (random.choice (integer_grid) +
random.choice (integer_grid) * I)
elif t == "clustered":
degree = 0
while degree < n:
step = random.choice (range (1, n - degree + 1))
polygen.add_cluster (random.choice (integer_grid) +
random.choice (integer_grid) * I, step)
degree += step
elif t == "mixed":
degree = 0
while degree < n:
if random.random() <= 0.7:
polygen.add_complex_root (random.choice (integer_grid) +
random.choice (integer_grid) * I)
degree += 1
else:
step = random.choice (range (1, n - degree + 1))
polygen.add_cluster (random.choice (integer_grid) +
random.choice (integer_grid) * I, step)
degree += step
elif t == "multiple":
pass
elif t == "all":
pass
else:
usage ()
print polygen.get_pol_file ()
|
