NZMATH  1.2.0 About: NZMATH is a Python based number theory oriented calculation system.   Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)
equation.py
Go to the documentation of this file.
1 """
2 equation -- methods to solve algebraic equations.
3
4 If you would like to solve an equation in the complex field and need
5 more precision, then you can import and use plugins.SETPRECISION.
6
7 e.g:
8 from nzmath.plugins import SETPRECISION
9 SETPRECISION(200)
10
11 Then, the following computations will be done in such precision, if
12 plugins.PRECISION_CHANGEABLE is true.
13 """
14
15 from __future__ import division
16 import logging
17
18 import nzmath.arith1 as arith1
19 import nzmath.bigrandom as bigrandom
20 import nzmath.bigrange as bigrange
21 import nzmath.finitefield as finitefield
22 import nzmath.imaginary as imaginary
23 import nzmath.poly.uniutil as uniutil
24
25 from nzmath.plugins import MATHMODULE as math, CMATHMODULE as cmath, \
26  FLOATTYPE as Float, COMPLEXTYPE as Complex
27
28 _log = logging.getLogger('nzmath.equation')
29 _log.setLevel(logging.DEBUG)
30
31 # x is (list,tuple)
32 # t is variable
33 def e1(x):
34  """
35  0 = x[0] + x[1]*t
36  """
37  if x[1] == 0:
38  raise ZeroDivisionError("No Solution")
39  else:
40  return -x[0]/x[1]
41
42 def e1_ZnZ(x, n):
43  """
44  Return the solution of x[0] + x[1]*t = 0 (mod n).
45  x[0], x[1] and n must be positive integers.
46  """
47  try:
48  return (-x[0] * arith1.inverse(x[1], n)) % n
49  except ZeroDivisionError:
50  raise ValueError("No Solution")
51
52 def e2(x):
53  """
54  0 = x[0] + x[1]*t + x[2]*t**2
55  """
56  c, b, a = x
57  d = b**2 - 4*a*c
58  if d >= 0:
59  sqrtd = math.sqrt(d)
60  else:
61  sqrtd = cmath.sqrt(d)
62  return ((-b + sqrtd)/(2*a), (-b - sqrtd)/(2*a))
63
64 def e2_Fp(x, p):
65  """
66  p is prime
67  f = x[0] + x[1]*t + x[2]*t**2
68  """
69  c, b, a = [_x % p for _x in x]
70  if a == 0:
71  return [e1_ZnZ([c, b], p)]
72  if p == 2:
73  solutions = []
74  if x[0] & 1 == 0:
75  solutions.append(0)
76  if (x[0] + x[1] + x[2]) & 1 == 0:
77  solutions.append(1)
78  if len(solutions) == 1:
79  return solutions * 2
80  return solutions
81  d = b**2 - 4*a*c
82  if arith1.legendre(d, p) == -1:
83  return []
84  sqrtd = arith1.modsqrt(d, p)
85  a = arith1.inverse(2*a, p)
86  return [((-b+sqrtd)*a)%p, ((-b-sqrtd)*a)%p]
87
88 def e3(x):
89  """
90  0 = x[0] + x[1]*t + x[2]*t**2 + x[3]*t**3
91  """
92  x3 = Float(x[3])
93  a = x[2] / x3
94  b = x[1] / x3
95  c = x[0] / x3
96  p = b - (a**2)/3
97  q = 2*(a**3)/27 - a*b/3 + c
98  w = (-1 + cmath.sqrt(-3)) / 2
99  W = (1, w, w.conjugate())
100  k = -q/2 + cmath.sqrt((q**2)/4 + (p**3)/27)
101  l = -q/2 - cmath.sqrt((q**2)/4 + (p**3)/27)
102  m = k ** (1/3)
103  n = l ** (1/3)
104
105  # choose n*W[i] by which m*n*W[i] be the nearest to -p/3
106  checker = [abs(3*m*n*z + p) for z in W]
107  n = n * W[checker.index(min(checker))]
108
109  sol = []
110  for i in range(3):
111  sol.append(W[i]*m + W[-i]*n - a/3)
112  return sol
113
114 def e3_Fp(x, p):
115  """
116  p is prime
117  0 = x[0] + x[1]*t + x[2]*t**2 + x[3]*t**3
118  """
119  x.reverse()
120  lc_inv = finitefield.FinitePrimeFieldElement(x[0], p).inverse()
121  coeff = []
122  for c in x[1:]:
123  coeff.append((c * lc_inv).n)
124  sol = []
125  for i in bigrange.range(p):
126  if (i**3 + coeff[0]*i**2 + coeff[1]*i + coeff[2]) % p == 0:
127  sol.append(i)
128  break
129  if len(sol) == 0:
130  return sol
131  X = e2_Fp([coeff[1] + (coeff[0] + sol[0])*sol[0], coeff[0] + sol[0], 1], p)
132  if len(X) != 0:
133  sol.extend(X)
134  return sol
135
136 def Newton(f, initial=1, repeat=250):
137  """
138  Compute x s.t. 0 = f[0] + f[1] * x + ... + f[n] * x ** n
139  """
140  length = len(f)
141  df = []
142  for i in range(1, length):
143  df.append(i * f[i])
144  l = initial
145  for k in range(repeat):
146  coeff = Float(0)
147  dfcoeff = Float(0)
148  for i in range(1, length):
149  coeff = coeff * l + f[-i]
150  dfcoeff += dfcoeff * l + df[-i]
151  coeff = coeff * l + f[0]
152  if coeff == 0:
153  return l
154  elif dfcoeff == 0: # Note coeff != 0
155  raise ValueError("There is not solution or Choose different initial")
156  else:
157  l = l - coeff / dfcoeff
158  return l
159
160
161 def SimMethod(f, NewtonInitial=1, repeat=250):
162  """
163  Return zeros of a polynomial given as a list.
164  """
165  if NewtonInitial != 1:
166  ni = NewtonInitial
167  else:
168  ni = None
169  return _SimMethod(f, newtoninitial=ni, repeat=repeat)
170
171
172 def _SimMethod(g, initials=None, newtoninitial=None, repeat=250):
173  """
174  Return zeros of a polynomial given as a list.
175
176  - g is the list of the polynomial coefficient in ascending order.
177  - initial (optional) is a list of initial approximations of zeros.
178  - newtoninitial (optional) is an initial value for Newton method to
179  obtain an initial approximations of zeros if 'initial' is not given.
180  - repeat (optional) is the number of iteration. The default is 250.
181  """
182  if initials is None:
183  if newtoninitial:
184  z = _initialize(g, newtoninitial)
185  else:
186  z = _initialize(g)
187  else:
188  z = initials
189
190  f = uniutil.polynomial(enumerate(g), imaginary.theComplexField)
191  deg = f.degree()
192  df = f.differentiate()
193
194  value_list = [f(z[i]) for i in range(deg)]
195  for loop in range(repeat):
196  sigma_list = [0] * deg
197  for i in range(deg):
198  if not value_list[i]:
199  continue
200  sigma = 0
201  for j in range(i):
202  sigma += 1 / (z[i] - z[j])
203  for j in range(i+1, deg):
204  sigma += 1 / (z[i] - z[j])
205  sigma_list[i] = sigma
206
207  for i in range(deg):
208  if not value_list[i]:
209  continue
210  ratio = value_list[i] / df(z[i])
211  z[i] -= ratio / (1 - ratio*sigma_list[i])
212  value_list[i] = f(z[i])
213
214  return z
215
216 def _initialize(g, newtoninitial=None):
217  """
218  create initial values of equation given as a list g.
219  """
220  g = [float(c) for c in g[:]]
221  q = [-abs(c) for c in g[:-1]]
222  q.append(abs(g[-1]))
223  if newtoninitial is None:
224  r = Newton(q, _upper_bound_of_roots(q))
225  else:
226  r = Newton(q, newtoninitial)
227
228  deg = len(g) - 1
229  center = -g[-2]/(deg*g[-1])
231  angular_step = cmath.exp(Complex(0, 1) * about_two_pi / deg)
232  angular_move = r
233  z = []
234  for i in range(deg):
235  z.append(center + angular_move)
236  angular_move *= angular_step
237
238  return z
239
241  """
242  Return an upper bound of roots.
243  """
244  weight = len(filter(None, g))
245  assert g[-1]
246  return max([pow(weight*abs(c/g[-1]), 1/len(g)) for c in g])
247
248 def root_Fp(g, p, flag=True):
249  """
250  Return a root over F_p of nonzero polynomial g.
251  p must be prime.
252  If flag = False, return a root randomly
253  """
254  if isinstance(g, list):
255  if not isinstance(g[0], tuple):
256  g = zip(range(len(g)), g)
258  g = uniutil.FinitePrimeFieldPolynomial(g, Fp)
259  h = uniutil.FinitePrimeFieldPolynomial({1:-1, p:1}, Fp)
260  g = g.gcd(h)
261  deg_g = g.degree()
262  if g[0] == 0:
263  deg_g = deg_g - 1
264  g = g.shift_degree_to(deg_g)
265  while True:
266  if deg_g == 0:
267  return None
268  if deg_g == 1:
269  return (-g[0]/g[1]).toInteger()
270  elif deg_g == 2:
271  d = g[1]*g[1] - 4*g[0]
272  e = arith1.modsqrt(d.toInteger(), p)
273  return ((-g[1]-e)/(2*g[2])).toInteger()
274  deg_h = 0
275  x = uniutil.FinitePrimeFieldPolynomial({0:-1, (p-1)>>1:1}, Fp)
276  if flag:
277  a = 0
278  while deg_h == 0 or deg_h == deg_g:
279  b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
280  v = g(b)
281  h = x.gcd(v)
282  a = a + 1
283  deg_h = h.degree()
284  b = uniutil.FinitePrimeFieldPolynomial({0:a-1, 1:1}, Fp)
285  else:
286  while deg_h == 0 or deg_h == deg_g:
287  a = bigrandom.randrange(p)
288  b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
289  v = g(b)
290  h = x.gcd(v)
291  deg_h = h.degree()
292  b = uniutil.FinitePrimeFieldPolynomial({0:a, 1:1}, Fp)
293  g = h(b)
294  deg_g = deg_h
295
296 def allroots_Fp(g, p):
297  """
298  Return roots over F_p of nonzero polynomial g.
299  p must be prime.
300  """
301  if isinstance(g, list):
302  if not isinstance(g[0], tuple):
303  g = zip(range(len(g)), g)
305  g = uniutil.FinitePrimeFieldPolynomial(g, Fp)
306  h = uniutil.FinitePrimeFieldPolynomial({1:-1, p:1}, Fp)
307  g = g.gcd(h)
308  deg_g = g.degree()
309  roots = []
310  if g[0] == 0:
311  roots.append(0)
312  deg_g = deg_g - 1
313  g = g.shift_degree_to(deg_g)
314  return roots + roots_loop(g, deg_g, p, Fp)
315
316 def roots_loop(g, deg_g, p, Fp):
317  if deg_g == 0:
318  return []
319  if deg_g == 1:
320  return [(-g[0]/g[1]).toInteger()]
321  elif deg_g == 2:
322  d = g[1]*g[1]-4*g[0]
323  e = arith1.modsqrt(d.toInteger(), p)
324  g1 = -g[1]
325  g2 = 2*g[2]
326  return [((g1 - e) / g2).toInteger(), ((g1 + e) / g2).toInteger()]
327  deg_h = 0
328  x = uniutil.FinitePrimeFieldPolynomial({0:-1, (p-1)>>1:1}, Fp)
329  a = 0
330  while deg_h == 0 or deg_h == deg_g:
331  b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
332  v = g(b)
333  h = x.gcd(v)
334  a = a + 1
335  deg_h = h.degree()
336  b = uniutil.FinitePrimeFieldPolynomial({0:a-1, 1:1}, Fp)
337  s = h(b)
338  deg_s = deg_h
339  t = g.exact_division(s)
340  deg_t = t.degree()
341  return roots_loop(s, deg_s, p, Fp) + roots_loop(t, deg_t, p, Fp)
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.equation.allroots_Fp
def allroots_Fp(g, p)
Definition: equation.py:296
nzmath.equation.root_Fp
def root_Fp(g, p, flag=True)
Definition: equation.py:248
nzmath.equation.e1_ZnZ
def e1_ZnZ(x, n)
Definition: equation.py:42
nzmath.equation.SimMethod
def SimMethod(f, NewtonInitial=1, repeat=250)
Definition: equation.py:161
nzmath.equation._SimMethod
def _SimMethod(g, initials=None, newtoninitial=None, repeat=250)
Definition: equation.py:172
nzmath.equation._upper_bound_of_roots
def _upper_bound_of_roots(g)
Definition: equation.py:240
nzmath.equation.e1
def e1(x)
Definition: equation.py:33
nzmath.equation.roots_loop
def roots_loop(g, deg_g, p, Fp)
Definition: equation.py:316
nzmath.equation.e3
def e3(x)
Definition: equation.py:88
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.uniutil
Definition: uniutil.py:1
nzmath.arith1.inverse
def inverse(x, n)
Definition: arith1.py:183
nzmath.plugins
Definition: plugins.py:1
nzmath.equation.Newton
def Newton(f, initial=1, repeat=250)
Definition: equation.py:136
nzmath.bigrange
Definition: bigrange.py:1
nzmath.equation.e2
def e2(x)
Definition: equation.py:52
nzmath.equation._initialize
def _initialize(g, newtoninitial=None)
Definition: equation.py:216
nzmath.imaginary
Definition: imaginary.py:1
nzmath.equation.e3_Fp
def e3_Fp(x, p)
Definition: equation.py:114
nzmath.equation.e2_Fp
def e2_Fp(x, p)
Definition: equation.py:64
nzmath.finitefield
Definition: finitefield.py:1