/*
  Author:  Pate Williams (c) 1997

  "Algorithm 1.2.3 (Left-Right Binary, Using Bits).
  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 ^ e <= |n| < 2 ^ (e + 1). We write 1 for the
  unit element of G." -Henri Cohen- See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen page 9.  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 binary(long n, int *bit, int *count)
{
  int found = 0, i;
  long t = n;

  /* find the highest order 1 bit */
  for (i = 0; !found && i < BITS_PER_LONG; i++) {
    found = t & 0x80000000;
    if (found) *count = BITS_PER_LONG - i;
    else t <<= 1;
  }
  for (i = 0; i < *count; i++)
    bit[i] = n & 1, n >>= 1;
}

void left_right_binary(verylong zg, long n,
                       verylong zp, verylong *zy)
{
  int bit[BITS_PER_LONG], count;
  long e, f, N;
  verylong za = 0, zz = 0;

  if (n == 0) {
    zone(zy);
    return;
  }
  e = floor(log(labs(n)) / log(2.0));
  f = e;
  if (n < 0) {
    N = - n;
    zinvmod(zg, zp, &zz);
  }
  else {
    N = n;
    zcopy(zg, &zz);
  }
  zcopy(zz, zy);
  binary(N, bit, &count);
  while (f != 0) {
    f--;
    zmulmod(*zy, *zy, zp, &za);
    zcopy(za, zy);
    if (bit[f]) {
      zmulmod(zz, *zy, zp, &za);
      zcopy(za, zy);
    }
  }
  zfree(&za);
  zfree(&zz);
}

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_binary(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;
}