## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/poly/array.py" (19 Nov 2012, 14884 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 "array.py" see the Fossies "Dox" file reference documentation.

```    1 import nzmath.arith1 as arith1
2
3 def check_zero_poly(coefficients):
4     """
5     This function checks whether all elements of coefficients equal
6     zero or not.  If all elements of coefficients equal zero, this
7     function returns True.  Else this function returns False.
8     """
9     for i in coefficients:
10         if i != 0:
11             return False
12     return True
13
14 def arrange_coefficients(coefficients):
15     """
16     This function arranges coefficient.
17     For example, [1,2,0,3,0] => [1,2,0,3].
18     """
19     if check_zero_poly(coefficients):
20         return [0]
21     while coefficients[-1] == 0:
22         coefficients = coefficients[:-1]
23     return coefficients
24
25 class ArrayPoly():
26     """
27     Polynomial with integer number coefficients.
28     Coefficients has to be a initializer for list.
29     """
30
31     def __init__(self, coefficients = [0]):
32         """
33         Initialize the polynomial.
34         coefficients:initializer for polynomial coefficients
35         """
36         self.coefficients = arrange_coefficients(coefficients)
37         self.degree = len(self.coefficients) - 1
38
39     def coefficients_to_dict(self):
40         """
41         Return coefficients as dict.
42         """
43         if self.coefficients == [0]:
44             return {0:0}
45         dict_coefficients = {}
46         for i in range(self.degree + 1):
47             if self.coefficients[i] != 0:
48                 dict_coefficients.update({i:self.coefficients[i]})
49         return dict_coefficients
50
51     def __repr__(self):
52         poly_repr = self.coefficients_to_dict()
53         return "ArrayPoly : %s" % poly_repr
54
55     def __str__(self):
56         poly_str = self.coefficients_to_dict()
57         return "polynomial : %s" % poly_str
58
59     def __add__(self, other):
60         """
61         self + other
62         """
63         add = []
64         deg1 = self.degree
65         deg2 = other.degree
66         if deg1 >= deg2:
67             long_coefficients = self.coefficients[:]
68             short_coefficients = other.coefficients[:]
69             deg = deg2 + 1
70         else:
71             long_coefficients = other.coefficients[:]
72             short_coefficients = self.coefficients[:]
73             deg = deg1 + 1
74         for i in range(deg):
75             add.append(long_coefficients[i] + short_coefficients[i])
76         add += long_coefficients[deg:]
77         return ArrayPoly(add)
78
79     def __sub__(self, other):
80         """
81         self - other
82         """
83         sub = []
84         deg1 = self.degree
85         deg2 = other.degree
86         if deg1 >= deg2:
87             long_coefficients = self.coefficients[:]
88             short_coefficients = other.coefficients[:]
89             deg = deg2 + 1
90         else:
91             long_coefficients = other.coefficients[:]
92             short_coefficients = self.coefficients[:]
93             deg = deg1 + 1
94         for i in range(deg):
95             sub.append(long_coefficients[i] - short_coefficients[i])
96         sub += long_coefficients[deg:]
97         return ArrayPoly(sub)
98
99     def scalar_mul(self, scalar):
100         """
101         Return the result of scalar multiplicaton.
102         """
103         scalar_mul = [i * scalar for i in self.coefficients]
104         return ArrayPoly(scalar_mul)
105
106     def upshift_degree(self, slide):
107         """
108         Return the polynomial obtained by shifting upward all terms
109         with degrees of 'slide'.
110         """
111         if slide == 0:
112             return ArrayPoly(self.coefficients[:])
113         up_degree = [0] * slide
114         return ArrayPoly(up_degree + self.coefficients[:])
115
116     def downshift_degree(self, slide):
117         """
118         Return the polynomial obtained by shifting downward all terms
119         with degrees of 'slide'.
120         """
121         if slide == 0:
122             return ArrayPoly(self.coefficients[:])
123         elif slide > self.degree:
124             return ArrayPoly()
125         down_degree= self.coefficients[slide:]
126         return ArrayPoly(down_degree)
127
128     def __eq__(self, other):
129         """
130         self == other
131         """
132         self_dict_coefficients = self.coefficients_to_dict()
133         other_dict_coefficients = other.coefficients_to_dict()
134         return self_dict_coefficients == other_dict_coefficients
135
136     def __ne__(self, other):
137         """
138         self != other
139         """
140         return not self.__eq__(other)
141
142     def __mul__(self, other):
143         """
144         self * other
145         """
146         total = [0] * (self.degree + other.degree + 1)
147         for i in range(self.degree + 1):
148             for j in range(other.degree + 1):
149                 total[i + j] += self.coefficients[i] * other.coefficients[j]
150         return ArrayPoly(total)
151
152     def power(self):
153         """
154         self * self
155         """
156         total = [0] * (self.degree + self.degree + 1)
157         for i in range(self.degree + 1):
158             total[i + i] += self.coefficients[i] * self.coefficients[i]
159             for j in range(i + 1, self.degree + 1):
160                 total[i + j] += ((self.coefficients[i] * self.coefficients[j]) << 1)
161         return ArrayPoly(total)
162
163     def split_at(self, split_point):
164         """
165         Return tuple of two polynomials, which are splitted at the
166         given degree.  The term of the given degree, if exists,
167         belongs to the lower degree polynomial.
168         """
169         split = self.coefficients[:]
170         if split_point == 0:
171             return (ArrayPoly(), ArrayPoly(split))
172         elif split_point >= self.degree:
173             return (ArrayPoly(split), ArrayPoly())
174         split1 = split[:split_point + 1]
175         split2 = split[:]
176         for i in range(split_point + 1):
177             split2[i] = 0
178         return (ArrayPoly(split1), ArrayPoly(split2))
179
180     def FFT_mul(self, other):
181         """
182         self * other by Fast Fourier Transform.
183         """
184         coefficients1 = [abs(i) for i in self.coefficients]
185         coefficients2 = [abs(i) for i in other.coefficients]
186         bound1 = max(coefficients1)
187         bound2 = max(coefficients2)
188         bound = bound1 * bound2 * (max(self.degree, other.degree) + 1)
189         bound = ceillog(bound, 2)
190         bound = 1 << bound
191         if bound < (self.degree + other.degree + 1):
192             bound = self.degree + other.degree + 1
193             bound = ceillog(bound, 2)
194             bound = 1 << bound
195         fft1 = ArrayPoly(self.coefficients[:])
196         fft2 = ArrayPoly(other.coefficients[:])
197         fft_mul1 = FFT(fft1, bound)
198         fft_mul2 = FFT(fft2, bound)
199         fft_mul = []
200         for i in range(bound):
201             fft_mul.append(fft_mul1[i] * fft_mul2[i])
202         #print fft_mul
203         total = reverse_FFT(fft_mul, bound)
204         #print total
205         return ArrayPoly(total)
206
207
208 class ArrayPolyMod(ArrayPoly):
209     """
210     Polynomial with modulo n coefficients, n is a nutural number.
211     Coefficients has to be a initializer for list.
212     """
213
214     def __init__(self, coefficients, mod):
215         """
216         Initialize the polynomial.
217         coefficients:initializer for polynomial coefficients
218         mod:initializer for coefficients modulo mod
219         """
220         if mod <= 0:
221             raise ValueError("Please input a positive integer in 'mod'.")
222         mod_coefficients = [i % mod for i in coefficients]
223         ArrayPoly.__init__(self, mod_coefficients)
224         self.mod = mod
225
226     def __repr__(self):
227         poly_repr = self.coefficients_to_dict()
228         return "polynomial_mod(%d):%s" % (self.mod, poly_repr)
229
230     def __str__(self):
231         poly_str = self.coefficients_to_dict()
232         return "polynomial_mod(%d):%s" % (self.mod, poly_str)
233
234     def __add__(self, other):
235         """
236         self + other
237         """
238         if self.mod != other.mod:
239             raise ValueError("mod mismatch")
240         add = ArrayPoly.__add__(self, other)
241         return ArrayPolyMod(add.coefficients, self.mod)
242
243     def __sub__(self, other):
244         """
245         self - other
246         """
247         if self.mod != other.mod:
248             raise ValueError("mod mismatch")
249         sub = ArrayPoly.__sub__(self, other)
250         return ArrayPolyMod(sub.coefficients, self.mod)
251
252     def scalar_mul(self, scalar):
253         """
254         Return the result of scalar multiplicaton.
255         """
256         scalar_mul = ArrayPoly.scalar_mul(self, scalar)
257         return ArrayPolyMod(scalar_mul.coefficients, self.mod)
258
259     def upshift_degree(self, slide):
260         """
261         Return the polynomial obtained by shifting upward all terms
262         with degrees of 'slide'.
263         """
264         up_degree = ArrayPoly.upshift_degree(self, slide)
265         return ArrayPolyMod(up_degree.coefficients, self.mod)
266
267     def downshift_degree(self, slide):
268         """
269         Return the polynomial obtained by shifting downward all terms
270         with degrees of 'slide'.
271         """
272         down_degree = ArrayPoly.downshift_degree(self, slide)
273         return ArrayPolyMod(down_degree.coefficients, self.mod)
274
275     def __eq__(self, other):
276         """
277         self == other
278         """
279         if self.mod != other.mod:
280             raise ValueError("mod mismatch")
281         eq = ArrayPoly.__eq__(self, other)
282         return eq
283
284     def __ne__(self, other):
285         """
286         self != other
287         """
288         if self.mod != other.mod:
289             raise ValueError("mod mismatch")
290         ne = ArrayPoly.__ne__(self, other)
291         return ne
292
293     def __mul__(self, other):
294         """
295         self * other
296         """
297         if self.mod != other.mod:
298             raise ValueError("mod mismatch")
299         total = [0] * (self.degree + other.degree + 1)
300         for i in range(self.degree + 1):
301             for j in range(other.degree + 1):
302                 total[i + j] = (total[i + j] + self.coefficients[i] * other.coefficients[j]) % self.mod
303         return ArrayPolyMod(total, self.mod)
304
305     def power(self):
306         """
307         self * self
308         """
309         total = [0] * (self.degree + self.degree + 1)
310         for i in range(self.degree + 1):
311             total[i + i] = (total[i + i] + self.coefficients[i] * self.coefficients[i]) % self.mod
312             for j in range(i + 1, self.degree + 1):
313                 total[i + j] = (total[i + j] + ((self.coefficients[i] * self.coefficients[j]) << 1)) % self.mod
314         return ArrayPolyMod(total, self.mod)
315
316     def split_at(self, split_point):
317         """
318         Return tuple of two polynomials, which are splitted at the
319         given degree.  The term of the given degree, if exists,
320         belongs to the lower degree polynomial.
321         """
322         if split_point == 0:
323             return (ArrayPolyMod([0], self.mod), ArrayPolyMod(self.coefficients, self.mod))
324         elif split_point >= self.degree:
325             return (ArrayPolyMod(self.coefficients, self.mod), ArrayPolyMod([0], self.mod))
326         split = ArrayPoly.split_at(self, split_point)
327         split1 = split[0].coefficients
328         split2 = split[1].coefficients
329         return (ArrayPolyMod(split1, self.mod), ArrayPolyMod(split2, self.mod))
330
331     def FFT_mul(self, other):
332         """
333         self * other by Fast Fourier Transform.
334         """
335         if self.mod != other.mod:
336             raise ValueError("mod mismatch")
337         bound1 = max(self.coefficients)
338         bound2 = max(other.coefficients)
339         bound = bound1 * bound2 * (max(self.degree, other.degree) + 1)
340         bound = ceillog(bound, 2)
341         bound = 1 << bound
342         if bound < (self.degree + other.degree + 1):
343             bound = self.degree + other.degree + 1
344             bound = ceillog(bound, 2)
345             bound = 1 << bound
346         fft1 = ArrayPolyMod(self.coefficients[:], self.mod)
347         fft2 = ArrayPolyMod(other.coefficients[:], other.mod)
348         fft_mul1 = FFT(fft1, bound)
349         fft_mul2 = FFT(fft2, bound)
350         fft_mul = []
351         for i in range(bound):
352             fft_mul.append(fft_mul1[i] * fft_mul2[i])
353         total = reverse_FFT(fft_mul, bound)
354         return ArrayPolyMod(total, self.mod)
355
356 #
357 #Some functions for FFT(Fast Fourier Transform).
358 #
359 def min_abs_mod(a, b):
360     """
361     This function returns absolute minimum modulo of a over b.
362     """
363     if a >= 0:
364         return a - b * ((a + a + b) // (b + b))
365     a = -a
366     return -(a - b * ((a + a + b) // (b + b)))
367
368 def bit_reverse(n, bound):
369     """
370     This function returns the result reversed bit of n.
371     bound:number of significant figures of bit.
372     """
373     bit = []
374     while n > 0:
375         bit.append(n & 1)
376         n = n >> 1
377     bound_bitlen = ceillog(bound, 2)
378     if bound_bitlen > len(bit):
379         bit.extend([0] * (bound_bitlen - len(bit)))
380     bit.reverse()
381     total = 0
382     count = 0
383     for i in bit:
384         if i != 0:
385             total += 1 << count
386         count += 1
387     return total
388
389 def ceillog(n, base=2):
390     """
391     Return ceiling of log(n, 2)
392     """
393     floor = arith1.log(n, base)
394     if base ** floor == n:
395         return floor
396     else:
397         return floor + 1
398
399 def perfect_shuffle(List):
400     """
401     This function returns list arranged by divide-and-conquer method.
402     """
403     length = len(List)
404     shuffle = [0] * length
405     for i in range(length):
406         rebit = bit_reverse(i, length)
407         shuffle[rebit] = List[i]
408     return shuffle
409
410 def FFT(f, bound):
411     """
412     Fast Fourier Transform.
413     This function returns the result of valuations of f by fast fourier transform
414     against number of bound different values.
415     """
416     count = 1 << (bound >> 1)
417     mod = count + 1
418     if isinstance(f, ArrayPoly):
419         coefficients = f.coefficients[:]
420     else:
421         coefficients = f[:]
422     coefficients.extend([0] * (bound - len(coefficients)))
423     shuf_coefficients = perfect_shuffle(coefficients)
424     shuf_coefficients = [min_abs_mod(i, mod) for i in shuf_coefficients]
425     bound_log = arith1.log(bound, 2)
426     for i in range(1, bound_log + 1):
427         m = 1 << i
428         wm = count
429         wm = wm % mod
430         w = 1
431         m1 = m >> 1
432         for j in range(m1):
433             for k in range(j, bound, m):
434                 m2 = k + m1
435                 plus = shuf_coefficients[k]
436                 minus = w * shuf_coefficients[m2]
437                 shuf_coefficients[k] = (plus + minus) % mod
438                 shuf_coefficients[m2] = (plus - minus) % mod
439             w = w * wm
440         shuf_coefficietns = [min_abs_mod(i, mod) for i in shuf_coefficients]
441         count = arith1.floorsqrt(count)
442     return shuf_coefficients
443
444 def reverse_FFT(values, bound):
445     """
446     Reverse Fast Fourier Tfransform.
447     """
448     mod = (1 << (bound >> 1)) + 1
449     shuf_values = values[:]
450     reverse_FFT_coeffficients = FFT(shuf_values, bound)
451     inverse = arith1.inverse(bound, mod)
452     reverse_FFT_coeffficients = [min_abs_mod(inverse * i, mod) for i in reverse_FFT_coeffficients]
453     reverse_FFT_coeffficients1 = reverse_FFT_coeffficients[:1]
454     reverse_FFT_coeffficients2 = reverse_FFT_coeffficients[1:]
455     reverse_FFT_coeffficients2.reverse()
456     reverse_FFT_coeffficients_total = reverse_FFT_coeffficients1 + reverse_FFT_coeffficients2
457     return reverse_FFT_coeffficients_total
```