-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathchebyshev_factorization_method.sf
114 lines (89 loc) · 3.28 KB
/
chebyshev_factorization_method.sf
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#!/usr/bin/ruby
# Daniel "Trizen" Șuteu
# Date: 11 December 2019
# https://github.com/trizen
# A simple factorization method, using the Chebyshev T_n(x) polynomials, based on the identity:
# T_{m n}(x) = T_m(T_n(x))
# where:
# T_n(x) = (1/2) * V_n(2x, 1)
# where V_n(P, Q) is the Lucas V sequence.
# See also:
# https://oeis.org/A001075
# https://en.wikipedia.org/wiki/Lucas_sequence
# https://en.wikipedia.org/wiki/Iterated_function
# https://en.wikipedia.org/wiki/Chebyshev_polynomials
func lucasVmod(P, n, m) { # assumes Q = 1
var(V1 = 2, V2 = P)
for bit in (n.bits) {
if (bit) {
V1 = submod(mulmod(V2, V1, m), P, m)
V2 = mulmod(V2, V2, m)-2
}
else {
V2 = submod(mulmod(V2, V1, m), P, m)
V1 = mulmod(V1, V1, m)-2
}
}
return V1
}
func chebyshev_factorization(n, B1 = n.ilog(6)**3, a = 3) {
var x = a
var G = B1*B1
var i = invmod(2, n)
func chebyshevTmod(a, x) {
mulmod(lucasVmod(2*x, a, n), i, n)
}
primes_each(2, B1.isqrt, {|p|
G.ilog(p).times {
x = chebyshevTmod(p, x) # T_k(x) (mod n)
}
})
primes_each(B1.isqrt+1, B1, {|p|
x = chebyshevTmod(p, x) # T_k(x) (mod n)
is_coprime(x-1, n) || return gcd(x-1, n)
})
return gcd(x-1, n)
}
say chebyshev_factorization(257221 * 470783, 1000); #=> 257221 (p+1 is 1000-smooth)
say chebyshev_factorization(1124075136413 * 3556516507813, 4000); #=> 1124075136413 (p+1 is 4000-smooth)
say chebyshev_factorization(7553377229 * 588103349, 800); #=> 7553377229 (p+1 is 800-smooth)
say chebyshev_factorization(333732865481 * 1632480277613, 3000); #=> 333732865481 (p-1 is 3000-smooth)
say ''
say chebyshev_factorization(15597344393 * 12388291753, 3000) #=> 15597344393 (p-1 is 3000-smooth)
say chebyshev_factorization(43759958467 * 59037829639, 3200) #=> 43759958467 (p+1 is 3200-smooth)
say chebyshev_factorization(112601635303 * 83979783007, 700) #=> 112601635303 (p-1 is 700-smooth)
say chebyshev_factorization(228640480273 * 224774973299, 2000) #=> 228640480273 (p-1 is 2000-smooth)
say ''
say chebyshev_factorization(5140059121 * 8382882743, 2500, 2) #=> 5140059121 (p-1 is 2500-smooth)
say chebyshev_factorization(18114813019 * 17402508649, 6000, 2) #=> 18114813019 (p+1 is 6000-smooth)
say chebyshev_factorization(533091092393 * 440050095029, 300, 2) #=> 533091092393 (p+1 is 300-smooth)
say ''
for k in (3..20) {
var n = 2.of { random_nbit_prime(k) }.prod
var p = chebyshev_factorization(n, 2 * n.ilog2**2)
if (p.is_prime) {
say "#{n} = #{p} * #{n/p}"
}
}
__END__
143 = 11 * 13
323 = 17 * 19
2773 = 59 * 47
12091 = 107 * 113
42781 = 179 * 239
181427 = 433 * 419
685603 = 967 * 709
3240379 = 1999 * 1621
12254063 = 3089 * 3967
38271817 = 6763 * 5659
119391583 = 9511 * 12553
367073731 = 20929 * 17539
3578349739 = 60961 * 58699
15888478703 = 121571 * 130693
36182994919 = 143593 * 251983
545782820497 = 668713 * 816169
3656558375893 = 1912639 * 1911787
12346497466267 = 3737821 * 3303127
31369209076123 = 5655031 * 5547133
148520973663353 = 13394431 * 11088263
412638129424819 = 17334017 * 23805107