"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
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)
```