1implement Mersenne; 2 3include "sys.m"; 4 sys : Sys; 5include "draw.m"; 6include "ipints.m"; 7 ipints: IPints; 8 IPint: import ipints; 9 10# Test primality of Mersenne numbers 11 12Mersenne: module 13{ 14 init: fn(nil: ref Draw->Context, argv: list of string); 15}; 16 17init(nil: ref Draw->Context, argv: list of string) 18{ 19 sys = load Sys Sys->PATH; 20 ipints = load IPints IPints->PATH; 21 p := 3; 22 if(tl argv != nil) 23 p = int hd tl argv; 24 if(isprime(p) && (p == 2 || lucas(p))) 25 s := ""; 26 else 27 s = "not "; 28 sys->print("2^%d-1 is %sprime\n", p, s); 29} 30 31# s such that s^2 <= n 32sqrt(n: int): int 33{ 34 v := n; 35 r := 0; 36 for(t := 1<<30; t; t >>= 2){ 37 if(t+r <= v){ 38 v -= t+r; 39 r = (r>>1)|t; 40 } 41 else 42 r = r>>1; 43 } 44 return r; 45} 46 47isprime(n: int): int 48{ 49 if(n < 2) 50 return 0; 51 if(n == 2) 52 return 1; 53 if((n&1) == 0) 54 return 0; 55 s := sqrt(n); 56 for(i := 3; i <= s; i += 2) 57 if(n%i == 0) 58 return 0; 59 return 1; 60} 61 62pow(b : ref IPint, n : int): ref IPint 63{ 64 zero := IPint.inttoip(0); 65 one := IPint.inttoip(1); 66 if((b.cmp(zero) == 0 && n != 0) || b.cmp(one) == 0 || n == 1) 67 return b; 68 if(n == 0) 69 return one; 70 c := b; 71 b = one; 72 while(n){ 73 while(!(n & 1)){ 74 n >>= 1; 75 c = c.mul(c); 76 } 77 n--; 78 b = c.mul(b); 79 } 80 return b; 81} 82 83lucas(p: int): int 84{ 85 zero := IPint.inttoip(0); 86 one := IPint.inttoip(1); 87 two := IPint.inttoip(2); 88 bigp := pow(two, p).sub(one); 89 u := IPint.inttoip(4); 90 for(i := 2; i < p; i++){ 91 u = u.mul(u); 92 if(u.cmp(two) <= 0) 93 u = two.sub(u); 94 else 95 u = u.sub(two).expmod(one, bigp); 96 } 97 return u.cmp(zero) == 0; 98} 99 100