functions.qcl

set allow-redefines 1;

// returns the smallest factor > 1 of n or 1 if n is prime 

int findfactor(int n) {
  int i;
  if n<=0 { exit "findfactor takes only positive args"; }
  for i=2 to floor(sqrt(n)) {
    if n mod i == 0 { return i; }
  }
  return 1;
}

// test if n is a prime number

boolean testprime(int n) {
  int i;
  if n<=1 { return false; }
  for i=2 to floor(sqrt(n)) {
    if n mod i == 0 { return false; }
  }
  return true;
}

// test if n is a prime power

boolean testprimepower(int n) {
  int i;
  int f;
  i=2;
  while i<=floor(sqrt(n)) and f==0 {
    if n mod i == 0 { f=i; }
    i=i+1;
  }
  for i=2 to floor(log(n,f)) {
    if f^i==n { return true; }
  }
  return false;
}

// returns x^a mod n

int powmod(int x,int a,int n) {
  int u=x;
  int y=1;
  int i;
  
  for i=0 to 30 {
    if a/2^i mod 2 == 1 { y=y*u mod n; }
    u=u^2 mod n;
  }
  return y;
}

// return the modular inverse to a mod n or 0 if gcd(a,n)>1

int invmod(int a,int n) {
  int b=a;
  int i;

  if gcd(a,n)>1 { return 0; }  
  for i=1 to n {
    if b*a mod n == 1 { return b; } 
    b=b*a mod n;
  }
  return 0;
}

// finds the denominator q of the best rational approximation p/q
// for x with q<qmax

int denominator(real x,int qmax) {
  real y=x;
  real z;
  int q0;
  int q1=1;
  int q2;
  
  while true {
    z=y-floor(y);
    if z<0.5/qmax^2 { return q1; }
    y=1/z;
    q2=floor(y)*q1+q0;
    if q2>=qmax { return q1; }
    q0=q1; q1=q2;
  }
}

set allow-redefines 0;