-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathis_almost_prime.sf
72 lines (54 loc) · 1.89 KB
/
is_almost_prime.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
#!/usr/bin/ruby
# A fast method for checking if a given integer n has exactly k prime factors counted with multiplicity (i.e.: bigomega(n) = k).
# Based on algorithm from the PARI program posted by Jinyuan Wang:
# https://oeis.org/history?seq=A219019
func my_is_almost_prime (n, k) {
if (k == 1) {
return n.is_prime
}
return false if (n <= 1)
var f = [n]
var r = 0
while (k > 1) {
var t = f.last.iroot(k)+1
t > r || return false
r = t
f = factor_upto(f.last, r)
if (f.len == 1) {
return false
}
if (f.last.is_prime) {
return (k == f.len)
}
k = (k - f.len + 1)
}
return false
}
for n in (1..100) {
var t = irand(1 << n)+1
say "Testing: #{t}"
for k in (1..20) {
if (t.is_almost_prime(k)) {
assert(my_is_almost_prime(t, k), [t,k])
}
elsif (my_is_almost_prime(t, k)) {
die "counter-example: #{[t, k]}"
}
}
}
# These tests may randomly fail
assert(my_is_almost_prime(1762610652661, 4))
assert(my_is_almost_prime(17295389739104958689918, 6))
assert(my_is_almost_prime(3194569977264671866214863610, 4))
assert(my_is_almost_prime(72960029911372484469592, 6))
assert(my_is_almost_prime(158441958340314419522552997965, 5))
assert(my_is_almost_prime(1038735417774685624754922, 8))
assert(my_is_almost_prime(323463936073052242257157, 2))
assert(my_is_almost_prime(63625815508153350513817311207, 4))
assert(my_is_almost_prime(2468336534484213369852218139, 4))
assert(my_is_almost_prime(307764215443043830937725274, 3))
assert(my_is_almost_prime(2068360226936926347380742889, 4))
assert(my_is_almost_prime(3997087877499983596169090, 6))
assert(my_is_almost_prime(51403227993662852987455497928, 7))
assert(my_is_almost_prime(2011545292777067446258687026, 5))
assert(my_is_almost_prime(4810554005148753711362794456, 6))