Math Solutions 1995 October

home *** CD-ROM | disk | FTP | other *** search

/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / integer.g < prev next >

Wrap

Text File | 1993-05-05 | 44.7 KB | 1,255 lines

############################################################################# ## #A integer.g GAP library Martin Schoenert #A & Alice Niemeyer #A & Werner Nickel #A & Alex Wegner ## #A @(#)$Id: integer.g,v 3.24 1993/02/15 09:21:28 fceller Rel $ ## #Y Copyright 1990-1992, Lehrstuhl D fuer Mathematik, RWTH Aachen, Germany ## ## This file contains those functions that mainly deal with integers. ## #H $Log: integer.g,v $ #H Revision 3.24 1993/02/15 09:21:28 fceller #H fixed 'Integers.units' and 'IntegerOps.IsIrreducible' #H #H Revision 3.23 1993/02/10 21:20:25 martin #H fixed the error message of 'ChineseRemainder' #H #H Revision 3.22 1992/12/16 19:47:27 martin #H replaced quoted record names with escaped ones #H #H Revision 3.21 1992/09/10 10:39:33 martin #H fixed the problem in 'IsPrimeInt' #H #H Revision 3.20 1992/06/04 07:41:10 fceller #H added some primes occuring in q^n-1 #H #H Revision 3.19 1992/03/19 18:16:39 martin #H added 'IntegersOps.Field' #H #H Revision 3.18 1992/02/06 11:47:31 martin #H removed 'InverseMod' and added 'QuotientMod' #H #H Revision 3.17 1991/12/27 15:06:20 martin #H moved 'IsPrimeInt', 'FactorsInt', etc. to 'integer.g' #H #H Revision 3.16 1991/12/27 15:01:19 martin #H moved 'Rationals' to 'rational.g' #H #H Revision 3.15 1991/12/04 10:57:26 martin #H added 'isCyclotomicField' to 'Rationals' #H #H Revision 3.14 1991/10/24 11:58:52 martin #H added 'LcmInt' again #H #H Revision 3.13 1991/10/24 11:33:29 martin #H changed package for new domain concept with inheritance #H #H Revision 3.12 1991/08/08 16:00:00 martin #H added functions for rationals #H #H Revision 3.11 1991/08/08 15:30:00 martin #H changed the integer package for domains #H #H Revision 3.10 1991/06/28 14:00:00 martin #H moved 'Maximum' and 'Minimum' to the list package #H #H Revision 3.9 1991/06/28 13:00:00 martin #H removed 'Random' #H #H Revision 3.8 1991/06/02 13:00:00 martin #H improved 'Gcdex' #H #H Revision 3.7 1991/06/02 12:00:00 martin #H removed 'GcdInt' since there is an internal function of that name #H #H Revision 3.6 1991/06/01 17:00:00 martin #H changed 'IsPrime' to 'IsPrimeInt', 'Factors' to 'FactorsInt', ... #H #H Revision 3.5 1991/06/01 16:00:00 martin #H changed 'Gcd' to 'GcdInt', 'Lcm' to 'LcmInt', ... #H #H Revision 3.4 1991/06/01 12:00:00 martin #H changed 'Quo' to 'QuoInt' #H #H Revision 3.3 1991/01/07 12:00:00 martin #H changed 'Sign' to refuse permutations #H #H Revision 3.2 1990/11/26 12:00:00 martin #H renamed 'Combination' to 'ChineseRem' #H #H Revision 3.1 1990/11/16 12:00:00 martin #H moved functions from integer to combinatorics package #H #H Revision 3.0 1990/01/01 12:00:00 martin #H initial revision under RCS #H ## ############################################################################# ## #V Integers . . . . . . . . . . . . . . . . . . . . . ring of the integers #V IntegersOps . . . . . . . . . . . . . . . . operation record for integers ## IntegersOps := Copy( RingOps ); Integers := rec( isDomain := true, isRing := true, generators := [ 1 ], zero := 0, one := 1, name := "Integers", size := "infinity", isFinite := false, isCommutativeRing := true, isIntegralRing := true, isUniqueFactorizationRing := true, isEuclideanRing := true, units := [ -1, 1 ], operations := IntegersOps ); ############################################################################# ## #F IntegersOps.Ring(<elms>) . . . . . . . . ring generated by some integers ## IntegersOps.Ring := function ( elms ) return Integers; end; ############################################################################# ## #F IntegersOps.DefaultRing(<elms>) . . . . . . default ring of some integers ## IntegersOps.DefaultRing := function ( elms ) return Integers; end; ############################################################################# ## #F IntegersOps.Field(<elms>) . . . . . . . field generated by some integers ## IntegersOps.Field := function ( elms ) return Rationals; end; ############################################################################# ## #F IntegersOps.DefaultField(<elms>) . . . . . default field of some integers ## IntegersOps.DefaultField := function ( elms ) return Rationals; end; ############################################################################# ## #F IntegersOps.\in(<n>,<Integers>) . . . . . . mebership test for integers ## IntegersOps.\in := function ( n, Integers ) return IsInt( n ); end; ############################################################################# ## #F IntegersOps.Random(<Integers>) . . . . . . . . . . . . . random integer ## NrBitsInt := function ( n ) local nr, nr64; nr64:=[0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5, 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6]; nr := 0; while 0 < n do nr := nr + nr64[ n mod 64 + 1 ]; n := QuoInt( n, 64 ); od; return nr; end; IntegersOps.Random := function ( Integers ) return NrBitsInt( RandomList( [0..2^20-1] ) ) - 10; end; ############################################################################# ## #F IntegersOps.Quotient(<Integers>,<n>,<m>) . . . quotient of two integers ## IntegersOps.Quotient := function ( Integers, n, m ) local q; q := QuoInt( n, m ); if n <> q * m then q := false; fi; return q; end; ############################################################################# ## #F IntegersOps.StandardAssociate(<Integers>,<n>) . . . . . . absolute value ## IntegersOps.StandardAssociate := function ( Integers, n ) if n < 0 then return -n; else return n; fi; end; ############################################################################# ## #F IntegersOps.IsPrime(<Integers>,<n>) . test whether an integer is a prime ## IntegersOps.IsPrime := function ( Integers, n ) return IsPrimeInt( n ); end; ############################################################################# ## #F IntegersOps.IsIrreducible(<Integers>,<n>) ## IntegersOps.IsIrreducible := IntegersOps.IsPrime; ############################################################################# ## #V Primes[] . . . . . . . . . . . . . . . . . . . . . . . . list of primes ## ## 'Primes' is a set, i.e., sorted list, of the 168 primes less than 1000. ## ## This is used in 'IsPrimeInt' and 'FactorsInt' to cast out small primes ## quickly. ## Primes := [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137,139, 149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239, 241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349, 353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457, 461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577, 587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683, 691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821, 823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941, 947,953,967,971,977,983,991,997 ]; IsSet( Primes ); ############################################################################# ## #V Primes2[] . . . . . . . . . . . . . . . . . . . . . additional prime list ## ## 'Primes2' contains those primes found by 'IsPrimeInt' that are not in ## 'Primes'. 'Primes2' is kept sorted, but may contain holes. ## ## 'IsPrimeInt' and 'FactorsInt' use this list to cast out already found ## primes quickly. If 'IsPrimeInt' is called only for random integers this ## list would be quite useless. However, users do not behave randomly. ## Instead, it is not uncommon to factor the same integer twice. Likewise, ## once we have tested that $2^{31}-1$ is prime, factoring $2^{62}-1$ is ## very cheap, because the former divides the latter. ## ## This list is initialized to contain all the prime factors of the integers ## $2^n-1$ with $n < 201$, $3^n-1$ with $n < 101$, $5^n-1$ with $n < 101$, ## $7^n-1$ with $n < 91$, $11^n-1$ with $n < 79$, and $13^n-1$ with $n < 37$ ## that are larger than $10^7$. ## Primes2 := [ 10047871, 10567201, 10746341, 12112549, 12128131, 12207031, 12323587, 12553493, 12865927, 13097927, 13264529, 13473433, 13821503, 13960201, 14092193, 14597959, 15216601, 15790321, 16018507, 18837001, 20381027, 20394401, 20515111, 20515909, 21207101, 21523361, 22253377, 22366891, 22996651, 23850061, 25781083, 26295457, 28325071, 28878847, 29010221, 29247661, 29423041, 29866451, 32234893, 32508061, 36855109, 41540861, 42521761, 43249589, 44975113, 47392381, 47763361, 48544121, 48912491, 49105547, 49892851, 51457561, 55527473, 56409643, 56737873, 59302051, 59361349, 59583967, 60816001, 62020897, 65628751, 69566521, 75068993, 76066181, 85280581, 93507247, 96656723, 97685839, 106431697, 107367629, 109688713, 110211473, 112901153, 119782433, 127540261, 134818753, 134927809, 136151713, 147300841, 160465489, 164511353, 177237331, 183794551, 184481113, 190295821, 190771747, 193707721, 195019441, 202029703, 206244761, 212601841, 212885833, 228511817, 231769777, 234750601, 272010961, 283763713, 297315901, 305175781, 308761441, 319020217, 359390389, 407865361, 420778751, 424256201, 432853009, 457315063, 466344409, 510810301, 515717329, 527093491, 529510939, 536903681, 540701761, 550413361, 603926681, 616318177, 632133361, 715827883, 724487149, 745988807, 815702161, 834019001, 852133201, 857643277, 879399649, 1001523179, 1036745531, 1065264019, 1106131489, 1169382127, 1390636259, 1503418321, 1527007411, 1636258751, 1644512641, 1743831169, 1824179209, 1824726041, 1826934301, 1866013003, 1990415149, 2127431041, 2147483647, 2238236249, 2316281689, 2413941289, 2481791513, 2550183799, 2576743207, 2664097031, 2767631689, 2903110321, 2931542417, 3158528101, 3173389601, 3357897971, 4011586307, 4058036683, 4278255361, 4375578271, 4562284561, 4649919401, 4698932281, 4795973261, 4885168129, 5960555749, 6809710909, 7068569257, 7151459701, 7484047069, 7685542369, 7830118297, 7866608083, 8209475377, 8831418697, 9598959833, 10879733611, 11898664849, 12447002677, 13455809771, 13564461457, 13841169553, 13971969971, 14425532687, 15085812853, 15768033143, 15888756269, 16148168401, 17154094481, 17189128703, 19707683773, 22434744889, 23140471537, 23535794707, 24127552321, 25480398173, 25829691707, 25994736109, 27669118297, 27989941729, 28086211607, 30327152671, 32952799801, 33057806959, 35532364099, 39940132241, 43872038849, 45076044553, 47072139617, 50150933101, 54410972897, 56625998353, 60726444167, 61070817601, 62983048367, 70845409351, 76831835389, 77158673929, 77192844961, 78009515593, 83960385389, 86950696619, 88959882481, 99810171997, 115868130379, 125096112091, 127522693159, 128011456717, 128653413121, 131105292137, 152587500001, 158822951431, 159248456569, 164504919713, 165768537521, 168749965921, 229890275929, 241931001601, 269089806001, 282429005041, 332207361361, 374857981681, 386478495679, 392038110671, 402011881627, 441019876741, 447600088289, 487824887233, 531968664833, 555915824341, 593554036769, 598761682261, 641625222857, 654652168021, 761838257287, 810221830361, 840139875599, 918585913061, 1030330938209, 1047623475541, 1113491139767, 1133836730401, 1273880539247, 1534179947851, 1628744948329, 1654058017289, 1759217765581, 1856458657451, 2098303812601, 2454335007529, 2481357870461, 2549755542947, 2663568851051, 2879347902817, 2932031007403, 3138426605161, 3203431780337, 3421169496361, 3740221981231, 4363953127297, 4432676798593, 4446437759531, 4534166740403, 4981857697937, 5625767248687, 6090817323763, 6493405343627, 6713103182899, 6740339310641, 7432339208719, 8090594434231, 8157179360521, 8737481256739, 8868050880709, 9361973132609, 9468940004449, 9857737155463, 10052678938039, 10979607179423, 13952598148481, 15798461357509, 18158209813151, 22125996444329, 22542470482159, 22735632934561, 23161037562937, 23792163643711, 24517014940753, 24587411156281, 28059810762433, 29078814248401, 31280679788951, 31479823396757, 33232924804801, 42272797713043, 44479210368001, 45920153384867, 49971617830801, 57583418699431, 62911130477521, 67280421310721, 70601370627701, 71316922984999, 83181652304609, 89620825374601, 110133112994711, 140737471578113, 145295143558111, 150224123975857, 204064664440913, 205367807127911, 242099935645987, 270547105429567, 303567967057423, 332584516519201, 434502978835771, 475384700124973, 520518327319589, 560088668384411, 608459012088799, 637265428480297, 643170158708221, 707179356161321, 926510094425921, 990643452963163, 1034150930241911, 1066818132868207, 1120648576818041, 1357105535093947, 1416258521793067, 1587855697992791, 1611479891519807, 1628413557556843, 1958423494433591, 2134387368610417, 2646507710984041, 2649263870814793, 2752135920929651, 2864226125209369, 4889988840047743, 5420506947192709, 6957533874046531, 9460375336977361, 9472026608675509, 12557612956332313, 13722816749522711, 14436295738510501, 18584774046020617, 18624275418445601, 20986207825565581, 21180247636732981, 22666879066355177, 27145365052629449, 46329453543600481, 50544702849929377, 59509429687890001, 60081451169922001, 70084436712553223, 76394148218203559, 77001139434480073, 79787519018560501, 96076791871613611, 133088039373662309, 144542918285300809, 145171177264407947, 153560376376050799, 166003607842448777, 177722253954175633, 196915704073465747, 316825425410373433, 341117531003194129, 380808546861411923, 489769993189671059, 538953023961943033, 581283643249112959, 617886851384381281, 625552508473588471, 645654335737185721, 646675035253258729, 658812288653553079, 768614336404564651, 862970652262943171, 909456847814334401, 1100876018364883721, 1195857367853217109, 1245576402371959291, 1795918038741070627, 2192537062271178641, 2305843009213693951, 2312581841562813841, 2461243576713869557, 2615418118891695851, 2691614274040036601, 3011347479614249131, 3358335487319458201, 3421093417510114543, 3602372010909260861, 3747607031112307667, 3999088279399464409, 4710883168879506001, 5079304643216687969, 5559917315850179173, 5782172113400990737, 6106505825833677713, 6115909044841454629, 9213624084535989031, 9520972806333758431, 10527743181888260981, 14808607715315782481, 18446744069414584321, 26831423036065352611, 32032215596496435569, 34563155350221618511, 36230454570129675721, 58523123221688392679, 82064241848634269407, 86656268566282183151, 87274497124602996457, 157571957584602258799, 162715052426691233701, 172827552198815888791, 195489390796456327201, 240031591394168814433, 344120456368919234899, 358475907408445923469, 846041103974872866961, 2519545342349331183143, 3658524738455131951223, 3976656429941438590393, 5439042183600204290159, 8198241112969626815581, 11600321878916922053491, 12812432238302009985937, 17551032119981679046729, 18489605314740987765913, 27665283091695977275201, 42437717969530394595211, 57912614113275649087721, 61654440233248340616559, 63681511996418550459487, 105293313660391861035901, 155285743288572277679887, 201487636602438195784363, 231669654363683130095909, 235169662395069356312233, 402488219476647465854701, 535347624791488552837151, 604088623657497125653141, 870035986098720987332873, 950996059627210897943351, 1412900479108654932024439, 1431185706701868962383741, 2047572230657338751575051, 2048568835297380486760231, 2741672362528725535068727, 3042645634792541312037847, 3745603812007166116831643, 4362139336229068656094783, 4805345109492315767981401, 5042939439565996049162197, 7289088383388253664437433, 8235109336690846723986161, 9680647790568589086355559, 9768997162071483134919121, 9842332430037465033595921, 11053036065049294753459639, 11735415506748076408140121, 13842607235828485645766393, 17499733663152976533452519, 26273701844015319144827917, 75582488424179347083438319, 88040095945103834627376781, 100641220283951395639601683, 140194179307171898833699259, 207617485544258392970753527, 291280009243618888211558641, 303309617049998388989376043, 354639323684545612988577649, 618970019642690137449562111, 913242407367610843676812931, 7222605228105536202757606969, 7248808599285760001152755641, 8170509011431363408568150369, 8206973609150536446402438593, 9080418348371887359375390001, 14732265321145317331353282383, 15403468930064931175264655869, 15572244900182528777225808449, 18806327041824690595747113889, 21283620033217629539178799361, 37201708625305146303973352041, 42534656091583268045915654719, 48845962828028421155731228333, 123876132205208335762278423601, 134304196845099262572814573351, 172974812463239310024750410929, 217648180992721729506406538251, 227376585863531112677002031251, 1786393878363164227858270210279, 2598696228942460402343442913969, 2643999917660728787808396988849, 3340762283952395329506327023033, 5465713352000770660547109750601, 28870194250662203210437116612769, 70722308812401674174993533367023, 78958087694609321439660131899631, 88262612316754526107621113329689, 162259276829213363391578010288127, 163537220852725398851434325720959, 177635683940025046467781066894531, 2679895157783862814690027494144991, 3754733257489862401973357979128773, 5283012903770196631383821046101707, 5457586804596062091175455674392801, 10052011757370829033540932021825161, 11419697846380955982026777206637491, 38904276017035188056372051839841219, 1914662449813727660680530326064591907, 7923871097285295625344647665764672671, 9519524151770349914726200576714027279, 10350794431055162386718619237468234569, 170141183460469231731687303715884105727, 1056836588644853738704557482552056406147, 6918082374901313855125397665325977135579, 235335702141939072378977155172505285655211, 360426336941693434048414944508078750920763, 1032670816743843860998850056278950666491537, 1461808298382111034194027645506019619578037, 79638304766856507377778616296087448490695649, 169002145064468556765676975247413756542145739, 8166146875847876762859119015147004762656450569, 18607929421228039083223253529869111644362732899, 33083146850190391025301565142735000331370209599, 138497973518827432485604572537024087153816681041, 673267426712748387612994804392183645147042355211, 1489459109360039866456940197095433721664951999121, 4884164093883941177660049098586324302977543600799, 466345922275629775763320748688970211803553256223529, 26828803997912886929710867041891989490486893845712448833, 153159805660301568024613754993807288151489686913246436306439, 1051153199500053598403188407217590190707671147285551702341089650185945215953 ]; IsSet( Primes2 ); ############################################################################# ## #F IsPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . . test for a prime ## ## 'IsPrimeInt' returns 'true' if the integer <n> is a prime and 'false' ## otherwise. By convention 'IsPrimeInt( 0 ) = IsPrimeInt( 1 ) = false' and ## we define 'IsPrimeInt( -<n> ) = IsPrimeInt( <n> )'. ## ## 'IsPrimeInt' does trial divisions by the primes less than 1000 to detect ## composites with a factor less than 1000 and primes less than 1000000. ## ## 'IsPrimeInt' then checks that $n$ is a strong pseudoprime to the base 2. ## This uses Fermats theorem which says $2^{n-1}=1$ mod $n$ for a prime $n$. ## If $2^{n-1}\<>1$ mod $n$, $n$ is composite, 'IsPrimeInt' returns 'false'. ## There are composite numbers for which $2^{n-1}=1$, but they are seldom. ## ## Then 'IsPrimeInt' checks that $n$ is a Lucas pseudoprime for $p$, choosen ## so that the discriminant $d=p^2/4-1$ is an quadratic nonresidue mod $n$. ## I.e., 'IsPrimeInt' takes the root $a = p/2+\sqrt{d}$ of $x^2 - px + 1$ in ## the ring $Z_n[\sqrt{d}] and computes the traces of $a^n$ and $a^{n+1}$. ## If $n$ is a prime, this ring is the field of order $n^2$ and raising to ## the $n$th power is conjugation, so $trace(a^n)=p$ and $trace(a^{n+1})=2$. ## However, these identities hold only for extremly few composite numbers. ## ## Note that this test for $trace(a^n) = p$ and $trace(a^{n+1}) = 2$ is ## usually formulated using the Lucas sequences $U_k = (a^k-b^k)/(a-b)$ and ## $V_k = (a^k+b^k)=trace(a^k)$, where one tests $U_{n+1} = 0, V_{n+1} = 2$. ## However, the trace test is equivalent and requires fewer multiplications. ## Thanks to Daniel R. Grayson (dan@symcom.math.uiuc.edu) for telling me. ## ## Better descriptions of the algorithm and related topics can be found in: ## G. Miller, cf. Algorithms and Complexity ed. Traub, AcademPr, 1976, 35-36 ## C. Pomerance et.al., Pseudoprimes to 25*10^9, MathComp 35 1980, 1003-1026 ## D. Knuth, Seminumerical Algorithms (TACP II), AddiWesl, 1973, 378-380 ## G. Gonnet, Heuristic Primality Testing, Maple Newsletter 4, 1989, 36-38 ## R. Baillie, S. Wagstaff, Lucas Pseudoprimes, MathComp 35 1980, 1391-1417 ## TraceModQF := function ( p, k, n ) local trc; if k = 1 then trc := [ p, 2 ]; elif k mod 2 = 0 then trc := TraceModQF( p, k/2, n ); trc := [ (trc[1]^2 - 2) mod n, (trc[1]*trc[2] - p) mod n ]; else trc := TraceModQF( p, (k+1)/2, n ); trc := [ (trc[1]*trc[2] - p) mod n, (trc[2]^2 - 2) mod n ]; fi; return trc; end; IsPrimeInt := function ( n ) local p, e, o, x, i, d; # make $n$ positive and handle trivial cases if n < 0 then n := -n; fi; if n in Primes then return true; fi; if n in Primes2 then return true; fi; if n <= 1000 then return false; fi; # do trial divisions by the primes less than 1000 # faster than anything fancier because $n$ mod <small int> is very fast for p in Primes do if n mod p = 0 then return false; fi; if n < (p+1)^2 then AddSet(Primes2,n); return true; fi; od; # do trial division by the other known primes for p in Primes2 do if n mod p = 0 then return false; fi; od; # find $e$ and $o$ odd such that $n-1 = 2^e * o$ e := 0; o := n-1; while o mod 2 = 0 do e := e+1; o := o/2; od; # look at the seq $2^o, 2^{2 o}, 2^{4 o}, .., 2^{2^e o}=2^{n-1}$ x := PowerModInt( 2, o, n ); i := 0; while i < e and x <> 1 and x <> n-1 do x := x * x mod n; i := i + 1; od; # if it is not of the form $.., -1, 1, 1, ..$ then $n$ is composite if not (x = n-1 or (i = 0 and x = 1)) then return false; fi; # there are no strong pseudo-primes to base 2 smaller than 2047 if n < 2047 then AddSet( Primes2, n ); return true; fi; # make sure that $n$ is not a perfect power (especially not a square) if SmallestRootInt(n) < n then return false; fi; # find a quadratic nonresidue $d = p^2/4-1$ mod $n$ p := 2; while Jacobi( p^2-4, n ) <> -1 do p := p+1; od; # for a prime $n$ the trace of $(p/2+\sqrt{d})^n$ must be $p$ # and the trace of $(p/2+\sqrt{d})^{n+1}$ must be 2 if TraceModQF( p, n+1, n ) = [ 2, p ] then AddSet( Primes2, n ); return true; fi; # $n$ is not a prime return false; end; ############################################################################# ## #F IsPrimePowerInt( <n> ) . . . . . . . . . . . test for a power of a prime ## ## 'IsPrimePowerInt' returns 'true' if the integer <n> is a prime power and ## 'false' otherwise. ## IsPrimePowerInt := function ( n ) return IsPrimeInt( SmallestRootInt( n ) ); end; ############################################################################# ## #F NextPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . next larger prime ## ## 'NextPrimeInt' returns the smallest prime which is strictly larger than ## the integer <n>. ## NextPrimeInt := function ( n ) if -3 = n then n := -2; elif -3 < n and n < 2 then n := 2; elif n mod 2 = 0 then n := n+1; else n := n+2; fi; while not IsPrimeInt(n) do if n mod 6 = 1 then n := n+4; else n := n+2; fi; od; return n; end; ############################################################################# ## #F PrevPrimeInt( <n> ) . . . . . . . . . . . . . . . previous smaller prime ## ## 'PrevPrimeInt' returns the largest prime which is strictly smaller than ## the integer <n>. ## PrevPrimeInt := function ( n ) if 3 = n then n := 2; elif -2 < n and n < 3 then n := -2; elif n mod 2 = 0 then n := n-1; else n := n-2; fi; while not IsPrimeInt(n) do if n mod 6 = 5 then n := n-4; else n := n-2; fi; od; return n; end; ############################################################################# ## #F IntegersOps.Factors(<Integers>,<n>) . . . . . factorization of an integer ## IntegersOps.Factors := function ( Integers, n ) return FactorsInt( n ); end; ############################################################################# ## #F FactorsInt( <n> ) . . . . . . . . . . . . . . prime factors of an integer #F FactorsRho(<n>,<inc>,<cluster>,<limit>) . . . Pollards rho factorization ## ## 'FactorsInt' returns a list of prime factors of the integer <n>. ## ## 'FactorsInt' does trial divisions by the primes less than 1000 to detect ## all composites with a factor less than 1000 and primes less than 1000000. ## After that it calls 'FactorsRho(<n>,1,16,8192)' to do the hard work. ## ## 'FactorsInt' and 'FactorsRho' return the same value for all integers. ## Usually 'FactorsInt' factors integers with prime factors $\<1000$ faster. ## However for integers with no factor $\<1000$ 'FactorsRho' will be faster. ## ## 'FactorsRho' uses Pollards $\rho$ method to factor the integer $n = p q$. ## For a small simple example lets assume we want to factor $667 = 23 * 29$. ## 'FactorsRho' first calls 'IsPrimeInt' to avoid trying to factor a prime. ## ## Then it uses the sequence defined by $x_0=1, x_{i+1}=(x_i^2+1)$ mod $n$. ## In our example this is $1, 2, 5, 26, 10, 101, 197, 124, 36, 630, .. $. ## ## Modulo $p$ it takes on at most $p-1$ different values, thus it eventually ## becomes recurrent, usually this happens after roughly $2 \sqrt{p}$ steps. ## In our example modulo 23 we get $1, 2, 5, 3, 10, 9, 13, 9, 13, 9, .. $. ## ## Thus there exist pairs $i, j$ such that $x_i = x_j$ mod $p$, i.e., such ## that $p$ divides $Gcd( n, x_j-x_i )$. With a bit of luck no other factor ## of $n$ divides $x_j - x_i$ so we find $p$ if we know such a pair. In our ## example $5, 7$ is the first pair, $x_7-x_5=23$, and $Gcd(667,23) = 23$. ## ## Now it is too expensive to check all pairs, but there also must be pairs ## of the form $2^i-1, j$ with $3*2^{i-1} <= j < 4*2^{i-1}$. In our example ## $7, 13$ is the first such pair, $x_13-x_7=506$, and $Gcd(667,506) = 23$. ## ## Thus by taking the gcds of $n$ and $x_j-x_i$ for such pairs, we will find ## the factor $p$ after approximately $2 \sqrt{p} \<= 2 \sqrt^4{n}$ steps. ## ## If $Gcd( n, x_j - x_i )$ is not a prime 'FactorsRho' will call itself ## recursivly with a different value for <inc>, i.e., it will try to factor ## the gcd using a different sequence $x_{i+1} = (x_i^2 + inc)$ mod $n$. ## ## Since the gcd computations are by far the most time consuming part of the ## algorithm one can save time by clustering differences and computing the ## gcd only every <cluster> iteration. This slightly increases the chance ## that a gcd is composite, but reduces the runtime by a large amount. ## ## Finally 'FactorsRho' accepts an argument <limit> which is the number of ## iterations performed by 'FactorsRho' before giving up. The default value ## is 8192 which corresponds to a few minutes while guaranteing that all ## prime factors less than $10^6$ and most less than $10^9$ are found. ## ## Better descriptions of the algorithm and related topics can be found in: ## J. Pollard, A Monte Carlo Method for Factorization, BIT 15, 1975, 331-334 ## R. Brent, An Improved Monte Carlo Method for Fact., BIT 20, 1980, 176-184 ## D. Knuth, Seminumerical Algorithms (TACP II), AddiWesl, 1973, 369-371 ## FactorsRho := function ( n, inc, cluster, limit ) local sign, factors, x, y, z, i, g, k; # make $n$ positive and handle trivial cases sign := 1; if n < 0 then sign := -sign; n := -n; fi; if n < 4 then return [ sign * n ]; fi; factors := []; while n mod 2 = 0 do Add( factors, 2 ); n := n / 2; od; while n mod 3 = 0 do Add( factors, 3 ); n := n / 3; od; if IsPrimeInt(n) then Add( factors, n ); n := 1; fi; # initialize $x_0$ x := 1; z := 1; i := 0; # loop until we have factored $n$ completely or run out of patience while 1 < n and 2^i <= limit do # $y = x_{2^i-1}$ y := x; i := i + 1; # $x_{2^i}, .., x_{3*2^{i-1}-1}$ need not be compared to $x_{2^i-1}$ for k in [1..2^(i-1)] do x := (x^2 + inc) mod n; od; # compare $x_{3*2^{i-1}}, .., x_{4*2^{i-1}-1}$ with $x_{2^i-1}$ for k in [1..2^(i-1)] do x := (x^2 + inc) mod n; z := z * (x - y) mod n; # from time to time compute the gcd if k mod cluster = 0 then g := GcdInt( n, z ); # if it is > 1 we have found a factor which need not be prime if g > 1 then factors := Concatenation( factors, FactorsRho(g,inc+1,QuoInt(cluster+1,2),limit)); n := n / g; if IsPrimeInt(n) then Add( factors, n ); n := 1; fi; fi; fi; od; od; # sort the list of factors and return it if 1 < n then Error("could not factor",n); fi; Sort( factors ); factors[1] := sign * factors[1]; return factors; end; FactorsInt := function ( n ) local sign, factors, p; # make $n$ positive and handle trivial cases sign := 1; if n < 0 then sign := -sign; n := -n; fi; if n < 4 then return [ sign * n ]; fi; factors := []; # do trial divisions by the primes less than 1000 # faster than anything fancier because $n$ mod <small int> is very fast for p in Primes do while n mod p = 0 do Add( factors, p ); n := n / p; od; if n < (p+1)^2 and 1 < n then Add(factors,n); n := 1; fi; if n = 1 then factors[1] := sign*factors[1]; return factors; fi; od; # do trial divisions by known factors for p in Primes2 do while n mod p = 0 do Add( factors, p ); n := n / p; od; if n = 1 then factors[1] := sign*factors[1]; return factors; fi; od; # handle perfect powers p := SmallestRootInt( n ); if p < n then while 1 < n do Append( factors, FactorsInt(p) ); n := n / p; od; Sort( factors ); factors[1] := sign * factors[1]; return factors; fi; # let 'FactorsRho' do the work factors := Concatenation( factors, FactorsRho( n, 1, 16, 8192 ) ); Sort( factors ); factors[1] := sign * factors[1]; return factors; end; ############################################################################# ## #F DivisorsInt( <n> ) . . . . . . . . . . . . . . . divisors of an integer ## ## 'DivisorsInt' returns a list of all divisors of the integer <n>. The ## list is sorted, so that it starts with 1 and ends with <n>. We define ## that 'Divisors( -<n> ) = Divisors( <n> )'. ## DivisorsSmall := [,[1],[1,2],[1,3],[1,2,4],[1,5],[1,2,3,6],[1,7]]; DivisorsInt := function ( n ) local divisors, factors, divs; # make <n> it nonnegative, handle trivial cases, and get prime factors if n < 0 then n := -n; fi; if n = 0 then Error("DivisorsInt: <n> must not be 0"); fi; if n < 8 then return DivisorsSmall[n+1]; fi; factors := FactorsInt( n ); # recursive function to compute the divisors divs := function ( i, m ) if Length(factors) < i then return [ m ]; elif m mod factors[i] = 0 then return divs(i+1,m*factors[i]); else return Concatenation( divs(i+1,m), divs(i+1,m*factors[i]) ); fi; end; divisors := divs( 1, 1 ); Sort( divisors ); return divisors; end; ############################################################################# ## #F Sigma( <n> ) . . . . . . . . . . . . . . . sum of divisors of an integer ## ## 'Sigma' returns the sum of the positive divisors of the integer <n>. ## ## 'Sigma' is a multiplicative arithmetic function, i.e., if $n$ and $m$ are ## relative prime we have $\sigma(n m) = \sigma(n) \sigma(m)$. ## Sigma := function ( n ) local sigma, p, q, k; # make <n> it nonnegative, handle trivial cases if n < 0 then n := -n; fi; if n = 0 then Error("Sigma: <n> must not be 0"); fi; if n < 8 then return Sum(DivisorsSmall[n+1]); fi; # loop over all prime $p$ factors of $n$ sigma := 1; for p in Set(FactorsInt(n)) do # compute $p^e$ and $k = 1+p+p^2+..p^e$ q := p; k := 1 + p; while n mod (q * p) = 0 do q := q * p; k := k + q; od; # combine with the value found so far sigma := sigma * k; od; return sigma; end; ############################################################################# ## #F Tau( <n> ) . . . . . . . . . . . . . . number of divisors of an integer ## ## 'Tau' returns the number of the positive divisors of the integer <n>. ## ## 'Tau' is a multiplicative arithmetic function, i.e., if $n$ and $m$ are ## relative prime we have $\tau(n m) = \tau(n) \tau(m)$. ## Tau := function ( n ) local tau, p, q, k; # make <n> it nonnegative, handle trivial cases if n < 0 then n := -n; fi; if n = 0 then Error("Tau: <n> must not be 0"); fi; if n < 8 then return Length(DivisorsSmall[n+1]); fi; # loop over all prime factors $p$ of $n$ tau := 1; for p in Set(FactorsInt(n)) do # compute $p^e$ and $k = e+1$ q := p; k := 2; while n mod (q * p) = 0 do q := q * p; k := k + 1; od; # combine with the value found so far tau := tau * k; od; return tau; end; ############################################################################# ## #F MoebiusMu( <n> ) . . . . . . . . . . . . . . Moebius inversion function ## ## 'MoebiusMu' computes the value of Moebius inversion function for the ## integer <n>. This is 0 for integers which are not squarefree, i.e., ## which are divided by a square $r^2$. Otherwise it is 1 if <n> has a even ## number and -1 if <n> has an odd number of prime factors. ## MoebiusMu := function ( n ) local factors; if n < 0 then n := -n; fi; if n = 0 then Error("MoebiusMu: <n> must be nonzero"); fi; if n = 1 then return 1; fi; factors := FactorsInt( n ); if factors <> Set( factors ) then return 0; fi; return (-1) ^ Length(factors); end; ############################################################################# ## #F IntegersOps.EuclideanDegree(<Integers>,<n>) . . . . . . . . absolut value ## IntegersOps.EuclideanDegree := function ( Integers, n ) if n < 0 then return -n; else return n; fi; end; ############################################################################# ## #F IntegersOps.Mod(<Integers>,<n>,<m>) . . . . . . remainder of two integers ## IntegersOps.Mod := function ( Integers, n, m ) return n mod m; end; ############################################################################# ## #F IntegersOps.QuotientMod(<Integers>,<r>,<s>,<m>) quotient of two integers #F modulo another ## IntegersOps.QuotientMod := function ( Integers, r, s, m ) if r mod GcdInt( s, m ) = 0 then return r/s mod m; else return false; fi; end; ############################################################################# ## #F IntegersOps.PowerMod(<Integers>,<r>,<e>,<m>) . . power of an integer mod #F another #F PowerModInt(<r>,<e>,<m>) . . . . . . power of one integer modulo another ## PowerModInt := function ( r, e, m ) local pow, f; # reduce r initially r := r mod m; # handle special case if e = 0 then return 1; fi; # if e is negative then invert n modulo m with Euclids algorithm if e < 0 then r := 1/r mod m; e := -e; fi; # now use the repeated squaring method (right-to-left) pow := 1; f := 2 ^ (LogInt( e, 2 ) + 1); while 1 < f do pow := (pow * pow) mod m; f := QuoInt( f, 2 ); if f <= e then pow := (pow * r) mod m; e := e - f; fi; od; # return the power return pow; end; IntegersOps.PowerMod := function ( Integers, r, e, m ) return PowerModInt( r, e, m ); end; ############################################################################# ## #F IntegersOps.Gcd(<Integers>,<n>,<m>) . . . . . . . . . gcd of two integers #F GcdInt(<n>,<m>) . . . . . . . . . . . . . . . . . . . gcd of two integers ## IntegersOps.Gcd := function ( Integers, n, m ) return GcdInt( n, m ); end; ############################################################################# ## #F IntegersOps.Lcm(<Integers>,<n>,<m>) . . least common multiple of integers #F LcmInt( <m>, <n> ) . . . . . . . . . . least common multiple of integers ## IntegersOps.Lcm := function ( Integers, n, m ) return LcmInt( n, m ); end; LcmInt := function ( n, m ) if m = 0 and n = 0 then return 0; else return AbsInt( m / GcdInt( m, n ) * n ); fi; end; ############################################################################# ## #F Gcdex( <m>, <n> ) . . . . . . . . . . greatest common divisor of integers ## Gcdex := function ( m, n ) local f, g, h, fm, gm, hm, q; if 0 <= m then f:=m; fm:=1; else f:=-m; fm:=-1; fi; if 0 <= n then g:=n; gm:=0; else g:=-n; gm:=0; fi; while g <> 0 do q := QuoInt( f, g ); h := g; hm := gm; g := f - q * g; gm := fm - q * gm; f := h; fm := hm; od; if n = 0 then return rec( gcd := f, coeff1 := fm, coeff2 := 0, coeff3 := gm, coeff4 := 1 ); else return rec( gcd := f, coeff1 := fm, coeff2 := (f - fm * m) / n, coeff3 := gm, coeff4 := (0 - gm * m) / n ); fi; end; ############################################################################# ## #F Int( <obj> ) . . . . . . . . . . . . . . . . . . . convert to an integer ## Int := function ( obj ) if IsInt( obj ) then return obj; elif IsRat( obj ) then return QuoInt( Numerator( obj ), Denominator( obj ) ); elif IsFFE( obj ) then return IntFFE( obj ); else Error("<obj> must be rational or a finite field element"); fi; end; ############################################################################# ## #F AbsInt( <n> ) . . . . . . . . . . . . . . . absolute value of an integer ## AbsInt := function ( n ) if 0 <= n then return n; else return -n; fi; end; ############################################################################# ## #F SignInt( <n> ) . . . . . . . . . . . . . . . . . . . sign of an integer ## SignInt := function ( n ) if 0 = n then return 0; elif 0 <= n then return 1; else return -1; fi; end; ############################################################################# ## #F ChineseRem( <moduli>, <residues> ) . . . . . . . . . . chinese remainder ## ChineseRem := function ( moduli, residues ) local i, c, l, g; # combine the residues modulo the moduli i := 1; c := residues[1]; l := moduli[1]; while i < Length(moduli) do i := i + 1; g := Gcdex( l, moduli[i] ); if g.gcd <> 1 and (residues[i]-c) mod g.gcd <> 0 then Error("the residues must be equal modulo ",g.gcd); fi; c := l * (((residues[i]-c) / g.gcd * g.coeff1) mod moduli[i]) + c; l := moduli[i] / g.gcd * l; od; # reduce c into the range [0..l-1] c := c mod l; return c; end; ############################################################################# ## #F LogInt( <n>, <base> ) . . . . . . . . . . . . . . logarithm of an integer ## LogInt := function ( n, base ) local log; # check arguments if n <= 0 then Error("<n> must be positive"); fi; if base <= 1 then Error("<base> must be greater than 1"); fi; # 'log(b)' returns $log_b(n)$ and divides 'n' by 'b^log(b)' log := function ( b ) local i; if b > n then return 0; fi; i := log( b^2 ); if b > n then return 2 * i; else n := QuoInt( n, b ); return 2 * i + 1; fi; end; return log( base ); end; ############################################################################# ## #F RootInt( <n>, <k> ) . . . . . . . . . . . . . . . . . root of an integer ## RootInt := function ( arg ) local n, k, r, s, t; # get the arguments if Length(arg) = 1 then n := arg[1]; k := 2; elif Length(arg) = 2 then n := arg[1]; k := arg[2]; else Error("usage: 'Root( <n> )' or 'Root( <n>, <k> )'"); fi; # check the arguments and handle trivial cases if k <= 0 then Error("<k> must be positive"); elif k = 1 then return n; elif n < 0 and k mod 2 = 0 then Error("<n> must be positive"); elif n < 0 and k mod 2 = 1 then return -RootInt( -n, k ); elif n = 0 then return 0; elif n <= k then return 1; fi; # r is the first approximation, s the second, we need: root <= s < r r := n; s := 2^( QuoInt( LogInt(n,2), k ) + 1 ) - 1; # do Newton iterations until the approximations stop decreasing while s < r do r := s; t := r^(k-1); s := QuoInt( n + (k-1)*r*t, k*t ); od; # and thats the integer part of the root return r; end; ############################################################################# ## #F SmallestRootInt( <n> ) . . . . . . . . . . . smallest root of an integer ## SmallestRootInt := function ( n ) local k, r, s, p, l, q; # check the argument if n > 0 then k := 2; s := 1; elif n < 0 then k := 3; s := -1; n := -n; else return 0; fi; # exclude small divisors, and thereby large exponents if n mod 2 = 0 then p := 2; else p := 3; while p < 100 and n mod p <> 0 do p := p+2; od; fi; l := LogInt( n, p ); # loop over the possible prime divisors of exponents # use Euler's criterion to cast out impossible ones while k <= l do q := 2*k+1; while not IsPrimeInt(q) do q := q+2*k; od; if PowerModInt( n, (q-1)/k, q ) <= 1 then r := RootInt( n, k ); if r ^ k = n then n := r; l := QuoInt( l, k ); else k := NextPrimeInt( k ); fi; else k := NextPrimeInt( k ); fi; od; return s * n; end; ############################################################################# ## #F IntegersOps.AsGroup(<Integers>) . . . . . . . . . . view the integers as #F multiplicative group ## IntegersOps.AsGroup := function ( Integers ) Error("sorry, Z is not finitely generated as multiplicative group"); end; ############################################################################# ## #F IntegersOps.AsAdditiveGroup(<Integers>) . . . . . . view the integers as #F additive group ## #N 14-Oct-91 martin this should be #N IntegersAsAddtiveGroupOps := Copy( AdditveGroupOps ); ## IntegersAsAdditiveGroupOps := Copy( DomainOps ); IntegersAsAdditiveGroupOps.\in := IntegersOps.\in; IntegersAsAdditiveGroupOps.Random := IntegersOps.Random; IntegersOps.AsAdditiveGroup := function ( Integers ) return rec( isDomain := true, isAdditiveGroup := true, generators := [ 1 ], zero := 0, isFinite := false, size := "infinity", isCyclic := true, operations := IntegersAsAdditiveGroupOps ); end; ############################################################################# ## #E Emacs . . . . . . . . . . . . . . . . . . . . . . . local emacs variables ## ## Local Variables: ## mode: outline ## outline-regexp: "#F\\|#V\\|#E" ## fill-column: 73 ## fill-prefix: "## " ## eval: (hide-body) ## End: ##