1set library 1;
2
3// returns the smallest factor > 1 of n or 1 if n is prime
4
5int findfactor(int n) {
6  int i;
7  if n<=0 { exit "findfactor takes only positive args"; }
8  for i=2 to floor(sqrt(n)) {
9    if n mod i == 0 { return i; }
10  }
11  return 1;
12}
13
14// test if n is a prime number
15
16boolean testprime(int n) {
17  int i;
18  if n<=1 { return false; }
19  if n==2 or n==3 { return true; }
20  if n mod 2 == 0 { return false; }
21  for i=3 to floor(sqrt(n)) step 2 {
22    if n mod i == 0 { return false; }
23  }
24  return true;
25}
26
27// test if n is a prime power
28
29boolean testprimepower(int n) {
30  int i;
31  int f;
32  i=2;
33  while i<=floor(sqrt(n)) and f==0 {
34    if n mod i == 0 { f=i; }
35    i=i+1;
36  }
37  for i=2 to floor(log(n,f)) {
38    if f^i==n { return true; }
39  }
40  return false;
41}
42
43// returns x^a mod n
44
45int powmod(int x,int a,int n) {
46  int u=x;
47  int y=1;
48  int i;
49
50  for i=0 to 30 {
51    if a/2^i mod 2 == 1 { y=y*u mod n; }
52    u=u^2 mod n;
53  }
54  return y;
55}
56
57// return the modular inverse to a mod n or 0 if gcd(a,n)>1
58
59int invmod(int a,int n) {
60  int b=a;
61  int i;
62
63  if gcd(a,n)>1 { return 0; }
64  for i=1 to n {
65    if b*a mod n == 1 { return b; }
66    b=b*a mod n;
67  }
68  return 0;
69}
70
71// finds the denominator q of the best rational approximation p/q
72// for x with q<qmax
73
74int denominator(real x,int qmax) {
75  real y=x;
76  real z;
77  int q0;
78  int q1=1;
79  int q2;
80
81  while true {
82    z=y-floor(y);
83    if z<0.5/qmax^2 { return q1; }
84    y=1/z;
85    q2=floor(y)*q1+q0;
86    if q2>=qmax { return q1; }
87    q0=q1; q1=q2;
88  }
89}
90
91set library 0;
92