/*
  Author:  Pate Williams (c) 1997

  "Algorithm 1.2.3 (Left-Right Base 2 ^ k).
  Given g in G and n in Z, this algorithm computes
  g ^ n in G. If n != 0, we assume also that we
  are given the unique integer e such that
  2 ^ ke <= |n| < 2 ^ (ke + 1). We write 1 for the
  unit element of G." -Henri Cohen- See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen page 10.  We specialize the code to
  the field Zp where p is of course prime.
*/

#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include "lip.h"

#define BITS_PER_LONG 32

void radix_representation(long b, long x,
                          long *a, long *t)
{
  long i = 0, q;

  q = x / b;
  a[i++] = x - q * b;
  while (q > 0) {
    x = q;
    q = x / b;
    a[i++] = x - q * b;
  }
  *t = i;
}

void left_right_base(verylong zg, long n,
                     verylong zp, verylong *zy)
{
  long a, b, e, f, i, k, ke, l, ln, N;
  long base, digit[BITS_PER_LONG], t;
  verylong za = 0, zz = 0, zx[16];

  if (n == 0) {
    zone(zy);
    return;
  }
  ln = labs(n);
  ke = floor(log(ln) / log(2.0));
  if (ln <= pow(2, 8)) k = 1;
  else if (ln <= pow(2, 24)) k = 2;
  else k = 3;
  e = ke / k;
  f = e;
  base = pow(2, k);
  l = base - 1;
  if (n < 0) {
    N = - n;
    zinvmod(zg, zp, &zz);
  }
  else {
    N = n;
    zcopy(zg, &zz);
  }
  zcopy(zz, zy);
  for (i = 1; i <= l; i += 2) zx[i] = 0;
  zcopy(zz, &zx[1]);
  for (i = 3; i <= l; i += 2)
    zsexpmod(zz, i, zp, &zx[i]);
  radix_representation(base, N, digit, &t);
  do {
    a = digit[f];
    if (a == 0) {
      for (i = 0; i < k; i++) {
        zmulmod(*zy, *zy, zp, &za);
        zcopy(za, zy);
      }
    }
    else {
      t = 0;
      b = a;
      while (!(b & 1)) b /= 2, t++;
      if (f != e) {
        for (i = 0; i < k - t; i++) {
          zmulmod(*zy, *zy, zp, &za);
          zcopy(za, zy);
        }
        zmulmod(zx[b], *zy, zp, &za);
        zcopy(za, zy);
      }
      else
        zcopy(zx[b], zy);
      for (i = 0; i < t; i++) {
        zmulmod(*zy, *zy, zp, &za);
        zcopy(za, zy);
      }
    }
    f--;
  } while (f >= 0);
  zfree(&za);
  zfree(&zz);
  for (i = 1; i <= l; i += 2) zfree(&zx[i]);
}

long get_number(void)
{
  int c, sign = 1;
  long number;

  while ((c = getchar()) <= ' ');
  if (c == '+') {
    sign = 1;
    c = getchar();
  }
  else if (c == '-') {
    sign = - 1;
    c = getchar();
  }
  number = c - '0';
  while ((c = getchar()) != '\n')
    number = 10 * number + c - '0';
  return sign * number;
}

int main(void)
{
  char answer[256];
  long n;
  verylong zg = 0, zp = 0, zy = 0, zz = 0;

  do {
    printf("g = "); zread(&zg);
    printf("n = "); n = get_number();
    do {
      printf("p = "); zread(&zp);
    } while (!zprobprime(zp, 5l));
    left_right_base(zg, n, zp, &zy);
    printf("y = g ^ n = "); zwriteln(zy);
    zsexpmod(zg, n, zp, &zz);
    printf("z = g ^ n = "); zwriteln(zz);
    if (zcompare(zy, zz) != 0)
      printf("error - in calculation\n");
    printf("another calculation (n or y) ? ");
    scanf("%s", answer);
  } while (tolower(answer[0]) == 'y');
  zfree(&zg);
  zfree(&zp);
  zfree(&zy);
  zfree(&zz);
  return 0;
}