/*

  Author:  Pate Williams (c) 1997



  Barrett modular reduction. See "Handbook of Applied

  Cryptography" by Alfred J. Menezes et al 14.3.3

  Section pages 603 - 604.

*/



#include <stdio.h>

#include "lip.h"



#define DEBUG



#ifdef DEBUG

void radix_representation(long b, long n, long *a, verylong za)

{

  long i = 0;

  verylong zb = 0, zq = 0, zr = 0, zx = 0;



  zintoz(b, &zb);

  zcopy(za, &zx);

  do {

    zdiv(zx, zb, &zq, &zr);

    a[i++] = ztoint(zr);

    zcopy(zq, &zx);

  } while (zscompare(zq, 0l) != 0);

  while (i < n) a[i++] = 0;

  zfree(&zb);

  zfree(&zq);

  zfree(&zr);

  zfree(&zx);

}

#endif



void Barrett_reduction(long b, long k,

                       verylong zm, verylong zmu,

                       verylong zx, verylong *zr)

{

  verylong zb = 0, zbk = 0, zq1 = 0, zq2 = 0, zq3 =0;

  verylong zr1 = 0, zr2 = 0;



  zintoz(b, &zb);

  zsexp(zb, k - 1, &zbk);

  zdiv(zx, zbk, &zq1, &zq2);

  zmul(zq1, zmu, &zq2);

  zsexp(zb, k + 1, &zbk);

  zdiv(zq2, zbk, &zq3, &zr1);

  zmod(zx, zbk, &zr1);

  zmulmod(zq3, zm, zbk, &zr2);

  zsub(zr1, zr2, zr);

  if (zscompare(*zr, 0l) < 0) {

    zadd(*zr, zbk, &zb);

    zcopy(zb, zr);

  }

  #ifdef DEBUG

  {

    long a[256], i, k2 = (k + 1) * 2;



    radix_representation(b, k2, a, zq1);

    printf("q1 = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(zq1);

    radix_representation(b, k2, a, zq2);

    printf("q2 = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(zq2);

    radix_representation(b, k2, a, zq3);

    printf("q3 = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(zq3);

    radix_representation(b, k2, a, zr1);

    printf("r1 = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(zr1);

    radix_representation(b, k2, a, zr2);

    printf("r2 = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(zr2);

    radix_representation(b, k2, a, *zr);

    printf("r  = ");

    for (i = k2 - 1; i >= 0; i--)

      printf("%ld", a[i]);

    printf(" = ");

    zwriteln(*zr);

  }

  #endif

  while (zcompare(*zr, zm) >= 0) {

    zsub(*zr, zm, &zb);

    zcopy(zb, zr);

  }

  zfree(&zb);

  zfree(&zbk);

  zfree(&zq1);

  zfree(&zq2);

  zfree(&zq3);

  zfree(&zr1);

  zfree(&zr2);

}



int main(void)

{

  verylong zm = 0, zmu = 0, zr = 0, zx = 0;



  zintoz(47l, &zm);

  zintoz(87l, &zmu);

  zintoz(3561l, &zx);

  Barrett_reduction(4, 3, zm, zmu, zx, &zr);

  zfree(&zm);

  zfree(&zmu);

  zfree(&zr);

  zfree(&zx);

  return 0;

}