Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
| Download
GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: 418346/****************************************************************************1**2*W cyclotom.c GAP source Martin Schönert3**4**5*Y Copyright (C) 1996, Lehrstuhl D für Mathematik, RWTH Aachen, Germany6*Y (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland7*Y Copyright (C) 2002 The GAP Group8**9** This file implements the arithmetic for elements from cyclotomic fields10** $Q(e^{{2 \pi i}/n}) = Q(e_n)$, which we call cyclotomics for short.11**12** The obvious way to represent cyclotomics is to write them as a polynom in13** $e_n$, the primitive <n>th root of unity. However, if we do this it14** happens that various polynomials actually represent the same cyclotomic,15** e.g., $2+e_3^2 = -2e_3-e_3^2$. This is because, if viewed as a vector16** space over the rationals, $Q(e_n)$ has dimension $\phi(n)$ and not $n$.17**18** This is solved by taking a system of $\phi(n)$ linear independent19** vectors, i.e., a base, instead of the $n$ linear dependent roots $e_n^i$,20** and writing cyclotomics as linear combinations in those vectors. A21** possible base would be the set of $e_n^i$ with $i=0..\phi(n)-1$. In this22** representation we have $2+e_3^2 = 1-e_3$.23**24** However we take a different base. We take the set of roots $e_n^i$ such25** that $i \notin (n/q)*[-(q/p-1)/2..(q/p-1)/2]$ mod $q$, for every odd26** prime divisor $p$ of $n$, where $q$ is the maximal power of $p$ in $n$,27** and $i \notin (n/q)*[q/2..q-1]$, if $q$ is the maximal power of 2 in $n$.28** It is not too difficult to see, that this gives in fact $\phi(n)$ roots.29**30** For example for $n = 45$ we take the roots $e_{45}^i$ such that $i$ is31** not congruent to $(45/5)*[-(5/5-1)/2..(5/5-1)/2]$ mod $5$, i.e., is not32** divisable by 5, and is not congruent to $(45/9)[-(9/3-1)/2 .. (9/3-1)/2]33** = [-5,0,5]$ mod $9$, i.e., $i \in [1,2,3,6,7,8,11,12,16,17,19,21,24,26,34** 28,29,33,34,37,38,39,42,43,44]$.35**36** This base has two properties, which make computing with this base easy.37** First we can convert an arbitrary polynom in $e_n$ into this base without38** doing polynom arithmetic. This is necessary for the base $e_n^i$ with39** $i=0..\phi(n)$, where we have to compute modulo the cyclotomic polynom.40** The algorithm for this is given in the description of 'ConvertToBase'.41**42** It follows from this algorithm that the set of roots is in fact a43** generating system, and because the set contains exactly $\phi(n)$ roots44** it is also linear independent system, so it is in fact a base. Actually45** it is even an integral base, but this is not so easy to prove.46**47** On the other hand we can test if a cyclotomic lies in fact in a proper48** cyclotomic subfield of $Q(e_n)$ and if so reduce it into this field. So49** each cyclotomic now has a unique representation in its minimal cyclotomic50** field, which makes testing for equality easy. Also a reduced cyclotomic51** has less terms than the unreduced cyclotomic, which makes arithmetic52** operations, whose effort depends on the number of terms, cheaper.53**54** For odd $n$ this base is also closed under complex conjugation, i.e.,55** complex conjugation just permutes the roots of the base in this case.56** This is not possible if $n$ is even for any base. This shows again that57** 2 is the oddest of all primes.58**59** Better descriptions of the base and related topics can be found in:60** Matthias Zumbroich,61** Grundlagen der Kreisteilungskoerper und deren Implementation in CAS,62** Diplomarbeit Mathematik, Lehrstuhl D für Mathematik, RWTH Aachen, 198963**64** We represent a cyclotomic with <d> terms, i.e., <d> nonzero coefficients65** in the linear combination, by a bag of type 'T_CYC' with <d>+1 subbags66** and <d>+1 unsigned short integers. All the bag identifiers are stored at67** the beginning of the bag and all unsigned short integers are stored at68** the end of the bag.69**70** +-------+-------+-------+-------+- - - -+----+----+----+----+- - -71** | order | coeff | coeff | coeff | | un | exp| exp| exp|72** | | 1 | 2 | 3 | |used| 1 | 2 | 3 |73** +-------+-------+-------+-------+- - - -+----+----+----+----+- - -74**75** The first subbag is the order of the primitive root of the cyclotomic76** field in which the cyclotomic lies. It is an immediate positive integer,77** therefore 'INT_INTOBJ( ADDR_OBJ(<cyc>)[ 0 ] )' gives you the order. The78** first unsigned short integer is unused (but reserved for future use :-).79**80** The other subbags and shorts are paired and each pair describes one term.81** The subbag is the coefficient and the unsigned short gives the exponent.82** The coefficient will usually be an immediate integer, but could as well83** be a large integer or even a rational.84**85** The terms are sorted with respect to the exponent. Note that none of the86** arithmetic functions need this, but it makes the equality test simpler.87**88** Chnaged the exponent size from 2 to 4 bytes to avoid overflows SL, 200889*/90#include "system.h" /* Ints, UInts */919293#include "gasman.h" /* garbage collector */94#include "objects.h" /* objects */95#include "scanner.h" /* scanner */9697#include "gap.h" /* error handling, initialisation */9899#include "gvars.h" /* global variables */100101#include "calls.h" /* generic call mechanism */102#include "opers.h" /* generic operations */103104#include "ariths.h" /* basic arithmetic */105106#include "bool.h" /* booleans */107108#include "integer.h" /* integers */109110#include "cyclotom.h" /* cyclotomics */111112#include "records.h" /* generic records */113#include "precord.h" /* plain records */114115#include "lists.h" /* generic lists */116#include "plist.h" /* plain lists */117#include "string.h" /* strings */118119#include "saveload.h" /* saving and loading */120121#include "code.h"122#include "tls.h"123124/****************************************************************************125**126*/127#define SIZE_CYC(cyc) (SIZE_OBJ(cyc) / (sizeof(Obj)+sizeof(UInt4)))128#define COEFS_CYC(cyc) (ADDR_OBJ(cyc))129#define EXPOS_CYC(cyc,len) ((UInt4*)(ADDR_OBJ(cyc)+(len)))130#define NOF_CYC(cyc) (COEFS_CYC(cyc)[0])131#define XXX_CYC(cyc,len) (EXPOS_CYC(cyc,len)[0])132133134/****************************************************************************135**136137*V ResultCyc . . . . . . . . . . . . temporary buffer for the result, local138**139** 'ResultCyc' is used by all the arithmetic functions as a buffer for the140** result. Unlike bags of type 'T_CYC' it stores the cyclotomics unpacked,141** i.e., 'ADDR_OBJ( ResultCyc )[<i+1>]' is the coefficient of $e_n^i$.142**143** It is created in 'InitCyc' with room for up to 1000 coefficients and is144** resized when need arises.145*/146Obj ResultCyc;147148void GrowResultCyc(UInt size) {149Obj *res;150UInt i;151if ( LEN_PLIST(TLS(ResultCyc)) < size ) {152GROW_PLIST( TLS(ResultCyc), size );153SET_LEN_PLIST( TLS(ResultCyc), size );154res = &(ELM_PLIST( TLS(ResultCyc), 1 ));155for ( i = 0; i < size; i++ ) { res[i] = INTOBJ_INT(0); }156}157}158159/****************************************************************************160**161*V LastECyc . . . . . . . . . . . . last constructed primitive root, local162*V LastNCyc . . . . . . . . order of last constructed primitive root, local163**164** 'LastECyc' remembers the primitive root that was last constructed by165** 'FunE'.166**167** 'LastNCyc' is the order of this primitive root.168**169** These values are used in 'FunE' to avoid constructing the same primitive170** root over and over again. This might be expensive, because $e_n$ need171** itself not belong to the base.172**173** Also these values are used in 'PowCyc' which thereby can recognize if it174** is called to compute $e_n^i$ and can then do this easier by just putting175** 1 at the <i>th place in 'ResultCyc' and then calling 'Cyclotomic'.176*/177Obj LastECyc;178UInt LastNCyc;179180181/****************************************************************************182**183*F TypeCyc( <cyc> ) . . . . . . . . . . . . . . . . . type of a cyclotomic184**185** 'TypeCyc' returns the type of a cyclotomic.186**187** 'TypeCyc' is the function in 'TypeObjFuncs' for cyclotomics.188*/189Obj TYPE_CYC;190191Obj TypeCyc (192Obj cyc )193{194return TYPE_CYC;195}196197198/****************************************************************************199**200*F PrintCyc( <cyc> ) . . . . . . . . . . . . . . . . . . print a cyclotomic201**202** 'PrintCyc' prints the cyclotomic <cyc> in the standard form.203**204** In principle this is very easy, but it is complicated because we do not205** want to print stuff like '+1*', '-1*', 'E(<n>)^0', 'E(<n>)^1, etc.206*/207void PrintCyc (208Obj cyc )209{210UInt n; /* order of the field */211UInt len; /* number of terms */212Obj * cfs; /* pointer to the coefficients */213UInt4 * exs; /* pointer to the exponents */214UInt i; /* loop variable */215216n = INT_INTOBJ( NOF_CYC(cyc) );217len = SIZE_CYC(cyc);218Pr("%>",0L,0L);219for ( i = 1; i < len; i++ ) {220/* get pointers, they can change during Pr */221cfs = COEFS_CYC(cyc);222exs = EXPOS_CYC(cyc,len);223if ( cfs[i]==INTOBJ_INT(1) && exs[i]==0 )224Pr("1",0L,0L);225else if ( cfs[i]==INTOBJ_INT(1) && exs[i]==1 && i==1 )226Pr("%>E(%d%<)",n,0L);227else if ( cfs[i]==INTOBJ_INT(1) && exs[i]==1 )228Pr("%>+E(%d%<)",n,0L);229else if ( cfs[i]==INTOBJ_INT(1) && i==1 )230Pr("%>E(%d)%>^%2<%d",n,(Int)exs[i]);231else if ( cfs[i]==INTOBJ_INT(1) )232Pr("%>+E(%d)%>^%2<%d",n,(Int)exs[i]);233else if ( LT(INTOBJ_INT(0),cfs[i]) && exs[i]==0 )234PrintObj(cfs[i]);235else if ( LT(INTOBJ_INT(0),cfs[i]) && exs[i]==1 && i==1 ) {236Pr("%>",0L,0L); PrintObj(cfs[i]); Pr("%>*%<E(%d%<)",n,0L); }237else if ( LT(INTOBJ_INT(0),cfs[i]) && exs[i]==1 ) {238Pr("%>+",0L,0L); PrintObj(cfs[i]); Pr("%>*%<E(%d%<)",n,0L); }239else if ( LT(INTOBJ_INT(0),cfs[i]) && i==1 ) {240Pr("%>",0L,0L); PrintObj(cfs[i]);241Pr("%>*%<E(%d)%>^%2<%d",n,(Int)exs[i]); }242else if ( LT(INTOBJ_INT(0),cfs[i]) ) {243Pr("%>+",0L,0L); PrintObj(cfs[i]);244Pr("%>*%<E(%d)%>^%2<%d",n,(Int)exs[i]); }245else if ( cfs[i]==INTOBJ_INT(-1) && exs[i]==0 )246Pr("%>-%<1",0L,0L);247else if ( cfs[i]==INTOBJ_INT(-1) && exs[i]==1 )248Pr("%>-E(%d%<)",n,0L);249else if ( cfs[i]==INTOBJ_INT(-1) )250Pr("%>-E(%d)%>^%2<%d",n,(Int)exs[i]);251else if ( exs[i]==0 )252PrintObj(cfs[i]);253else if ( exs[i]==1 ) {254Pr("%>",0L,0L); PrintObj(cfs[i]); Pr("%>*%<E(%d%<)",n,0L); }255else {256Pr("%>",0L,0L); PrintObj(cfs[i]);257Pr("%>*%<E(%d)%>^%2<%d",n,(Int)exs[i]); }258}259Pr("%<",0L,0L);260}261262263/****************************************************************************264**265*F EqCyc( <opL>, <opR> ) . . . . . . . . . test if two cyclotomics are equal266**267** 'EqCyc' returns 'true' if the two cyclotomics <opL> and <opR> are equal268** and 'false' otherwise.269**270** 'EqCyc' is pretty simple because every cyclotomic has a unique271** representation, so we just have to compare the terms.272*/273Int EqCyc (274Obj opL,275Obj opR )276{277UInt len; /* number of terms */278Obj * cfl; /* ptr to coeffs of left operand */279UInt4 * exl; /* ptr to expnts of left operand */280Obj * cfr; /* ptr to coeffs of right operand */281UInt4 * exr; /* ptr to expnts of right operand */282UInt i; /* loop variable */283284/* compare the order of both fields */285if ( NOF_CYC(opL) != NOF_CYC(opR) )286return 0L;287288/* compare the number of terms */289if ( SIZE_CYC(opL) != SIZE_CYC(opR) )290return 0L;291292/* compare the cyclotomics termwise */293len = SIZE_CYC(opL);294cfl = COEFS_CYC(opL);295cfr = COEFS_CYC(opR);296exl = EXPOS_CYC(opL,len);297exr = EXPOS_CYC(opR,len);298for ( i = 1; i < len; i++ ) {299if ( exl[i] != exr[i] )300return 0L;301else if ( ! EQ(cfl[i],cfr[i]) )302return 0L;303}304305/* all terms are equal */306return 1L;307}308309310/****************************************************************************311**312*F LtCyc( <opL>, <opR> ) . . . . test if one cyclotomic is less than another313**314** 'LtCyc' returns 'true' if the cyclotomic <opL> is less than the315** cyclotomic <opR> and 'false' otherwise.316**317** Cyclotomics are first sorted according to the order of the primitive root318** they are written in. That means that the rationals are smallest, then319** come cyclotomics from $Q(e_3)$ followed by cyclotomics from $Q(e_4)$ etc.320** Cyclotomics from the same field are sorted lexicographicaly with respect321** to their representation in the base of this field. That means that the322** cyclotomic with smaller coefficient for the first base root is smaller,323** for cyclotomics with the same first coefficient the second decides which324** is smaller, etc.325**326** 'LtCyc' is pretty simple because every cyclotomic has a unique327** representation, so we just have to compare the terms.328*/329Int LtCyc (330Obj opL,331Obj opR )332{333UInt lel; /* nr of terms of left operand */334Obj * cfl; /* ptr to coeffs of left operand */335UInt4 * exl; /* ptr to expnts of left operand */336UInt ler; /* nr of terms of right operand */337Obj * cfr; /* ptr to coeffs of right operand */338UInt4 * exr; /* ptr to expnts of right operand */339UInt i; /* loop variable */340341/* compare the order of both fields */342if ( NOF_CYC(opL) != NOF_CYC(opR) ) {343if ( INT_INTOBJ( NOF_CYC(opL) ) < INT_INTOBJ( NOF_CYC(opR) ) )344return 1L;345else346return 0L;347}348349/* compare the cyclotomics termwise */350lel = SIZE_CYC(opL);351ler = SIZE_CYC(opR);352cfl = COEFS_CYC(opL);353cfr = COEFS_CYC(opR);354exl = EXPOS_CYC(opL,lel);355exr = EXPOS_CYC(opR,ler);356for ( i = 1; i < lel && i < ler; i++ ) {357if ( exl[i] != exr[i] )358if ( exl[i] < exr[i] )359return LT( cfl[i], INTOBJ_INT(0) );360else361return LT( INTOBJ_INT(0), cfr[i] );362else if ( ! EQ(cfl[i],cfr[i]) )363return LT( cfl[i], cfr[i] );364}365366/* if one cyclotomic has more terms than the other compare it agains 0 */367if ( lel < ler )368return LT( INTOBJ_INT(0), cfr[i] );369else if ( ler < lel )370return LT( cfl[i], INTOBJ_INT(0) );371else372return 0L;373}374375Int LtCycYes (376Obj opL,377Obj opR )378{379return 1L;380}381382Int LtCycNot (383Obj opL,384Obj opR )385{386return 0L;387}388389390/****************************************************************************391**392*F ConvertToBase(<n>) . . . . . . convert a cyclotomic into the base, local393**394** 'ConvertToBase' converts the cyclotomic 'TLS(ResultCyc)' from the cyclotomic395** field of <n>th roots of unity, into the base form. This means that it396** replaces every root $e_n^i$ that does not belong to the base by a sum of397** other roots that do.398**399** Suppose that $c*e_n^i$ appears in 'TLS(ResultCyc)' but $e_n^i$ does not lie in400** the base. This happens because, for some prime $p$ dividing $n$, with401** maximal power $q$, $i \in (n/q)*[-(q/p-1)/2..(q/p-1)/2]$ mod $q$.402**403** We take the identity $1+e_p+e_p^2+..+e_p^{p-1}=0$, write it using $n$th404** roots of unity, $0=1+e_n^{n/p}+e_n^{2n/p}+..+e_n^{(p-1)n/p}$ and multiply405** 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}$.406** Now we subtract $c$ times the left hand side from 'TLS(ResultCyc)'.407**408** If $p^2$ does not divide $n$ then the roots that are not in the base409** because of $p$ are those whose exponent is divisable by $p$. But $n/p$410** is not divisable by $p$, so neither of the exponent $k*n/p+i, k=1..p-1$411** is divisable by $p$, so those new roots are acceptable w.r.t. $p$.412**413** A similar argument shows that the new roots are also acceptable w.r.t.414** $p$ even if $p^2$ divides $n$...415**416** Note that the new roots might still not lie in the case because of some417** other prime $p2$. However, because $i = k*n/p+i$ mod $p2$, this can only418** happen if $e_n^i$ did also not lie in the base because of $p2$. So if we419** remove all roots that lie in the base because of $p$, the later steps,420** which remove the roots that are not in the base because of larger primes,421** will not add new roots that do not lie in the base because of $p$ again.422**423** For an example, suppose 'TLS(ResultCyc)' is $e_{45}+e_{45}^5 =: e+e^5$. $e^5$424** does not lie in the base because $5 \in 5*[-1,0,1]$ mod $9$ and also425** because it is divisable by 5. After subtracting $e^5*(1+e_3+e_3^2) =426** e^5+e^{20}+e^{35}$ from 'TLS(ResultCyc)' we get $e-e^{20}-e^{35}$. Those two427** roots are still not in the base because of 5. But after subtracting428** $-e^{20}*(1+e_5+e_5^2+e_5^3+e_5^4)=-e^{20}-e^{29}-e^{38}-e^2-e^{11}$ and429** $-e^{35}*(1+e_5+e_5^2+e_5^3+e_5^4)=-e^{35}-e^{44}-e^8-e^{17}-e^{26}$ we430** get $e+e^{20}+e^{29}+e^{38}+e^2+e^{11}+e^{35}+e^{44}+e^8+e^{17}+e^{26}$,431** which contains only roots that lie in the base.432**433** 'ConvertToBase' and 'Cyclotomic' are the functions that know about the434** structure of the base. 'EqCyc' and 'LtCyc' only need the property that435** the representation of all cyclotomic integers is unique. All other436** functions dont even require that cyclotomics are written as a linear437** combination of linear independent roots, they would work also if438** cyclotomic integers were written as polynomials in $e_n$.439**440** The inner loops in this function have been duplicated to avoid using the441** modulo ('%') operator to reduce the exponents into the range $0..n-1$.442** Those divisions are quite expensive on some processors, e.g., MIPS and443** SPARC, and they may singlehanded account for 20 percent of the runtime.444*/445void ConvertToBase (446UInt n )447{448Obj * res; /* pointer to the result */449UInt nn; /* copy of n to factorize */450UInt p, q; /* prime and prime power */451UInt i, k, l; /* loop variables */452UInt t; /* temporary holds n+i+(n/p-n/q)/2 */453Obj sum; /* sum of two coefficients */454455/* get a pointer to the cyclotomic and a copy of n to factor */456res = &(ELM_PLIST( TLS(ResultCyc), 1 ));457nn = n;458459/* first handle 2 */460if ( nn % 2 == 0 ) {461q = 2; while ( nn % (2*q) == 0 ) q = 2*q;462nn = nn / q;463464/* get rid of all terms e^{a*q+b*(n/q)} a=0..(n/q)-1 b=q/2..q-1 */465for ( i = 0; i < n; i += q ) {466t = i + (n/q)*(q-1) + n/q; /* end (n <= t < 2n) */467k = i + (n/q)*(q/2); /* start (0 <= k <= t) */468for ( ; k < n; k += n/q ) {469if ( res[k] != INTOBJ_INT(0) ) {470l = (k + n/2) % n;471if ( ! ARE_INTOBJS( res[l], res[k] )472|| ! DIFF_INTOBJS( sum, res[l], res[k] ) ) {473CHANGED_BAG( TLS(ResultCyc) );474sum = DIFF( res[l], res[k] );475res = &(ELM_PLIST( TLS(ResultCyc), 1 ));476}477res[l] = sum;478res[k] = INTOBJ_INT(0);479}480}481t = t - n; /* end (0 <= t < n) */482k = k - n; /* cont. (0 <= k ) */483for ( ; k < t; k += n/q ) {484if ( res[k] != INTOBJ_INT(0) ) {485l = (k + n/2) % n;486if ( ! ARE_INTOBJS( res[l], res[k] )487|| ! DIFF_INTOBJS( sum, res[l], res[k] ) ) {488CHANGED_BAG( TLS(ResultCyc) );489sum = DIFF( res[l], res[k] );490res = &(ELM_PLIST( TLS(ResultCyc), 1 ));491}492res[l] = sum;493res[k] = INTOBJ_INT(0);494}495}496}497}498499/* now handle the odd primes */500for ( p = 3; p <= nn; p += 2 ) {501if ( nn % p != 0 ) continue;502q = p; while ( nn % (p*q) == 0 ) q = p*q;503nn = nn / q;504505/* get rid of e^{a*q+b*(n/q)} a=0..(n/q)-1 b=-(q/p-1)/2..(q/p-1)/2 */506for ( i = 0; i < n; i += q ) {507if ( n <= i+(n/p-n/q)/2 ) {508t = i + (n/p-n/q)/2; /* end (n <= t < 2n) */509k = i - (n/p-n/q)/2; /* start (t-n <= k <= t) */510}511else {512t = i + (n/p-n/q)/2+n; /* end (n <= t < 2n) */513k = i - (n/p-n/q)/2+n; /* start (t-n <= k <= t) */514}515for ( ; k < n; k += n/q ) {516if ( res[k] != INTOBJ_INT(0) ) {517for ( l = k+n/p; l < k+n; l += n/p ) {518if ( ! ARE_INTOBJS( res[l%n], res[k] )519|| ! DIFF_INTOBJS( sum, res[l%n], res[k] ) ) {520CHANGED_BAG( TLS(ResultCyc) );521sum = DIFF( res[l%n], res[k] );522res = &(ELM_PLIST( TLS(ResultCyc), 1 ));523}524res[l%n] = sum;525}526res[k] = INTOBJ_INT(0);527}528}529t = t - n; /* end (0 <= t < n) */530k = k - n; /* start (0 <= k ) */531for ( ; k <= t; k += n/q ) {532if ( res[k] != INTOBJ_INT(0) ) {533for ( l = k+n/p; l < k+n; l += n/p ) {534if ( ! ARE_INTOBJS( res[l%n], res[k] )535|| ! DIFF_INTOBJS( sum, res[l%n], res[k] ) ) {536CHANGED_BAG( TLS(ResultCyc) );537sum = DIFF( res[l%n], res[k] );538res = &(ELM_PLIST( TLS(ResultCyc), 1 ));539}540res[l%n] = sum;541}542res[k] = INTOBJ_INT(0);543}544}545}546}547548/* notify Gasman */549CHANGED_BAG( TLS(ResultCyc) );550}551552553/****************************************************************************554**555*F Cyclotomic(<n>,<m>) . . . . . . . . . . create a packed cyclotomic, local556**557** 'Cyclotomic' reduces the cyclotomic 'TLS(ResultCyc)' into the smallest558** possible cyclotomic subfield and returns it in packed form.559**560** 'TLS(ResultCyc)' must also be already converted into the base by561** 'ConvertToBase'. <n> must be the order of the primitive root in which562** written.563**564** <m> must be a divisor of $n$ and gives a hint about possible subfields.565** If a prime $p$ divides <m> then no reduction into a subfield whose order566** is $n / p$ is possible. In the arithmetic functions you can take567** $lcm(n_l,n_r) / gcd(n_l,n_r) = n / gcd(n_l,n_r)$. If you can not provide568** such a hint just pass 1.569**570** A special case of the reduction is the case that the cyclotomic is a571** rational. If this is the case 'Cyclotomic' reduces it into the rationals572** and returns it as a rational.573**574** After 'Cyclotomic' has done its work it clears the 'TLS(ResultCyc)' bag, so575** that it only contains 'INTOBJ_INT(0)'. Thus the arithmetic functions can576** use this buffer without clearing it first.577**578** 'ConvertToBase' and 'Cyclotomic' are the functions that know about the579** structure of the base. 'EqCyc' and 'LtCyc' only need the property that580** the representation of all cyclotomic integers is unique. All other581** functions dont even require that cyclotomics are written as a linear582** combination of linear independent roots, they would work also if583** cyclotomic integers were written as polynomials in $e_n$.584*/585Obj Cyclotomic (586UInt n,587UInt m )588{589Obj cyc; /* cyclotomic, result */590UInt len; /* number of terms */591Obj * cfs; /* pointer to the coefficients */592UInt4 * exs; /* pointer to the exponents */593Obj * res; /* pointer to the result */594UInt gcd, s, t; /* gcd of the exponents, temporary */595UInt eql; /* are all coefficients equal? */596Obj cof; /* if so this is the coefficient */597UInt i, k; /* loop variables */598UInt nn; /* copy of n to factorize */599UInt p; /* prime factor */600static UInt lastN; /* rember last n, dont recompute: */601static UInt phi; /* Euler phi(n) */602static UInt isSqfree; /* is n squarefree? */603static UInt nrp; /* number of its prime factors */604605/* get a pointer to the cyclotomic and a copy of n to factor */606res = &(ELM_PLIST( TLS(ResultCyc), 1 ));607608/* count the terms and compute the gcd of the exponents with n */609len = 0;610gcd = n;611eql = 1;612cof = 0;613for ( i = 0; i < n; i++ ) {614if ( res[i] != INTOBJ_INT(0) ) {615len++;616if ( gcd != 1 ) {617s = i; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }618}619if ( eql && cof == 0 )620cof = res[i];621else if ( eql && ! EQ(cof,res[i]) )622eql = 0;623}624}625626/* if all exps are divisable 1 < k replace $e_n^i$ by $e_{n/k}^{i/k}$ */627/* this is the only way a prime whose square divides $n$ could reduce */628if ( 1 < gcd ) {629for ( i = 1; i < n/gcd; i++ ) {630res[i] = res[i*gcd];631res[i*gcd] = INTOBJ_INT(0);632}633n = n / gcd;634}635636/* compute $phi(n)$, test if n is squarefree, compute number of primes */637if ( n != lastN ) {638lastN = n;639phi = n; k = n;640isSqfree = 1;641nrp = 0;642for ( p = 2; p <= k; p++ ) {643if ( k % p == 0 ) {644phi = phi * (p-1) / p;645if ( k % (p*p) == 0 ) isSqfree = 0;646nrp++;647while ( k % p == 0 ) k = k / p;648}649}650}651652/* if possible reduce into the rationals, clear buffer bag */653if ( len == phi && eql && isSqfree ) {654for ( i = 0; i < n; i++ )655res[i] = INTOBJ_INT(0);656/* return as rational $(-1)^{number primes}*{common coefficient}$ */657if ( nrp % 2 == 0 )658res[0] = cof;659else {660CHANGED_BAG( TLS(ResultCyc) );661res[0] = DIFF( INTOBJ_INT(0), cof );662}663n = 1;664}665CHANGED_BAG( TLS(ResultCyc) );666667/* for all primes $p$ try to reduce from $Q(e_n)$ into $Q(e_{n/p})$ */668gcd = phi; s = len; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }669nn = n;670for ( p = 3; p <= nn && p-1 <= gcd; p += 2 ) {671if ( nn % p != 0 ) continue;672nn = nn / p; while ( nn % p == 0 ) nn = nn / p;673674/* if $p$ is not quadratic and the number of terms is divisiable */675/* $p-1$ and $p$ divides $m$ not then a reduction is possible */676if ( n % (p*p) != 0 && len % (p-1) == 0 && m % p != 0 ) {677678/* test that coeffs for expnts congruent mod $n/p$ are equal */679eql = 1;680for ( i = 0; i < n && eql; i += p ) {681cof = res[(i+n/p)%n];682for ( k = i+2*n/p; k < i+n && eql; k += n/p )683if ( ! EQ(res[k%n],cof) )684eql = 0;685}686687/* if all coeffs for expnts in all classes are equal reduce */688if ( eql ) {689690/* replace every sum of $p-1$ terms with expnts congruent */691/* to $i*p$ mod $n/p$ by the term with exponent $i*p$ */692/* is just the inverse transformation of 'ConvertToBase' */693for ( i = 0; i < n; i += p ) {694cof = res[(i+n/p)%n];695res[i] = INTOBJ_INT( - INT_INTOBJ(cof) );696if ( ! IS_INTOBJ(cof)697|| (cof == INTOBJ_INT(-(1L<<NR_SMALL_INT_BITS))) ) {698CHANGED_BAG( TLS(ResultCyc) );699cof = DIFF( INTOBJ_INT(0), cof );700res = &(ELM_PLIST( TLS(ResultCyc), 1 ));701res[i] = cof;702}703for ( k = i+n/p; k < i+n && eql; k += n/p )704res[k%n] = INTOBJ_INT(0);705}706len = len / (p-1);707CHANGED_BAG( TLS(ResultCyc) );708709/* now replace $e_n^{i*p}$ by $e_{n/p}^{i}$ */710for ( i = 1; i < n/p; i++ ) {711res[i] = res[i*p];712res[i*p] = INTOBJ_INT(0);713}714n = n / p;715716}717718}719720}721722/* if the cyclotomic is a rational return it as a rational */723if ( n == 1 ) {724cyc = res[0];725res[0] = INTOBJ_INT(0);726}727728/* otherwise copy terms into a new 'T_CYC' bag and clear 'TLS(ResultCyc)' */729else {730cyc = NewBag( T_CYC, (len+1)*(sizeof(Obj)+sizeof(UInt4)) );731cfs = COEFS_CYC(cyc);732exs = EXPOS_CYC(cyc,len+1);733cfs[0] = INTOBJ_INT(n);734exs[0] = 0;735k = 1;736res = &(ELM_PLIST( TLS(ResultCyc), 1 ));737for ( i = 0; i < n; i++ ) {738if ( res[i] != INTOBJ_INT(0) ) {739cfs[k] = res[i];740exs[k] = i;741k++;742res[i] = INTOBJ_INT(0);743}744}745/* 'CHANGED_BAG' not needed for last bag */746}747748/* return the result */749return cyc;750}751752/****************************************************************************753**754*F find smallest field size containing CF(nl) and CF(nr)755** Also adjusts the results bag to ensure that it is big enough756** returns the field size n, and sets *ml and *mr to n/nl and n/nr757** respectively.758*/759760UInt4 CyclotomicsLimit = 1000000;761762static UInt FindCommonField(UInt nl, UInt nr, UInt *ml, UInt *mr)763{764UInt n,a,b,c;765UInt8 n8;766767/* get the smallest field that contains both cyclotomics */768/* First Euclid's Algorithm for gcd */769if (nl > nr) {770a = nl;771b = nr;772} else {773a = nr;774b = nl;775}776while (b > 0) {777c = a % b;778a = b;779b = c;780}781*ml = nr/a;782/* Compute the result (lcm) in 64 bit */783n8 = (UInt8)nl * ((UInt8)*ml);784/* Check if it is too large for a small int */785if (n8 > ((UInt8)(1) << NR_SMALL_INT_BITS))786ErrorMayQuit("This computation would require a cyclotomic field too large to be handled",0L, 0L);787788/* Switch to UInt now we know we can*/789n = (UInt)n8;790791/* Handle the soft limit */792while (n > CyclotomicsLimit) {793ErrorReturnVoid("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");794}795796/* Finish up */797*mr = n/nr;798799/* make sure that the result bag is large enough */800GrowResultCyc(n);801return n;802}803804Obj FuncSetCyclotomicsLimit(Obj self, Obj NewLimit) {805UInt ok;806Int limit;807UInt ulimit;808do {809ok = 1;810if (!IS_INTOBJ(NewLimit)) {811ok = 0;812NewLimit = ErrorReturnObj("Cyclotomic Field size limit must be a small integer, not a %s ",(Int)TNAM_OBJ(NewLimit), 0L, "You can return a new value");813} else {814limit = INT_INTOBJ(NewLimit);815if (limit <= 0) {816ok = 0;817NewLimit = ErrorReturnObj("Cyclotomic Field size limit must be positive",0L, 0L, "You can return a new value");818} else {819ulimit = limit;820if (ulimit < CyclotomicsLimit) {821ok = 0;822NewLimit = ErrorReturnObj("Cyclotomic Field size limit must not be less than old limit of %d",CyclotomicsLimit, 0L, "You can return a new value");823}824#ifdef SYS_IS_64_BIT825else if (ulimit > (1L << 32)) {826ok = 0;827NewLimit = ErrorReturnObj("Cyclotomic field size limit must be less than 2^32", 0L, 0L, "You can return a new value");828}829#endif830831}832}833} while (!ok);834CyclotomicsLimit = ulimit;835return (Obj) 0L;836}837838Obj FuncGetCyclotomicsLimit( Obj self) {839return INTOBJ_INT(CyclotomicsLimit);840}841842/****************************************************************************843**844*F SumCyc( <opL>, <opR> ) . . . . . . . . . . . . . sum of two cyclotomics845**846** 'SumCyc' returns the sum of the two cyclotomics <opL> and <opR>.847** Either operand may also be an integer or a rational.848**849** This function is lengthy because we try to use immediate integer850** arithmetic if possible to avoid the function call overhead.851*/852Obj SumCyc (853Obj opL,854Obj opR )855{856UInt nl, nr; /* order of left and right field */857UInt n; /* order of smallest superfield */858UInt ml, mr; /* cofactors into the superfield */859UInt len; /* number of terms */860Obj * cfs; /* pointer to the coefficients */861UInt4 * exs; /* pointer to the exponents */862Obj * res; /* pointer to the result */863Obj sum; /* sum of two coefficients */864UInt i; /* loop variable */865866/* take the cyclotomic with less terms as the right operand */867if ( TNUM_OBJ(opL) != T_CYC868|| (TNUM_OBJ(opR) == T_CYC && SIZE_CYC(opL) < SIZE_CYC(opR)) ) {869sum = opL; opL = opR; opR = sum;870}871872nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));873nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));874875n = FindCommonField(nl, nr, &ml, &mr);876877/* Copy the left operand into the result */878if ( TNUM_OBJ(opL) != T_CYC ) {879res = &(ELM_PLIST( TLS(ResultCyc), 1 ));880res[0] = opL;881CHANGED_BAG( TLS(ResultCyc) );882}883else {884len = SIZE_CYC(opL);885cfs = COEFS_CYC(opL);886exs = EXPOS_CYC(opL,len);887res = &(ELM_PLIST( TLS(ResultCyc), 1 ));888if ( ml == 1 ) {889for ( i = 1; i < len; i++ )890res[exs[i]] = cfs[i];891}892else {893for ( i = 1; i < len; i++ )894res[exs[i]*ml] = cfs[i];895}896CHANGED_BAG( TLS(ResultCyc) );897}898899/* add the right operand to the result */900if ( TNUM_OBJ(opR) != T_CYC ) {901res = &(ELM_PLIST( TLS(ResultCyc), 1 ));902sum = SUM( res[0], opR );903res = &(ELM_PLIST( TLS(ResultCyc), 1 ));904res[0] = sum;905CHANGED_BAG( TLS(ResultCyc) );906}907else {908len = SIZE_CYC(opR);909cfs = COEFS_CYC(opR);910exs = EXPOS_CYC(opR,len);911res = &(ELM_PLIST( TLS(ResultCyc), 1 ));912for ( i = 1; i < len; i++ ) {913if ( ! ARE_INTOBJS( res[exs[i]*mr], cfs[i] )914|| ! SUM_INTOBJS( sum, res[exs[i]*mr], cfs[i] ) ) {915CHANGED_BAG( TLS(ResultCyc) );916sum = SUM( res[exs[i]*mr], cfs[i] );917cfs = COEFS_CYC(opR);918exs = EXPOS_CYC(opR,len);919res = &(ELM_PLIST( TLS(ResultCyc), 1 ));920}921res[exs[i]*mr] = sum;922}923CHANGED_BAG( TLS(ResultCyc) );924}925926/* return the base reduced packed cyclotomic */927if ( nl % ml != 0 || nr % mr != 0 ) ConvertToBase( n );928return Cyclotomic( n, ml * mr );929}930931932/****************************************************************************933**934*F ZeroCyc( <op> ) . . . . . . . . . . . . . . . . . . zero of a cyclotomic935**936** 'ZeroCyc' returns the additive neutral element of the cyclotomic <op>.937*/938Obj ZeroCyc (939Obj op )940{941return INTOBJ_INT( 0L );942}943944945/****************************************************************************946**947*F AInvCyc( <op> ) . . . . . . . . . . . . additive inverse of a cyclotomic948**949** 'AInvCyc' returns the additive inverse element of the cyclotomic <op>.950*/951Obj AInvCyc (952Obj op )953{954Obj res; /* inverse, result */955UInt len; /* number of terms */956Obj * cfs; /* ptr to coeffs of left operand */957UInt4 * exs; /* ptr to expnts of left operand */958Obj * cfp; /* ptr to coeffs of product */959UInt4 * exp; /* ptr to expnts of product */960UInt i; /* loop variable */961Obj prd; /* product of two coefficients */962963/* simply invert the coefficients */964res = NewBag( T_CYC, SIZE_CYC(op) * (sizeof(Obj)+sizeof(UInt4)) );965NOF_CYC(res) = NOF_CYC(op);966len = SIZE_CYC(op);967cfs = COEFS_CYC(op);968cfp = COEFS_CYC(res);969exs = EXPOS_CYC(op,len);970exp = EXPOS_CYC(res,len);971for ( i = 1; i < len; i++ ) {972prd = INTOBJ_INT( - INT_INTOBJ(cfs[i]) );973if ( ! IS_INTOBJ( cfs[i] ) ||974cfs[i] == INTOBJ_INT(-(1L<<NR_SMALL_INT_BITS)) ) {975CHANGED_BAG( res );976prd = AINV( cfs[i] );977cfs = COEFS_CYC(op);978cfp = COEFS_CYC(res);979exs = EXPOS_CYC(op,len);980exp = EXPOS_CYC(res,len);981}982cfp[i] = prd;983exp[i] = exs[i];984}985CHANGED_BAG( res );986987/* return the result */988return res;989}990991992/****************************************************************************993**994*F DiffCyc( <opL>, <opR> ) . . . . . . . . . . difference of two cyclotomics995**996** 'DiffCyc' returns the difference of the two cyclotomic <opL> and <opR>.997** Either operand may also be an integer or a rational.998**999** This function is lengthy because we try to use immediate integer1000** arithmetic if possible to avoid the function call overhead.1001*/1002Obj DiffCyc (1003Obj opL,1004Obj opR )1005{1006UInt nl, nr; /* order of left and right field */1007UInt n; /* order of smallest superfield */1008UInt ml, mr; /* cofactors into the superfield */1009UInt len; /* number of terms */1010Obj * cfs; /* pointer to the coefficients */1011UInt4 * exs; /* pointer to the exponents */1012Obj * res; /* pointer to the result */1013Obj sum; /* difference of two coefficients */1014UInt i; /* loop variable */10151016/* get the smallest field that contains both cyclotomics */1017nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));1018nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));1019n = FindCommonField(nl, nr, &ml, &mr);10201021/* copy the left operand into the result */1022if ( TNUM_OBJ(opL) != T_CYC ) {1023res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1024res[0] = opL;1025CHANGED_BAG( TLS(ResultCyc) );1026}1027else {1028len = SIZE_CYC(opL);1029cfs = COEFS_CYC(opL);1030exs = EXPOS_CYC(opL,len);1031res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1032if ( ml == 1 ) {1033for ( i = 1; i < len; i++ )1034res[exs[i]] = cfs[i];1035}1036else {1037for ( i = 1; i < len; i++ )1038res[exs[i]*ml] = cfs[i];1039}1040CHANGED_BAG( TLS(ResultCyc) );1041}10421043/* subtract the right operand from the result */1044if ( TNUM_OBJ(opR) != T_CYC ) {1045res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1046sum = DIFF( res[0], opR );1047res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1048res[0] = sum;1049CHANGED_BAG( TLS(ResultCyc) );1050}1051else {1052len = SIZE_CYC(opR);1053cfs = COEFS_CYC(opR);1054exs = EXPOS_CYC(opR,len);1055res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1056for ( i = 1; i < len; i++ ) {1057if ( ! ARE_INTOBJS( res[exs[i]*mr], cfs[i] )1058|| ! DIFF_INTOBJS( sum, res[exs[i]*mr], cfs[i] ) ) {1059CHANGED_BAG( TLS(ResultCyc) );1060sum = DIFF( res[exs[i]*mr], cfs[i] );1061cfs = COEFS_CYC(opR);1062exs = EXPOS_CYC(opR,len);1063res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1064}1065res[exs[i]*mr] = sum;1066}1067CHANGED_BAG( TLS(ResultCyc) );1068}10691070/* return the base reduced packed cyclotomic */1071if ( nl % ml != 0 || nr % mr != 0 ) ConvertToBase( n );1072return Cyclotomic( n, ml * mr );1073}107410751076/****************************************************************************1077**1078*F ProdCycInt( <opL>, <opR> ) . . . product of a cyclotomic and an integer1079**1080** 'ProdCycInt' returns the product of a cyclotomic and a integer or1081** rational. Which operand is the cyclotomic and wich the integer does not1082** matter.1083**1084** This is a special case, because if the integer is not 0, the product will1085** automatically be base reduced. So we dont need to call 'ConvertToBase'1086** or 'Reduce' and directly write into a result bag.1087**1088** This function is lengthy because we try to use immediate integer1089** arithmetic if possible to avoid the function call overhead.1090*/1091Obj ProdCycInt (1092Obj opL,1093Obj opR )1094{1095Obj hdP; /* product, result */1096UInt len; /* number of terms */1097Obj * cfs; /* ptr to coeffs of left operand */1098UInt4 * exs; /* ptr to expnts of left operand */1099Obj * cfp; /* ptr to coeffs of product */1100UInt4 * exp; /* ptr to expnts of product */1101UInt i; /* loop variable */1102Obj prd; /* product of two coefficients */11031104/* for $rat * rat$ delegate */1105if ( TNUM_OBJ(opL) != T_CYC && TNUM_OBJ(opR) != T_CYC ) {1106return PROD( opL, opR );1107}11081109/* make the right operand the non cyclotomic */1110if ( TNUM_OBJ(opL) != T_CYC ) { hdP = opL; opL = opR; opR = hdP; }11111112/* for $cyc * 0$ return 0 and for $cyc * 1$ return $cyc$ */1113if ( opR == INTOBJ_INT(0) ) {1114hdP = INTOBJ_INT(0);1115}1116else if ( opR == INTOBJ_INT(1) ) {1117hdP = opL;1118}11191120/* for $cyc * -1$ need no multiplication or division */1121else if ( opR == INTOBJ_INT(-1) ) {1122return AInvCyc( opL );1123}11241125/* for $cyc * small$ use immediate multiplication if possible */1126else if ( TNUM_OBJ(opR) == T_INT ) {1127hdP = NewBag( T_CYC, SIZE_CYC(opL) * (sizeof(Obj)+sizeof(UInt4)) );1128NOF_CYC(hdP) = NOF_CYC(opL);1129len = SIZE_CYC(opL);1130cfs = COEFS_CYC(opL);1131cfp = COEFS_CYC(hdP);1132exs = EXPOS_CYC(opL,len);1133exp = EXPOS_CYC(hdP,len);1134for ( i = 1; i < len; i++ ) {1135if ( ! IS_INTOBJ( cfs[i] )1136|| ! PROD_INTOBJS( prd, cfs[i], opR ) ) {1137CHANGED_BAG( hdP );1138prd = PROD( cfs[i], opR );1139cfs = COEFS_CYC(opL);1140cfp = COEFS_CYC(hdP);1141exs = EXPOS_CYC(opL,len);1142exp = EXPOS_CYC(hdP,len);1143}1144cfp[i] = prd;1145exp[i] = exs[i];1146}1147CHANGED_BAG( hdP );1148}11491150/* otherwise multiply every coefficent */1151else {1152hdP = NewBag( T_CYC, SIZE_CYC(opL) * (sizeof(Obj)+sizeof(UInt4)) );1153NOF_CYC(hdP) = NOF_CYC(opL);1154len = SIZE_CYC(opL);1155cfs = COEFS_CYC(opL);1156cfp = COEFS_CYC(hdP);1157exs = EXPOS_CYC(opL,len);1158exp = EXPOS_CYC(hdP,len);1159for ( i = 1; i < len; i++ ) {1160CHANGED_BAG( hdP );1161prd = PROD( cfs[i], opR );1162cfs = COEFS_CYC(opL);1163cfp = COEFS_CYC(hdP);1164exs = EXPOS_CYC(opL,len);1165exp = EXPOS_CYC(hdP,len);1166cfp[i] = prd;1167exp[i] = exs[i];1168}1169CHANGED_BAG( hdP );1170}11711172/* return the result */1173return hdP;1174}117511761177/****************************************************************************1178**1179*F ProdCyc( <opL>, <opR> ) . . . . . . . . . . . product of two cyclotomics1180**1181** 'ProdCyc' returns the product of the two cyclotomics <opL> and <opR>.1182** Either operand may also be an integer or a rational.1183**1184** This function is lengthy because we try to use immediate integer1185** arithmetic if possible to avoid the function call overhead.1186*/1187Obj ProdCyc (1188Obj opL,1189Obj opR )1190{1191UInt nl, nr; /* order of left and right field */1192UInt n; /* order of smallest superfield */1193UInt ml, mr; /* cofactors into the superfield */1194Obj c; /* one coefficient of the left op */1195UInt e; /* one exponent of the left op */1196UInt len; /* number of terms */1197Obj * cfs; /* pointer to the coefficients */1198UInt4 * exs; /* pointer to the exponents */1199Obj * res; /* pointer to the result */1200Obj sum; /* sum of two coefficients */1201Obj prd; /* product of two coefficients */1202UInt i, k; /* loop variable */12031204/* for $rat * cyc$ and $cyc * rat$ delegate */1205if ( TNUM_OBJ(opL) != T_CYC || TNUM_OBJ(opR) != T_CYC ) {1206return ProdCycInt( opL, opR );1207}12081209/* take the cyclotomic with less terms as the right operand */1210if ( SIZE_CYC(opL) < SIZE_CYC(opR) ) {1211prd = opL; opL = opR; opR = prd;1212}12131214/* get the smallest field that contains both cyclotomics */1215nl = (TNUM_OBJ(opL) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opL) ));1216nr = (TNUM_OBJ(opR) != T_CYC ? 1 : INT_INTOBJ( NOF_CYC(opR) ));1217n = FindCommonField(nl, nr, &ml, &mr);12181219/* loop over the terms of the right operand */1220for ( k = 1; k < SIZE_CYC(opR); k++ ) {1221c = COEFS_CYC(opR)[k];1222e = (mr * EXPOS_CYC( opR, SIZE_CYC(opR) )[k]) % n;12231224/* if the coefficient is 1 just add */1225if ( c == INTOBJ_INT(1) ) {1226len = SIZE_CYC(opL);1227cfs = COEFS_CYC(opL);1228exs = EXPOS_CYC(opL,len);1229res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1230for ( i = 1; i < len; i++ ) {1231if ( ! ARE_INTOBJS( res[(e+exs[i]*ml)%n], cfs[i] )1232|| ! SUM_INTOBJS( sum, res[(e+exs[i]*ml)%n], cfs[i] ) ) {1233CHANGED_BAG( TLS(ResultCyc) );1234sum = SUM( res[(e+exs[i]*ml)%n], cfs[i] );1235cfs = COEFS_CYC(opL);1236exs = EXPOS_CYC(opL,len);1237res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1238}1239res[(e+exs[i]*ml)%n] = sum;1240}1241CHANGED_BAG( TLS(ResultCyc) );1242}12431244/* if the coefficient is -1 just subtract */1245else if ( c == INTOBJ_INT(-1) ) {1246len = SIZE_CYC(opL);1247cfs = COEFS_CYC(opL);1248exs = EXPOS_CYC(opL,len);1249res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1250for ( i = 1; i < len; i++ ) {1251if ( ! ARE_INTOBJS( res[(e+exs[i]*ml)%n], cfs[i] )1252|| ! DIFF_INTOBJS( sum, res[(e+exs[i]*ml)%n], cfs[i] ) ) {1253CHANGED_BAG( TLS(ResultCyc) );1254sum = DIFF( res[(e+exs[i]*ml)%n], cfs[i] );1255cfs = COEFS_CYC(opL);1256exs = EXPOS_CYC(opL,len);1257res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1258}1259res[(e+exs[i]*ml)%n] = sum;1260}1261CHANGED_BAG( TLS(ResultCyc) );1262}12631264/* if the coefficient is a small integer use immediate operations */1265else if ( TNUM_OBJ(c) == T_INT ) {1266len = SIZE_CYC(opL);1267cfs = COEFS_CYC(opL);1268exs = EXPOS_CYC(opL,len);1269res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1270for ( i = 1; i < len; i++ ) {1271if ( ! ARE_INTOBJS( cfs[i], res[(e+exs[i]*ml)%n] )1272|| ! PROD_INTOBJS( prd, cfs[i], c )1273|| ! SUM_INTOBJS( sum, res[(e+exs[i]*ml)%n], prd ) ) {1274CHANGED_BAG( TLS(ResultCyc) );1275prd = PROD( cfs[i], c );1276exs = EXPOS_CYC(opL,len);1277res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1278sum = SUM( res[(e+exs[i]*ml)%n], prd );1279cfs = COEFS_CYC(opL);1280exs = EXPOS_CYC(opL,len);1281res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1282}1283res[(e+exs[i]*ml)%n] = sum;1284}1285CHANGED_BAG( TLS(ResultCyc) );1286}12871288/* otherwise do it the normal way */1289else {1290len = SIZE_CYC(opL);1291for ( i = 1; i < len; i++ ) {1292CHANGED_BAG( TLS(ResultCyc) );1293cfs = COEFS_CYC(opL);1294prd = PROD( cfs[i], c );1295exs = EXPOS_CYC(opL,len);1296res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1297sum = SUM( res[(e+exs[i]*ml)%n], prd );1298exs = EXPOS_CYC(opL,len);1299res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1300res[(e+exs[i]*ml)%n] = sum;1301}1302CHANGED_BAG( TLS(ResultCyc) );1303}13041305}13061307/* return the base reduced packed cyclotomic */1308ConvertToBase( n );1309return Cyclotomic( n, ml * mr );1310}131113121313/****************************************************************************1314**1315*F OneCyc( <op> ) . . . . . . . . . . . . . . . . . . . one of a cyclotomic1316**1317** 'OneCyc' returns the multiplicative neutral element of the cyclotomic1318** <op>.1319*/1320Obj OneCyc (1321Obj op )1322{1323return INTOBJ_INT( 1L );1324}132513261327/****************************************************************************1328**1329*F InvCyc( <op> ) . . . . . . . . . . . . . . . . . inverse of a cyclotomic1330**1331** 'InvCyc' returns the multiplicative inverse element of the cyclotomic1332** <op>.1333**1334** 'InvCyc' computes the inverse of <op> by computing the product $prd$ of1335** nontrivial Galois conjugates of <op>. Then $op * (prd / (op * prd)) = 1$1336** so $prd / (op * prd)$ is the inverse of $op$. Because the denominator1337** $op * prd$ is the norm of $op$ over the rationals it is rational so we1338** can compute the quotient $prd / (op * prd)$ with 'ProdCycInt'.1339*T better multiply only the *different* conjugates?1340*/1341Obj InvCyc (1342Obj op )1343{1344Obj prd; /* product of conjugates */1345UInt n; /* order of the field */1346UInt sqr; /* if n < sqr*sqr n is squarefree */1347UInt len; /* number of terms */1348Obj * cfs; /* pointer to the coefficients */1349UInt4 * exs; /* pointer to the exponents */1350Obj * res; /* pointer to the result */1351UInt i, k; /* loop variable */1352UInt gcd, s, t; /* gcd of i and n, temporaries */13531354/* get the order of the field, test if it is squarefree */1355n = INT_INTOBJ( NOF_CYC(op) );1356for ( sqr = 2; sqr*sqr <= n && n % (sqr*sqr) != 0; sqr++ )1357;13581359/* compute the product of all nontrivial galois conjugates of <opL> */1360len = SIZE_CYC(op);1361prd = INTOBJ_INT(1);1362for ( i = 2; i < n; i++ ) {13631364/* if i gives a galois automorphism apply it */1365gcd = n; s = i; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }1366if ( gcd == 1 ) {13671368/* permute the terms */1369cfs = COEFS_CYC(op);1370exs = EXPOS_CYC(op,len);1371res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1372for ( k = 1; k < len; k++ )1373res[(i*exs[k])%n] = cfs[k];1374CHANGED_BAG( TLS(ResultCyc) );13751376/* if n is squarefree conversion and reduction are unnecessary */1377if ( n < sqr*sqr ) {1378prd = ProdCyc( prd, Cyclotomic( n, n ) );1379}1380else {1381ConvertToBase( n );1382prd = ProdCyc( prd, Cyclotomic( n, 1 ) );1383}13841385}13861387}13881389/* the inverse is the product divided by the norm */1390return ProdCycInt( prd, INV( ProdCyc( op, prd ) ) );1391}139213931394/****************************************************************************1395**1396*F PowCyc( <opL>, <opR> ) . . . . . . . . . . . . . . power of a cyclotomic1397**1398** 'PowCyc' returns the <opR>th, which must be an integer, power of the1399** cyclotomic <opL>. The left operand may also be an integer or a rational.1400*/1401Obj PowCyc (1402Obj opL,1403Obj opR )1404{1405Obj pow; /* power (result) */1406Int exp; /* exponent (right operand) */1407UInt n; /* order of the field */1408Obj * res; /* pointer into the result */1409UInt i; /* exponent of left operand */14101411/* get the exponent */1412exp = INT_INTOBJ(opR);14131414/* for $cyc^0$ return 1, for $cyc^1$ return cyc, for $rat^exp$ delegate*/1415if ( exp == 0 ) {1416pow = INTOBJ_INT(1);1417}1418else if ( exp == 1 ) {1419pow = opL;1420}1421else if ( TNUM_OBJ(opL) != T_CYC ) {1422pow = PowInt( opL, opR );1423}14241425/* for $e_n^exp$ just put a 1 at the <exp>th position and convert */1426else if ( opL == TLS(LastECyc) ) {1427exp = (exp % TLS(LastNCyc) + TLS(LastNCyc)) % TLS(LastNCyc);1428SET_ELM_PLIST( TLS(ResultCyc), exp, INTOBJ_INT(1) );1429CHANGED_BAG( TLS(ResultCyc) );1430ConvertToBase( TLS(LastNCyc) );1431pow = Cyclotomic( TLS(LastNCyc), 1 );1432}14331434/* for $(c*e_n^i)^exp$ if $e_n^i$ belongs to the base put 1 at $i*exp$ */1435else if ( SIZE_CYC(opL) == 2 ) {1436n = INT_INTOBJ( NOF_CYC(opL) );1437pow = POW( COEFS_CYC(opL)[1], opR );1438i = EXPOS_CYC(opL,2)[1];1439res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1440res[((exp*i)%n+n)%n] = pow;1441CHANGED_BAG( TLS(ResultCyc) );1442ConvertToBase( n );1443pow = Cyclotomic( n, 1 );1444}14451446/* otherwise compute the power with repeated squaring */1447else {14481449/* if neccessary invert the cyclotomic */1450if ( exp < 0 ) {1451opL = InvCyc( opL );1452exp = -exp;1453}14541455/* compute the power using repeated squaring */1456pow = INTOBJ_INT(1);1457while ( exp != 0 ) {1458if ( exp % 2 == 1 ) pow = ProdCyc( pow, opL );1459if ( exp > 1 ) opL = ProdCyc( opL, opL );1460exp = exp / 2;1461}14621463}14641465/* return the result */1466return pow;1467}146814691470/****************************************************************************1471**1472*F FuncE( <self>, <n> ) . . . . . . . . . . . . create a new primitive root1473**1474** 'FuncE' implements the internal function 'E'.1475**1476** 'E( <n> )'1477**1478** 'E' returns a primitive root of order <n>, which must be a positive1479** integer, represented as cyclotomic.1480*/1481Obj EOper;14821483Obj FuncE (1484Obj self,1485Obj n )1486{1487Obj * res; /* pointer into result bag */14881489/* do full operation */1490if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(n) ) {1491return DoOperation1Args( self, n );1492}14931494/* get and check the argument */1495while ( TNUM_OBJ(n) != T_INT || INT_INTOBJ(n) <= 0 ) {1496n = ErrorReturnObj(1497"E: <n> must be a positive integer (not a %s)",1498(Int)TNAM_OBJ(n), 0L,1499"you can replace <n> via 'return <n>;'" );1500}15011502/* for $e_1$ return 1 and for $e_2$ return -1 */1503if ( n == INTOBJ_INT(1) )1504return INTOBJ_INT(1);1505else if ( n == INTOBJ_INT(2) )1506return INTOBJ_INT(-1);15071508/* if the root is not known already construct it */1509if ( TLS(LastNCyc) != INT_INTOBJ(n) ) {1510TLS(LastNCyc) = INT_INTOBJ(n);1511GrowResultCyc(TLS(LastNCyc));1512res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1513res[1] = INTOBJ_INT(1);1514CHANGED_BAG( TLS(ResultCyc) );1515ConvertToBase( TLS(LastNCyc) );1516TLS(LastECyc) = Cyclotomic( TLS(LastNCyc), 1 );1517}15181519/* return the root */1520return TLS(LastECyc);1521}152215231524/****************************************************************************1525**1526*F FuncIS_CYC( <self>, <val> ) . . . . . . test if an object is a cyclomtic1527**1528** 'FuncIS_CYC' implements the internal function 'IsCyc'.1529**1530** 'IsCyc( <val> )'1531**1532** 'IsCyc' returns 'true' if the value <val> is a cyclotomic and 'false'1533** otherwise.1534*/1535Obj IsCycFilt;15361537Obj FuncIS_CYC (1538Obj self,1539Obj val )1540{1541/* return 'true' if <obj> is a cyclotomic and 'false' otherwise */1542if ( TNUM_OBJ(val) == T_CYC1543|| TNUM_OBJ(val) == T_INT || TNUM_OBJ(val) == T_RAT1544|| TNUM_OBJ(val) == T_INTPOS || TNUM_OBJ(val) == T_INTNEG )1545return True;1546else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {1547return False;1548}1549else {1550return DoFilter( self, val );1551}1552}155315541555/****************************************************************************1556**1557*F FuncIS_CYC_INT( <self>, <val> ) test if an object is a cyclomtic integer1558**1559** 'FuncIS_CYC_INT' implements the internal function 'IsCycInt'.1560**1561** 'IsCycInt( <val> )'1562**1563** 'IsCycInt' returns 'true' if the value <val> is a cyclotomic integer and1564** 'false' otherwise.1565**1566** 'IsCycInt' relies on the fact that the base is an integral base.1567*/1568Obj IsCycIntOper;15691570Obj FuncIS_CYC_INT (1571Obj self,1572Obj val )1573{1574UInt len; /* number of terms */1575Obj * cfs; /* pointer to the coefficients */1576UInt i; /* loop variable */15771578/* return 'true' if <obj> is a cyclotomic integer and 'false' otherwise*/1579if ( TNUM_OBJ(val) == T_INT1580|| TNUM_OBJ(val) == T_INTPOS || TNUM_OBJ(val) == T_INTNEG ) {1581return True;1582}1583else if ( TNUM_OBJ(val) == T_RAT ) {1584return False;1585}1586else if ( TNUM_OBJ(val) == T_CYC ) {1587len = SIZE_CYC(val);1588cfs = COEFS_CYC(val);1589for ( i = 1; i < len; i++ ) {1590if ( TNUM_OBJ(cfs[i]) == T_RAT )1591return False;1592}1593return True;1594}1595else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {1596return False;1597}1598else {1599return DoOperation1Args( self, val );1600}1601}160216031604/****************************************************************************1605**1606*F FuncCONDUCTOR( <self>, <cyc> ) . . . . . . . . . . . . N of a cyclotomic1607**1608** 'FuncCONDUCTOR' implements the internal function 'Conductor'.1609**1610** 'Conductor( <cyc> )'1611**1612** 'Conductor' returns the N of the cyclotomic <cyc>, i.e., the order of the1613** roots of which <cyc> is written as a linear combination.1614*/1615Obj ConductorAttr;16161617Obj FuncCONDUCTOR (1618Obj self,1619Obj cyc )1620{1621UInt n; /* N of the cyclotomic, result */1622UInt m; /* N of element of the list */1623UInt gcd, s, t; /* gcd of n and m, temporaries */1624Obj list; /* list of cyclotomics */1625UInt i; /* loop variable */16261627/* do full operation */1628if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(cyc) ) {1629return DoAttribute( ConductorAttr, cyc );1630}16311632/* check the argument */1633while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT1634&& TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG1635&& TNUM_OBJ(cyc) != T_CYC1636&& ! IS_SMALL_LIST(cyc) ) {1637cyc = ErrorReturnObj(1638"Conductor: <cyc> must be a cyclotomic or a small list (not a %s)",1639(Int)TNAM_OBJ(cyc), 0L,1640"you can replace <cyc> via 'return <cyc>;'" );1641}16421643/* handle cyclotomics */1644if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT1645|| TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) {1646n = 1;1647}1648else if ( TNUM_OBJ(cyc) == T_CYC ) {1649n = INT_INTOBJ( NOF_CYC(cyc) );1650}16511652/* handle a list by computing the lcm of the entries */1653else {1654list = cyc;1655n = 1;1656for ( i = 1; i <= LEN_LIST( list ); i++ ) {1657cyc = ELMV_LIST( list, i );1658while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT1659&& TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG1660&& TNUM_OBJ(cyc) != T_CYC ) {1661cyc = ErrorReturnObj(1662"Conductor: <list>[%d] must be a cyclotomic (not a %s)",1663(Int)i, (Int)TNAM_OBJ(cyc),1664"you can replace the list element with <cyc> via 'return <cyc>;'" );1665}1666if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT1667|| TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) {1668m = 1;1669}1670else /* if ( TNUM_OBJ(cyc) == T_CYC ) */ {1671m = INT_INTOBJ( NOF_CYC(cyc) );1672}1673gcd = n; s = m; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }1674n = n / gcd * m;1675}1676}16771678/* return the N of the cyclotomic */1679return INTOBJ_INT( n );1680}168116821683/****************************************************************************1684**1685*F FuncCOEFFS_CYC( <self>, <cyc> ) . . . . . . coefficients of a cyclotomic1686**1687** 'FuncCOEFFS_CYC' implements the internal function 'COEFFSCYC'.1688**1689** 'COEFFSCYC( <cyc> )'1690**1691** 'COEFFSCYC' returns a list of the coefficients of the cyclotomic <cyc>.1692** The list has length <n> if <n> is the order of the primitive root $e_n$1693** of which <cyc> is written as a linear combination. The <i>th element of1694** the list is the coefficient of $e_l^{i-1}$.1695*/1696Obj CoeffsCycOper;16971698Obj FuncCOEFFS_CYC (1699Obj self,1700Obj cyc )1701{1702Obj list; /* list of coefficients, result */1703UInt n; /* order of field */1704UInt len; /* number of terms */1705Obj * cfs; /* pointer to the coefficients */1706UInt4 * exs; /* pointer to the exponents */1707UInt i; /* loop variable */17081709/* do full operation */1710if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(cyc) ) {1711return DoOperation1Args( self, cyc );1712}17131714/* check the argument */1715while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT1716&& TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG1717&& TNUM_OBJ(cyc) != T_CYC ) {1718cyc = ErrorReturnObj(1719"COEFFSCYC: <cyc> must be a cyclotomic (not a %s)",1720(Int)TNAM_OBJ(cyc), 0L,1721"you can replace <cyc> via 'return <cyc>;'" );1722}17231724/* if <cyc> is rational just put it in a list of length 1 */1725if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT1726|| TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) {1727list = NEW_PLIST( T_PLIST, 1 );1728SET_LEN_PLIST( list, 1 );1729SET_ELM_PLIST( list, 1, cyc );1730/* 'CHANGED_BAG' not needed for last bag */1731}17321733/* otherwise make a list and fill it with zeroes and copy the coeffs */1734else {1735n = INT_INTOBJ( NOF_CYC(cyc) );1736list = NEW_PLIST( T_PLIST, n );1737SET_LEN_PLIST( list, n );1738len = SIZE_CYC(cyc);1739cfs = COEFS_CYC(cyc);1740exs = EXPOS_CYC(cyc,len);1741for ( i = 1; i <= n; i++ )1742SET_ELM_PLIST( list, i, INTOBJ_INT(0) );1743for ( i = 1; i < len; i++ )1744SET_ELM_PLIST( list, exs[i]+1, cfs[i] );1745/* 'CHANGED_BAG' not needed for last bag */1746}17471748/* return the result */1749return list;1750}175117521753/****************************************************************************1754**1755*F FuncGALOIS_CYC( <self>, <cyc>, <ord> ) image of a cyc. under galois aut1756**1757** 'FuncGALOIS_CYC' implements the internal function 'GaloisCyc'.1758**1759** 'GaloisCyc( <cyc>, <ord> )'1760**1761** 'GaloisCyc' computes the image of the cyclotomic <cyc> under the galois1762** automorphism given by <ord>, which must be an integer.1763**1764** The galois automorphism is the mapping that takes $e_n$ to $e_n^ord$.1765** <ord> may be any integer, of course if it is not relative prime to $ord$1766** the mapping will not be an automorphism, though still an endomorphism.1767*/1768Obj GaloisCycOper;17691770Obj FuncGALOIS_CYC (1771Obj self,1772Obj cyc,1773Obj ord )1774{1775Obj gal; /* galois conjugate, result */1776Obj sum; /* sum of two coefficients */1777Int n; /* order of the field */1778UInt sqr; /* if n < sqr*sqr n is squarefree */1779Int o; /* galois automorphism */1780UInt gcd, s, t; /* gcd of n and ord, temporaries */1781UInt len; /* number of terms */1782Obj * cfs; /* pointer to the coefficients */1783UInt4 * exs; /* pointer to the exponents */1784Obj * res; /* pointer to the result */1785UInt i; /* loop variable */1786UInt tnumord, tnumcyc;17871788/* do full operation for any but standard arguments */17891790tnumord = TNUM_OBJ(ord);1791tnumcyc = TNUM_OBJ(cyc);1792if ( FIRST_EXTERNAL_TNUM <= tnumcyc1793|| FIRST_EXTERNAL_TNUM <= tnumord1794|| (tnumord != T_INT && tnumord != T_INTNEG && tnumord != T_INTPOS)1795|| ( tnumcyc != T_INT && tnumcyc != T_RAT1796&& tnumcyc != T_INTPOS && tnumcyc != T_INTNEG1797&& tnumcyc != T_CYC )1798)1799{1800return DoOperation2Args( self, cyc, ord );1801}18021803/* get and check <ord> */1804if ( TNAM_OBJ(ord) == T_INT ) {1805o = INT_INTOBJ(ord);1806}1807else {1808ord = MOD( ord, FuncCONDUCTOR( 0, cyc ) );1809o = INT_INTOBJ(ord);1810}18111812/* every galois automorphism fixes the rationals */1813if ( tnumcyc == T_INT || tnumcyc == T_RAT1814|| tnumcyc == T_INTPOS || tnumcyc == T_INTNEG ) {1815return cyc;1816}18171818/* get the order of the field, test if it squarefree */1819n = INT_INTOBJ( NOF_CYC(cyc) );1820for ( sqr = 2; sqr*sqr <= n && n % (sqr*sqr) != 0; sqr++ )1821;18221823/* force <ord> into the range 0..n-1, compute the gcd of <ord> and <n> */1824o = (o % n + n) % n;1825gcd = n; s = o; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; }18261827/* if <ord> = 1 just return <cyc> */1828if ( o == 1 ) {1829gal = cyc;1830}18311832/* if <ord> == 0 compute the sum of the entries */1833else if ( o == 0 ) {1834len = SIZE_CYC(cyc);1835cfs = COEFS_CYC(cyc);1836gal = INTOBJ_INT(0);1837for ( i = 1; i < len; i++ ) {1838if ( ! ARE_INTOBJS( gal, cfs[i] )1839|| ! SUM_INTOBJS( sum, gal, cfs[i] ) ) {1840sum = SUM( gal, cfs[i] );1841cfs = COEFS_CYC( cyc );1842}1843gal = sum;1844}1845}18461847/* if <ord> == n/2 compute alternating sum since $(e_n^i)^ord = -1^i$ */1848else if ( n % 2 == 0 && o == n/2 ) {1849gal = INTOBJ_INT(0);1850len = SIZE_CYC(cyc);1851cfs = COEFS_CYC(cyc);1852exs = EXPOS_CYC(cyc,len);1853for ( i = 1; i < len; i++ ) {1854if ( exs[i] % 2 == 1 ) {1855if ( ! ARE_INTOBJS( gal, cfs[i] )1856|| ! DIFF_INTOBJS( sum, gal, cfs[i] ) ) {1857sum = DIFF( gal, cfs[i] );1858cfs = COEFS_CYC(cyc);1859exs = EXPOS_CYC(cyc,len);1860}1861gal = sum;1862}1863else {1864if ( ! ARE_INTOBJS( gal, cfs[i] )1865|| ! SUM_INTOBJS( sum, gal, cfs[i] ) ) {1866sum = SUM( gal, cfs[i] );1867cfs = COEFS_CYC(cyc);1868exs = EXPOS_CYC(cyc,len);1869}1870gal = sum;1871}1872}1873}18741875/* if <ord> is prime to <n> (automorphism) permute the coefficients */1876else if ( gcd == 1 ) {18771878/* permute the coefficients */1879len = SIZE_CYC(cyc);1880cfs = COEFS_CYC(cyc);1881exs = EXPOS_CYC(cyc,len);1882res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1883for ( i = 1; i < len; i++ ) {1884res[(UInt8)exs[i]*(UInt8)o%(UInt8)n] = cfs[i];1885}1886CHANGED_BAG( TLS(ResultCyc) );18871888/* if n is squarefree conversion and reduction are unnecessary */1889if ( n < sqr*sqr || (o == n-1 && n % 2 != 0) ) {1890gal = Cyclotomic( n, n );1891}1892else {1893ConvertToBase( n );1894gal = Cyclotomic( n, 1 );1895}18961897}18981899/* if <ord> is not prime to <n> (endomorphism) compute it the hard way */1900else {19011902/* multiple roots may be mapped to the same root, add the coeffs */1903len = SIZE_CYC(cyc);1904cfs = COEFS_CYC(cyc);1905exs = EXPOS_CYC(cyc,len);1906res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1907for ( i = 1; i < len; i++ ) {1908if ( ! ARE_INTOBJS( res[(UInt8)exs[i]*(UInt8)o%(UInt8)n], cfs[i] )1909|| ! SUM_INTOBJS( sum, res[(UInt8)exs[i]*(UInt8)o%(UInt8)n], cfs[i] ) ) {1910CHANGED_BAG( TLS(ResultCyc) );1911sum = SUM( res[(UInt8)exs[i]*(UInt8)o%(UInt8)n], cfs[i] );1912cfs = COEFS_CYC(cyc);1913exs = EXPOS_CYC(cyc,len);1914res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1915}1916res[exs[i]*o%n] = sum;1917}1918CHANGED_BAG( TLS(ResultCyc) );19191920/* if n is squarefree conversion and reduction are unnecessary */1921if ( n < sqr*sqr ) {1922gal = Cyclotomic( n, 1 ); /*N?*/1923}1924else {1925ConvertToBase( n );1926gal = Cyclotomic( n, 1 );1927}19281929}19301931/* return the result */1932return gal;1933}193419351936/****************************************************************************1937**1938*F FuncCycList( <self>, <list> ) . . . . . . . . . . . . create a cyclotomic1939**1940** 'FuncCycList' implements the internal function 'CycList'.1941**1942** 'CycList( <list> )'1943**1944** 'CycList' returns the cyclotomic described by the list <list>1945** of rationals.1946*/1947Obj CycListOper;19481949Obj FuncCycList (1950Obj self,1951Obj list )1952{1953UInt i; /* loop variable */1954Obj * res; /* pointer into result bag */1955Obj val; /* one list entry */1956UInt n; /* length of the given list */19571958/* do full operation */1959if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ( list ) ) {1960return DoOperation1Args( self, list );1961}19621963/* get and check the argument */1964if ( ! IS_PLIST( list ) || ! IS_DENSE_LIST( list ) ) {1965ErrorQuit( "CycList: <list> must be a dense plain list (not a %s)",1966(Int)TNAM_OBJ( list ), 0L );1967}19681969/* enlarge the buffer if necessary */1970n = LEN_PLIST( list );1971GrowResultCyc(n);19721973/* transfer the coefficients into the buffer */1974res = &(ELM_PLIST( TLS(ResultCyc), 1 ));1975for ( i = 0; i < n; i++ ) {1976val = ELM_PLIST( list, i+1 );1977if ( ! ( TNUM_OBJ(val) == T_INT ||1978TNUM_OBJ(val) == T_RAT ||1979TNUM_OBJ(val) == T_INTPOS ||1980TNUM_OBJ(val) == T_INTNEG ) ) {1981ErrorQuit( "CycList: each entry must be a rational (not a %s)",1982(Int)TNAM_OBJ( val ), 0L );1983}1984res[i] = val;1985}19861987/* return the base reduced packed cyclotomic */1988CHANGED_BAG( TLS(ResultCyc) );1989ConvertToBase( n );1990return Cyclotomic( n, 1 );1991}199219931994/****************************************************************************1995**1996*F MarkCycSubBags( <bag> ) . . . . . . . . marking function for cyclotomics1997**1998** 'MarkCycSubBags' is the marking function for bags of type 'T_CYC'.1999*/2000void MarkCycSubBags (2001Obj cyc )2002{2003Obj * cfs; /* pointer to coeffs */2004UInt i; /* loop variable */20052006/* mark everything */2007cfs = COEFS_CYC( cyc );2008for ( i = SIZE_CYC(cyc); 0 < i; i-- ) {2009MARK_BAG( cfs[i-1] );2010}20112012}201320142015/****************************************************************************2016**20172018*F SaveCyc() . . . . . . . . . . . . . . . . . . . . . . save a cyclotyomic2019**2020** We do not save the XXX_CYC field, since it is not used.2021*/2022void SaveCyc ( Obj cyc )2023{2024UInt len, i;2025Obj *coefs;2026UInt4 *expos;2027len = SIZE_CYC(cyc);2028coefs = COEFS_CYC(cyc);2029for (i = 0; i < len; i++)2030SaveSubObj(*coefs++);2031expos = EXPOS_CYC(cyc,len);2032expos++; /*Skip past the XXX */2033for (i = 1; i < len; i++)2034SaveUInt4(*expos++);2035return;2036}20372038/****************************************************************************2039**2040*F LoadCyc() . . . . . . . . . . . . . . . . . . . . . . load a cyclotyomic2041**2042** We do not load the XXX_CYC field, since it is not used.2043*/2044void LoadCyc ( Obj cyc )2045{2046UInt len, i;2047Obj *coefs;2048UInt4 *expos;2049len = SIZE_CYC(cyc);2050coefs = COEFS_CYC(cyc);2051for (i = 0; i < len; i++)2052*coefs++ = LoadSubObj();2053expos = EXPOS_CYC(cyc,len);2054expos++; /*Skip past the XXX */2055for (i = 1; i < len; i++)2056*expos++ = LoadUInt4();2057return;2058}205920602061/****************************************************************************2062**2063**20642065*F * * * * * * * * * * * * * initialize package * * * * * * * * * * * * * * *2066*/206720682069/****************************************************************************2070**20712072*V GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export2073*/2074static StructGVarFilt GVarFilts [] = {20752076{ "IS_CYC", "obj", &IsCycFilt,2077FuncIS_CYC, "src/cyclotom.c:IS_CYC" },20782079{ 0 }20802081};208220832084/****************************************************************************2085**2086*V GVarAttrs . . . . . . . . . . . . . . . . . list of attributes to export2087*/2088static StructGVarAttr GVarAttrs [] = {20892090{ "CONDUCTOR", "cyc", &ConductorAttr,2091FuncCONDUCTOR, "src/cyclotom.c:CONDUCTOR" },20922093{ 0 }20942095};209620972098/****************************************************************************2099**2100*V GVarOpers . . . . . . . . . . . . . . . . . list of operations to export2101*/2102static StructGVarOper GVarOpers [] = {21032104{ "E", 1, "n", &EOper,2105FuncE, "src/cyclotom.c:E" },21062107{ "IS_CYC_INT", 1, "obj", &IsCycIntOper,2108FuncIS_CYC_INT, "src/cyclotom.c:IS_CYC_INT" },21092110{ "COEFFS_CYC", 1, "cyc", &CoeffsCycOper,2111FuncCOEFFS_CYC, "src/cyclotom.c:COEFFS_CYC" },21122113{ "GALOIS_CYC", 2, "cyc, n", &GaloisCycOper,2114FuncGALOIS_CYC, "src/cyclotom.c:GALOIS_CYC" },21152116{ "CycList", 1, "list", &CycListOper,2117FuncCycList, "src/cyclotom.c:CycList" },211821192120{ 0 }21212122};212321242125/****************************************************************************2126**2127*V GVarFuncs . . . . . . . . . . . . . . . . . list of functions to export2128*/2129static StructGVarFunc GVarFuncs [] = {2130{ "SetCyclotomicsLimit", 1, "newlimit",2131FuncSetCyclotomicsLimit, "src/cyclotomics.c:SetCyclotomicsLimit" },21322133{ "GetCyclotomicsLimit", 0, "",2134FuncGetCyclotomicsLimit, "src/cyclotomics.c:GetCyclotomicsLimit" },21352136{0}2137};21382139/****************************************************************************2140**21412142*F InitKernel( <module> ) . . . . . . . . initialise kernel data structures2143*/2144static Int InitKernel (2145StructInitInfo * module )2146{2147TLS(LastECyc) = (Obj)0;2148TLS(LastNCyc) = 0;21492150/* install the marking function */2151InfoBags[ T_CYC ].name = "cyclotomic";2152InitMarkFuncBags( T_CYC, MarkCycSubBags );21532154/* create the result buffer */2155InitGlobalBag( &ResultCyc , "src/cyclotom.c:ResultCyc" );21562157/* tell Gasman about the place were we remember the primitive root */2158InitGlobalBag( &LastECyc, "src/cyclotom.c:LastECyc" );21592160/* install the type function */2161ImportGVarFromLibrary( "TYPE_CYC", &TYPE_CYC );2162TypeObjFuncs[ T_CYC ] = TypeCyc;21632164/* init filters and functions */2165InitHdlrFiltsFromTable( GVarFilts );2166InitHdlrAttrsFromTable( GVarAttrs );2167InitHdlrOpersFromTable( GVarOpers );2168InitHdlrFuncsFromTable( GVarFuncs );21692170/* and the saving function */2171SaveObjFuncs[ T_CYC ] = SaveCyc;2172LoadObjFuncs[ T_CYC ] = LoadCyc;21732174/* install the evaluation and print function */2175PrintObjFuncs[ T_CYC ] = PrintCyc;21762177/* install the comparison methods */2178EqFuncs[ T_CYC ][ T_CYC ] = EqCyc;2179LtFuncs[ T_CYC ][ T_CYC ] = LtCyc;2180LtFuncs[ T_INT ][ T_CYC ] = LtCycYes;2181LtFuncs[ T_INTPOS ][ T_CYC ] = LtCycYes;2182LtFuncs[ T_INTNEG ][ T_CYC ] = LtCycYes;2183LtFuncs[ T_RAT ][ T_CYC ] = LtCycYes;2184LtFuncs[ T_CYC ][ T_INT ] = LtCycNot;2185LtFuncs[ T_CYC ][ T_INTPOS ] = LtCycNot;2186LtFuncs[ T_CYC ][ T_INTNEG ] = LtCycNot;2187LtFuncs[ T_CYC ][ T_RAT ] = LtCycNot;21882189/* install the unary arithmetic methods */2190ZeroFuncs[ T_CYC ] = ZeroCyc;2191ZeroMutFuncs[ T_CYC ] = ZeroCyc;2192AInvFuncs[ T_CYC ] = AInvCyc;2193AInvMutFuncs[ T_CYC ] = AInvCyc;2194OneFuncs [ T_CYC ] = OneCyc;2195OneMutFuncs [ T_CYC ] = OneCyc;2196InvFuncs [ T_CYC ] = InvCyc;2197InvMutFuncs [ T_CYC ] = InvCyc;21982199/* install the arithmetic methods */2200SumFuncs[ T_CYC ][ T_CYC ] = SumCyc;2201SumFuncs[ T_INT ][ T_CYC ] = SumCyc;2202SumFuncs[ T_INTPOS ][ T_CYC ] = SumCyc;2203SumFuncs[ T_INTNEG ][ T_CYC ] = SumCyc;2204SumFuncs[ T_RAT ][ T_CYC ] = SumCyc;2205SumFuncs[ T_CYC ][ T_INT ] = SumCyc;2206SumFuncs[ T_CYC ][ T_INTPOS ] = SumCyc;2207SumFuncs[ T_CYC ][ T_INTNEG ] = SumCyc;2208SumFuncs[ T_CYC ][ T_RAT ] = SumCyc;2209DiffFuncs[ T_CYC ][ T_CYC ] = DiffCyc;2210DiffFuncs[ T_INT ][ T_CYC ] = DiffCyc;2211DiffFuncs[ T_INTPOS ][ T_CYC ] = DiffCyc;2212DiffFuncs[ T_INTNEG ][ T_CYC ] = DiffCyc;2213DiffFuncs[ T_RAT ][ T_CYC ] = DiffCyc;2214DiffFuncs[ T_CYC ][ T_INT ] = DiffCyc;2215DiffFuncs[ T_CYC ][ T_INTPOS ] = DiffCyc;2216DiffFuncs[ T_CYC ][ T_INTNEG ] = DiffCyc;2217DiffFuncs[ T_CYC ][ T_RAT ] = DiffCyc;2218ProdFuncs[ T_CYC ][ T_CYC ] = ProdCyc;2219ProdFuncs[ T_INT ][ T_CYC ] = ProdCycInt;2220ProdFuncs[ T_INTPOS ][ T_CYC ] = ProdCycInt;2221ProdFuncs[ T_INTNEG ][ T_CYC ] = ProdCycInt;2222ProdFuncs[ T_RAT ][ T_CYC ] = ProdCycInt;2223ProdFuncs[ T_CYC ][ T_INT ] = ProdCycInt;2224ProdFuncs[ T_CYC ][ T_INTPOS ] = ProdCycInt;2225ProdFuncs[ T_CYC ][ T_INTNEG ] = ProdCycInt;2226ProdFuncs[ T_CYC ][ T_RAT ] = ProdCycInt;22272228/* return success */2229return 0;2230}223122322233/****************************************************************************2234**2235*F InitLibrary( <module> ) . . . . . . . initialise library data structures2236*/2237static Int InitLibrary (2238StructInitInfo * module )2239{2240Obj * res; /* pointer into the result */2241UInt i; /* loop variable */22422243/* create the result buffer */2244TLS(ResultCyc) = NEW_PLIST( T_PLIST, 1024 );2245SET_LEN_PLIST( TLS(ResultCyc), 1024 );2246res = &(ELM_PLIST( TLS(ResultCyc), 1 ));2247for ( i = 0; i < 1024; i++ ) { res[i] = INTOBJ_INT(0); }22482249/* init filters and functions */2250InitGVarFiltsFromTable( GVarFilts );2251InitGVarAttrsFromTable( GVarAttrs );2252InitGVarOpersFromTable( GVarOpers );2253InitGVarFuncsFromTable( GVarFuncs );22542255/* return success */2256return 0;2257}225822592260/****************************************************************************2261**2262*F InitInfoCyc() . . . . . . . . . . . . . . . . . . table of init functions2263*/2264static StructInitInfo module = {2265MODULE_BUILTIN, /* type */2266"cyclotom", /* name */22670, /* revision entry of c file */22680, /* revision entry of h file */22690, /* version */22700, /* crc */2271InitKernel, /* initKernel */2272InitLibrary, /* initLibrary */22730, /* checkInit */22740, /* preSave */22750, /* postSave */22760 /* postRestore */2277};22782279StructInitInfo * InitInfoCyc ( void )2280{2281return &module;2282}228322842285/****************************************************************************2286**22872288*E cyclotom.c . . . . . . . . . . . . . . . . . . . . . . . . . . ends here2289*/229022912292