Quelle  power.c   Sprache: C

/****************************************************************************
**
*A  power.c                     ANUPQ source                   Eamonn O'Brien
**
*Y  Copyright 1995-2001,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
*Y  Copyright 1995-2001,  School of Mathematical Sciences, ANU,     Australia
**
*/

#include "pq_defs.h"
#include "pcp_vars.h"

static void zero_array(int ptr, int length);
static void copy_array(int old, int length, int new);

/* power routine - written by M J Smith, May 1991.

   raise exponent vector with base address cp to power exp;
   the method used depends on the value of the prime P;

   P = 2, 3: For each factor of P in the power the word is raised
   to the power P by multiplication P - 1 times. Any factor of the
   power remaining is done using a P-ary decomposition.

   P > 3: A P-ary decomposition is used. Calculation of the powers
   up to P at each step are done using a binary expansion.

   Storage:
   X - exponent vect contains word on entry, will contain
       on exit of routine. Always accumulating answer.
   A - exponent string - used for multiplying a word a number of
       times in a loop (small primes) or for squaring words in binary.
   Z - In binary method, accumulates the word to the prime p
       exponent. This is used as the base word in the next iteration.
   B - Is used only in binary expansion to square the string A.
       A is unpacked into B, then collected onto B, then packed from B. */

void power(int exp, int cp, struct pcp_vars *pcp)
{
   register int *y = y_address;

   register int p = pcp->p;
   register int lastg = pcp->lastg;
   register int x = cp;
   register int a = pcp->submlg - (lastg + 1);
   register int b = a - (lastg + 1);
   register int z = b - (lastg + 1);
   register int q, r, pp, nn;
   register int i;

   if (exp == 1)
      return;

   /* nn is the exponent requested */
   nn = exp;

   /* first consider small primes */
   if (p == 2 || p == 3) {

      /* extract all powers of the prime from the exponent */
      while (MOD(nn, p) == 0) {
         nn /= p;

         /* pack word in X into string for multiplication p - 1 times */
         vector_to_string(x, a, pcp);
         if (y[a + 1] == 0)
            return;

         /* now multiply p - 1 times to get X^p */
         for (i = 1; i <= p - 1; ++i) {
            collect(-a, x, pcp);
         }
      }

      if (nn == 1)
         return;

      /* have extracted all powers of p from exponent -
         now do rest using prime p expansion */

      /* move X into Z, set X to 1 */
      copy_array(x, lastg, z);
      zero_array(x, lastg);

      while (nn > 0) {
         r = MOD(nn, p);
         nn /= p;

         /* move Z into A to multiply onto Z p - 1 times and
            onto X r times */
         vector_to_string(z, a, pcp);

         /* now calculate Z = Z^p and X = X * Z^r */
         if (y[a + 1] != 0) {
            for (i = 1; i <= p - 1; ++i) {
               if (i <= r)
                  collect(-a, x, pcp);
               collect(-a, z, pcp);
            }
         }
      }
   }

   /* for larger primes, use prime p decomposition and subsequent
      binary expansion */
   else {
      /* move X into Z and set X to 1 */
      vector_to_string(x, z, pcp);
      zero_array(x, lastg);

      while (nn > 0) {

         /* move word w in Z into A, and set Z to 1; A will square each
            iteration, and Z will accumulate some of these powers to
            end up with w^p at end of while loop */

         string_to_vector(z, a, pcp);
         zero_array(z, lastg);

         q = nn / p;
         r = MOD(nn, p);
         pp = p;

         /* Now use binary expansion of both PP (ie p) and remainder R
            to accumulate w^p in Z and w^R onto X from squaring of w.
            Must continue until we have last w^R on X or until we get
            w^p in Z if there is any remaining exponent (ie Q > 0) */

         while (r > 0 || (pp > 0 && q > 0)) {

            /* collect onto answer if needed (ie R = 1) */
            if (MOD(r, 2) == 1) {
               copy_array(a, lastg, b);
               if (y[x + 1] > 0) {
                  collect(-x, b, pcp);
               }
               vector_to_string(b, x, pcp);
            }

            /* collect onto Z for next power of w if next iteration reqd */
            if (MOD(pp, 2) == 1 && q > 0) {
               copy_array(a, lastg, b);
               if (y[z + 1] > 0) {
                  collect(-z, b, pcp);
               }
               vector_to_string(b, z, pcp);
            }

            r = r >> 1;
            pp = pp >> 1;

            /* if powers still needed for answer X or for w^p in Z for
               another iteration (ie Q > 0) then square A by unpacking into
               exponent vector B, collecting A and B, then repacking into A */

            if (r > 0 || (pp > 0 && q > 0)) {
               /* square A */
               vector_to_string(a, b, pcp);
               if (y[b + 1] > 0) {
                  collect(-b, a, pcp);
               }
            }
         }
         nn = q;
      }

      /* now X is the answer as a string, so convert to exponent vector */
      string_to_vector(x, b, pcp);
      copy_array(b, lastg, x);
   }
}

/* zero a section of the array, y */

static void zero_array(int ptr, int length)
{
   register int *y = y_address;

   register int i;

   for (i = 1; i <= length; ++i)
      y[ptr + i] = 0;
}

/* copy a section of the array, y */

static void copy_array(int old, int length, int new)
{
   register int *y = y_address;

   register int i;

   for (i = 1; i <= length; ++i)
      y[new + i] = y[old + i];
}