"Fossies" - the Fresh Open Source Software Archive

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

    1 """
    2 Abstract (higher order) functions
    3 """
    4 
    5 import functools
    6 
    7 def digital_method(coefficients, val, add, mul, act, power, zero, one):
    8     """
    9     Evaluate a univariate polynomial at val.
   10 
   11     The polynomial is given as an iterator 'coefficients' which yields
   12     (degree, coefficient) pair in descending order of degree.
   13 
   14     If the polynomial is of R-coefficients, then val should be in an
   15     R-algebra D.
   16 
   17     All operations 'add', 'mul', 'act', 'power', 'zero', 'one' should
   18     be explicitly given, where: 
   19       'add' means addition (D x D -> D),
   20       'mul' multiplication (D x D -> D),
   21       'act' action of R (R x D -> D),
   22       'power' powering (D x Z -> D) and
   23       'zero', 'one' constants in D.
   24     """
   25     result = zero
   26     val_pow = {0:one, 1:val}
   27     d = None
   28     for d_next, c_next in coefficients:
   29         if d:
   30             diff = d - d_next
   31             if diff not in val_pow:
   32                 val_pow[diff] = power(val, diff)
   33             result = add(mul(result, val_pow[diff]), act(c_next, one))
   34         else:
   35             result = act(c_next, one)
   36         d = d_next
   37     if d:
   38         return mul(result, val_pow.get(d, power(val, d)))
   39     else:
   40         return result
   41 
   42 def digital_method_func(add, mul, act, power, zero, one):
   43     """
   44     Return a function which evaluates polynomial (binary function)
   45     constructed from given functions.
   46     """
   47     result = functools.partial(digital_method, add=add, mul=mul,
   48                                act=act, power=power, zero=zero, one=one)
   49     result.___doc__ = """
   50     Evaluate a univariate polynomial at val.
   51 
   52     The polynomial is given as an iterator 'coefficients' which yields
   53     (degree, coefficient) pair in descending order of degree.
   54 
   55     If the polynomial is of R-coefficients, then val should be in an
   56     R-algebra D.
   57     """
   58     return result
   59 
   60 def rl_binary_powering(element, index, mul, square=None, one=None):
   61     """
   62      powering by using right-left binary method
   63      (assume that index >= 0)
   64     """
   65     square = _set_square(mul, square)
   66     sol = one
   67     mul_part = element
   68     while True:
   69         if index & 1:
   70             try:
   71                 sol = mul(sol, mul_part)
   72             except: # one is None!
   73                 sol = mul_part
   74         index >>= 1
   75         if not index:
   76             return sol
   77         else:
   78             mul_part = square(mul_part)
   79 
   80 def lr_binary_powering(element, index, mul, square=None, one=None):
   81     """
   82     powering by using left-right binary method.
   83     (assume that index >= 0)
   84     """
   85     import math
   86     square = _set_square(mul, square)
   87     if not index:
   88         return one
   89     spot = 1 << long(math.log(index, 2))
   90     sol = one
   91     while True:
   92         if spot & index:
   93             try:
   94                 sol = mul(sol, element)
   95             except: # one is None!
   96                 sol = element
   97         spot >>= 1
   98         if not spot:
   99             return sol
  100         sol = square(sol)
  101 
  102 def window_powering(element, index, mul, square=None, one=None):
  103     """
  104     powering by using small-window method
  105     window size is selected by average analystic optimization
  106     (assume that index >= 0)
  107     """
  108     import math
  109     square = _set_square(mul, square)
  110     if not index:
  111         return one
  112     log_n = long(math.log(index, 2))
  113     # Find the proper window size
  114     size = 2
  115     pow_size = 2
  116     while log_n > (size + 1) * (size + 2) * (pow_size - 1):
  117         pow_size <<= 1
  118         size += 1
  119     # Precomputation
  120     sqr = square(element)
  121     pre_table = [element]
  122     pow_size -= 1
  123     for i in range(pow_size):
  124         pre_table.append(mul(pre_table[-1], sqr))
  125 
  126     near_n = 1 << log_n
  127     spot = near_n
  128     while spot:
  129         if not(spot & index):
  130             sol = square(sol)
  131             spot >>= 1
  132         else:
  133             # Find the window
  134             f_spot, e_spot = spot, spot >> (size - 1)
  135             t_size = size
  136             if not(e_spot):
  137                 e_spot = 1
  138                 t_size = int(math.log(f_spot, 2)) + 1
  139             while True:
  140                 if e_spot & index:
  141                     spot = (e_spot >> 1)
  142                     window = 0
  143                     sqr = 1
  144                     while f_spot != e_spot: # Compute value of window
  145                         if index & e_spot:
  146                             window += sqr
  147                         sqr <<= 1
  148                         e_spot <<= 1
  149                     window += sqr
  150                     break
  151                 t_size -= 1
  152                 e_spot <<= 1
  153             # Compute a part of powering
  154             if f_spot == near_n:
  155                 sol = pre_table[(window - 1) >> 1]
  156             else:
  157                 for i in range(t_size):
  158                     sol = square(sol)
  159                 sol = mul(sol, pre_table[(window - 1) >> 1])
  160     return sol
  161 
  162 def _set_square(mul, square):
  163     """
  164     Set square if square is None
  165     (for the initialization of methods for powering)
  166     """
  167     if square == None:
  168         return lambda x : mul(x, x)
  169     else:
  170         return square
  171 
  172 def powering_func(mul, square=None, one=None, type=0):
  173     """
  174     Return a function which computes powering (binary function)
  175     constructed from given functions.
  176     """
  177     if type == 0:
  178         result = functools.partial(
  179             rl_binary_powering, mul=mul, one=one, square=square)
  180         result.___doc__ = """
  181         compute powering, that is, element ** index.
  182 
  183         This function uses right-left binary method.
  184         (assume that index >= 0)
  185         """
  186     elif type == 1:
  187         result = functools.partial(
  188             lr_binary_powering, mul=mul, one=one, square=square)
  189         result.___doc__ = """
  190         compute powering, that is, element ** index.
  191 
  192         This function uses left-right binary method.
  193         (assume that index >= 0)
  194         """
  195     elif type == 2:
  196         result = functools.partial(
  197             window_powering, mul=mul, one=one, square=square)
  198         result.___doc__ = """
  199         compute powering, that is, element ** index.
  200 
  201         This function uses window method.
  202         (assume that index >= 0)
  203         """
  204     return result