blob: 7efc4a97a750738a13f6e956869e8f5d765dc8be (
plain)
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
|
implement Mersenne;
include "sys.m";
sys : Sys;
include "draw.m";
include "ipints.m";
ipints: IPints;
IPint: import ipints;
# Test primality of Mersenne numbers
Mersenne: module
{
init: fn(nil: ref Draw->Context, argv: list of string);
};
init(nil: ref Draw->Context, argv: list of string)
{
sys = load Sys Sys->PATH;
ipints = load IPints IPints->PATH;
p := 3;
if(tl argv != nil)
p = int hd tl argv;
if(isprime(p) && (p == 2 || lucas(p)))
s := "";
else
s = "not ";
sys->print("2^%d-1 is %sprime\n", p, s);
}
# s such that s^2 <= n
sqrt(n: int): int
{
v := n;
r := 0;
for(t := 1<<30; t; t >>= 2){
if(t+r <= v){
v -= t+r;
r = (r>>1)|t;
}
else
r = r>>1;
}
return r;
}
isprime(n: int): int
{
if(n < 2)
return 0;
if(n == 2)
return 1;
if((n&1) == 0)
return 0;
s := sqrt(n);
for(i := 3; i <= s; i += 2)
if(n%i == 0)
return 0;
return 1;
}
pow(b : ref IPint, n : int): ref IPint
{
zero := IPint.inttoip(0);
one := IPint.inttoip(1);
if((b.cmp(zero) == 0 && n != 0) || b.cmp(one) == 0 || n == 1)
return b;
if(n == 0)
return one;
c := b;
b = one;
while(n){
while(!(n & 1)){
n >>= 1;
c = c.mul(c);
}
n--;
b = c.mul(b);
}
return b;
}
lucas(p: int): int
{
zero := IPint.inttoip(0);
one := IPint.inttoip(1);
two := IPint.inttoip(2);
bigp := pow(two, p).sub(one);
u := IPint.inttoip(4);
for(i := 2; i < p; i++){
u = u.mul(u);
if(u.cmp(two) <= 0)
u = two.sub(u);
else
u = u.sub(two).expmod(one, bigp);
}
return u.cmp(zero) == 0;
}
|
