/*
  Author:  Pate Williams (c) 1997

  "Exercise 4.2. Fibonacci numbers.
  Use (4.45)-(4.46) and the binary representation
  of m to deduce a set of formulas for the
  computation of U_m and V_m. Write a computer
  program Fibonacci which reads P, Q, m, and N."
  -Hans Riesel-
  See "Prime Numbers and Computer Methods for
  Factorization" by Hans Riesel second edition
  page 111.
*/

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

void mcopy(verylong **zb, verylong **za)
/* a = b */
{
  zcopy(zb[0][0], &za[0][0]);
  zcopy(zb[0][1], &za[0][1]);
  zcopy(zb[1][0], &za[1][0]);
  zcopy(zb[1][1], &za[1][1]);
}

verylong **mcreate(void)
{
  verylong **zm = calloc(2, sizeof(verylong *));

  zm[0] = calloc(2, sizeof(verylong));
  zm[1] = calloc(2, sizeof(verylong));
  return zm;
}

void mfree(verylong **za)
{
  zfree(&za[0][0]);
  zfree(&za[0][1]);
  zfree(&za[1][0]);
  zfree(&za[1][1]);
  free(za[0]);
  free(za[1]);
  free(za);
}

void mmultiply(verylong **za, verylong **zb,
               verylong **zc)
/* c = a * b */
{
  verylong zd = 0, ze = 0;

  zmul(za[0][0], zb[0][0], &zd);
  zmul(za[0][1], zb[1][0], &ze);
  zadd(zd, ze, &zc[0][0]);
  zmul(za[0][0], zb[0][1], &zd);
  zmul(za[0][1], zb[1][1], &ze);
  zadd(zd, ze, &zc[0][1]);
  zmul(za[1][0], zb[0][0], &zd);
  zmul(za[1][1], zb[1][0], &ze);
  zadd(zd, ze, &zc[1][0]);
  zmul(za[1][0], zb[0][1], &zd);
  zmul(za[1][1], zb[1][1], &ze);
  zadd(zd, ze, &zc[1][1]);
  zfree(&zd);
  zfree(&ze);
}

void msexp(verylong **za, long m, verylong **zx)
/* returns x = a ^ m */
{
  verylong **zb = mcreate();
  verylong **zs = mcreate();

  mcopy(za, zs);
  zone(&zx[0][0]);
  zzero(&zx[0][1]);
  zzero(&zx[1][0]);
  zone(&zx[1][1]);
  while (m != 0) {
    if (m & 1) {
      mmultiply(zx, zs, zb);
      mcopy(zb, zx);
    }
    m >>= 1;
    if (m != 0) {
      mmultiply(zs, zs, zb);
      mcopy(zb, zs);
    }
  }
  mfree(zb);
  mfree(zs);
}

void mwrite(verylong **za)
{
  zwrite(za[0][0]);
  printf(" ");
  zwriteln(za[0][1]);
  zwrite(za[1][0]);
  printf(" ");
  zwriteln(za[1][1]);
}

int main(void)
{
  long P, Q, m;
  verylong **za = mcreate();
  verylong **zx = mcreate();
  verylong **zy = mcreate();
  verylong **zz = mcreate();

  printf("P = ");
  scanf("%ld", &P);
  printf("Q = ");
  scanf("%ld", &Q);
  printf("m = ");
  scanf("%ld", &m);
  zintoz(P, &za[0][0]);
  zintoz(- Q, &za[0][1]);
  zone(&za[1][0]);
  zzero(&za[1][1]);
  msexp(za, m, zx);
  mwrite(zx);
  zone(&zy[0][0]);
  zone(&zy[0][1]);
  zzero(&zy[1][0]);
  zintoz(2, &zy[1][1]);
  mmultiply(zx, zy, zz);
  printf("U_%ld = ", m + 1);
  zwriteln(zz[0][0]);
  printf("V_%ld = ", m + 1);
  zwriteln(zz[0][1]);
  printf("U_%ld = ", m);
  zwriteln(zz[1][0]);
  printf("V_%ld = ", m);
  zwriteln(zz[1][1]);
  mfree(za);
  mfree(zx);
  mfree(zy);
  mfree(zz);
  return 0;
}