176 lines
2.7 KiB
CoffeeScript
176 lines
2.7 KiB
CoffeeScript
# Factor using the Pollard rho method
|
|
|
|
|
|
|
|
n_factor_number = 0
|
|
|
|
factor_number = ->
|
|
h = 0
|
|
|
|
save()
|
|
|
|
p1 = pop()
|
|
|
|
# 0 or 1?
|
|
|
|
if (equaln(p1, 0) || equaln(p1, 1) || equaln(p1, -1))
|
|
push(p1)
|
|
restore()
|
|
return
|
|
|
|
n_factor_number = p1.q.a
|
|
|
|
h = tos
|
|
|
|
factor_a()
|
|
|
|
if (tos - h > 1)
|
|
list(tos - h)
|
|
push_symbol(MULTIPLY)
|
|
swap()
|
|
cons()
|
|
|
|
restore()
|
|
|
|
# factor using table look-up, then switch to rho method if necessary
|
|
|
|
# From TAOCP Vol. 2 by Knuth, p. 380 (Algorithm A)
|
|
|
|
factor_a = ->
|
|
k = 0
|
|
|
|
if (n_factor_number.isNegative())
|
|
n_factor_number = setSignTo(n_factor_number, 1)
|
|
push_integer(-1)
|
|
|
|
for k in [0...10000]
|
|
|
|
try_kth_prime(k)
|
|
|
|
# if n_factor_number is 1 then we're done
|
|
|
|
if (n_factor_number.compare(1) == 0)
|
|
return
|
|
|
|
factor_b()
|
|
|
|
try_kth_prime = (k) ->
|
|
count = 0
|
|
|
|
d = mint(primetab[k])
|
|
|
|
count = 0
|
|
|
|
while (1)
|
|
|
|
# if n_factor_number is 1 then we're done
|
|
|
|
if (n_factor_number.compare(1) == 0)
|
|
if (count)
|
|
push_factor(d, count)
|
|
return
|
|
|
|
[q,r] = mdivrem(n_factor_number, d)
|
|
|
|
# continue looping while remainder is zero
|
|
|
|
if (r.isZero())
|
|
count++
|
|
n_factor_number = q
|
|
else
|
|
break
|
|
|
|
if (count)
|
|
push_factor(d, count)
|
|
|
|
# q = n_factor_number/d, hence if q < d then
|
|
# n_factor_number < d^2 so n_factor_number is prime
|
|
|
|
if (mcmp(q, d) == -1)
|
|
push_factor(n_factor_number, 1)
|
|
n_factor_number = mint(1)
|
|
|
|
|
|
# From TAOCP Vol. 2 by Knuth, p. 385 (Algorithm B)
|
|
|
|
factor_b = ->
|
|
k = 0
|
|
l = 0
|
|
|
|
bigint_one = mint(1)
|
|
x = mint(5)
|
|
xprime = mint(2)
|
|
|
|
k = 1
|
|
l = 1
|
|
|
|
while (1)
|
|
|
|
if (mprime(n_factor_number))
|
|
push_factor(n_factor_number, 1)
|
|
return 0
|
|
|
|
while (1)
|
|
|
|
if (esc_flag)
|
|
stop("esc")
|
|
|
|
# g = gcd(x' - x, n_factor_number)
|
|
|
|
t = msub(xprime, x)
|
|
t = setSignTo(t,1)
|
|
g = mgcd(t, n_factor_number)
|
|
|
|
if (MEQUAL(g, 1))
|
|
if (--k == 0)
|
|
xprime = x
|
|
l *= 2
|
|
k = l
|
|
|
|
# x = (x ^ 2 + 1) mod n_factor_number
|
|
|
|
t = mmul(x, x)
|
|
x = madd(t, bigint_one)
|
|
t = mmod(x, n_factor_number)
|
|
x = t
|
|
|
|
continue
|
|
|
|
push_factor(g, 1)
|
|
|
|
if (mcmp(g, n_factor_number) == 0)
|
|
return -1
|
|
|
|
# n_factor_number = n_factor_number / g
|
|
|
|
t = mdiv(n_factor_number, g)
|
|
n_factor_number = t
|
|
|
|
# x = x mod n_factor_number
|
|
|
|
t = mmod(x, n_factor_number)
|
|
x = t
|
|
|
|
# xprime = xprime mod n_factor_number
|
|
|
|
t = mmod(xprime, n_factor_number)
|
|
xprime = t
|
|
|
|
break
|
|
|
|
push_factor = (d, count) ->
|
|
p1 = new U()
|
|
p1.k = NUM
|
|
p1.q.a = d
|
|
p1.q.b = mint(1)
|
|
push(p1)
|
|
if (count > 1)
|
|
push_symbol(POWER)
|
|
swap()
|
|
p1 = new U()
|
|
p1.k = NUM
|
|
p1.q.a = mint(count)
|
|
p1.q.b = mint(1)
|
|
push(p1)
|
|
list(3)
|