Quelle collect.c Sprache: C

/****************************************************************************
**
*A  collect.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"
#include "constants.h"
#define MAXEXP (1 << 30)

void stack_overflow(void);
void integer_overflow(void);
void add_string(int string,
                int length,
                int exponent,
                int collected_part,
                struct pcp_vars *pcp);

/* an exponent vector with base address collected_part
   is multiplied on the right by a string with base address
   pointer; the string is assumed to be normal

   this routine follows the algorithm devised by
   M.R. Vaughan-Lee for collection from the left;
   step numbers correspond to step numbers in
   "Collection from the Left" by M.R. Vaughan-Lee,
   J. Symbolic Computation (1990) 9, 725-733

   this version embodies a simpler form of combinatorial
   collection as described on page 733, which lessens the
   chance of integer overflow; it also incorporates recursion
   when adding strings and generators to the exponent vector

   during the collection algorithm, pointers to a number of strings
   are stored on a stack; theoretically, the stack depth could get
   as large as (p - 1) * pcp->lastg * pcp->cc, but, in practice, depths
   this large do not arise; however, it could be necessary to increase
   the dimensions of the stack arrays for some calculations with large
   groups; STACK_SIZE is the maximum depth of the stack; this constant
   is declared in the constants header file; overflow is tested in
   this algorithm

   integer overflow is also tested for, although it can only arise
   if p^3 > 2^31; if p^2 > 2^31, integer overflow can arise which
   is not tested for

   if either overflow arises, a message is printed out and
   the program terminates;

   #####################################################################

   note the sections in this code enclosed within

   #ifndef EXPONENT_P
   #endif

   if you insert as part of the header material of this file

   #define EXPONENT_P

   or supply EXPONENT_P as a -D flag to the compiler, then the
   resulting collector assumes that all generators of the pcp
   have order p and does not execute these portions of the code */

/* stack size on Apollo causes problems */

#ifdef APOLLO
static t_int strstk[STACK_SIZE]; /* string stack */
static t_int lenstk[STACK_SIZE]; /* length stack */
static t_int expstk[STACK_SIZE]; /* exponent stack */
#endif

/* if Lie program, use different collect routine */
#if defined(GROUP)

void collect(int pointer, int collected_part, struct pcp_vars *pcp)
{
   register int *y = y_address;

   register int p1;          /* string pointer */
   register int ce;          /* collected exponent */
   register int cg;          /* collected generator */
   register int ug;          /* uncollected generator */
   register int ue;          /* uncollected exponent */
   register int sp = 0;      /* stack pointer */
   register int exp;         /* exponent */
   register int len = 1;     /* length */
   register int str;         /* string pointer */
   register int halfwt;      /* last generator with weight <= cc/2 */
   register int thirdwt;     /* last generator with weight <= cc/3 */
   register int weight_diff; /* current class - weight of ug */
   register int entry;       /* value to be inserted in collected part */
   register int firstcg;     /* first collected generator for loop counter */
   register int lastcg;      /* last collected generator for loop counter */
   register int maxexp;      /* max exponent allowed in call to add_string */

   register int cp = collected_part;
   register int class_end = pcp->clend;
   register int current_class = pcp->cc;
   register int prime = pcp->p;
   register int pm1 = pcp->pm1;
   register int p_pcomm = pcp->ppcomm;
   register int p_power = pcp->ppower;
   register int structure = pcp->structure;

#ifndef APOLLO
   int strstk[STACK_SIZE]; /* string stack */
   int lenstk[STACK_SIZE]; /* length stack */
   int expstk[STACK_SIZE]; /* exponent stack */
#endif

   register int i;

#include "access.h"

   /* if prime is 2, use special collector */
   if (prime == 2) {
      collectp2(pointer, collected_part, pcp);
      return;
   }

   /* Step (0) --
      initialize collector */

   if (pointer < 0)
      lenstk[0] = y[-pointer + 1];
   else if (pointer == 0)
      return;

   strstk[0] = pointer;
   expstk[0] = 1;

   maxexp = (MAXEXP / prime) * 2;
   halfwt = y[class_end + current_class / 2];
   thirdwt = y[class_end + current_class / 3];

   /* Step (1) --
      process next word on stack */

   while (sp >= 0) {

      str = strstk[sp];
      exp = expstk[sp];

      /* check if str is a string or a generator */

      if (str < 0) {
         /* we have a genuine string */
         len = lenstk[sp];
         sp--;

         /* get first generator exponent pair from string */
         i = y[-str + 2];
         ug = FIELD2(i);

         /* if ug > halfwt, the string can be added to the
            collected part without creating any commutators */
         if (ug > halfwt) {
            add_string(str, len, exp, cp, pcp);
            continue;
         }

         /* ug <= halfwt and so exp must equal 1; stack remainder of string */
         ue = FIELD1(i);
         if (len != 1) {
            strstk[++sp] = str - 1;
            lenstk[sp] = len - 1;
         }
      } else {
         /* str is a generator */
         ug = str;
         ue = exp;
         sp--;
         /* if ug > halfwt, ug commutes with all higher generators */
         if (ug > halfwt) {
            add_string(ug, 1, ue, cp, pcp);
            continue;
         }
      }

      /* ug <= halfwt; if ug > thirdwt, any commutators arising in
         collecting ug commute with all generators after ug, so ug
         can be collected without stacking up collected part */

      if (ug <= thirdwt) {

         /* Step (2) --
            combinatorial collection;
            scan collected part towards the left;
            bypass generators we know must commute with ug;
            when 2 * WT(cg) + WT(ug) > current_class, all generators
            occurring in [cg, ug] commute with each other;
            [cg^ce, ug] = [cg, ug]^ce;
            if cg1,..., cgk all satisfy this weight condition then
            [cg1 * ... * cgk, ug] = [cg1, ug] ... [cgk, ug] */

         if (ue != 1) {
            /* we only move one ug at a time; stack ug^(ue - 1) */
            if (++sp >= STACK_SIZE)
               stack_overflow();
            strstk[sp] = ug;
            expstk[sp] = ue - 1;
            ue = 1;
         }

         weight_diff = current_class - WT(y[structure + ug]);
         firstcg = y[class_end + weight_diff];
         lastcg = y[class_end + weight_diff / 2];

         /* scan collected part to the left, bypassing generators
            which must commute with ug; the collected part between
            lastcg and firstcg contains a word w; we add in [w, ug] */

         for (cg = firstcg; cg > ug; cg--) {
            ce = y[cp + cg];
            if (ce != 0) {
               /* add [cg, ug]^ce to the collected part */
               p1 = y[p_pcomm + cg] + ug;
               p1 = y[p1];
               if (p1 != 0) {
                  if (cg <= lastcg)
                     break;
                  if (p1 < 0)
                     len = y[-p1 + 1];
                  add_string(p1, len, ce, cp, pcp);
               }
            }
         }

         if (cg == ug) {
            /* we have reached ug position during combinatorial
               collection; add 1 to ug entry of collected part
               without stacking any entries if appropriate */
            if (y[cp + ug] == pm1) {
               if (y[p_power + ug] == 0) {
                  y[cp + ug] = 0;
                  continue;
               }
            } else {
               ++y[cp + ug];
               continue;
            }
         }

         /* we have now added in [w, ug]; stack up
            collected part between firstcg and cg + 1 */
         for (i = firstcg; i > cg; i--) {
            ce = y[cp + i];
            if (ce != 0) {
               y[cp + i] = 0;
               if (++sp >= STACK_SIZE)
                  stack_overflow();
               strstk[sp] = i;
               expstk[sp] = ce;
            }
         }

         /* Step (3) --
            ordinary collection; we have moved ug to the cg position;
            continue scanning to the left */
         for (; cg > ug; cg--) {
            ce = y[cp + cg];
            if (ce != 0) {
               /* zero the cg entry of collected part */
               y[cp + cg] = 0;

               /* get [cg, ug] */
               p1 = y[p_pcomm + cg] + ug;
               p1 = y[p1];
               if (p1 == 0) {
                  /* cg commutes with ug so stack cg^ce */
                  if (++sp >= STACK_SIZE)
                     stack_overflow();
                  strstk[sp] = cg;
                  expstk[sp] = ce;
               } else {
                  /* cg does not commute with ug;
                     we can only move ug past one cg at a time;
                     stack [cg, ug] and then cg a total of ce times */

                  if (sp + ce + ce >= STACK_SIZE)
                     stack_overflow();
                  if (p1 < 0)
                     len = y[-p1 + 1];

                  for (; ce > 0; --ce) {
                     strstk[++sp] = p1;
                     lenstk[sp] = len;
                     expstk[sp] = 1;
                     strstk[++sp] = cg;
                     expstk[sp] = 1;
                  }
               }
            }
         }

         /* we have moved ug to the correct position;
            add 1 to the ug entry of collected part */
         add_string(ug, 1, 1, cp, pcp);
         continue;
      } /* ug <= thirdwt */

      /* Step (6) -- ug > thirdwt;
         move ug^ue past entries in the collected part, adding
         commutators directly to the collected part;
         scan collected part towards the left, bypassing generators
         we know must commute with ug; if the cg position in
         the collected part contains an entry ce then this
         represents cg^ce;  [cg^ce, ug^ue] = [cg, ug]^(ce * ue),
         and we add the ce * ue power of [cg, ug] directly to
         the collected part */

      weight_diff = current_class - WT(y[structure + ug]);
      firstcg = y[class_end + weight_diff];

      for (cg = firstcg; cg > ug; cg--) {
         ce = y[cp + cg];
         if (ce != 0) {
            /* add [cg, ug]^(ce * ue) directly to the collected part */
            p1 = y[p_pcomm + cg];
            p1 = y[p1 + ug];
            if (p1 != 0) {
               exp = ce * ue;
               if (exp > maxexp)
                  integer_overflow();
               if (p1 < 0)
                  len = y[-p1 + 1];
               add_string(p1, len, exp, cp, pcp);
            }
         }
      }

      /* add ue to the ug entry of collected part */
      entry = y[cp + ug] + ue;
      if (entry < prime) {
         y[cp + ug] = entry;
         continue;
      } else {
         y[cp + ug] = entry - prime;
         p1 = y[p_power + ug];
         if (p1 == 0)
            continue;
      }

      /* adding ue to the ug entry has created an entry >= prime;
         we have to stack some of collected part */
      for (cg = firstcg; cg > ug; cg--) {
         ce = y[cp + cg];
         if (ce != 0) {
            /* set entry to zero and stack cg^ce */
            y[cp + cg] = 0;
            if (++sp >= STACK_SIZE)
               stack_overflow();
            strstk[sp] = cg;
            expstk[sp] = ce;
         }
      }

      /* add in ug^p; p1 is a pointer to ug^p */
      if (p1 < 0)
         len = y[-p1 + 1];
      add_string(p1, len, 1, cp, pcp);
      continue;
   }
}

#endif /* GROUP*/

/* add exponent times the string with address string and length length
   directly to the collected part with base address collected_part,
   recursively adding powers as required */

void add_string(int string,
                int length,
                int exponent,
                int collected_part,
                struct pcp_vars *pcp)
{
   register int *y = y_address;

   register int cp = collected_part;
   register int exp = exponent;
   register int len = length;
   register int str = string;
   register int entry;
   register int ug;
   register int ue;
   register int power;

   register int class_begin = pcp->ccbeg;
   register int prime = pcp->p;
   register int p_power = pcp->ppower;

   register int i;
   int lower, upper;
#include "access.h"

   if (str > 0) {
      /* Step (4) --
         we have moved generator str to the correct position;
         add exp to the str entry of the collected part;
         reduce entry modulo p and add a power of str^p
         to the collected part if necessary */

      entry = y[cp + str] + exp;
      y[cp + str] = entry % prime;
      if (str < class_begin) {
         exp = entry / prime;
#ifndef EXPONENT_P
         if (exp != 0) {
            /* we need to recursively add in (str^p)^exp */
            str = y[p_power + str];
            if (str != 0) {
               if (str < 0)
                  len = y[-str + 1];
               add_string(str, len, exp, cp, pcp);
            }
         }
#endif
      }
   } else {
      /* Step (5) --
         add string with base address -str and length len
         directly to the collected part exp times; if this
         creates an entry >= prime we reduce the entry modulo
         prime and add in the appropriate power */

      lower = -str + 2;
      upper = -str + len + 1;

      /* get one generator exponent pair at a time from string */

      for (i = lower; i <= upper; i++) {
         ug = FIELD2(y[i]);
         ue = FIELD1(y[i]) * exp;
         entry = y[cp + ug] + ue;

         /* add ue to ug entry of the collected part and reduce mod p */
         y[cp + ug] = entry % prime;

#ifndef EXPONENT_P
         /* we need to recursively add in (ug^p)^power */
         if (ug < class_begin) {
            power = entry / prime;
            if (power != 0) {
               str = y[p_power + ug];
               if (str != 0) {
                  if (str < 0)
                     len = y[-str + 1];
                  add_string(str, len, power, cp, pcp);
               }
            }
         }
#endif
      }
   }
}

/* stack is not big enough */

void stack_overflow(void)
{
   printf("Stack overflow in collection routine; you should increase\n");
   printf("value of STACK_SIZE in constants.h and recompile.\n");
   exit(FAILURE);
}

/* arithmetic overflow */

void integer_overflow(void)
{
   printf("Arithmetic overflow may occur in collection ");
   printf("routine. Results may be invalid.\n");
   exit(FAILURE);
}

quality97%

¤ Dauer der Verarbeitung: 0.9 Sekunden ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.