"Fossies" - the Fresh Open Source Software Archive

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

    1 import random
    2 import nzmath.arith1 as arith1
    3 import nzmath.arygcd as arygcd
    4 
    5 def c_root_p(a, p):
    6     """
    7     Return the cubic root modulo p.
    8     (i.e. a^3 = x (mod p))
    9     """
   10     if (a % p) == 0:
   11         return [0]
   12     if p == 2 or p == 3 :
   13         return [a % p]
   14     if (p % 3) == 2:
   15         return [pow(a, (((2 * p) - 1) // 3), p)]
   16     p_div_3, p_mod_3 = divmod((p - 1), 3)
   17     # Compute e,q
   18     e = 0
   19     temp = p_div_3
   20     tempmod = p_mod_3
   21 
   22     while tempmod == 0:
   23         e += 1
   24         temp, tempmod = divmod(temp, 3)
   25     q = (p - 1) // (3 ** e)
   26     search_range = (p - 1) >> 1
   27     h = 2
   28     while pow(h, p_div_3, p) == 1:
   29         h = random.randrange(2, search_range)
   30     sym = pow(h, p_div_3, p)
   31     g = pow(h, q, p)
   32     # Initialize
   33     y = g
   34     r = e
   35     if q % 3 == 2:
   36         x = pow(a, (q - 2) // 3, p)
   37     else:
   38         x = pow(a, (((2 * q) - 2) // 3), p)
   39 
   40     b = (pow(a, 2, p) * pow(x, 3, p)) % p
   41     x = (a * x) % p
   42     while (b % p) != 1:
   43         # Find exponent
   44         b_pow = pow(b, 3, p)
   45         m = 1
   46         while b_pow != 1:
   47             b_pow = (b_pow ** 3) % p
   48             m += 1
   49         if m == r:
   50             raise ValueError("there is no cubic root mod p")
   51         # Reduce exponent
   52         if sym == pow(b, pow(3, m - 1, p), p):
   53             t = pow(y, 2, p)
   54             sym = pow(sym, 2, p)
   55         else:
   56             t = y
   57         t = pow(t, pow(3, r - m - 1, p), p)
   58         y = pow(t, 3, p)
   59         r = m
   60         x = (x * t) % p
   61         b = (b * y) % p
   62     return [ x, (x * sym) % p, (x * pow(sym, 2, p)) % p ]
   63 
   64 def c_residue(a, p):
   65     """
   66     Check whether rational integer a is cubic residue mod prime p.
   67     If a % p == 0, then return 0,
   68     elif cubic residue, then return 1,
   69     elif cubic non-residue, then return -1.
   70     """
   71     if p == 3:
   72         if a % p == 0:
   73             return 0
   74         else:
   75             return 1
   76     elif p % 3 == 1:
   77         b1, b2 = decomposite_p(p)
   78         return c_symbol(a, 0, b1, b2)
   79     else:
   80         return c_symbol(a, 0, p, 0)
   81 
   82 def c_symbol(a1, a2, b1, b2):
   83     """
   84     Return the (Jacobi) cubic residue symbol of 2 Eisenstein Integers ((a1+a2*w)/(b1+b2*w)).
   85     assume that b1+b2*w is not divisable 1-w.
   86     """
   87     r1, r2, j1 = _divides(a1, a2)
   88     r1, r2, i1 = _FormAdj_w(r1, r2)
   89     d1, d2, i2 = _FormAdj_w(b1, b2)
   90     m = (d1 - 1) / 3
   91     n = d2 / 3
   92     t = ((m * j1) + (-(m + n) * i1)) % 3
   93     a1, a2 = r1, r2
   94     b1, b2 = d1, d2
   95     nrm_a, nrm_b = arygcd._ap_norm_w(a1, a2, b1, b2)
   96     if nrm_a < nrm_b:
   97         tmpre, tmpim = a1, a2
   98         a1, a2 = b1, b2
   99         b1, b2 = tmpre, tmpim
  100         m = (b1 - 1) / 3
  101         n = b2 / 3
  102     while not((a1 == b1) and (a2 == b2)):
  103         sub_re_ab, sub_im_ab = a1 - b1, a2 - b2
  104         tmp_sub_re_ab, tmp_sub_im_ab, j = _divides(sub_re_ab, sub_im_ab)
  105         r1, r2, i = _FormAdj_w(tmp_sub_re_ab, tmp_sub_im_ab)
  106         t = (t + (m * j) + (-(m + n) * i)) % 3
  107         a1, a2 = r1, r2
  108         nrm_a, nrm_b = arygcd._ap_norm_w(a1, a2, b1, b2)
  109         if nrm_a < nrm_b:
  110             tmpre, tmpim = a1, a2
  111             a1, a2 = b1, b2
  112             b1, b2 = tmpre, tmpim
  113             m = (b1 - 1) / 3
  114             n = b2 / 3
  115     if not ((a1 == 1) and (a2 == 0)) :
  116         return 0
  117     elif t  == 0:
  118         return 1
  119     else:
  120         return -1
  121 
  122 def _FormAdj_w(a1, a2):
  123     """
  124     Transform Eisenstein Integer a1+a2*w  ->  (-w)^i * (x+y*w)
  125     x+y*w is primary element
  126     assume that a1+a2*w is not divisible 1-w
  127     """
  128     if a1 % 3 == 0:
  129         if ((a2 % 3) == 2) or ((a2 % 3) == -1):
  130             return a1 - a2, a1, 1
  131         elif ((a2 % 3) == 1) or ((a2 % 3) == -2):
  132             return a2 - a1, -a1, 4
  133     elif ((a1 % 3) == 1) or ((a1 % 3) == -2):
  134         if ((a2 % 3) == 1) or ((a2 % 3) == -2):
  135             return a2, a2 - a1, 5
  136         elif a2 % 3 == 0:
  137             return a1, a2, 0
  138     else:
  139         if ((a2 % 3) == 2) or ((a2 % 3) == -1):
  140             return -a2, a1 - a2, 2
  141         elif a2 % 3 == 0:
  142             return -a1, -a2, 3
  143 
  144 def _divides(a1, a2):
  145     """
  146     divide 1-w
  147     (i.e. a1+a2*w -> (1-w)^k * (x+y*w))
  148     """
  149     if (a1 == 0) and (a2 == 0):
  150         return a1, a2, 0
  151     j = 0
  152     while (a1 + a2) % 3 == 0:
  153         tmpa1 = a1
  154         a1 = ((a1 + a1) - a2) / 3
  155         a2 = (tmpa1 + a2) / 3
  156         j += 1
  157     return a1, a2, j
  158 
  159 def decomposite_p(p):
  160     """
  161     Decomposite p = 1(mod 3) into its primary prime factors in z[w]
  162     """
  163     s, t = cornacchia(3, p)
  164     if t % 3 == 0:
  165         a = s + t
  166         b = 2 * t
  167     if (t % 3) != 0:
  168         a = 2 * t
  169         b = s + t
  170     return a, b
  171 
  172 def cornacchia(d, p):
  173     """
  174     Return the solution of x^2 + d * y^2 = p .
  175     p be a prime and d be an integer such that 0 < d < p.
  176     """
  177     if (d <= 0) or (d >= p):
  178         raise ValueError("invalid input")
  179     k = arith1.legendre(-d, p)
  180     if k == -1:
  181         raise ValueError("no solution")
  182     x0 = arith1.modsqrt(-d, p)
  183     if x0 < (p / 2):
  184         x0 = p - x0
  185     a = p
  186     b = x0
  187     l = arith1.floorsqrt(p)
  188     while b > l:
  189         a, b = b, a % b
  190     c, r = divmod(p - b * b, d)
  191     if r:
  192         raise ValueError("no solution")
  193     t = arith1.issquare(c)
  194     if t == 0:
  195         raise ValueError("no solution")
  196     else:
  197         return (b, t)