UNPKG

2.56 kBtext/coffeescriptView Raw
1# Factor using the Pollard rho method
2
3
4
5n_factor_number = 0
6
7factor_number = ->
8 h = 0
9
10 save()
11
12 p1 = pop()
13
14 # 0 or 1?
15
16 if (equaln(p1, 0) || equaln(p1, 1) || equaln(p1, -1))
17 push(p1)
18 restore()
19 return
20
21 n_factor_number = p1.q.a
22
23 h = tos
24
25 factor_a()
26
27 if (tos - h > 1)
28 list(tos - h)
29 push_symbol(MULTIPLY)
30 swap()
31 cons()
32
33 restore()
34
35# factor using table look-up, then switch to rho method if necessary
36
37# From TAOCP Vol. 2 by Knuth, p. 380 (Algorithm A)
38
39factor_a = ->
40 k = 0
41
42 if (n_factor_number.isNegative())
43 n_factor_number = setSignTo(n_factor_number, 1)
44 push_integer(-1)
45
46 for k in [0...10000]
47
48 try_kth_prime(k)
49
50 # if n_factor_number is 1 then we're done
51
52 if (n_factor_number.compare(1) == 0)
53 return
54
55 factor_b()
56
57try_kth_prime = (k) ->
58 count = 0
59
60 d = mint(primetab[k])
61
62 count = 0
63
64 while (1)
65
66 # if n_factor_number is 1 then we're done
67
68 if (n_factor_number.compare(1) == 0)
69 if (count)
70 push_factor(d, count)
71 return
72
73 [q,r] = mdivrem(n_factor_number, d)
74
75 # continue looping while remainder is zero
76
77 if (r.isZero())
78 count++
79 n_factor_number = q
80 else
81 break
82
83 if (count)
84 push_factor(d, count)
85
86 # q = n_factor_number/d, hence if q < d then
87 # n_factor_number < d^2 so n_factor_number is prime
88
89 if (mcmp(q, d) == -1)
90 push_factor(n_factor_number, 1)
91 n_factor_number = mint(1)
92
93
94# From TAOCP Vol. 2 by Knuth, p. 385 (Algorithm B)
95
96factor_b = ->
97 k = 0
98 l = 0
99
100 bigint_one = mint(1)
101 x = mint(5)
102 xprime = mint(2)
103
104 k = 1
105 l = 1
106
107 while (1)
108
109 if (mprime(n_factor_number))
110 push_factor(n_factor_number, 1)
111 return 0
112
113 while (1)
114
115 if (esc_flag)
116 stop("esc")
117
118 # g = gcd(x' - x, n_factor_number)
119
120 t = msub(xprime, x)
121 t = setSignTo(t,1)
122 g = mgcd(t, n_factor_number)
123
124 if (MEQUAL(g, 1))
125 if (--k == 0)
126 xprime = x
127 l *= 2
128 k = l
129
130 # x = (x ^ 2 + 1) mod n_factor_number
131
132 t = mmul(x, x)
133 x = madd(t, bigint_one)
134 t = mmod(x, n_factor_number)
135 x = t
136
137 continue
138
139 push_factor(g, 1)
140
141 if (mcmp(g, n_factor_number) == 0)
142 return -1
143
144 # n_factor_number = n_factor_number / g
145
146 t = mdiv(n_factor_number, g)
147 n_factor_number = t
148
149 # x = x mod n_factor_number
150
151 t = mmod(x, n_factor_number)
152 x = t
153
154 # xprime = xprime mod n_factor_number
155
156 t = mmod(xprime, n_factor_number)
157 xprime = t
158
159 break
160
161push_factor = (d, count) ->
162 p1 = new U()
163 p1.k = NUM
164 p1.q.a = d
165 p1.q.b = mint(1)
166 push(p1)
167 if (count > 1)
168 push_symbol(POWER)
169 swap()
170 p1 = new U()
171 p1.k = NUM
172 p1.q.a = mint(count)
173 p1.q.b = mint(1)
174 push(p1)
175 list(3)