#include <stdint.h>
#include <stdio.h>
#include <inttypes.h>

uint32_t modexp(uint32_t a, unsigned char e[32], uint32_t m)
{
  /* TODO: implement */
  int i, j;
  uint64_t precomputed[16];
  precomputed[0] = 1;
  for(i=1; i<15;i++)
  {
      precomputed[i] = precomputed[i-1]*a;
      precomputed[i] = precomputed [i] % m;
  }

  uint64_t r = 1;
  for(i=31; i>=0; i--)
  {
      for(j=3; j>=0; j--)
      {
          r = r*r;
          r = r % m;
      }
      unsigned char x = e[i]<<4;
      x = x>>4;
      r = r * precomputed[x];
      r = r % m;
      for(j=3; j>=0; j--)
      {
          r = r*r;
          r = r % m;
      }
      x = e[i]>>4;
      r = r * precomputed[x];
      r = r % m;
  }
  return r;
}

/* Pipe output through sage */
#define NTESTS 20
int main(void)
{
  //FILE *urandom = fopen("/dev/urandom", "r");
  uint32_t a,r,ei, m;
  unsigned char e[32];
  int i,j;

  for(i=0;i<NTESTS;i++)
  {
    a = (uint32_t)rand();
    m = (uint32_t)rand();
    for(j=0;j<32;j++)
        e[j] = (char) rand();
    r = modexp(a,e,m);
    ei = 0;
    for(j=0;j<4;j++)
      ei |= (uint32_t)e[j] << 8*j;

    printf("power_mod(%u,%u,%u) - %u\n", a,ei,m,r);
  }

  //fclose(urandom);
  return 0;
}