/* arith-update.c:  MORE Number-Theoretic arithmetic routines
** Written by Michael J. Gourlay
** Fall 1996, University of Colorado at Boulder
** Algorithms
*/

/*

Top level description:
----------------------

These are supplementary routines to the previous "arith.c" package
for the NumberT abstract data type.  Two new operations are
implemented here:  Exponentiation and "Inversion", where "Inversion"
is a by-product of a more general routine which finds the Greatest
Common Divisor of two numbers.  In addition to the above operations,
an auxiliary routine is introduced which tests whether a NumberT is
zero.

The exponentiation routine is implemented using a repeated squaring
algorithm.

The greatest common divisor routine is implemented using the extended
euclid algorithm.

Both of these algorithms are described in detail in the textbook
"Introduction to Algorithms" by .., .., Rivest.


*/




/* numCompareZero: see if a big number is zero
**
** Return 0 if number is zero
** Return nonzero if number is not zero
*/
int
numCompareZero(const NumberT *num)
{
  int ni;

  for(ni=num->len - 1; ni >= 0; ni--) {
    if(num->units[ni])
      return ni+1;
  }
  return 0;
}




/* Forward declarations.  Would not be needed if this file could be
** merged into to the original arith.c file
*/
void numSubtract(const NumberT *n1,const NumberT *n2,NumberT *dif,const NumberT *mod);
void numMultiply(const NumberT *n1,const NumberT *n2,NumberT *prod,const NumberT *mod);

void numLongDiv(const NumberT *num, const NumberT *den, NumberT *quo, NumberT *rem);




/* numExponentiate: Exponentiate a big number, MODULUS mod
**
** Return result = num ^ exp MODULUS mod
**
** Assumes result->len >= num->len.
** Other assumptions are those which apply to numMultiply,
** in particular any assumptions about num->base ?= mod->base
** and num->len ?= mod->len.
**
** However, this routine NOT assume that num->base == exp->base
** or that num->len == exp->len.
**
** result->base is set to num->base
**
** Algorithm is a version of "raising to powers with repeated squaring"
**
** Algorithm description:
**   The original number is squared a number of times until the
**   effective power-of-two would exceed the original exponent.  Then,
**   the effective power-of-two is subtracted from the original
**   exponent.  The repeatedly squared number is multiplied by the
**   result and stored in the result.  This procedure loops until the
**   exponent is zero.
**
** Algorithm pseudo-code:
**   Set exp_rem <- (original exponent)
**   Set result <- 1
**   Loop: While exp_rem is nonzero,
**     pow_2 <- 2
**     super_sqr <- (original number)
**
**     ; The following loop is the repeated squaring
**     Loop: While pow_2 <= exp_rem
**       super_sqr <- super_sqr * super_sqr MODULUS mod
**       pow_2 <- pow_2 + pow_2
**
**     result *= super_sqr
**     exp_rem <- exp_rem - pow_2 / 2
*/
void
numExponentiate(const NumberT *num, const NumberT *exp,
                const NumberT *mod, NumberT *result)
{
  NumberT exp_rem;    /* exponent remaining */
  NumberT pow_2;      /* Current power of 2 being operated with */
  NumberT super_sqr;  /* Current number raised to power of 2 */
  NumberT quo;        /* temp space */
  NumberT two;        /* the value "2" as a big number */

#if defined(VERBOSE) || defined(DEBUG_EXP)
  printf("num   = "); numPrintContiguous(num);
  printf("exp   = "); numPrintContiguous(exp);
  printf("mod   = "); numPrintContiguous(mod);
#endif

  numAlloc(&exp_rem, exp->len, exp->base);
  numAlloc(&pow_2, exp->len + 1, exp->base);
  numAlloc(&quo, exp->len + 1, exp->base);
  numAlloc(&super_sqr, num->len, num->base);
  numAlloc(&two, 1, exp->base);

  /* Initial values */
  two.units[0] = 2;        /* two <- 2 */
  result->base = num->base;
  numZero(result);
  result->units[0] = 1;    /* result <- 1 */
  numCopy(&exp_rem, exp);  /* exp_rem <- exp */

  while(numCompareZero(&exp_rem)) {
    numZero(&pow_2);
    pow_2.units[0] = 2;        /* pow_2 <- 2 */
    numCopy(&super_sqr, num);  /* super_sqr <- num */

    /* Repeat squaring */
    while(numCompare(&pow_2, &exp_rem) <= 0) {
      debug = &exp_rem;
      numMultiply(&super_sqr, &super_sqr, &super_sqr, mod);
      numAddPrim(&pow_2, &pow_2, &pow_2);
    }

    /* Compute pow /= 2 , exp_rem -= pow */
    numLongDiv(&pow_2, &two, &quo, NULL);
    numCopy(&pow_2, &quo);
    numSubPrim(&exp_rem, &pow_2, &exp_rem);

    /* Compute result *= super_sqr */
    numMultiply(result, &super_sqr, result, mod);
  }

  numFree(&exp_rem);
  numFree(&pow_2);
  numFree(&quo);
  numFree(&super_sqr);
  numFree(&two);
}




/* numEuclidGCD: Find Greatest Common Denominator using extended Euclid
**
** Returns dn, xn, yn such that dn == gcd(an, bn) == an * xn + bn * yn
**
** Incidentally finds the "inverse" of a number, in the modulus "mod" Number
** system, assuming that bn == mod at the top level call.
**
** mod is passed in at the top level as the modulus of the Number
** system.  This value is passed to each recursive call of this routine,
** unchanged.  The purpose is in the use of modular arithmetic routines
** within this routine which would otherwise produce negative
** intermediate results, which are not allowed in this implementation
** of my Number system.
**
** Assumptions are as per the calls to numLongDiv, numMultPrim below.
**
** If dn == 1 (i.e. if the GCD==1), this routine gives the solution to
**   an * xn == 1 (mod b)
** which means that the value (an^(-1) mod bn) exists and is given by xn.
**
** (Note that the number "yn" which is returned might also be an
** inverse in a sense, although it is unlikely to have any utility.
** "yn" is the inverse of "bn" when the modulus is "an".  This routine
** is reciprocal in an, bn in that sense.  However, it is expected that
** the number passed in as "mod" is equal to the modulus of the Number
** system.)
**
** Algorithm:
** This is a variation of the "Extended Euclid" algorith in the textbook
** "Introduction to Algorithms" by .., .., Rivest.
*/
void
numEuclidGCD(const NumberT *an, const NumberT *bn, NumberT *dn,
             NumberT *xn, NumberT *yn,
             const NumberT *mod)
{
  NumberT dp;
  NumberT xp;
  NumberT yp;
  NumberT a_over_b;
  NumberT y_times_a_div_b;

  dn->base = an->base;
  xn->base = an->base;
  yn->base = an->base;

  /* EE1 */
  if(!numCompareZero(bn)) {
    /* EE2 */
    numCopy(dn, an);   /* dn <- an */
    numZero(xn);
    xn->units[0] = 1;  /* xn <- 1 */
    numZero(yn);       /* yn <- 0 */

    return;
  }

  numAlloc(&dp, dn->len, dn->base);
  numAlloc(&xp, xn->len, xn->base);
  numAlloc(&yp, yn->len, yn->base);
  numAlloc(&a_over_b, an->len, an->base);
  numAlloc(&y_times_a_div_b, an->len + yn->len, an->base);

  /* EE3 */
  numLongDiv(an, bn, NULL, &a_over_b);  /* a_over_b <- an MOD bn */
  numEuclidGCD(bn, &a_over_b, &dp, &xp, &yp, mod);

  /* EE4 */
  numCopy(dn, &dp);  /* dn <- dp */
  numCopy(xn, &yp);  /* xn <- yp */
  /* Compute  yn <- (xp - [an/bn] * yp) */
  numLongDiv(an, bn, &a_over_b, NULL);  /* a_over_b <- [an / bn] */
  numMultiply(&a_over_b, &yp, &y_times_a_div_b, mod);
  numSubtract(&xp, &y_times_a_div_b, yn, mod);

  /* Make sure results are modular */
  numLongDiv(dn, mod, NULL, &dp);  /* dp <- dn MOD mod */
  numCopy(dn, &dp);                /* dn <- dp */

  numFree(&dp);
  numFree(&xp);
  numFree(&yp);
  numFree(&a_over_b);
  numFree(&y_times_a_div_b);
}