#include <stdio.h>
#include <stdlib.h>

const int N = 19;
const int k = 5;
int R;
int invN;

/* Fonction pour reduire modulo R, en extrayant les k derniers bits a
l'aide du "et" bit-a-bit */
int moduloR (int n) {
 return n & (R-1);
}

/* Inverse de N modulo R */
/* On programme l'algo d'Euclide etendu... 
   Renvoie l'inverse de a modulo b */
int euclEt (int a, int b) {
  int r = a; int r2 = b; int u = 1; int v = 0; int u2 = 0; int v2 = 1;
  int q, rs, us, vs;
  while (r2 != 0) {
   q = r/r2;
   rs = r; us = u; vs = v;
   r = r2; u = u2; v = v2;
   r2 = rs -q *r2; u2 = us - q*u2; v2 = vs -q*v2;
  }
 if (u < 0) u += R;
 return u;
}

/* reduction de Montgomery, renvoie x*R^(-1) modulo N */
int reduction (x) {
 int q = moduloR(moduloR(x)*(-invN));
 int y = (x+q*N) >> k; // decalage de digits : revient a diviser par R
 if (y < N) return y; else return y-N;
}

/* exponentiation rapide de Montgomery */
/* renvoie y^n*R */
int RmodN;
int expoM(y,n) {
  //  int y = x << k;
  if (n==0) return RmodN;
  else {
    int aux = expoM(y,n/2);
    if (n%2) return reduction(reduction(y*aux)*aux);
    else return reduction(aux*aux);
  }
}

/* renvoie x^n */
int expo(x,n) {
  return reduction(expoM(x<<k,n));
}

int main () {
 R = 1<<k; // ou pow(2,k) avec math.h
 invN = euclEt(N,R);
 RmodN = R%N;

 printf("%d\n",reduction(8));
 printf("%d\n",expo(5,6));
 return 0;
}