"Fossies" - the Fresh Open Source Software Archive

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

    1 """
    2 binary-like gcd algorithms for rational, gauss, and eisenstein integer
    3 """
    4 
    5 import math
    6 
    7 
    8 def bit_num(a):
    9     """
   10     return the number of bits
   11     """
   12     if a == 0:
   13         return 1
   14     else:
   15         return int(math.log(a, 2)) + 1
   16 
   17 
   18 def binarygcd(a, b):
   19     """
   20     Return the greatest common divisor of 2 integers a and b
   21     by binary gcd algorithm.
   22     """
   23     if a < b:
   24         a, b = b, a
   25     if b == 0:
   26         return a
   27     a, b = b, a % b
   28     if b == 0:
   29         return a
   30     k = 0
   31     while not a&1 and not b&1:
   32         k += 1
   33         a >>= 1
   34         b >>= 1
   35     while not a&1:
   36         a >>= 1
   37     while not b&1:
   38         b >>= 1
   39     while a:
   40         while not a & 1:
   41             a >>= 1
   42         if abs(a) < abs(b):
   43             a, b = b, a
   44         a = (a - b) >> 1
   45     return b << k
   46 
   47 
   48 def arygcd_i(a1, a2, b1, b2):
   49     """
   50     Return the greatest common divisor of 2 gauss-integers a1+a2*i and b1+b2*i
   51     by (1+i)-ary gcd algorithm.
   52     """
   53     if a1 == 0 and a2 == 0:
   54        return b1, b2
   55     elif b1 == 0 and b2 == 0:
   56        return a1, a2
   57     ap, bp = 0, 0
   58     while (a1-a2)&1 == 0:
   59         a1, a2 = (a1+a2) >> 1, (a2-a1) >> 1
   60         ap += 1
   61     while (b1-b2)&1 == 0:
   62         b1, b2 = (b1+b2) >> 1, (b2-b1) >> 1
   63         bp += 1
   64     k = min(ap, bp)
   65     while a1 != 0 or a2 != 0:
   66         while (a1-a2) & 1 == 0:
   67             a1, a2 = (a1+a2) >> 1, (a2-a1) >> 1
   68         norm_a, norm_b = _ap_norm_i(a1, a2, b1, b2)
   69         if norm_a < norm_b:
   70             a1, a2, b1, b2 = b1, b2, a1, a2
   71         a1, a2 = _FormAdj_i(a1, a2)
   72         b1, b2 = _FormAdj_i(b1, b2)
   73         a1, a2 = a1-b1, a2-b2
   74     if 0 in (ap, bp):
   75         return b1, b2
   76     else:
   77         s, t = _n_pow_i(1, 1, k)
   78         return (b1*s)-(b2*t), (b1*t)+(b2*s)
   79 
   80 def _ap_norm_i(a, b, c, d):
   81     """
   82     Return approximately norm of 2 gaussInteger a+b*i and c+d*i  
   83     """
   84     a, b, c, d = abs(a), abs(b), abs(c), abs(d)
   85     max_dig = max(bit_num(a), bit_num(b), bit_num(c), bit_num(d))
   86     if max_dig > 6:
   87         sft_num = max_dig - 6
   88         a, b, c, d = a >> sft_num, b >> sft_num, c >> sft_num, d >> sft_num
   89         return a*a + b*b, c*c + d*d
   90     else:
   91         return a*a + b*b, c*c + d*d
   92 
   93 def _n_pow_i(a, b, n):
   94     """
   95     return (1+i)**k 
   96     """
   97     x = a
   98     y = b
   99 
  100     for i in range(1, n):
  101         x1 = (x*a) - (y*b)
  102         y1 = (y*a) + (x*b)
  103         x = x1
  104         y = y1
  105     return x, y
  106 
  107 def _FormAdj_i(a, b):
  108     """
  109     transform gaussInteger a+b*i ->  form 1+2(1+i)*(x+y*i)
  110     """
  111     if a & 1 == 0:
  112         a, b = -b, a
  113     if (b - a + 1) & 3 == 0:
  114         return a, b
  115     else:
  116         return -a, -b
  117 
  118 
  119 def arygcd_w(a1, a2, b1, b2):
  120     """
  121     Return the greatest common divisor of 2 eisensteinIntegers a1+a2*w and b1+b2*w
  122     by (1-w)-ary gcd algorithm.
  123     """
  124     if a1 == 0 and a2 == 0:
  125        return b1, b2
  126     elif b1 == 0 and b2 == 0:
  127        return a1, a2
  128     ap, bp = 0, 0
  129     while (a1 + a2) % 3 == 0:
  130         a1, a2 = (a1 + a1 - a2) / 3, (a1 + a2) / 3
  131         ap += 1
  132     while (b1 + b2) % 3 == 0:
  133         b1, b2 = (b1 + b1 - b2) / 3, (b1 + b2) / 3
  134         bp += 1
  135     k = min(ap, bp)
  136     while a1 != 0 or a2 != 0:
  137         while (a1 + a2) % 3 == 0:
  138             a1, a2 = (a1 + a1 - a2) / 3, (a1 + a2) / 3
  139         nrm_a, nrm_b = _ap_norm_w(a1, a2, b1, b2)
  140         if nrm_a < nrm_b:
  141             a1, a2, b1, b2 = b1, b2, a1, a2
  142         a1, a2 = _FormAdj_w(a1, a2)
  143         b1, b2 = _FormAdj_w(b1, b2)
  144         a1, a2 = a1 - b1, a2 - b2
  145     k1, k2 = _n_pow_w(k)
  146     return k1*b1 - k2*b2, k1*b2 + k2*b1 - k2*b2
  147 
  148 def _ap_norm_w(a, b, c, d):
  149     """
  150     Return approximately norm of 2 eisensteinInteger a+b*w and c+d*w  
  151     """
  152     a, b, c, d = abs(a), abs(b), abs(c), abs(d)
  153     max_dig = max(bit_num(a), bit_num(b), bit_num(c), bit_num(d))
  154     if max_dig > 6:
  155         sft_num = max_dig - 6
  156         a, b, c, d = a >> sft_num, b >> sft_num, c >> sft_num, d >> sft_num
  157         subst1, subst2 = (a - b) >> sft_num, (c - d) >> sft_num
  158         return a*a + b*b + subst1*subst1, c*c + d*d + subst2*subst2
  159     else:
  160         return a*a + b*b + (a - b)*(a - b), c*c + d*d + (c -d)*(c -d)
  161 
  162 def _n_pow_w(n):
  163     """
  164     return (1-w)**k
  165     """
  166     x, y = divmod(n, 2)
  167     k1 = 3**x
  168     if y == 1:
  169         k2 = -k1
  170     else:
  171         k2 = 0
  172     return k1, k2
  173 
  174 def _FormAdj_w(a1, a2):
  175     """
  176     transform eisensteinInteger a1+a2*w ->  form 1+3*(x+y*w)
  177     """
  178     if a1 % 3 == 0:
  179         if a2 % 3 == -1 or a2 % 3 == 2:
  180             return a1 - a2, a1
  181         else:
  182             return a2 - a1, -a1
  183     elif a1 % 3 == 1:
  184         if a2 % 3 == 1:
  185             return a2, a2 -a1
  186         else:
  187             return a1, a2
  188     else:
  189         if a2 % 3 == -1 or a2 % 3 == 2:
  190             return -a2, a1 -a2
  191         else:
  192             return -a1, -a2