File: cyclotom.c

package info (click to toggle)
gap 4r4p12-2
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217 /**************************************************************************** ** *W cyclotom.c GAP source Martin Schoenert ** *H @(#)$Id: cyclotom.c,v 4.38.6.2 2005/04/27 15:40:04 gap Exp$ ** *Y Copyright (C) 1996, Lehrstuhl D fuer Mathematik, RWTH Aachen, Germany *Y (C) 1998 School Math and Comp. Sci., University of St. Andrews, Scotland *Y Copyright (C) 2002 The GAP Group ** ** 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 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 ** divisable 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 fuer Mathematik, RWTH Aachen, 1989 ** ** We represent a cyclotomic with terms, i.e., nonzero coefficients ** in the linear combination, by a bag of type 'T_CYC' with +1 subbags ** and +1 unsigned short integers. All the bag identifiers are stored at ** the beginning of the bag and all unsigned short 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()[ 0 ] )' gives you the order. The ** first unsigned short integer is unused (but reserved for future use :-). ** ** The other subbags and shorts 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 "system.h" /* Ints, UInts */ const char * Revision_cyclotom_c = "@(#)$Id: cyclotom.c,v 4.38.6.2 2005/04/27 15:40:04 gap Exp$"; #include "gasman.h" /* garbage collector */ #include "objects.h" /* objects */ #include "scanner.h" /* scanner */ #include "gap.h" /* error handling, initialisation */ #include "gvars.h" /* global variables */ #include "calls.h" /* generic call mechanism */ #include "opers.h" /* generic operations */ #include "ariths.h" /* basic arithmetic */ #include "bool.h" /* booleans */ #include "integer.h" /* integers */ #define INCLUDE_DECLARATION_PART #include "cyclotom.h" /* cyclotomics */ #undef INCLUDE_DECLARATION_PART #include "records.h" /* generic records */ #include "precord.h" /* plain records */ #include "lists.h" /* generic lists */ #include "plist.h" /* plain lists */ #include "string.h" /* strings */ #include "saveload.h" /* saving and loading */ /**************************************************************************** ** */ #define SIZE_CYC(cyc) (SIZE_OBJ(cyc) / (sizeof(Obj)+sizeof(UInt2))) #define COEFS_CYC(cyc) (ADDR_OBJ(cyc)) #define EXPOS_CYC(cyc,len) ((UInt2*)(ADDR_OBJ(cyc)+(len))) #define NOF_CYC(cyc) (COEFS_CYC(cyc)[0]) #define XXX_CYC(cyc,len) (EXPOS_CYC(cyc,len)[0]) /**************************************************************************** ** *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 )[]' 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 ** 'FunE'. ** ** 'LastNCyc' is the order of this primitive root. ** ** These values are used in 'FunE' 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 th place in 'ResultCyc' and then calling 'Cyclotomic'. */ Obj LastECyc; UInt LastNCyc; /**************************************************************************** ** *F TypeCyc( ) . . . . . . . . . . . . . . . . . type of a cyclotomic ** ** 'TypeCyc' returns the type of a cyclotomic. ** ** 'TypeCyc' is the function in 'TypeObjFuncs' for cyclotomics. */ Obj TYPE_CYC; Obj TypeCyc ( Obj cyc ) { return TYPE_CYC; } /**************************************************************************** ** *F PrintCyc( ) . . . . . . . . . . . . . . . . . . print a cyclotomic ** ** 'PrintCyc' prints the cyclotomic 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()^0', 'E()^1, etc. */ void PrintCyc ( Obj cyc ) { UInt n; /* order of the field */ UInt len; /* number of terms */ Obj * cfs; /* pointer to the coefficients */ UInt2 * exs; /* pointer to the exponents */ UInt i; /* loop variable */ n = INT_INTOBJ( NOF_CYC(cyc) ); len = SIZE_CYC(cyc); Pr("%>",0L,0L); for ( i = 1; i < len; i++ ) { /* get pointers, they can change during Pr */ cfs = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len); if ( cfs[i]==INTOBJ_INT(1) && exs[i]==0 ) Pr("1",0L,0L); else if ( cfs[i]==INTOBJ_INT(1) && exs[i]==1 && i==1 ) Pr("%>E(%d%<)",n,0L); else if ( cfs[i]==INTOBJ_INT(1) && exs[i]==1 ) Pr("%>+E(%d%<)",n,0L); else if ( cfs[i]==INTOBJ_INT(1) && i==1 ) Pr("%>E(%d)%>^%2<%d",n,(Int)exs[i]); else if ( cfs[i]==INTOBJ_INT(1) ) Pr("%>+E(%d)%>^%2<%d",n,(Int)exs[i]); else if ( LT(INTOBJ_INT(0),cfs[i]) && exs[i]==0 ) PrintObj(cfs[i]); else if ( LT(INTOBJ_INT(0),cfs[i]) && exs[i]==1 && i==1 ) { Pr("%>",0L,0L); PrintObj(cfs[i]); Pr("%>*%+",0L,0L); PrintObj(cfs[i]); Pr("%>*%",0L,0L); PrintObj(cfs[i]); Pr("%>*%^%2<%d",n,(Int)exs[i]); } else if ( LT(INTOBJ_INT(0),cfs[i]) ) { Pr("%>+",0L,0L); PrintObj(cfs[i]); Pr("%>*%^%2<%d",n,(Int)exs[i]); } else if ( cfs[i]==INTOBJ_INT(-1) && exs[i]==0 ) Pr("%>-%<1",0L,0L); else if ( cfs[i]==INTOBJ_INT(-1) && exs[i]==1 ) Pr("%>-E(%d%<)",n,0L); else if ( cfs[i]==INTOBJ_INT(-1) ) Pr("%>-E(%d)%>^%2<%d",n,(Int)exs[i]); else if ( exs[i]==0 ) PrintObj(cfs[i]); else if ( exs[i]==1 ) { Pr("%>",0L,0L); PrintObj(cfs[i]); Pr("%>*%",0L,0L); PrintObj(cfs[i]); Pr("%>*%^%2<%d",n,(Int)exs[i]); } } Pr("%<",0L,0L); } /**************************************************************************** ** *F EqCyc( , ) . . . . . . . . . test if two cyclotomics are equal ** ** 'EqCyc' returns 'true' if the two cyclotomics and are equal ** and 'false' otherwise. ** ** 'EqCyc' is pretty simple because every cyclotomic has a unique ** representation, so we just have to compare the terms. */ Int EqCyc ( Obj opL, Obj opR ) { UInt len; /* number of terms */ Obj * cfl; /* ptr to coeffs of left operand */ UInt2 * exl; /* ptr to expnts of left operand */ Obj * cfr; /* ptr to coeffs of right operand */ UInt2 * 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 0L; /* compare the number of terms */ if ( SIZE_CYC(opL) != SIZE_CYC(opR) ) return 0L; /* compare the cyclotomics termwise */ len = SIZE_CYC(opL); cfl = COEFS_CYC(opL); cfr = COEFS_CYC(opR); exl = EXPOS_CYC(opL,len); exr = EXPOS_CYC(opR,len); for ( i = 1; i < len; i++ ) { if ( exl[i] != exr[i] ) return 0L; else if ( ! EQ(cfl[i],cfr[i]) ) return 0L; } /* all terms are equal */ return 1L; } /**************************************************************************** ** *F LtCyc( , ) . . . . test if one cyclotomic is less than another ** ** 'LtCyc' returns 'true' if the cyclotomic is less than the ** cyclotomic 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. */ Int LtCyc ( Obj opL, Obj opR ) { UInt lel; /* nr of terms of left operand */ Obj * cfl; /* ptr to coeffs of left operand */ UInt2 * exl; /* ptr to expnts of left operand */ UInt ler; /* nr of terms of right operand */ Obj * cfr; /* ptr to coeffs of right operand */ UInt2 * 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 1L; else return 0L; } /* compare the cyclotomics termwise */ lel = SIZE_CYC(opL); ler = SIZE_CYC(opR); cfl = COEFS_CYC(opL); cfr = COEFS_CYC(opR); exl = EXPOS_CYC(opL,lel); exr = 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 agains 0 */ if ( lel < ler ) return LT( INTOBJ_INT(0), cfr[i] ); else if ( ler < lel ) return LT( cfl[i], INTOBJ_INT(0) ); else return 0L; } Int LtCycYes ( Obj opL, Obj opR ) { return 1L; } Int LtCycNot ( Obj opL, Obj opR ) { return 0L; } /**************************************************************************** ** *F ConvertToBase() . . . . . . convert a cyclotomic into the base, local ** ** 'ConvertToBase' converts the cyclotomic 'ResultCyc' from the cyclotomic ** field of 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 divisable by $p$. But $n/p$ ** is not divisable by $p$, so neither of the exponent $k*n/p+i, k=1..p-1$ ** is divisable 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 divisable 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. */ 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 = &(ELM_PLIST( ResultCyc, 1 )); 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 = &(ELM_PLIST( ResultCyc, 1 )); } 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 = &(ELM_PLIST( ResultCyc, 1 )); } 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 = &(ELM_PLIST( ResultCyc, 1 )); } 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 = &(ELM_PLIST( ResultCyc, 1 )); } res[l%n] = sum; } res[k] = INTOBJ_INT(0); } } } } /* notify Gasman */ CHANGED_BAG( ResultCyc ); } /**************************************************************************** ** *F Cyclotomic(,) . . . . . . . . . . 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'. must be the order of the primitive root in which ** written. ** ** must be a divisor of $n$ and gives a hint about possible subfields. ** If a prime $p$ divides 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 can not 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$. */ Obj Cyclotomic ( UInt n, UInt m ) { Obj cyc; /* cyclotomic, result */ UInt len; /* number of terms */ Obj * cfs; /* pointer to the coefficients */ UInt2 * 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; /* rember last n, dont recompute: */ static UInt phi; /* Euler phi(n) */ static UInt 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 = &(ELM_PLIST( ResultCyc, 1 )); /* 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 divisable 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 = 1; nrp = 0; for ( p = 2; p <= k; p++ ) { if ( k % p == 0 ) { phi = phi * (p-1) / p; if ( k % (p*p) == 0 ) isSqfree = 0; 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 ); res[0] = DIFF( INTOBJ_INT(0), cof ); } 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 divisiable */ /* $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]; res[i] = INTOBJ_INT( - INT_INTOBJ(cof) ); if ( ! IS_INTOBJ(cof) || (cof == INTOBJ_INT(-(1L<<28))) ) { CHANGED_BAG( ResultCyc ); cof = DIFF( INTOBJ_INT(0), cof ); res = &(ELM_PLIST( ResultCyc, 1 )); res[i] = 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(UInt2)) ); cfs = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len+1); cfs[0] = INTOBJ_INT(n); exs[0] = 0; k = 1; res = &(ELM_PLIST( ResultCyc, 1 )); 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 the result */ return cyc; } /**************************************************************************** ** *F SumCyc( , ) . . . . . . . . . . . . . sum of two cyclotomics ** ** 'SumCyc' returns the sum of the two cyclotomics and . ** 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. */ 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 */ Obj * cfs; /* pointer to the coefficients */ UInt2 * 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; } /* 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 = nl; while ( n % nr != 0 ) n += nl; ml = n / nl; mr = n / nr; /* make sure that the result back is large enough */ if ( LEN_PLIST(ResultCyc) < n ) { GROW_PLIST( ResultCyc, n ); SET_LEN_PLIST( ResultCyc, n ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < n; i++ ) { res[i] = INTOBJ_INT(0); } } /* copy the left operand into the result */ if ( TNUM_OBJ(opL) != T_CYC ) { res = &(ELM_PLIST( ResultCyc, 1 )); res[0] = opL; CHANGED_BAG( ResultCyc ); } else { len = SIZE_CYC(opL); cfs = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = &(ELM_PLIST( ResultCyc, 1 )); sum = SUM( res[0], opR ); res = &(ELM_PLIST( ResultCyc, 1 )); res[0] = sum; CHANGED_BAG( ResultCyc ); } else { len = SIZE_CYC(opR); cfs = COEFS_CYC(opR); exs = EXPOS_CYC(opR,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = COEFS_CYC(opR); exs = EXPOS_CYC(opR,len); res = &(ELM_PLIST( ResultCyc, 1 )); } 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( ) . . . . . . . . . . . . . . . . . . zero of a cyclotomic ** ** 'ZeroCyc' returns the additive neutral element of the cyclotomic . */ Obj ZeroCyc ( Obj op ) { return INTOBJ_INT( 0L ); } /**************************************************************************** ** *F AInvCyc( ) . . . . . . . . . . . . additive inverse of a cyclotomic ** ** 'AInvCyc' returns the additive inverse element of the cyclotomic . */ Obj AInvCyc ( Obj op ) { Obj res; /* inverse, result */ UInt len; /* number of terms */ Obj * cfs; /* ptr to coeffs of left operand */ UInt2 * exs; /* ptr to expnts of left operand */ Obj * cfp; /* ptr to coeffs of product */ UInt2 * 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(UInt2)) ); NOF_CYC(res) = NOF_CYC(op); len = SIZE_CYC(op); cfs = COEFS_CYC(op); cfp = COEFS_CYC(res); exs = EXPOS_CYC(op,len); exp = EXPOS_CYC(res,len); for ( i = 1; i < len; i++ ) { prd = INTOBJ_INT( - INT_INTOBJ(cfs[i]) ); if ( ! IS_INTOBJ( cfs[i] ) || cfs[i] == INTOBJ_INT(-(1L<<28)) ) { CHANGED_BAG( res ); prd = AINV( cfs[i] ); cfs = COEFS_CYC(op); cfp = COEFS_CYC(res); exs = EXPOS_CYC(op,len); exp = EXPOS_CYC(res,len); } cfp[i] = prd; exp[i] = exs[i]; } CHANGED_BAG( res ); /* return the result */ return res; } /**************************************************************************** ** *F DiffCyc( , ) . . . . . . . . . . difference of two cyclotomics ** ** 'DiffCyc' returns the difference of the two cyclotomic and . ** 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. */ 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 */ Obj * cfs; /* pointer to the coefficients */ UInt2 * 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 = nl; while ( n % nr != 0 ) n += nl; ml = n / nl; mr = n / nr; /* make sure that the result back is large enough */ if ( LEN_PLIST(ResultCyc) < n ) { GROW_PLIST( ResultCyc, n ); SET_LEN_PLIST( ResultCyc, n ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < n; i++ ) { res[i] = INTOBJ_INT(0); } } /* copy the left operand into the result */ if ( TNUM_OBJ(opL) != T_CYC ) { res = &(ELM_PLIST( ResultCyc, 1 )); res[0] = opL; CHANGED_BAG( ResultCyc ); } else { len = SIZE_CYC(opL); cfs = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = &(ELM_PLIST( ResultCyc, 1 )); sum = DIFF( res[0], opR ); res = &(ELM_PLIST( ResultCyc, 1 )); res[0] = sum; CHANGED_BAG( ResultCyc ); } else { len = SIZE_CYC(opR); cfs = COEFS_CYC(opR); exs = EXPOS_CYC(opR,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = COEFS_CYC(opR); exs = EXPOS_CYC(opR,len); res = &(ELM_PLIST( ResultCyc, 1 )); } 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( , ) . . . 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 wich 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. */ Obj ProdCycInt ( Obj opL, Obj opR ) { Obj hdP; /* product, result */ UInt len; /* number of terms */ Obj * cfs; /* ptr to coeffs of left operand */ UInt2 * exs; /* ptr to expnts of left operand */ Obj * cfp; /* ptr to coeffs of product */ UInt2 * 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 ( TNUM_OBJ(opR) == T_INT ) { hdP = NewBag( T_CYC, SIZE_CYC(opL) * (sizeof(Obj)+sizeof(UInt2)) ); NOF_CYC(hdP) = NOF_CYC(opL); len = SIZE_CYC(opL); cfs = COEFS_CYC(opL); cfp = COEFS_CYC(hdP); exs = EXPOS_CYC(opL,len); 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 = COEFS_CYC(opL); cfp = COEFS_CYC(hdP); exs = EXPOS_CYC(opL,len); exp = EXPOS_CYC(hdP,len); } cfp[i] = prd; exp[i] = exs[i]; } CHANGED_BAG( hdP ); } /* otherwise multiply every coefficent */ else { hdP = NewBag( T_CYC, SIZE_CYC(opL) * (sizeof(Obj)+sizeof(UInt2)) ); NOF_CYC(hdP) = NOF_CYC(opL); len = SIZE_CYC(opL); cfs = COEFS_CYC(opL); cfp = COEFS_CYC(hdP); exs = EXPOS_CYC(opL,len); exp = EXPOS_CYC(hdP,len); for ( i = 1; i < len; i++ ) { CHANGED_BAG( hdP ); prd = PROD( cfs[i], opR ); cfs = COEFS_CYC(opL); cfp = COEFS_CYC(hdP); exs = EXPOS_CYC(opL,len); exp = EXPOS_CYC(hdP,len); cfp[i] = prd; exp[i] = exs[i]; } CHANGED_BAG( hdP ); } /* return the result */ return hdP; } /**************************************************************************** ** *F ProdCyc( , ) . . . . . . . . . . . product of two cyclotomics ** ** 'ProdCyc' returns the product of the two cyclotomics and . ** 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. */ 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 */ Obj * cfs; /* pointer to the coefficients */ UInt2 * 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 = nl; while ( n % nr != 0 ) n += nl; ml = n / nl; mr = n / nr; /* make sure that the result bag is large enough */ if ( LEN_PLIST(ResultCyc) < n ) { GROW_PLIST( ResultCyc, n ); SET_LEN_PLIST( ResultCyc, n ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < n; i++ ) { res[i] = INTOBJ_INT(0); } } /* loop over the terms of the right operand */ for ( k = 1; k < SIZE_CYC(opR); k++ ) { c = COEFS_CYC(opR)[k]; e = (mr * 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 = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); } 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 = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); } res[(e+exs[i]*ml)%n] = sum; } CHANGED_BAG( ResultCyc ); } /* if the coefficient is a small integer use immediate operations */ else if ( TNUM_OBJ(c) == T_INT ) { len = SIZE_CYC(opL); cfs = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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 = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); sum = SUM( res[(e+exs[i]*ml)%n], prd ); cfs = COEFS_CYC(opL); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); } 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 = COEFS_CYC(opL); prd = PROD( cfs[i], c ); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); sum = SUM( res[(e+exs[i]*ml)%n], prd ); exs = EXPOS_CYC(opL,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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( ) . . . . . . . . . . . . . . . . . . . one of a cyclotomic ** ** 'OneCyc' returns the multiplicative neutral element of the cyclotomic ** . */ Obj OneCyc ( Obj op ) { return INTOBJ_INT( 1L ); } /**************************************************************************** ** *F InvCyc( ) . . . . . . . . . . . . . . . . . inverse of a cyclotomic ** ** 'InvCyc' returns the multiplicative inverse element of the cyclotomic ** . ** ** 'InvCyc' computes the inverse of by computing the product $prd$ of ** nontrivial Galois conjugates of . 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? */ 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 */ Obj * cfs; /* pointer to the coefficients */ UInt2 * 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 */ 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 = COEFS_CYC(op); exs = EXPOS_CYC(op,len); res = &(ELM_PLIST( ResultCyc, 1 )); 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( , ) . . . . . . . . . . . . . . power of a cyclotomic ** ** 'PowCyc' returns the th, which must be an integer, power of the ** cyclotomic . The left operand may also be an integer or a rational. */ Obj PowCyc ( Obj opL, Obj opR ) { Obj pow; /* power (result) */ Int exp; /* exponent (right operand) */ UInt n; /* order of the field */ Obj * res; /* pointer into the result */ 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 th position and convert */ else if ( opL == LastECyc ) { exp = (exp % LastNCyc + LastNCyc) % LastNCyc; SET_ELM_PLIST( ResultCyc, exp, 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( COEFS_CYC(opL)[1], opR ); i = EXPOS_CYC(opL,2)[1]; res = &(ELM_PLIST( ResultCyc, 1 )); res[((exp*i)%n+n)%n] = pow; CHANGED_BAG( ResultCyc ); ConvertToBase( n ); pow = Cyclotomic( n, 1 ); } /* otherwise compute the power with repeated squaring */ else { /* if neccessary 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 the result */ return pow; } /**************************************************************************** ** *F FuncE( , ) . . . . . . . . . . . . create a new primitive root ** ** 'FuncE' implements the internal function 'E'. ** ** 'E( )' ** ** 'E' returns a primitive root of order , which must be a positive ** integer, represented as cyclotomic. */ Obj EOper; Obj FuncE ( Obj self, Obj n ) { UInt i; /* loop variable */ Obj * res; /* pointer into result bag */ /* do full operation */ if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(n) ) { return DoOperation1Args( self, n ); } /* get and check the argument */ while ( TNUM_OBJ(n) != T_INT || INT_INTOBJ(n) <= 0 ) { n = ErrorReturnObj( "E: must be a positive integer (not a %s)", (Int)TNAM_OBJ(n), 0L, "you can replace via 'return ;'" ); } /* 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); if ( LEN_PLIST(ResultCyc) < LastNCyc ) { GROW_PLIST( ResultCyc, LastNCyc ); SET_LEN_PLIST( ResultCyc, LastNCyc ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < LastNCyc; i++ ) { res[i] = INTOBJ_INT(0); } } res = &(ELM_PLIST( ResultCyc, 1 )); res[1] = INTOBJ_INT(1); CHANGED_BAG( ResultCyc ); ConvertToBase( LastNCyc ); LastECyc = Cyclotomic( LastNCyc, 1 ); } /* return the root */ return LastECyc; } /**************************************************************************** ** *F FuncIS_CYC( , ) . . . . . . test if an object is a cyclomtic ** ** 'FuncIS_CYC' implements the internal function 'IsCyc'. ** ** 'IsCyc( )' ** ** 'IsCyc' returns 'true' if the value is a cyclotomic and 'false' ** otherwise. */ Obj IsCycFilt; Obj FuncIS_CYC ( Obj self, Obj val ) { /* return 'true' if is a cyclotomic and 'false' otherwise */ if ( TNUM_OBJ(val) == T_CYC || TNUM_OBJ(val) == T_INT || TNUM_OBJ(val) == T_RAT || TNUM_OBJ(val) == T_INTPOS || TNUM_OBJ(val) == T_INTNEG ) return True; else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) { return False; } else { return DoFilter( self, val ); } } /**************************************************************************** ** *F FuncIS_CYC_INT( , ) test if an object is a cyclomtic integer ** ** 'FuncIS_CYC_INT' implements the internal function 'IsCycInt'. ** ** 'IsCycInt( )' ** ** 'IsCycInt' returns 'true' if the value is a cyclotomic integer and ** 'false' otherwise. ** ** 'IsCycInt' relies on the fact that the base is an integral base. */ Obj IsCycIntOper; Obj FuncIS_CYC_INT ( Obj self, Obj val ) { UInt len; /* number of terms */ Obj * cfs; /* pointer to the coefficients */ UInt i; /* loop variable */ /* return 'true' if is a cyclotomic integer and 'false' otherwise*/ if ( TNUM_OBJ(val) == T_INT || TNUM_OBJ(val) == T_INTPOS || TNUM_OBJ(val) == T_INTNEG ) { 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 FuncCONDUCTOR( , ) . . . . . . . . . . . . N of a cyclotomic ** ** 'FuncCONDUCTOR' implements the internal function 'Conductor'. ** ** 'Conductor( )' ** ** 'Conductor' returns the N of the cyclotomic , i.e., the order of the ** roots of which is written as a linear combination. */ Obj ConductorAttr; Obj FuncCONDUCTOR ( 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 ); } /* check the argument */ while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT && TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG && TNUM_OBJ(cyc) != T_CYC && ! IS_SMALL_LIST(cyc) ) { cyc = ErrorReturnObj( "Conductor: must be a cyclotomic or a small list (not a %s)", (Int)TNAM_OBJ(cyc), 0L, "you can replace via 'return ;'" ); } /* handle cyclotomics */ if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT || TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) { 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 ); while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT && TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG && TNUM_OBJ(cyc) != T_CYC ) { cyc = ErrorReturnObj( "Conductor: [%d] must be a cyclotomic (not a %s)", (Int)i, (Int)TNAM_OBJ(cyc), "you can replace the list element with via 'return ;'" ); } if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT || TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) { 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( , ) . . . . . . coefficients of a cyclotomic ** ** 'FuncCOEFFS_CYC' implements the internal function 'COEFFSCYC'. ** ** 'COEFFSCYC( )' ** ** 'COEFFSCYC' returns a list of the coefficients of the cyclotomic . ** The list has length if is the order of the primitive root $e_n$ ** of which is written as a linear combination. The th element of ** the list is the coefficient of $e_l^{i-1}$. */ Obj CoeffsCycOper; Obj FuncCOEFFS_CYC ( Obj self, Obj cyc ) { Obj list; /* list of coefficients, result */ UInt n; /* order of field */ UInt len; /* number of terms */ Obj * cfs; /* pointer to the coefficients */ UInt2 * exs; /* pointer to the exponents */ UInt i; /* loop variable */ /* do full operation */ if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ(cyc) ) { return DoOperation1Args( self, cyc ); } /* check the argument */ while ( TNUM_OBJ(cyc) != T_INT && TNUM_OBJ(cyc) != T_RAT && TNUM_OBJ(cyc) != T_INTPOS && TNUM_OBJ(cyc) != T_INTNEG && TNUM_OBJ(cyc) != T_CYC ) { cyc = ErrorReturnObj( "COEFFSCYC: must be a cyclotomic (not a %s)", (Int)TNAM_OBJ(cyc), 0L, "you can replace via 'return ;'" ); } /* if is rational just put it in a list of length 1 */ if ( TNUM_OBJ(cyc) == T_INT || TNUM_OBJ(cyc) == T_RAT || TNUM_OBJ(cyc) == T_INTPOS || TNUM_OBJ(cyc) == T_INTNEG ) { list = NEW_PLIST( T_PLIST, 1 ); SET_LEN_PLIST( list, 1 ); SET_ELM_PLIST( list, 1, 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 = COEFS_CYC(cyc); exs = 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 the result */ return list; } /**************************************************************************** ** *F FuncGALOIS_CYC( , , ) image of a cyc. under galois aut ** ** 'FuncGALOIS_CYC' implements the internal function 'GaloisCyc'. ** ** 'GaloisCyc( , )' ** ** 'GaloisCyc' computes the image of the cyclotomic under the galois ** automorphism given by , which must be an integer. ** ** The galois automorphism is the mapping that takes $e_n$ to $e_n^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. */ Obj GaloisCycOper; 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 */ Obj * cfs; /* pointer to the coefficients */ UInt2 * exs; /* pointer to the exponents */ Obj * res; /* pointer to the result */ UInt i; /* loop variable */ UInt tnumord, tnumcyc; /* do full operation for any but standard arguments */ tnumord = TNUM_OBJ(ord); tnumcyc = TNUM_OBJ(cyc); if ( FIRST_EXTERNAL_TNUM <= tnumcyc || FIRST_EXTERNAL_TNUM <= tnumord || (tnumord != T_INT && tnumord != T_INTNEG && tnumord != T_INTPOS) || ( tnumcyc != T_INT && tnumcyc != T_RAT && tnumcyc != T_INTPOS && tnumcyc != T_INTNEG && tnumcyc != T_CYC ) ) { return DoOperation2Args( self, cyc, ord ); } /* get and check */ if ( TNAM_OBJ(ord) == T_INT ) { o = INT_INTOBJ(ord); } else { ord = MOD( ord, FuncCONDUCTOR( 0, cyc ) ); o = INT_INTOBJ(ord); } /* every galois automorphism fixes the rationals */ if ( tnumcyc == T_INT || tnumcyc == T_RAT || tnumcyc == T_INTPOS || tnumcyc == T_INTNEG ) { 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 into the range 0..n-1, compute the gcd of and */ o = (o % n + n) % n; gcd = n; s = o; while ( s != 0 ) { t = s; s = gcd % s; gcd = t; } /* if = 1 just return */ if ( o == 1 ) { gal = cyc; } /* if == 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 == 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 = COEFS_CYC(cyc); exs = 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 = COEFS_CYC(cyc); exs = 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 = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len); } gal = sum; } } } /* if is prime to (automorphism) permute the coefficients */ else if ( gcd == 1 ) { /* permute the coefficients */ len = SIZE_CYC(cyc); cfs = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 1; i < len; i++ ) { res[exs[i]*o%n] = cfs[i]; } CHANGED_BAG( ResultCyc ); /* if n is squarefree conversion and reduction are unnecessary */ if ( n < sqr*sqr || (o == n-1 && n % 2 != 0) ) { gal = Cyclotomic( n, n ); } else { ConvertToBase( n ); gal = Cyclotomic( n, 1 ); } } /* if is not prime to (endomorphism) compute it the hard way */ else { /* multiple roots may be mapped to the same root, add the coeffs */ len = SIZE_CYC(cyc); cfs = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 1; i < len; i++ ) { if ( ! ARE_INTOBJS( res[exs[i]*o%n], cfs[i] ) || ! SUM_INTOBJS( sum, res[exs[i]*o%n], cfs[i] ) ) { CHANGED_BAG( ResultCyc ); sum = SUM( res[exs[i]*o%n], cfs[i] ); cfs = COEFS_CYC(cyc); exs = EXPOS_CYC(cyc,len); res = &(ELM_PLIST( ResultCyc, 1 )); } res[exs[i]*o%n] = sum; } CHANGED_BAG( ResultCyc ); /* if n is squarefree conversion and reduction are unnecessary */ if ( n < sqr*sqr ) { gal = Cyclotomic( n, 1 ); /*N?*/ } else { ConvertToBase( n ); gal = Cyclotomic( n, 1 ); } } /* return the result */ return gal; } /**************************************************************************** ** *F FuncCycList( , ) . . . . . . . . . . . . create a cyclotomic ** ** 'FuncCycList' implements the internal function 'CycList'. ** ** 'CycList( )' ** ** 'CycList' returns the cyclotomic described by the list ** of rationals. */ Obj CycListOper; Obj FuncCycList ( Obj self, Obj list ) { UInt i; /* loop variable */ Obj * res; /* pointer into result bag */ Obj val; /* one list entry */ UInt n; /* length of the given list */ /* do full operation */ if ( FIRST_EXTERNAL_TNUM <= TNUM_OBJ( list ) ) { return DoOperation1Args( self, list ); } /* get and check the argument */ if ( ! ( T_PLIST <= TNUM_OBJ( list ) && TNUM_OBJ( list ) <= LAST_PLIST_TNUM ) || ! IS_DENSE_LIST( list ) ) { ErrorQuit( "CycList: must be a dense plain list (not a %s)", (Int)TNAM_OBJ( list ), 0L ); } /* enlarge the buffer if necessary */ n = LEN_PLIST( list ); if ( LEN_PLIST(ResultCyc) < n ) { GROW_PLIST( ResultCyc, n ); SET_LEN_PLIST( ResultCyc, n ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < n; i++ ) { res[i] = INTOBJ_INT(0); } } /* transfer the coefficients into the buffer */ res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < n; i++ ) { val = ELM_PLIST( list, i+1 ); if ( ! ( TNUM_OBJ(val) == T_INT || TNUM_OBJ(val) == T_RAT || TNUM_OBJ(val) == T_INTPOS || TNUM_OBJ(val) == T_INTNEG ) ) { ErrorQuit( "CycList: each entry must be a rational (not a %s)", (Int)TNAM_OBJ( val ), 0L ); } res[i] = val; } /* return the base reduced packed cyclotomic */ CHANGED_BAG( ResultCyc ); ConvertToBase( n ); return Cyclotomic( n, 1 ); } /**************************************************************************** ** *F MarkCycSubBags( ) . . . . . . . . marking function for cyclotomics ** ** 'MarkCycSubBags' is the marking function for bags of type 'T_CYC'. */ void MarkCycSubBags ( Obj cyc ) { Obj * cfs; /* pointer to coeffs */ UInt i; /* loop variable */ /* mark everything */ cfs = COEFS_CYC( cyc ); for ( i = SIZE_CYC(cyc); 0 < i; i-- ) { MARK_BAG( cfs[i-1] ); } } /**************************************************************************** ** *F SaveCyc() . . . . . . . . . . . . . . . . . . . . . . save a cyclotyomic ** ** We do not save the XXX_CYC field, since it is not used. */ void SaveCyc ( Obj cyc ) { UInt len, i; Obj *coefs; UInt2 *expos; len = SIZE_CYC(cyc); coefs = COEFS_CYC(cyc); for (i = 0; i < len; i++) SaveSubObj(*coefs++); expos = EXPOS_CYC(cyc,len); expos++; /*Skip past the XXX */ for (i = 1; i < len; i++) SaveUInt2(*expos++); return; } /**************************************************************************** ** *F LoadCyc() . . . . . . . . . . . . . . . . . . . . . . load a cyclotyomic ** ** We do not load the XXX_CYC field, since it is not used. */ void LoadCyc ( Obj cyc ) { UInt len, i; Obj *coefs; UInt2 *expos; len = SIZE_CYC(cyc); coefs = COEFS_CYC(cyc); for (i = 0; i < len; i++) *coefs++ = LoadSubObj(); expos = EXPOS_CYC(cyc,len); expos++; /*Skip past the XXX */ for (i = 1; i < len; i++) *expos++ = LoadUInt2(); return; } /**************************************************************************** ** ** *F * * * * * * * * * * * * * initialize package * * * * * * * * * * * * * * * */ /**************************************************************************** ** *V GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export */ static StructGVarFilt GVarFilts [] = { { "IS_CYC", "obj", &IsCycFilt, FuncIS_CYC, "src/cyclotom.c:IS_CYC" }, { 0 } }; /**************************************************************************** ** *V GVarAttrs . . . . . . . . . . . . . . . . . list of attributes to export */ static StructGVarAttr GVarAttrs [] = { { "CONDUCTOR", "cyc", &ConductorAttr, FuncCONDUCTOR, "src/cyclotom.c:CONDUCTOR" }, { 0 } }; /**************************************************************************** ** *V GVarOpers . . . . . . . . . . . . . . . . . list of operations to export */ static StructGVarOper GVarOpers [] = { { "E", 1, "n", &EOper, FuncE, "src/cyclotom.c:E" }, { "IS_CYC_INT", 1, "obj", &IsCycIntOper, FuncIS_CYC_INT, "src/cyclotom.c:IS_CYC_INT" }, { "COEFFS_CYC", 1, "cyc", &CoeffsCycOper, FuncCOEFFS_CYC, "src/cyclotom.c:COEFFS_CYC" }, { "GALOIS_CYC", 2, "cyc, n", &GaloisCycOper, FuncGALOIS_CYC, "src/cyclotom.c:GALOIS_CYC" }, { "CYC_LIST", 1, "list", &CycListOper, FuncCycList, "src/cyclotom.c:CYC_LIST" }, { 0 } }; /**************************************************************************** ** *F InitKernel( ) . . . . . . . . initialise kernel data structures */ static Int InitKernel ( StructInitInfo * module ) { /* install the marking function */ InfoBags[ T_CYC ].name = "cyclotomic"; InitMarkFuncBags( T_CYC, MarkCycSubBags ); /* create the result buffer */ InitGlobalBag( &ResultCyc , "src/cyclotom.c:ResultCyc" ); /* tell Gasman about the place were we remember the primitive root */ InitGlobalBag( &LastECyc, "src/cyclotom.c:LastECyc" ); /* install the kind function */ ImportGVarFromLibrary( "TYPE_CYC", &TYPE_CYC ); TypeObjFuncs[ T_CYC ] = TypeCyc; /* init filters and functions */ InitHdlrFiltsFromTable( GVarFilts ); InitHdlrAttrsFromTable( GVarAttrs ); InitHdlrOpersFromTable( GVarOpers ); /* and the saving function */ SaveObjFuncs[ T_CYC ] = SaveCyc; LoadObjFuncs[ T_CYC ] = LoadCyc; /* install the evaluation and print function */ PrintObjFuncs[ T_CYC ] = PrintCyc; /* install the comparison methods */ EqFuncs[ T_CYC ][ T_CYC ] = EqCyc; LtFuncs[ T_CYC ][ T_CYC ] = LtCyc; LtFuncs[ T_INT ][ T_CYC ] = LtCycYes; LtFuncs[ T_INTPOS ][ T_CYC ] = LtCycYes; LtFuncs[ T_INTNEG ][ T_CYC ] = LtCycYes; LtFuncs[ T_RAT ][ T_CYC ] = LtCycYes; LtFuncs[ T_CYC ][ T_INT ] = LtCycNot; LtFuncs[ T_CYC ][ T_INTPOS ] = LtCycNot; LtFuncs[ T_CYC ][ T_INTNEG ] = LtCycNot; LtFuncs[ T_CYC ][ T_RAT ] = LtCycNot; /* install the unary arithmetic methods */ ZeroFuncs[ T_CYC ] = ZeroCyc; ZeroMutFuncs[ T_CYC ] = ZeroCyc; AInvFuncs[ T_CYC ] = AInvCyc; AInvMutFuncs[ T_CYC ] = AInvCyc; OneFuncs [ T_CYC ] = OneCyc; OneMutFuncs [ T_CYC ] = OneCyc; InvFuncs [ T_CYC ] = InvCyc; InvMutFuncs [ T_CYC ] = InvCyc; /* install the arithmetic methods */ SumFuncs[ T_CYC ][ T_CYC ] = SumCyc; SumFuncs[ T_INT ][ T_CYC ] = SumCyc; SumFuncs[ T_INTPOS ][ T_CYC ] = SumCyc; SumFuncs[ T_INTNEG ][ T_CYC ] = SumCyc; SumFuncs[ T_RAT ][ T_CYC ] = SumCyc; SumFuncs[ T_CYC ][ T_INT ] = SumCyc; SumFuncs[ T_CYC ][ T_INTPOS ] = SumCyc; SumFuncs[ T_CYC ][ T_INTNEG ] = SumCyc; SumFuncs[ T_CYC ][ T_RAT ] = SumCyc; DiffFuncs[ T_CYC ][ T_CYC ] = DiffCyc; DiffFuncs[ T_INT ][ T_CYC ] = DiffCyc; DiffFuncs[ T_INTPOS ][ T_CYC ] = DiffCyc; DiffFuncs[ T_INTNEG ][ T_CYC ] = DiffCyc; DiffFuncs[ T_RAT ][ T_CYC ] = DiffCyc; DiffFuncs[ T_CYC ][ T_INT ] = DiffCyc; DiffFuncs[ T_CYC ][ T_INTPOS ] = DiffCyc; DiffFuncs[ T_CYC ][ T_INTNEG ] = DiffCyc; DiffFuncs[ T_CYC ][ T_RAT ] = DiffCyc; ProdFuncs[ T_CYC ][ T_CYC ] = ProdCyc; ProdFuncs[ T_INT ][ T_CYC ] = ProdCycInt; ProdFuncs[ T_INTPOS ][ T_CYC ] = ProdCycInt; ProdFuncs[ T_INTNEG ][ T_CYC ] = ProdCycInt; ProdFuncs[ T_RAT ][ T_CYC ] = ProdCycInt; ProdFuncs[ T_CYC ][ T_INT ] = ProdCycInt; ProdFuncs[ T_CYC ][ T_INTPOS ] = ProdCycInt; ProdFuncs[ T_CYC ][ T_INTNEG ] = ProdCycInt; ProdFuncs[ T_CYC ][ T_RAT ] = ProdCycInt; /* return success */ return 0; } /**************************************************************************** ** *F InitLibrary( ) . . . . . . . initialise library data structures */ static Int InitLibrary ( StructInitInfo * module ) { Obj * res; /* pointer into the result */ UInt i; /* loop variable */ /* create the result buffer */ ResultCyc = NEW_PLIST( T_PLIST, 1024 ); SET_LEN_PLIST( ResultCyc, 1024 ); res = &(ELM_PLIST( ResultCyc, 1 )); for ( i = 0; i < 1024; i++ ) { res[i] = INTOBJ_INT(0); } /* init filters and functions */ InitGVarFiltsFromTable( GVarFilts ); InitGVarAttrsFromTable( GVarAttrs ); InitGVarOpersFromTable( GVarOpers ); /* return success */ return 0; } /**************************************************************************** ** *F InitInfoCyc() . . . . . . . . . . . . . . . . . . table of init functions */ static StructInitInfo module = { MODULE_BUILTIN, /* type */ "cyclotom", /* name */ 0, /* revision entry of c file */ 0, /* revision entry of h file */ 0, /* version */ 0, /* crc */ InitKernel, /* initKernel */ InitLibrary, /* initLibrary */ 0, /* checkInit */ 0, /* preSave */ 0, /* postSave */ 0 /* postRestore */ }; StructInitInfo * InitInfoCyc ( void ) { module.revision_c = Revision_cyclotom_c; module.revision_h = Revision_cyclotom_h; FillInVersion( &module ); return &module; } /**************************************************************************** ** *E cyclotom.c . . . . . . . . . . . . . . . . . . . . . . . . . . ends here */