"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/equation.py" (19 Nov 2012, 9601 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "equation.py" see the Fossies "Dox" file reference documentation.

```    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])
230     about_two_pi = 6
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
240 def _upper_bound_of_roots(g):
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)
257     Fp = finitefield.FinitePrimeField(p)
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)
304     Fp = finitefield.FinitePrimeField(p)
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)
```