Quelle cyclotom.c Sprache: C

/****************************************************************************
**
**  This file is part of GAP, a system for computational discrete algebra.
**
**  Copyright of GAP belongs to its developers, whose names are too numerous
**  to list here. Please refer to the COPYRIGHT file for details.
**
**  SPDX-License-Identifier: GPL-2.0-or-later
**
**  This file implements the arithmetic for elements from  cyclotomic  fields
**  $Q(e^{{2 \pi i}/n}) = Q(e_n)$,  which  we  call  cyclotomics  for  short.
**
**  The obvious way to represent cyclotomics is to write them as a polynom in
**  $e_n$, the  primitive <n>th root  of unity.  However,  if we  do  this it
**  happens that various  polynomials actually represent the same cyclotomic,
**  e.g., $2+e_3^2 = -2e_3-e_3^2$.  This  is because, if  viewed  as a vector
**  space over the rationals, $Q(e_n)$ has dimension $\phi(n)$ and not $n$.
**
**  This  is  solved by   taking  a  system of $\phi(n)$  linear  independent
**  vectors, i.e., a base, instead of the $n$ linear dependent roots $e_n^i$,
**  and writing  cyclotomics as linear combinations  in  those  vectors.    A
**  possible base would be the set of $e_n^i$ with $i=0..\phi(n)-1$.  In this
**  representation we have $2+e_3^2 = 1-e_3$.
**
**  However we take a different base.  We take  the set of roots $e_n^i$ such
**  that $i \notin (n/q)*[-(q/p-1)/2..(q/p-1)/2]$    mod $q$, for every   odd
**  prime divisor $p$ of $n$, where $q$ is the maximal power  of  $p$ in $n$,
**  and $i \notin (n/q)*[q/2..q-1]$, if $q$ is the maximal power of 2 in $n$.
**  It is not too difficult to see, that this gives in fact $\phi(n)$ roots.
**
**  For example for $n = 45$ we take the roots $e_{45}^i$ such  that  $i$  is
**  not congruent to $(45/5)*[-(5/5-1)/2..(5/5-1)/2]$ mod $5$, i.e.,   is not
**  divisible by 5, and is not congruent  to $(45/9)[-(9/3-1)/2 .. (9/3-1)/2]
**  = [-5,0,5]$ mod $9$,  i.e.,  $i \in [1,2,3,6,7,8,11,12,16,17,19,21,24,26,
**  28,29,33,34,37,38,39,42,43,44]$.
**
**  This base  has two properties, which make  computing with this base easy.
**  First we can convert an arbitrary polynom in $e_n$ into this base without
**  doing  polynom arithmetic.  This is necessary   for the base $e_n^i$ with
**  $i=0..\phi(n)$, where we have to compute modulo the  cyclotomic  polynom.
**  The algorithm for this is given in the description of 'ConvertToBase'.
**
**  It  follows  from this  algorithm that the  set   of roots is  in  fact a
**  generating system,  and because the set  contains exactly $\phi(n)$ roots
**  it is also linear independent system, so it is in fact a  base.  Actually
**  it is even an integral base, but this is not so easy to prove.
**
**  On the other hand we can test  if  a cyclotomic lies  in fact in a proper
**  cyclotomic subfield of $Q(e_n)$ and if so reduce it into  this field.  So
**  each cyclotomic now has a unique representation in its minimal cyclotomic
**  field, which makes testing for equality  easy.  Also a reduced cyclotomic
**  has less terms  than  the unreduced  cyclotomic,  which makes  arithmetic
**  operations, whose effort depends on the number of terms, cheaper.
**
**  For  odd $n$ this base  is also closed  under complex  conjugation, i.e.,
**  complex conjugation  just permutes the roots of the base  in  this  case.
**  This is not possible if $n$ is even for any  base.  This shows again that
**  2 is the oddest of all primes.
**
**  Better descriptions of the base and  related  topics  can  be  found  in:
**  Matthias Zumbroich,
**  Grundlagen  der  Kreisteilungskoerper  und deren  Implementation in  CAS,
**  Diplomarbeit Mathematik,  Lehrstuhl D für Mathematik, RWTH Aachen,  1989
**
**  We represent a cyclotomic with <d>  terms, i.e., <d> nonzero coefficients
**  in the linear  combination, by a bag  of type 'T_CYC'  with <d>+1 subbags
**  and <d>+1 unsigned integers.  All the bag identifiers are stored at the
**  beginning of the  bag and all unsigned integers are stored  at the end of
**  the bag.
**
**      +-------+-------+-------+-------+- - - -+----+----+----+----+- - -
**      | order | coeff | coeff | coeff |       | un | exp| exp| exp|
**      |       |   1   |   2   |   3   |       |used|  1 |  2 |  3 |
**      +-------+-------+-------+-------+- - - -+----+----+----+----+- - -
**
**  The first subbag is  the order  of  the primitive root of  the cyclotomic
**  field in which the cyclotomic lies.  It is an immediate positive integer,
**  therefore 'INT_INTOBJ( ADDR_OBJ(<cyc>)[ 0 ] )'  gives you the order.  The
**  first unsigned integer is unused (but reserved for future use :-).
**
**  The other subbags and exponents are paired and each pair describes one term.
**  The subbag is the coefficient and the  unsigned short gives the exponent.
**  The coefficient will usually be  an immediate integer,  but could as well
**  be a large integer or even a rational.
**
**  The terms are sorted with respect to the exponent.  Note that none of the
**  arithmetic functions need this, but it makes the equality test simpler.
*/

#include "cyclotom.h"

#include "ariths.h"
#include "bool.h"
#include "error.h"
#include "gapstate.h"
#include "integer.h"
#include "io.h"
#include "lists.h"
#include "modules.h"
#include "opers.h"
#include "plist.h"
#include "saveload.h"

/****************************************************************************
**
*/
static inline UInt SIZE_CYC(Obj cyc)
{
    return SIZE_OBJ(cyc) / (sizeof(Obj)+sizeof(UInt4));
}

static inline Obj * COEFS_CYC(Obj cyc)
{
    return ADDR_OBJ(cyc);
}

static inline const Obj * CONST_COEFS_CYC(Obj cyc)
{
    return CONST_ADDR_OBJ(cyc);
}

static inline UInt4 * EXPOS_CYC(Obj cyc, UInt len)
{
    return (UInt4 *)(ADDR_OBJ(cyc)+(len));
}

static inline const UInt4 * CONST_EXPOS_CYC(Obj cyc, UInt len)
{
    return (const UInt4 *)(CONST_ADDR_OBJ(cyc)+(len));
}

static inline Obj NOF_CYC(Obj cyc)
{
    return CONST_COEFS_CYC(cyc)[0];
}

static inline void SET_NOF_CYC(Obj cyc, Obj val)
{
    COEFS_CYC(cyc)[0] = val;
}

// #define XXX_CYC(cyc,len)        (EXPOS_CYC(cyc,len)[0])

static ModuleStateOffset CycStateOffset = -1;

struct CycModuleState {

/****************************************************************************
**
*V  ResultCyc . . . . . . . . . . . .  temporary buffer for the result, local
**
**  'ResultCyc' is used  by all the arithmetic functions  as a buffer for the
**  result.  Unlike bags of type 'T_CYC' it  stores the cyclotomics unpacked,
**  i.e., 'ADDR_OBJ( ResultCyc )[<i+1>]' is the coefficient of $e_n^i$.
**
**  It is created in 'InitCyc' with room for up to 1000 coefficients  and  is
**  resized when need arises.
*/
Obj ResultCyc;

/****************************************************************************
**
*V  LastECyc  . . . . . . . . . . . .  last constructed primitive root, local
*V  LastNCyc  . . . . . . . . order of last constructed primitive root, local
**
**  'LastECyc'  remembers  the primitive  root that  was last  constructed by
**  'FuncE'.
**
**  'LastNCyc' is the order of this primitive root.
**
**  These values are used in 'FuncE' to avoid constructing the same primitive
**  root over and over again.  This might be expensive,  because  $e_n$  need
**  itself not belong to the base.
**
**  Also these values are used in 'PowCyc' which thereby can recognize if  it
**  is called to compute $e_n^i$ and can then do this easier by just  putting
**  1 at the <i>th place in 'ResultCyc' and then calling 'Cyclotomic'.
*/
Obj  LastECyc;
UInt LastNCyc;

}; // end of struct CycModuleState

extern inline struct CycModuleState *CycState(void)
{
    return (struct CycModuleState *)StateSlotsAtOffset(CycStateOffset);
}

// For convenience and readability
#define ResultCyc   CycState()->ResultCyc
#define LastECyc    CycState()->LastECyc
#define LastNCyc    CycState()->LastNCyc

static void GrowResultCyc(UInt size)
{
    Obj *res;
    UInt i;
    if (ResultCyc == 0) {
        ResultCyc = NEW_PLIST( T_PLIST, size );
        res = BASE_PTR_PLIST(ResultCyc);
        for ( i = 0; i < size; i++ ) { res[i] = INTOBJ_INT(0); }
    } else if ( LEN_PLIST(ResultCyc) < size ) {
        GROW_PLIST( ResultCyc, size );
        SET_LEN_PLIST( ResultCyc, size );
        res = BASE_PTR_PLIST(ResultCyc);
        for ( i = 0; i < size; i++ ) { res[i] = INTOBJ_INT(0); }
    }
}

/****************************************************************************
**
*F  TypeCyc( <cyc> )  . . . . . . . . . . . . . . . . .  type of a cyclotomic
**
**  'TypeCyc' returns the type of a cyclotomic.
**
**  'TypeCyc' is the function in 'TypeObjFuncs' for cyclotomics.
*/
static Obj TYPE_CYC;

static Obj TypeCyc(Obj cyc)
{
    return TYPE_CYC;
}

/****************************************************************************
**
*F  PrintCyc( <cyc> ) . . . . . . . . . . . . . . . . . .  print a cyclotomic
**
**  'PrintCyc' prints the cyclotomic <cyc> in the standard form.
**
**  In principle this is very easy, but it is complicated because we  do  not
**  want to print stuff like '+1*', '-1*', 'E(<n>)^0', 'E(<n>)^1, etc.
*/
static void PrintCyc(Obj cyc)
{
    UInt                n;              // order of the field
    UInt                len;            // number of terms
    UInt                i;              // loop variable

    n   = INT_INTOBJ( NOF_CYC(cyc) );
    len = SIZE_CYC(cyc);
    Pr("%>", 0, 0);
    for ( i = 1; i < len; i++ ) {
        // Store value in local variable, as they can change during Pr
        const Obj *   cfs = CONST_COEFS_CYC(cyc);
        const UInt4 * exs = CONST_EXPOS_CYC(cyc, len);
        Obj           cfsi = cfs[i];
        UInt4         exsi = exs[i];

        if (cfsi == INTOBJ_INT(1) && exsi == 0)
            Pr("1", 0, 0);
        else if (cfsi == INTOBJ_INT(1) && exsi == 1 && i == 1)
            Pr("%>E(%d%<)", n, 0);
        else if (cfsi == INTOBJ_INT(1) && exsi == 1)
            Pr("%>+E(%d%<)", n, 0);
        else if (cfsi == INTOBJ_INT(1) && i == 1)
            Pr("%>E(%d)%>^%2<%d", n, (Int)exsi);
        else if (cfsi == INTOBJ_INT(1))
            Pr("%>+E(%d)%>^%2<%d", n, (Int)exsi);
        else if (LT(INTOBJ_INT(0), cfsi) && exsi == 0)
            PrintObj(cfsi);
        else if (LT(INTOBJ_INT(0), cfsi) && exsi == 1 && i == 1) {
            Pr("%>", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%, n, 0);
        }
        else if (LT(INTOBJ_INT(0), cfsi) && exsi == 1) {
            Pr("%>+", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%, n, 0);
        }
        else if (LT(INTOBJ_INT(0), cfsi) && i == 1) {
            Pr("%>", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%^%2<%d", n, (Int)exsi);
        }
        else if (LT(INTOBJ_INT(0), cfsi)) {
            Pr("%>+", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%^%2<%d", n, (Int)exsi);
        }
        else if (cfsi == INTOBJ_INT(-1) && exsi == 0)
            Pr("%>-%<1", 0, 0);
        else if (cfsi == INTOBJ_INT(-1) && exsi == 1)
            Pr("%>-E(%d%<)", n, 0);
        else if (cfsi == INTOBJ_INT(-1))
            Pr("%>-E(%d)%>^%2<%d", n, (Int)exsi);
        else if (exsi == 0)
            PrintObj(cfsi);
        else if (exsi == 1) {
            Pr("%>", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%, n, 0);
        }
        else {
            Pr("%>", 0, 0);
            PrintObj(cfsi);
            Pr("%>*%^%2<%d", n, (Int)exsi);
        }
    }
    Pr("%<", 0, 0);
}

/****************************************************************************
**
*F  EqCyc( <opL>, <opR> ) . . . . . . . . . test if two cyclotomics are equal
**
**  'EqCyc' returns 'true' if the two cyclotomics <opL>  and <opR>  are equal
**  and 'false' otherwise.
**
**  'EqCyc'  is  pretty  simple because   every    cyclotomic  has a   unique
**  representation, so we just have to compare the terms.
*/
static Int EqCyc(Obj opL, Obj opR)
{
    UInt                len;            // number of terms
    const Obj *         cfl;            // ptr to coeffs of left operand
    const UInt4 *       exl;            // ptr to expnts of left operand
    const Obj *         cfr;            // ptr to coeffs of right operand
    const UInt4 *       exr;            // ptr to expnts of right operand
    UInt                i;              // loop variable

    // compare the order of both fields
    if ( NOF_CYC(opL) != NOF_CYC(opR) )
        return 0;

    // compare the number of terms
    if ( SIZE_CYC(opL) != SIZE_CYC(opR) )
        return 0;

    // compare the cyclotomics termwise
    len = SIZE_CYC(opL);
    cfl = CONST_COEFS_CYC(opL);
    cfr = CONST_COEFS_CYC(opR);
    exl = CONST_EXPOS_CYC(opL,len);
    exr = CONST_EXPOS_CYC(opR,len);
    for ( i = 1; i < len; i++ ) {
        if ( exl[i] != exr[i] )
            return 0;
        else if ( ! EQ(cfl[i],cfr[i]) )
            return 0;
    }

    // all terms are equal
    return 1;
}

/****************************************************************************
**
*F  LtCyc( <opL>, <opR> ) . . . . test if one cyclotomic is less than another
**
**  'LtCyc'  returns  'true'  if  the  cyclotomic  <opL>  is  less  than  the
**  cyclotomic <opR> and 'false' otherwise.
**
**  Cyclotomics are first sorted according to the order of the primitive root
**  they are written in.  That means that the rationals  are  smallest,  then
**  come cyclotomics from $Q(e_3)$ followed by cyclotomics from $Q(e_4)$ etc.
**  Cyclotomics from the same field are sorted lexicographicaly with  respect
**  to their representation in the base of this field.  That means  that  the
**  cyclotomic with smaller coefficient for the first base root  is  smaller,
**  for cyclotomics with the same first coefficient the second decides  which
**  is smaller, etc.
**
**  'LtCyc'  is  pretty  simple because   every    cyclotomic  has a   unique
**  representation, so we just have to compare the terms.
*/
static Int LtCyc(Obj opL, Obj opR)
{
    UInt                lel;            // nr of terms of left operand
    const Obj *         cfl;            // ptr to coeffs of left operand
    const UInt4 *       exl;            // ptr to expnts of left operand
    UInt                ler;            // nr of terms of right operand
    const Obj *         cfr;            // ptr to coeffs of right operand
    const UInt4 *       exr;            // ptr to expnts of right operand
    UInt                i;              // loop variable

    // compare the order of both fields
    if ( NOF_CYC(opL) != NOF_CYC(opR) ) {
        if ( INT_INTOBJ( NOF_CYC(opL) ) < INT_INTOBJ( NOF_CYC(opR) ) )
            return 1;
        else
            return 0;
    }

    // compare the cyclotomics termwise
    lel = SIZE_CYC(opL);
    ler = SIZE_CYC(opR);
    cfl = CONST_COEFS_CYC(opL);
    cfr = CONST_COEFS_CYC(opR);
    exl = CONST_EXPOS_CYC(opL,lel);
    exr = CONST_EXPOS_CYC(opR,ler);
    for ( i = 1; i < lel && i < ler; i++ ) {
        if ( exl[i] != exr[i] )
            if ( exl[i] < exr[i] )
                return LT( cfl[i], INTOBJ_INT(0) );
            else
                return LT( INTOBJ_INT(0), cfr[i] );
        else if ( ! EQ(cfl[i],cfr[i]) )
            return LT( cfl[i], cfr[i] );
    }

    // if one cyclotomic has more terms than the other compare it against 0
    if ( lel < ler )
        return LT( INTOBJ_INT(0), cfr[i] );
    else if ( ler < lel )
        return LT( cfl[i], INTOBJ_INT(0) );
    else
        return 0;
}

static Int LtCycYes(Obj opL, Obj opR)
{
    return 1;
}

static Int LtCycNot(Obj opL, Obj opR)
{
    return 0;
}

/****************************************************************************
**
*F  ConvertToBase(<n>)  . . . . . . convert a cyclotomic into the base, local
**
**  'ConvertToBase'  converts the cyclotomic  'ResultCyc' from the cyclotomic
**  field  of <n>th roots of  unity, into the base  form.  This means that it
**  replaces every root $e_n^i$ that does not belong to the  base by a sum of
**  other roots that do.
**
**  Suppose that $c*e_n^i$ appears in 'ResultCyc' but $e_n^i$ does not lie in
**  the base.  This happens  because, for some  prime $p$ dividing $n$,  with
**  maximal power $q$, $i \in (n/q)*[-(q/p-1)/2..(q/p-1)/2]$ mod $q$.
**
**  We take the identity  $1+e_p+e_p^2+..+e_p^{p-1}=0$, write it  using $n$th
**  roots of unity, $0=1+e_n^{n/p}+e_n^{2n/p}+..+e_n^{(p-1)n/p}$ and multiply
**  it  by $e_n^i$,   $0=e_n^i+e_n^{n/p+i}+e_n^{2n/p+i}+..+e_n^{(p-1)n/p+i}$.
**  Now we subtract $c$ times the left hand side from 'ResultCyc'.
**
**  If $p^2$  does not divide  $n$ then the roots  that are  not in the  base
**  because of $p$ are those  whose exponent is divisible  by $p$.  But $n/p$
**  is not  divisible by $p$, so  neither of the exponent $k*n/p+i, k=1..p-1$
**  is divisible by $p$, so those new roots are acceptable w.r.t. $p$.
**
**  A similar argument shows that  the new  roots  are also acceptable w.r.t.
**  $p$ even if $p^2$ divides $n$...
**
**  Note that the new roots might still not lie  in the  case because of some
**  other prime $p2$.  However, because $i = k*n/p+i$ mod $p2$, this can only
**  happen if $e_n^i$ did also not lie in the base because of $p2$.  So if we
**  remove all  roots that lie in  the base because  of $p$, the later steps,
**  which remove the roots that are not in the base because of larger primes,
**  will not add new roots that do not lie in the base because of $p$ again.
**
**  For an example, suppose 'ResultCyc' is $e_{45}+e_{45}^5 =: e+e^5$.  $e^5$
**  does  not lie in the  base  because $5  \in 5*[-1,0,1]$  mod $9$ and also
**  because it is  divisible  by 5.  After  subtracting  $e^5*(1+e_3+e_3^2) =
**  e^5+e^{20}+e^{35}$ from  'ResultCyc' we get $e-e^{20}-e^{35}$.  Those two
**  roots are  still not  in the  base because of  5.  But  after subtracting
**  $-e^{20}*(1+e_5+e_5^2+e_5^3+e_5^4)=-e^{20}-e^{29}-e^{38}-e^2-e^{11}$  and
**  $-e^{35}*(1+e_5+e_5^2+e_5^3+e_5^4)=-e^{35}-e^{44}-e^8-e^{17}-e^{26}$   we
**  get  $e+e^{20}+e^{29}+e^{38}+e^2+e^{11}+e^{35}+e^{44}+e^8+e^{17}+e^{26}$,
**  which contains only roots that lie in the base.
**
**  'ConvertToBase' and 'Cyclotomic' are the functions that  know  about  the
**  structure of the base.  'EqCyc' and 'LtCyc' only need the  property  that
**  the representation of  all  cyclotomic  integers  is  unique.  All  other
**  functions dont even require that cyclotomics  are  written  as  a  linear
**  combination of   linear  independent  roots,  they  would  work  also  if
**  cyclotomic integers were written as polynomials in $e_n$.
**
**  The inner loops in this function have been duplicated to avoid using  the
**  modulo ('%') operator to reduce the exponents into  the  range  $0..n-1$.
**  Those divisions are quite expensive  on  some  processors, e.g., MIPS and
**  SPARC, and they may singlehanded account for 20 percent of the runtime.
*/
static void ConvertToBase(UInt n)
{
    Obj *               res;            // pointer to the result
    UInt                nn;             // copy of n to factorize
    UInt                p, q;           // prime and prime power
    UInt                i, k, l;        // loop variables
    UInt                t;              // temporary holds n+i+(n/p-n/q)/2
    Obj                 sum;            // sum of two coefficients

    // get a pointer to the cyclotomic and a copy of n to factor
    res = BASE_PTR_PLIST(ResultCyc);
    nn  = n;

    // first handle 2
    if ( nn % 2 == 0 ) {
        q = 2;  while ( nn % (2*q) == 0 )  q = 2*q;
        nn = nn / q;

        // get rid of all terms e^{a*q+b*(n/q)} a=0..(n/q)-1 b=q/2..q-1
        for ( i = 0; i < n; i += q ) {
            t = i + (n/q)*(q-1) + n/q;          // end   (n <= t < 2n)
            k = i + (n/q)*(q/2);                // start (0 <= k <= t)
            for ( ; k < n; k += n/q ) {
                if ( res[k] != INTOBJ_INT(0) ) {
                    l = (k + n/2) % n;
                    if ( ! ARE_INTOBJS( res[l], res[k] )
                      || ! DIFF_INTOBJS( sum, res[l], res[k] ) ) {
                        CHANGED_BAG( ResultCyc );
                        sum = DIFF( res[l], res[k] );
                        res = BASE_PTR_PLIST(ResultCyc);
                    }
                    res[l] = sum;
                    res[k] = INTOBJ_INT(0);
                }
            }
            t = t - n;                          // end   (0 <= t <  n)
            k = k - n;                          // cont. (0 <= k     )
            for ( ; k < t; k += n/q ) {
                if ( res[k] != INTOBJ_INT(0) ) {
                    l = (k + n/2) % n;
                    if ( ! ARE_INTOBJS( res[l], res[k] )
                      || ! DIFF_INTOBJS( sum, res[l], res[k] ) ) {
                        CHANGED_BAG( ResultCyc );
                        sum = DIFF( res[l], res[k] );
                        res = BASE_PTR_PLIST(ResultCyc);
                    }
                    res[l] = sum;
                    res[k] = INTOBJ_INT(0);
                }
            }
        }
    }

    // now handle the odd primes
    for ( p = 3; p <= nn; p += 2 ) {
        if ( nn % p != 0 )  continue;
        q = p;  while ( nn % (p*q) == 0 )  q = p*q;
        nn = nn / q;

        // get rid of e^{a*q+b*(n/q)} a=0..(n/q)-1 b=-(q/p-1)/2..(q/p-1)/2
        for ( i = 0; i < n; i += q ) {
            if ( n <= i+(n/p-n/q)/2 ) {
                t = i + (n/p-n/q)/2;    // end   (n   <= t < 2n)
                k = i - (n/p-n/q)/2;    // start (t-n <= k <= t)
            }
            else {
                t = i + (n/p-n/q)/2+n;  // end   (n   <= t < 2n)
                k = i - (n/p-n/q)/2+n;  // start (t-n <= k <= t)
            }
            for ( ; k < n; k += n/q ) {
                if ( res[k] != INTOBJ_INT(0) ) {
                    for ( l = k+n/p; l < k+n; l += n/p ) {
                        if ( ! ARE_INTOBJS( res[l%n], res[k] )
                          || ! DIFF_INTOBJS( sum, res[l%n], res[k] ) ) {
                            CHANGED_BAG( ResultCyc );
                            sum = DIFF( res[l%n], res[k] );
                            res = BASE_PTR_PLIST(ResultCyc);
                        }
                        res[l%n] = sum;
                    }
                    res[k] = INTOBJ_INT(0);
                }
            }
            t = t - n;                  // end   (0   <= t <  n)
            k = k - n;                  // start (0   <= k     )
            for ( ; k <= t; k += n/q ) {
                if ( res[k] != INTOBJ_INT(0) ) {
                    for ( l = k+n/p; l < k+n; l += n/p ) {
                        if ( ! ARE_INTOBJS( res[l%n], res[k] )
                          || ! DIFF_INTOBJS( sum, res[l%n], res[k] ) ) {
                            CHANGED_BAG( ResultCyc );
                            sum = DIFF( res[l%n], res[k] );
                            res = BASE_PTR_PLIST(ResultCyc);
                        }
                        res[l%n] = sum;
                    }
                    res[k] = INTOBJ_INT(0);
                }
            }
        }
    }

    // notify Gasman
    CHANGED_BAG( ResultCyc );
}

/****************************************************************************
**
*F  Cyclotomic(<n>,<m>) . . . . . . . . . . create a packed cyclotomic, local
**
**  'Cyclotomic'    reduces  the cyclotomic   'ResultCyc'   into the smallest
**  possible cyclotomic subfield and returns it in packed form.
**
**  'ResultCyc'  must   also    be already converted      into  the base   by
**  'ConvertToBase'.   <n> must be  the order of the  primitive root in which
**  written.
**
**  <m> must be a divisor of $n$ and  gives a  hint about possible subfields.
**  If a prime $p$ divides <m> then no  reduction into a subfield whose order
**  is $n /  p$  is possible.   In the  arithmetic   functions  you can  take
**  $lcm(n_l,n_r) / gcd(n_l,n_r) = n / gcd(n_l,n_r)$.  If you cannot provide
**  such a hint just pass 1.
**
**  A special case of the  reduction is the case that  the  cyclotomic  is  a
**  rational.  If this is the case 'Cyclotomic' reduces it into the rationals
**  and returns it as a rational.
**
**  After 'Cyclotomic' has  done its work it clears  the 'ResultCyc'  bag, so
**  that it only contains 'INTOBJ_INT(0)'.  Thus the arithmetic functions can
**  use this buffer without clearing it first.
**
**  'ConvertToBase' and 'Cyclotomic' are the functions that  know  about  the
**  structure of the base.  'EqCyc' and 'LtCyc' only need the  property  that
**  the representation of  all  cyclotomic  integers  is  unique.  All  other
**  functions dont even require that cyclotomics  are  written  as  a  linear
**  combination of   linear  independent  roots,  they  would  work  also  if
**  cyclotomic integers were written as polynomials in $e_n$.
*/
static Obj Cyclotomic(UInt n, UInt m)
{
    Obj                 cyc;            // cyclotomic, result
    UInt                len;            // number of terms
    Obj *               cfs;            // pointer to the coefficients
    UInt4 *             exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    UInt                gcd, s, t;      // gcd of the exponents, temporary
    UInt                eql;            // are all coefficients equal?
    Obj                 cof;            // if so this is the coefficient
    UInt                i, k;           // loop variables
    UInt                nn;             // copy of n to factorize
    UInt                p;              // prime factor
    static UInt         lastN;          // remember last n, dont recompute
    static UInt         phi;            // Euler phi(n)
    static BOOL         isSqfree;       // is n squarefree?
    static UInt         nrp;            // number of its prime factors

    // get a pointer to the cyclotomic and a copy of n to factor
    res = BASE_PTR_PLIST(ResultCyc);

    // count the terms and compute the gcd of the exponents with n
    len = 0;
    gcd = n;
    eql = 1;
    cof = 0;
    for ( i = 0; i < n; i++ ) {
        if ( res[i] != INTOBJ_INT(0) ) {
            len++;
            if ( gcd != 1 ) {
                s = i; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }
            }
            if ( eql && cof == 0 )
                cof = res[i];
            else if ( eql && ! EQ(cof,res[i]) )
                eql = 0;
        }
    }

    // if all exps are divisible 1 < k replace $e_n^i$ by $e_{n/k}^{i/k}$
    // this is the only way a prime whose square divides $n$ could reduce
    if ( 1 < gcd ) {
        for ( i = 1; i < n/gcd; i++ ) {
            res[i]     = res[i*gcd];
            res[i*gcd] = INTOBJ_INT(0);
        }
        n = n / gcd;
    }

    // compute $phi(n)$, test if n is squarefree, compute number of primes
    if ( n != lastN ) {
        lastN = n;
        phi = n;  k = n;
        isSqfree = TRUE;
        nrp = 0;
        for ( p = 2; p <= k; p++ ) {
            if ( k % p == 0 ) {
                phi = phi * (p-1) / p;
                if (k % (p * p) == 0)
                    isSqfree = FALSE;
                nrp++;
                while ( k % p == 0 )  k = k / p;
            }
        }
    }

    // if possible reduce into the rationals, clear buffer bag
    if ( len == phi && eql && isSqfree ) {
        for ( i = 0; i < n; i++ )
            res[i] = INTOBJ_INT(0);
        // return as rational $(-1)^{number primes}*{common coefficient}$
        if ( nrp % 2 == 0 )
            res[0] = cof;
        else {
            CHANGED_BAG( ResultCyc );
            Obj negcof = DIFF( INTOBJ_INT(0), cof );
            res = BASE_PTR_PLIST(ResultCyc);
            res[0] = negcof;
        }
        n = 1;
    }
    CHANGED_BAG( ResultCyc );

    // for all primes $p$ try to reduce from $Q(e_n)$ into $Q(e_{n/p})$
    gcd = phi; s = len; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }
    nn = n;
    for ( p = 3; p <= nn && p-1 <= gcd; p += 2 ) {
        if ( nn % p != 0 )  continue;
        nn = nn / p;  while ( nn % p == 0 )  nn = nn / p;

        // if $p$ is not quadratic and the number of terms is divisible
        // $p-1$ and $p$ divides $m$ not then a reduction is possible
        if ( n % (p*p) != 0 && len % (p-1) == 0 && m % p != 0 ) {

            // test that coeffs for expnts congruent mod $n/p$ are equal
            eql = 1;
            for ( i = 0; i < n && eql; i += p ) {
                cof = res[(i+n/p)%n];
                for ( k = i+2*n/p; k < i+n && eql; k += n/p )
                    if ( ! EQ(res[k%n],cof) )
                        eql = 0;
            }

            // if all coeffs for expnts in all classes are equal reduce
            if ( eql ) {

                // replace every sum of $p-1$ terms with expnts congruent
                // to $i*p$ mod $n/p$ by the term with exponent $i*p$
                // is just the inverse transformation of 'ConvertToBase'
                for ( i = 0; i < n; i += p ) {
                    cof = res[(i+n/p)%n];
                    if ( ! IS_INTOBJ(cof)
                      || (cof == INTOBJ_MIN) ) {
                        CHANGED_BAG( ResultCyc );
                        cof = DIFF( INTOBJ_INT(0), cof );
                        res = BASE_PTR_PLIST(ResultCyc);
                        res[i] = cof;
                    }
                    else {
                        res[i] = INTOBJ_INT( - INT_INTOBJ(cof) );
                    }
                    for ( k = i+n/p; k < i+n && eql; k += n/p )
                        res[k%n] = INTOBJ_INT(0);
                }
                len = len / (p-1);
                CHANGED_BAG( ResultCyc );

                // now replace $e_n^{i*p}$ by $e_{n/p}^{i}$
                for ( i = 1; i < n/p; i++ ) {
                    res[i]   = res[i*p];
                    res[i*p] = INTOBJ_INT(0);
                }
                n = n / p;

            }

        }

    }

    // if the cyclotomic is a rational return it as a rational
    if ( n == 1 ) {
        cyc  = res[0];
        res[0] = INTOBJ_INT(0);
    }

    // otherwise copy terms into a new 'T_CYC' bag and clear 'ResultCyc'
    else {
        cyc = NewBag( T_CYC, (len+1)*(sizeof(Obj)+sizeof(UInt4)) );
        cfs = COEFS_CYC(cyc);
        exs = EXPOS_CYC(cyc,len+1);
        cfs[0] = INTOBJ_INT(n);
        exs[0] = 0;
        k = 1;
        res = BASE_PTR_PLIST(ResultCyc);
        for ( i = 0; i < n; i++ ) {
            if ( res[i] != INTOBJ_INT(0) ) {
                cfs[k] = res[i];
                exs[k] = i;
                k++;
                res[i] = INTOBJ_INT(0);
            }
        }
        // 'CHANGED_BAG' not needed for last bag
    }

    return cyc;
}

/****************************************************************************
**
*F  find smallest field size containing CF(nl) and CF(nr)
** Also adjusts the results bag to ensure that it is big enough
** returns the field size n, and sets *ml and *mr to n/nl and n/nr
** respectively.
*/

static UInt4 CyclotomicsLimit = 1000000;

static UInt FindCommonField(UInt nl, UInt nr, UInt *ml, UInt *mr)
{
  UInt n,a,b,c;
  UInt8 n8;

  // get the smallest field that contains both cyclotomics
  // First Euclid's Algorithm for gcd
  if (nl > nr) {
    a = nl;
    b = nr;
  } else {
    a = nr;
    b = nl;
  }
  while (b > 0) {
    c = a % b;
    a = b;
    b = c;
  }
  *ml = nr/a;
  // Compute the result (lcm) in 64 bit
  n8 = (UInt8)nl * ((UInt8)*ml);
  // Check if it is too large for a small int
  if (n8 > INT_INTOBJ_MAX)
    ErrorMayQuit("This computation would require a cyclotomic field too large to be handled", 0, 0);

  // Switch to UInt now we know we can
  n = (UInt)n8;

  // Handle the soft limit
  while (n > CyclotomicsLimit) {
      ErrorReturnVoid(
          "This computation requires a cyclotomic field of degree %d, larger "
          "than the current limit of %d",
          n, (Int)CyclotomicsLimit,
          "You may return after raising the limit with SetCyclotomicsLimit");
  }

  // Finish up
  *mr = n/nr;

  // make sure that the result bag is large enough
  GrowResultCyc(n);
  return n;
}

static Obj FuncSetCyclotomicsLimit(Obj self, Obj newlimit)
{
    UInt ulimit = GetPositiveSmallInt(SELF_NAME, newlimit);

    if (ulimit < CyclotomicsLimit) {
        ErrorMayQuit("SetCyclotomicsLimit: must not be less than "
                     "old limit of %d",
                     CyclotomicsLimit, 0);
    }
#ifdef SYS_IS_64_BIT
    if (ulimit >= ((UInt)1 << 32)) {
        ErrorMayQuit("Cyclotomic field size limit must be less than 2^32", 0,
                     0);
    }
#endif

    CyclotomicsLimit = ulimit;
    return 0;
}

static Obj FuncGetCyclotomicsLimit(Obj self)
{
    return INTOBJ_INT(CyclotomicsLimit);
}

/****************************************************************************
**
*F  SumCyc( <opL>, <opR> )  . . . . . . . . . . . . .  sum of two cyclotomics
**
**  'SumCyc' returns  the  sum  of  the  two  cyclotomics  <opL>  and  <opR>.
**  Either operand may also be an integer or a rational.
**
**  This   function  is lengthy  because  we  try to  use immediate   integer
**  arithmetic if possible to avoid the function call overhead.
*/
static Obj SumCyc(Obj opL, Obj opR)
{
    UInt                nl, nr;         // order of left and right field
    UInt                n;              // order of smallest superfield
    UInt                ml, mr;         // cofactors into the superfield
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    Obj                 sum;            // sum of two coefficients
    UInt                i;              // loop variable

    // take the cyclotomic with less terms as the right operand
    if ( TNUM_OBJ(opL) != T_CYC
      || (TNUM_OBJ(opR) == T_CYC && SIZE_CYC(opL) < SIZE_CYC(opR)) ) {
        sum = opL;  opL = opR;  opR = sum;
    }

    nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));
    nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));

    n = FindCommonField(nl, nr, &ml, &mr);

    // Copy the left operand into the result
    if ( TNUM_OBJ(opL) != T_CYC ) {
        res = BASE_PTR_PLIST(ResultCyc);
        res[0] = opL;
        CHANGED_BAG( ResultCyc );
    }
    else {
        len = SIZE_CYC(opL);
        cfs = CONST_COEFS_CYC(opL);
        exs = CONST_EXPOS_CYC(opL,len);
        res = BASE_PTR_PLIST(ResultCyc);
        if ( ml == 1 ) {
            for ( i = 1; i < len; i++ )
                res[exs[i]] = cfs[i];
        }
        else {
            for ( i = 1; i < len; i++ )
                res[exs[i]*ml] = cfs[i];
        }
        CHANGED_BAG( ResultCyc );
    }

    // add the right operand to the result
    if ( TNUM_OBJ(opR) != T_CYC ) {
        res = BASE_PTR_PLIST(ResultCyc);
        sum = SUM( res[0], opR );
        res = BASE_PTR_PLIST(ResultCyc);
        res[0] = sum;
        CHANGED_BAG( ResultCyc );
    }
    else {
        len = SIZE_CYC(opR);
        cfs = CONST_COEFS_CYC(opR);
        exs = CONST_EXPOS_CYC(opR,len);
        res = BASE_PTR_PLIST(ResultCyc);
        for ( i = 1; i < len; i++ ) {
            if ( ! ARE_INTOBJS( res[exs[i]*mr], cfs[i] )
              || ! SUM_INTOBJS( sum, res[exs[i]*mr], cfs[i] ) ) {
                CHANGED_BAG( ResultCyc );
                sum = SUM( res[exs[i]*mr], cfs[i] );
                cfs = CONST_COEFS_CYC(opR);
                exs = CONST_EXPOS_CYC(opR,len);
                res = BASE_PTR_PLIST(ResultCyc);
            }
            res[exs[i]*mr] = sum;
        }
        CHANGED_BAG( ResultCyc );
    }

    // return the base reduced packed cyclotomic
    if ( nl % ml != 0 || nr % mr != 0 )  ConvertToBase( n );
    return Cyclotomic( n, ml * mr );
}

/****************************************************************************
**
*F  ZeroCyc( <op> ) . . . . . . . . . . . . . . . . . .  zero of a cyclotomic
**
**  'ZeroCyc' returns the additive neutral element of the cyclotomic <op>.
*/
static Obj ZeroCyc(Obj op)
{
    return INTOBJ_INT(0);
}

/****************************************************************************
**
*F  AInvCyc( <op> ) . . . . . . . . . . . .  additive inverse of a cyclotomic
**
**  'AInvCyc' returns the additive inverse element of the cyclotomic <op>.
*/
static Obj AInvCyc(Obj op)
{
    Obj                 res;            // inverse, result
    UInt                len;            // number of terms
    const Obj *         cfs;            // ptr to coeffs of left operand
    const UInt4 *       exs;            // ptr to expnts of left operand
    Obj *               cfp;            // ptr to coeffs of product
    UInt4 *             exp;            // ptr to expnts of product
    UInt                i;              // loop variable
    Obj                 prd;            // product of two coefficients

    // simply invert the coefficients
    res = NewBag( T_CYC, SIZE_CYC(op) * (sizeof(Obj)+sizeof(UInt4)) );
    SET_NOF_CYC(res, NOF_CYC(op));
    len = SIZE_CYC(op);
    cfs = CONST_COEFS_CYC(op);
    exs = CONST_EXPOS_CYC(op,len);
    cfp = COEFS_CYC(res);
    exp = EXPOS_CYC(res,len);
    for ( i = 1; i < len; i++ ) {
        if ( ! IS_INTOBJ( cfs[i] ) || cfs[i] == INTOBJ_MIN ) {
            CHANGED_BAG( res );
            prd = AINV_SAMEMUT(cfs[i]);
            cfs = CONST_COEFS_CYC(op);
            exs = CONST_EXPOS_CYC(op,len);
            cfp = COEFS_CYC(res);
            exp = EXPOS_CYC(res,len);
        }
        else {
            prd = INTOBJ_INT( - INT_INTOBJ(cfs[i]) );
        }
        cfp[i] = prd;
        exp[i] = exs[i];
    }
    CHANGED_BAG( res );

    return res;
}

/****************************************************************************
**
*F  DiffCyc( <opL>, <opR> ) . . . . . . . . . . difference of two cyclotomics
**
**  'DiffCyc' returns the difference of the two cyclotomic <opL>  and  <opR>.
**  Either operand may also be an integer or a rational.
**
**  This   function  is lengthy  because  we  try to  use immediate   integer
**  arithmetic if possible to avoid the function call overhead.
*/
static Obj DiffCyc(Obj opL, Obj opR)
{
    UInt                nl, nr;         // order of left and right field
    UInt                n;              // order of smallest superfield
    UInt                ml, mr;         // cofactors into the superfield
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    Obj                 sum;            // difference of two coefficients
    UInt                i;              // loop variable

    // get the smallest field that contains both cyclotomics
    nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));
    nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));
    n = FindCommonField(nl, nr, &ml, &mr);

    // copy the left operand into the result
    if ( TNUM_OBJ(opL) != T_CYC ) {
        res = BASE_PTR_PLIST(ResultCyc);
        res[0] = opL;
        CHANGED_BAG( ResultCyc );
    }
    else {
        len = SIZE_CYC(opL);
        cfs = CONST_COEFS_CYC(opL);
        exs = CONST_EXPOS_CYC(opL,len);
        res = BASE_PTR_PLIST(ResultCyc);
        if ( ml == 1 ) {
            for ( i = 1; i < len; i++ )
                res[exs[i]] = cfs[i];
        }
        else {
            for ( i = 1; i < len; i++ )
                res[exs[i]*ml] = cfs[i];
        }
        CHANGED_BAG( ResultCyc );
    }

    // subtract the right operand from the result
    if ( TNUM_OBJ(opR) != T_CYC ) {
        res = BASE_PTR_PLIST(ResultCyc);
        sum = DIFF( res[0], opR );
        res = BASE_PTR_PLIST(ResultCyc);
        res[0] = sum;
        CHANGED_BAG( ResultCyc );
    }
    else {
        len = SIZE_CYC(opR);
        cfs = CONST_COEFS_CYC(opR);
        exs = CONST_EXPOS_CYC(opR,len);
        res = BASE_PTR_PLIST(ResultCyc);
        for ( i = 1; i < len; i++ ) {
            if ( ! ARE_INTOBJS( res[exs[i]*mr], cfs[i] )
              || ! DIFF_INTOBJS( sum, res[exs[i]*mr], cfs[i] ) ) {
                CHANGED_BAG( ResultCyc );
                sum = DIFF( res[exs[i]*mr], cfs[i] );
                cfs = CONST_COEFS_CYC(opR);
                exs = CONST_EXPOS_CYC(opR,len);
                res = BASE_PTR_PLIST(ResultCyc);
            }
            res[exs[i]*mr] = sum;
        }
        CHANGED_BAG( ResultCyc );
    }

    // return the base reduced packed cyclotomic
    if ( nl % ml != 0 || nr % mr != 0 )  ConvertToBase( n );
    return Cyclotomic( n, ml * mr );
}

/****************************************************************************
**
*F  ProdCycInt( <opL>, <opR> )  . . .  product of a cyclotomic and an integer
**
**  'ProdCycInt'    returns the product  of a    cyclotomic and  a integer or
**  rational.  Which operand is the cyclotomic and which the integer does not
**  matter.
**
**  This is a special case, because if the integer is not 0, the product will
**  automatically be base reduced.  So we dont need to  call  'ConvertToBase'
**  or 'Reduce' and directly write into a result bag.
**
**  This   function  is lengthy  because  we  try to  use immediate   integer
**  arithmetic if possible to avoid the function call overhead.
*/
static Obj ProdCycInt(Obj opL, Obj opR)
{
    Obj                 hdP;            // product, result
    UInt                len;            // number of terms
    const Obj *         cfs;            // ptr to coeffs of left operand
    const UInt4 *       exs;            // ptr to expnts of left operand
    Obj *               cfp;            // ptr to coeffs of product
    UInt4 *             exp;            // ptr to expnts of product
    UInt                i;              // loop variable
    Obj                 prd;            // product of two coefficients

    // for $rat * rat$ delegate
    if ( TNUM_OBJ(opL) != T_CYC && TNUM_OBJ(opR) != T_CYC ) {
        return PROD( opL, opR );
    }

    // make the right operand the non cyclotomic
    if ( TNUM_OBJ(opL) != T_CYC ) { hdP = opL;  opL = opR;  opR = hdP; }

    // for $cyc * 0$ return 0 and for $cyc * 1$ return $cyc$
    if ( opR == INTOBJ_INT(0) ) {
        hdP = INTOBJ_INT(0);
    }
    else if ( opR == INTOBJ_INT(1) ) {
        hdP = opL;
    }

    // for $cyc * -1$ need no multiplication or division
    else if ( opR == INTOBJ_INT(-1) ) {
        return AInvCyc( opL );
    }

    // for $cyc * small$ use immediate multiplication if possible
    else if ( IS_INTOBJ(opR) ) {
        hdP = NewBag( T_CYC, SIZE_CYC(opL) * (sizeof(Obj)+sizeof(UInt4)) );
        SET_NOF_CYC(hdP, NOF_CYC(opL));
        len = SIZE_CYC(opL);
        cfs = CONST_COEFS_CYC(opL);
        exs = CONST_EXPOS_CYC(opL,len);
        cfp = COEFS_CYC(hdP);
        exp = EXPOS_CYC(hdP,len);
        for ( i = 1; i < len; i++ ) {
            if ( ! IS_INTOBJ( cfs[i] )
              || ! PROD_INTOBJS( prd, cfs[i], opR ) ) {
                CHANGED_BAG( hdP );
                prd = PROD( cfs[i], opR );
                cfs = CONST_COEFS_CYC(opL);
                exs = CONST_EXPOS_CYC(opL,len);
                cfp = COEFS_CYC(hdP);
                exp = EXPOS_CYC(hdP,len);
            }
            cfp[i] = prd;
            exp[i] = exs[i];
        }
        CHANGED_BAG( hdP );
    }

    // otherwise multiply every coefficient
    else {
        len = SIZE_OBJ(opL);
        hdP = NewBag(T_CYC, len);
        memcpy( ADDR_OBJ(hdP), CONST_ADDR_OBJ(opL), len);
        len = SIZE_CYC(opL);
        cfp = COEFS_CYC(hdP);
        for ( i = 1; i < len; i++ ) {
            prd = PROD( cfp[i], opR );
            cfp = COEFS_CYC(hdP);
            cfp[i] = prd;
            CHANGED_BAG( hdP );
        }
    }

    return hdP;
}

/****************************************************************************
**
*F  ProdCyc( <opL>, <opR> ) . . . . . . . . . . .  product of two cyclotomics
**
**  'ProdCyc' returns the product of the two  cyclotomics  <opL>  and  <opR>.
**  Either operand may also be an integer or a rational.
**
**  This   function  is lengthy  because  we  try to  use immediate   integer
**  arithmetic if possible to avoid the function call overhead.
*/
static Obj ProdCyc(Obj opL, Obj opR)
{
    UInt                nl, nr;         // order of left and right field
    UInt                n;              // order of smallest superfield
    UInt                ml, mr;         // cofactors into the superfield
    Obj                 c;              // one coefficient of the left op
    UInt                e;              // one exponent of the left op
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    Obj                 sum;            // sum of two coefficients
    Obj                 prd;            // product of two coefficients
    UInt                i, k;           // loop variable

    // for $rat * cyc$ and $cyc * rat$ delegate
    if ( TNUM_OBJ(opL) != T_CYC || TNUM_OBJ(opR) != T_CYC ) {
        return ProdCycInt( opL, opR );
    }

    // take the cyclotomic with less terms as the right operand
    if ( SIZE_CYC(opL) < SIZE_CYC(opR) ) {
        prd = opL;  opL = opR;  opR = prd;
    }

    // get the smallest field that contains both cyclotomics
    nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));
    nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));
    n = FindCommonField(nl, nr, &ml, &mr);

    // loop over the terms of the right operand
    for ( k = 1; k < SIZE_CYC(opR); k++ ) {
        c = COEFS_CYC(opR)[k];
        e = (mr * CONST_EXPOS_CYC( opR, SIZE_CYC(opR) )[k]) % n;

        // if the coefficient is 1 just add
        if ( c == INTOBJ_INT(1) ) {
            len = SIZE_CYC(opL);
            cfs = CONST_COEFS_CYC(opL);
            exs = CONST_EXPOS_CYC(opL,len);
            res = BASE_PTR_PLIST(ResultCyc);
            for ( i = 1; i < len; i++ ) {
                if ( ! ARE_INTOBJS( res[(e+exs[i]*ml)%n], cfs[i] )
                  || ! SUM_INTOBJS( sum, res[(e+exs[i]*ml)%n], cfs[i] ) ) {
                    CHANGED_BAG( ResultCyc );
                    sum = SUM( res[(e+exs[i]*ml)%n], cfs[i] );
                    cfs = CONST_COEFS_CYC(opL);
                    exs = CONST_EXPOS_CYC(opL,len);
                    res = BASE_PTR_PLIST(ResultCyc);
                }
                res[(e+exs[i]*ml)%n] = sum;
            }
            CHANGED_BAG( ResultCyc );
        }

        // if the coefficient is -1 just subtract
        else if ( c == INTOBJ_INT(-1) ) {
            len = SIZE_CYC(opL);
            cfs = CONST_COEFS_CYC(opL);
            exs = CONST_EXPOS_CYC(opL,len);
            res = BASE_PTR_PLIST(ResultCyc);
            for ( i = 1; i < len; i++ ) {
                if ( ! ARE_INTOBJS( res[(e+exs[i]*ml)%n], cfs[i] )
                  || ! DIFF_INTOBJS( sum, res[(e+exs[i]*ml)%n], cfs[i] ) ) {
                    CHANGED_BAG( ResultCyc );
                    sum = DIFF( res[(e+exs[i]*ml)%n], cfs[i] );
                    cfs = CONST_COEFS_CYC(opL);
                    exs = CONST_EXPOS_CYC(opL,len);
                    res = BASE_PTR_PLIST(ResultCyc);
                }
                res[(e+exs[i]*ml)%n] = sum;
            }
            CHANGED_BAG( ResultCyc );
        }

        // if the coefficient is a small integer use immediate operations
        else if ( IS_INTOBJ(c) ) {
            len = SIZE_CYC(opL);
            cfs = CONST_COEFS_CYC(opL);
            exs = CONST_EXPOS_CYC(opL,len);
            res = BASE_PTR_PLIST(ResultCyc);
            for ( i = 1; i < len; i++ ) {
                if ( ! ARE_INTOBJS( cfs[i], res[(e+exs[i]*ml)%n] )
                  || ! PROD_INTOBJS( prd, cfs[i], c )
                  || ! SUM_INTOBJS( sum, res[(e+exs[i]*ml)%n], prd ) ) {
                    CHANGED_BAG( ResultCyc );
                    prd = PROD( cfs[i], c );
                    exs = CONST_EXPOS_CYC(opL,len);
                    res = BASE_PTR_PLIST(ResultCyc);
                    sum = SUM( res[(e+exs[i]*ml)%n], prd );
                    cfs = CONST_COEFS_CYC(opL);
                    exs = CONST_EXPOS_CYC(opL,len);
                    res = BASE_PTR_PLIST(ResultCyc);
                }
                res[(e+exs[i]*ml)%n] = sum;
            }
            CHANGED_BAG( ResultCyc );
        }

        // otherwise do it the normal way
        else {
            len = SIZE_CYC(opL);
            for ( i = 1; i < len; i++ ) {
                CHANGED_BAG( ResultCyc );
                cfs = CONST_COEFS_CYC(opL);
                prd = PROD( cfs[i], c );
                exs = CONST_EXPOS_CYC(opL,len);
                res = BASE_PTR_PLIST(ResultCyc);
                sum = SUM( res[(e+exs[i]*ml)%n], prd );
                exs = CONST_EXPOS_CYC(opL,len);
                res = BASE_PTR_PLIST(ResultCyc);
                res[(e+exs[i]*ml)%n] = sum;
            }
            CHANGED_BAG( ResultCyc );
        }

    }

    // return the base reduced packed cyclotomic
    ConvertToBase( n );
    return Cyclotomic( n, ml * mr );
}

/****************************************************************************
**
*F  OneCyc( <op> )  . . . . . . . . . . . . . . . . . . . one of a cyclotomic
**
**  'OneCyc'  returns  the multiplicative neutral  element  of the cyclotomic
**  <op>.
*/
static Obj OneCyc(Obj op)
{
    return INTOBJ_INT(1);
}

/****************************************************************************
**
*F  InvCyc( <op> )  . . . . . . . . . . . . . . . . . inverse of a cyclotomic
**
**  'InvCyc' returns the  multiplicative  inverse element  of  the cyclotomic
**  <op>.
**
**  'InvCyc' computes the  inverse of <op> by computing  the product $prd$ of
**  nontrivial Galois conjugates of <op>.  Then $op * (prd / (op * prd)) = 1$
**  so $prd / (op  * prd)$ is the  inverse of $op$.  Because the  denominator
**  $op *  prd$ is the norm  of $op$ over the rationals  it is rational so we
**  can compute the quotient $prd / (op * prd)$ with 'ProdCycInt'.
*T better multiply only the *different* conjugates?
*/
static Obj InvCyc(Obj op)
{
    Obj                 prd;            // product of conjugates
    UInt                n;              // order of the field
    UInt                sqr;            // if n < sqr*sqr n is squarefree
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    UInt                i, k;           // loop variable
    UInt                gcd, s, t;      // gcd of i and n, temporaries

    // get the order of the field, test if it is squarefree
    n = INT_INTOBJ( NOF_CYC(op) );
    for ( sqr = 2; sqr*sqr <= n && n % (sqr*sqr) != 0; sqr++ )
        ;

    // compute the product of all nontrivial galois conjugates of <opL>
    len = SIZE_CYC(op);
    prd = INTOBJ_INT(1);
    for ( i = 2; i < n; i++ ) {

        // if i gives a galois automorphism apply it
        gcd = n; s = i; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }
        if ( gcd == 1 ) {

            // permute the terms
            cfs = CONST_COEFS_CYC(op);
            exs = CONST_EXPOS_CYC(op,len);
            res = BASE_PTR_PLIST(ResultCyc);
            for ( k = 1; k < len; k++ )
                res[(i*exs[k])%n] = cfs[k];
            CHANGED_BAG( ResultCyc );

            // if n is squarefree conversion and reduction are unnecessary
            if ( n < sqr*sqr ) {
                prd = ProdCyc( prd, Cyclotomic( n, n ) );
            }
            else {
                ConvertToBase( n );
                prd = ProdCyc( prd, Cyclotomic( n, 1 ) );
            }

        }

    }

    // the inverse is the product divided by the norm
    return ProdCycInt( prd, INV( ProdCyc( op, prd ) ) );
}

/****************************************************************************
**
*F  PowCyc( <opL>, <opR> )  . . . . . . . . . . . . . . power of a cyclotomic
**
**  'PowCyc' returns the <opR>th, which must be  an  integer,  power  of  the
**  cyclotomic <opL>.  The left operand may also be an integer or a rational.
*/
static Obj PowCyc(Obj opL, Obj opR)
{
    Obj                 pow;            // power (result)
    Int                 exp;            // exponent (right operand)
    Int                 n;              // order of the field
    UInt                i;              // exponent of left operand

    // get the exponent
    exp = INT_INTOBJ(opR);

    // for $cyc^0$ return 1, for $cyc^1$ return cyc, for $rat^exp$ delegate
    if ( exp == 0 ) {
        pow = INTOBJ_INT(1);
    }
    else if ( exp == 1 ) {
        pow = opL;
    }
    else if ( TNUM_OBJ(opL) != T_CYC ) {
        pow = PowInt( opL, opR );
    }

    // for $e_n^exp$ just put a 1 at the <exp>th position and convert
    else if ( opL == LastECyc ) {
        n = LastNCyc;
        exp = (exp % n + n) % n;
        SET_ELM_PLIST( ResultCyc, exp + 1, INTOBJ_INT(1) );
        CHANGED_BAG( ResultCyc );
        ConvertToBase( LastNCyc );
        pow = Cyclotomic( LastNCyc, 1 );
    }

    // for $(c*e_n^i)^exp$ if $e_n^i$ belongs to the base put 1 at $i*exp$
    else if ( SIZE_CYC(opL) == 2 ) {
        n = INT_INTOBJ( NOF_CYC(opL) );
        pow = POW( CONST_COEFS_CYC(opL)[1], opR );
        i = CONST_EXPOS_CYC(opL,2)[1];
        exp = ((exp*(Int)i) % n + n) % n;
        SET_ELM_PLIST( ResultCyc, exp + 1, pow );
        CHANGED_BAG( ResultCyc );
        ConvertToBase( n );
        pow = Cyclotomic( n, 1 );
    }

    // otherwise compute the power with repeated squaring
    else {

        // if necessary invert the cyclotomic
        if ( exp < 0 ) {
            opL = InvCyc( opL );
            exp = -exp;
        }

        // compute the power using repeated squaring
        pow = INTOBJ_INT(1);
        while ( exp != 0 ) {
            if ( exp % 2 == 1 )  pow = ProdCyc( pow, opL );
            if ( exp     >  1 )  opL = ProdCyc( opL, opL );
            exp = exp / 2;
        }

    }

    return pow;
}

/****************************************************************************
**
*F  FuncE( <self>, <n> )  . . . . . . . . . . . . create a new primitive root
**
**  'FuncE' implements the internal function 'E'.
**
**  'E( <n> )'
**
**  'E'  returns a  primitive  root of order <n>, which must  be  a  positive
**  integer, represented as cyclotomic.
*/
static Obj EOper;

static Obj FuncE(Obj self, Obj n)
{
    Obj *               res;            // pointer into result bag

    // do full operation
    if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(n) ) {
        return DoOperation1Args( self, n );
    }

    GetPositiveSmallInt("E", n);

    // for $e_1$ return 1 and for $e_2$ return -1
    if ( n == INTOBJ_INT(1) )
        return INTOBJ_INT(1);
    else if ( n == INTOBJ_INT(2) )
        return INTOBJ_INT(-1);

    // if the root is not known already construct it
    if ( LastNCyc != INT_INTOBJ(n) ) {
        LastNCyc = INT_INTOBJ(n);
        GrowResultCyc(LastNCyc);
        res = BASE_PTR_PLIST(ResultCyc);
        res[1] = INTOBJ_INT(1);
        CHANGED_BAG( ResultCyc );
        ConvertToBase( LastNCyc );
        LastECyc = Cyclotomic( LastNCyc, 1 );
    }

    // return the root
    return LastECyc;
}

/****************************************************************************
**
*F  FiltIS_CYC( <self>, <val> ) . . . . . .  test if an object is a cyclomtic
**
**  'FiltIS_CYC' implements the internal function 'IsCyc'.
**
**  'IsCyc( <val> )'
**
**  'IsCyc' returns 'true'  if the value <val> is   a cyclotomic and  'false'
**  otherwise.
*/
static Obj IsCycFilt;

static Obj FiltIS_CYC(Obj self, Obj val)
{
    // return 'true' if <obj> is a cyclotomic and 'false' otherwise
    if (IS_CYC(val))
        return True;
    else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
        return False;
    }
    else {
        return DoFilter( self, val );
    }
}

/****************************************************************************
**
*F  FuncIS_CYC_INT( <self>, <val> )  test if an object is a cyclomtic integer
**
**  'FuncIS_CYC_INT' implements the internal function 'IsCycInt'.
**
**  'IsCycInt( <val> )'
**
**  'IsCycInt' returns 'true' if the value  <val> is a cyclotomic integer and
**  'false' otherwise.
**
**  'IsCycInt' relies on the fact that the base is an integral base.
*/
static Obj IsCycIntOper;

static Obj FuncIS_CYC_INT(Obj self, Obj val)
{
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    UInt                i;              // loop variable

    // return 'true' if <obj> is a cyclotomic integer and 'false' otherwise
    if ( IS_INT(val) ) {
        return True;
    }
    else if ( TNUM_OBJ(val) == T_RAT ) {
        return False;
    }
    else if ( TNUM_OBJ(val) == T_CYC ) {
        len = SIZE_CYC(val);
        cfs = COEFS_CYC(val);
        for ( i = 1; i < len; i++ ) {
            if ( TNUM_OBJ(cfs[i]) == T_RAT )
                return False;
        }
        return True;
    }
    else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
        return False;
    }
    else {
        return DoOperation1Args( self, val );
    }
}

/****************************************************************************
**
*F  AttrCONDUCTOR( <self>, <cyc> )  . . . . . . . . . . . . N of a cyclotomic
**
**  'AttrCONDUCTOR' implements the internal function 'Conductor'.
**
**  'Conductor( <cyc> )'
**
**  'Conductor' returns the N of the cyclotomic <cyc>, i.e., the order of the
**  roots of which <cyc> is written as a linear combination.
*/
static Obj ConductorAttr;

static Obj AttrCONDUCTOR(Obj self, Obj cyc)
{
    UInt                n;              // N of the cyclotomic, result
    UInt                m;              // N of element of the list
    UInt                gcd, s, t;      // gcd of n and m, temporaries
    Obj                 list;           // list of cyclotomics
    UInt                i;              // loop variable

    // do full operation
    if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(cyc) ) {
        return DoAttribute( ConductorAttr, cyc );
    }

    if (!IS_CYC(cyc) && !IS_SMALL_LIST(cyc)) {
        RequireArgument(SELF_NAME, cyc,
                        "must be a cyclotomic or a small list");
    }

    // handle cyclotomics
    if ( IS_INT(cyc) || TNUM_OBJ(cyc) == T_RAT ) {
        n = 1;
    }
    else if ( TNUM_OBJ(cyc) == T_CYC ) {
        n = INT_INTOBJ( NOF_CYC(cyc) );
    }

    // handle a list by computing the lcm of the entries
    else {
        list = cyc;
        n = 1;
        for ( i = 1; i <= LEN_LIST( list ); i++ ) {
            cyc = ELMV_LIST( list, i );
            if (!IS_INT(cyc) && TNUM_OBJ(cyc) != T_RAT &&
                TNUM_OBJ(cyc) != T_CYC) {
                ErrorMayQuit(
                    "Conductor: [%d] must be a cyclotomic (not a %s)",
                    (Int)i, (Int)TNAM_OBJ(cyc));
            }
            if ( IS_INT(cyc) || TNUM_OBJ(cyc) == T_RAT ) {
                m = 1;
            }
            else /* if ( TNUM_OBJ(cyc) == T_CYC ) */ {
                m = INT_INTOBJ( NOF_CYC(cyc) );
            }
            gcd = n; s = m; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }
            n = n / gcd * m;
        }
    }

    // return the N of the cyclotomic
    return INTOBJ_INT( n );
}

/****************************************************************************
**
*F  FuncCOEFFS_CYC( <self>, <cyc> ) . . . . . .  coefficients of a cyclotomic
**
**  'FuncCOEFFS_CYC' implements the internal function 'COEFFSCYC'.
**
**  'COEFFSCYC( <cyc> )'
**
**  'COEFFSCYC' returns a list of the coefficients of the  cyclotomic  <cyc>.
**  The list has length <n> if <n> is the order of the primitive  root  $e_n$
**  of which <cyc> is written as a linear combination.  The <i>th element  of
**  the list is the coefficient of $e_l^{i-1}$.
*/
static Obj CoeffsCycOper;

static Obj FuncCOEFFS_CYC(Obj self, Obj cyc)
{
    Obj                 list;           // list of coefficients, result
    UInt                n;              // order of field
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    UInt                i;              // loop variable

    // do full operation
    if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(cyc) ) {
        return DoOperation1Args( self, cyc );
    }

    if (!IS_CYC(cyc)) {
        RequireArgument(SELF_NAME, cyc, "must be a cyclotomic");
    }

    // if <cyc> is rational just put it in a list of length 1
    if ( IS_INT(cyc) || TNUM_OBJ(cyc) == T_RAT ) {
        list = NewPlistFromArgs(cyc);
        // 'CHANGED_BAG' not needed for last bag
    }

    // otherwise make a list and fill it with zeroes and copy the coeffs
    else {
        n = INT_INTOBJ( NOF_CYC(cyc) );
        list = NEW_PLIST( T_PLIST, n );
        SET_LEN_PLIST( list, n );
        len = SIZE_CYC(cyc);
        cfs = CONST_COEFS_CYC(cyc);
        exs = CONST_EXPOS_CYC(cyc,len);
        for ( i = 1; i <= n; i++ )
            SET_ELM_PLIST( list, i, INTOBJ_INT(0) );
        for ( i = 1; i < len; i++ )
            SET_ELM_PLIST( list, exs[i]+1, cfs[i] );
        // 'CHANGED_BAG' not needed for last bag
    }

    return list;
}

/****************************************************************************
**
*F  FuncGALOIS_CYC( <self>, <cyc>, <ord> )   image of a cyc. under galois aut
**
**  'FuncGALOIS_CYC' implements the internal function 'GaloisCyc'.
**
**  'GaloisCyc( <cyc>, <ord> )'
**
**  'GaloisCyc' computes the image of the cyclotomic <cyc>  under  the galois
**  automorphism given by <ord>, which must be an integer.
**
**  The galois automorphism is the mapping that  takes  $e_n$  to  $e_n^ord$.
**  <ord> may be any integer, of course if it is not relative prime to  $ord$
**  the mapping will not be an automorphism, though still an endomorphism.
*/
static Obj GaloisCycOper;

static Obj FuncGALOIS_CYC(Obj self, Obj cyc, Obj ord)
{
    Obj                 gal;            // galois conjugate, result
    Obj                 sum;            // sum of two coefficients
    Int                 n;              // order of the field
    UInt                sqr;            // if n < sqr*sqr n is squarefree
    Int                 o;              // galois automorphism
    UInt                gcd, s, t;      // gcd of n and ord, temporaries
    UInt                len;            // number of terms
    const Obj *         cfs;            // pointer to the coefficients
    const UInt4 *       exs;            // pointer to the exponents
    Obj *               res;            // pointer to the result
    UInt                i;              // loop variable

    // do full operation for any but standard arguments
    if (!IS_INT(ord) || !IS_CYC(cyc)) {
        return DoOperation2Args( self, cyc, ord );
    }

    // get and check <ord>
    if ( ! IS_INTOBJ(ord) ) {
        ord = MOD( ord, AttrCONDUCTOR( 0, cyc ) );
    }
    o = INT_INTOBJ(ord);

    // every galois automorphism fixes the rationals
    if (TNUM_OBJ(cyc) != T_CYC) {
        return cyc;
    }

    // get the order of the field, test if it squarefree
    n = INT_INTOBJ( NOF_CYC(cyc) );
    for ( sqr = 2; sqr*sqr <= n && n % (sqr*sqr) != 0; sqr++ )
        ;

    // force <ord> into the range 0..n-1, compute the gcd of <ord> and <n>
    o = (o % n + n) % n;
    gcd = n; s = o;  while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }

    // if <ord> = 1 just return <cyc>
    if ( o == 1 ) {
        gal = cyc;
    }

    // if <ord> == 0 compute the sum of the entries
    else if ( o == 0 ) {
        len = SIZE_CYC(cyc);
        cfs = COEFS_CYC(cyc);
        gal = INTOBJ_INT(0);
        for ( i = 1; i < len; i++ ) {
            if ( ! ARE_INTOBJS( gal, cfs[i] )
              || ! SUM_INTOBJS( sum, gal, cfs[i] ) ) {
                sum = SUM( gal, cfs[i] );
                cfs = COEFS_CYC( cyc );
            }
            gal = sum;
        }
    }

    // if <ord> == n/2 compute alternating sum since $(e_n^i)^ord = -1^i$
    else if ( n % 2 == 0  && o == n/2 ) {
        gal = INTOBJ_INT(0);
        len = SIZE_CYC(cyc);
        cfs = CONST_COEFS_CYC(cyc);
        exs = CONST_EXPOS_CYC(cyc,len);
        for ( i = 1; i < len; i++ ) {
            if ( exs[i] % 2 == 1 ) {
                if ( ! ARE_INTOBJS( gal, cfs[i] )
                  || ! DIFF_INTOBJS( sum, gal, cfs[i] ) ) {
                    sum = DIFF( gal, cfs[i] );
                    cfs = CONST_COEFS_CYC(cyc);
                    exs = CONST_EXPOS_CYC(cyc,len);
                }
                gal = sum;
            }
            else {
                if ( ! ARE_INTOBJS( gal, cfs[i] )
                  || ! SUM_INTOBJS( sum, gal, cfs[i] ) ) {
                    sum = SUM( gal, cfs[i] );
                    cfs = CONST_COEFS_CYC(cyc);
                    exs = CONST_EXPOS_CYC(cyc,len);
                }
                gal = sum;
            }
        }
--> --------------------

--> maximum size reached

--> --------------------

quality92%

¤ Dauer der Verarbeitung: 0.15 Sekunden (vorverarbeitet) ¤

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.