"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