> load "tut7data.txt";
Loading "L:\win\Magma\MATH2068\tut7data.txt"
> k:=4096;
> time x:=Factorial(k);
Time: 0.000
> time x:=Factorial(2*k);
Time: 0.031
> time x:=Factorial(4*k);
Time: 0.141
> // Roughly 4 times the previous one.
> time x:=Factorial(8*k);
Time: 0.672
> // Again, roughly 4 times the previous one.
> time x:=Factorial(16*k);
Time: 5.156
> // More than 4 times the previous one. We probably crossed some
> // threshold here, and from now on the computer will have to do 
> // more storing values in memory and retrieving them again.
> // But the amount of this will probably be proportional to the number
> // of bit operations it has to do -- so we will probably return to
> // multiplying the time by 4 with each doubling of the argument.
> time x:=Factorial(32*k);
Time: 24.250
> // About 4 times the previous one.
> time x:=Factorial(64*k);
Time: 106.047
> // About 4 times the previous one.
> time x:=Factorial(128*k);
Time: 453.000
> // Again, about 4 times the previous one.
> time x:=8765^8192;
Time: 0.000
> // Certainly faster than Factorial(8192). 
> Log(2,512);
9.00000000000000000000000000000
> // 512 is 2 to the power 9
>
> // The trouble with the command Log(2,711^16000000) is that magma
> // would have to compute the enormous number 711^16000000 first and
> // then find the log of this enormous number.
> // However, let us do it, and see how long it takes.
> time Log(2,711^16000000);
151579291.993910644597356227677
Time: 8.875
> // It is quicker to compute Log(2,711) and multiply the answer
> // by 16000000 -- using the fact that log(a^b) equals b times log(a).
> time 16000000*Log(2,711);
151579291.993910644597356227677
Time: 0.000
> s:=82349;
> IsPrime(s);
true
> time (711^16000000) mod s;
30638
Time: 8.906
> time Modexp(711,16000000,s);
30638
Time: 0.000
> /*
> The first method requires magma to work out and store 711^16000000
> and then find the remainder when this number is divided by s.
> The second method uses "residue arithmetic" modulo p: this means
> that at each intermediate step in the calculation the answer is reduced
> modulo s, so that at the next step the numbers that have to be
> multiplied are less than s. To compute any number raised to the 
> power 16000000 requires performing 30 multiplications. If the numbers 
> are always less than 82349 this will be very quick, but with enormous
> numbers it is significantly slower.
> */
> a:=Random(10^1000,10^2000);
> b:=Random(10^1000,10^2000);
> c:=Random(10^1000,10^2000);
> a;
89114819549254512808642384030079203078874464796753579025664491766222299402229985838052340\
46395962261778465159358946802717994338382812760167477188675200098943259780434581960422624\
77729033664034670761691607792863519893259980476401186782393869781429119403539591402560329\
19649472409943799031721799913577593089327897888204164362488527812522681670566999906088381\
75447459591544038460725242801568997403544461341304642759492993335821984292793871769667648\
01522975162996711105190017138575237700867356003033509737164659642335991398554664077835193\
86552459487050094874230276399344496897939263991894460558866071645440111487255194541044711\
69348043523603843237515515089198447401133282901680649934877395843120944000811683789234512\
44220169931152498999604163407526677905092424562437786628783507104824615718023552977743727\
92712008390458857205036811367145202475936221230859938425810819763186153414331968750719391\
08986303808436889691744683920984177519613581891734585033706030708088007360239599058176073\
73354372437027565193343695130265956905815541774508064964299117028177506583871037528615798\
76686540602968068577811633710061850848754265899047115582087756023875273760342692369210949\
92843114593774881882362851648728581541268015363281909570013072753995375563717798586951174\
24476896763641052088773110257828958467460822492832448394761356769284443752163558997281303\
10160956178996536803444301927448323932981085681788048541071442356971648924373861169860089\
23792511471184973341143884832129597954058714620720988918233781397942686829784226903605073\
47790351516887415811903183456838652965274590445777788998801416194861666678053554128336829\
92423965293650429443636032465012583139974804368978793836040730628405639814578770878112971\
90886799288958293106880081154268987490600386496941692729982871658198243456426198926518967\
80536939967809873087812945949670799335950123440178562342803868728838776434694816710541037\
95770455872673833078343856936489131203266829991065962474129285094138381977093474955500052\
531361405745407515587411028853591451550153
> // That is a big number. And b and c are also big.
> time Modexp(a,b,c);
10341106544300802926641454734266141873988373546605453224886279485456599514585927696394380\
27068291559089425777965630305884522996784808653098724144428084559153107891116029199983272\
13132359585941139666176295187290931527683065602331469522163281774273806051621813080408836\
07451559878922399014082101359885561576633580855317290616090589236359962576531421544885817\
26329670194874519721218440690188165476285122145138866181036110774467733603218839841485114\
12433067747891158283177683715483972758760082150596504809162683162167454195172975749659048\
14751137215145040191923063010336102760601936976779853820684466310252686473118976942286038\
70954385223083644738683465119161129396871251555372662243474864436737875274477488583460905\
84736203944392419143196834717324124143208592447544523478254860038318519721442648827419732\
86580000297672031856354466230337347306455319922979441234025658078156377073639548352900379\
92881409731030592565428598091499262228404283942120247970201712172099180489850032419338550\
18845892337230958530268212694104744366305740532722135273071708968920861490251713057016002\
99009096552812729886139669334095490463593571302005484287629878036341844808358575525874031\
11405082474846149695934331222676348819562721287604302038175484152729821201400680732604031\
47355563617526049542418802907008166440623480546256579527619318584165268357630390336718938\
88270659981699601593034647723976057600626991483842770694807436654970068080696327936282277\
41295308875888541853571852798944840217277193628237178977380680414817401586360623503964558\
39413584307493741256019841951060285815772795885170834495888227814217543310536373985168290\
40052910954213307340768584421026102687269065063800464668689148121710211241089111210691584\
24243270257517771936124135072834697502706019055809083993768282031450303235341481919655904\
20052591131786251229038915362846685504961955462553284186454849364676526202618796316588070\
91010138553485658163404262465560477307200769372567061513344220888071390827722544691162474\
2710144496284012156486560733523712412617
Time: 1.109
> // Impressively fast, considering the sizes of the numbers.
> r:=RandomPrime(500);
> a:=2+Random(r-3);
> r;
27318611366797047982018293936724122189963166063869149228516541514577958203354495785685776\
54840382287631779270701745390265958116894143367433695711690483
> a;
20086264819101267823959497602912364073436282770974016705963004700783815047282133908522898\
36568076830191925660122214798975446290108979200019995444574989
> Modexp(a,r-1,r);
1
> // Fermat would be pleased! His Little Theorem is verified again.
> r:=RandomPrime(500);
> a:=2+Random(r-3);
> Modexp(a,r-1,r);
1
> r;
22280432505955451078204267739514253482097559839392890680750590148023608426497111408781004\
33142815072195871794375072531223202424459175166306531959904647
> a;
53106532539956911133440723995035535712040983777790202889978453796410640193300084589898459\
0732632085441868137254169814180671419978102714934656499668728
> a:=2+Random(r-3);
> a;
56786820342039204776679457624040137282417220702138582717114793030721114549166674781151243\
3247653077697386358534632487506581437466724293083333309953290
> Modexp(a,r-1,r);
1
> a:=2+Random(r-3);
> a;
64183522020288492413088390459171355217102189735824276088345712765338947466975017318157693\
5411874595788216680422047407464799337682359395215916039469459
> Modexp(a,r-1,r);
1
> p;
13559124823046793778874877583461476339003730746686494352905506298024201955142764868369379\
58029578386351546627119214396469
> q;
50248191594291484564740960681639530146235395963601585562284859415426419457705468645716112\
398542382996392421480295251767252359709783
> n;
68132150195936891622562626877462147079749459583334673791656693585971848879807618343190418\
27427977532721209021029265203683020452541780452131844351982875944483709821435049993754378\
3440768421588830132176464233743456342151973035679362308273158872339956227
> n eq p*q;
true
> e;
21344994621692026374746166212755925147861060304319836290537059630761753379082601340297938\
71895427914147302199282723617177988317305977153798606954048708391694732678629081257502307\
1391256027402621715749313049812318617534422848299055630791112277490871149
> // The decryption exponent is the inverse of e modulo (p-1)*(q-1)
> d:=InverseMod(e,(p-1)*(q-1));
> e*d mod ((p-1)*(q-1));
1
> NaiveDecoding([Modexp(m,d,n): m in encipheredphi]);
phivalue:=&*[(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(123456789)];
phivalue eq EulerPhi(123456789);
/*
Please check for other numbers besides 123456789.
*/
> phivalue:=&*[(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(123456789)];
> phivalue eq EulerPhi(123456789);
true
> phivalue:=&*[(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(12121212)];
> phivalue eq EulerPhi(12121212);
true
> phivalue:=&*[(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(10800)];
> phivalue eq EulerPhi(10800);
true
> /*
> OK, it works. Let us see why.
> Since 10800 is 2^4 times 3^3 times 5^2, magma will tell us that
> Factorization(10800) is [<2,4>,<3,3>,<5,2>]. This is a sequence
> with three terms, each of which is an ordered pair.
> Let t be the first term of this sequence. That is, t = <2,4>.
> This means that t[1] (the first component of t) is 2, and t[2] (the
> second component of t) is 4. 
> So for this value of t, (t[1]-1)*t[1]^(t[2]-1) is equal to
> 1*2^3.
> Now let t be the second term in Factorization(10800), i.e. t = <3,3>.
> Then t[1] and t[2] are both 3, and (t[1]-1)*t[1]^(t[2]-1) equals
> 2*3^2.
> Finally, let t be the third term of Factorization(10800), i.e. t = <5,2>.
> Then t[1] is 5 and t[2] is 2; so (t[1]-1)*t[1]^(t[2]-1) is equal to
> 4*5^1.
> So the sequence [(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(10800)]
> is [1*2^3, 2*3^2, 4*5^1], that is, [8, 18, 20].
> The three terms are EulerPhi(2^4), EulerPhi(3^3) and EulerPhi(5^2).
> And &*[8, 18, 20] = 8*18*20 = EulerPhi(2^4)*EulerPhi(3^3)*EulerPhi(5^2),
> which equals EulerPhi(2^4*3^3*5^2) since EulerPhi is multiplicative.
> */
> Factorization(10800);
[ <2, 4>, <3, 3>, <5, 2> ]
> [(t[1]-1)*t[1]^(t[2]-1) : t in Factorization(10800)];
[ 8, 18, 20 ]
> [EulerPhi(2^4), EulerPhi(3^3), EulerPhi(5^2)];
[ 8, 18, 20 ]
> &*[ 8, 18, 20 ];
2880
> EulerPhi(10800);
2880
> NaiveDecoding([Modexp(m,d,n): m in encipheredtau]);
tauvalue:=&*[t[2]+1 : t in Factorization(123456789)];
tauvalue eq NumberOfDivisors(123456789);
> tauvalue:=&*[t[2]+1 : t in Factorization(123456789)];
> tauvalue eq NumberOfDivisors(123456789);
true
> tauvalue eq NumberOfDivisors(2*3*5*7*11);
false
> tauvalue:=&*[t[2]+1 : t in Factorization(2*3*5*7*11)];
> tauvalue;
32
> tauvalue eq NumberOfDivisors(2*3*5*7*11);
true
> Factorization(2*3*5*7*11);
[ <2, 1>, <3, 1>, <5, 1>, <7, 1>, <11, 1> ]
> [t[2]+1 : t in Factorization(2*3*5*7*11)];
[ 2, 2, 2, 2, 2 ]
> &*[ 2, 2, 2, 2, 2 ];
32
> tauvalue:=&*[t[2]+1 : t in Factorization(2*3^2*5^3*7^4*11^5)];
> [t[2]+1 : t in Factorization(2*3^2*5^3*7^4*11^5)];
[ 2, 3, 4, 5, 6 ]
> tauvalue;
720
> 2*3*4*5*6;
720
> NumberOfDivisors(2*3^2*5^3*7^4*11^5);
720
> NaiveDecoding([Modexp(m,d,n): m in encipheredsigma]);
sigmavalue:=&*[(t[1]^(t[2]+1)-1)/(t[1]-1) : t in Factorization(987654321)];
sigmavalue eq SumOfDivisors(987654321);
> sigmavalue:=&*[(t[1]^(t[2]+1)-1)/(t[1]-1) : t in Factorization(987654321)];
> sigmavalue eq SumOfDivisors(987654321);
true
> [(t[1]^(t[2]+1)-1)/(t[1]-1) : t in Factorization(10800)];
[ 31, 40, 31 ]
> [(2^5-1)/(2-1), (3^4-1)/(3-1), (5^3-1)/(5-1)];
[ 31, 40, 31 ]
> &*[(t[1]^(t[2]+1)-1)/(t[1]-1) : t in Factorization(10800)];
38440
> SumOfDivisors(10800);
38440
> NaiveDecoding([Modexp(m,d,n): m in encipheredmu]);
muvalue:=&*[-1*Floor(1/t[2]): t in Factorization(11111111111111)];
muvalue eq MoebiusMu(11111111111111);
/*
If x is a real number, Floor(x) is the integer part of x, which
means the largest integer less than or equal to x.
*/
> muvalue:=&*[-1*Floor(1/t[2]): t in Factorization(11111111111111)];
> muvalue eq MoebiusMu(11111111111111);
true
> Factorization(11111111111111);
[ <11, 1>, <239, 1>, <4649, 1>, <909091, 1> ]
> [-1*Floor(1/t[2]): t in Factorization(11111111111111)];
[ -1, -1, -1, -1 ]
> &*[-1*Floor(1/t[2]): t in Factorization(11111111111111)];
1
> MoebiusMu(11111111111111);
1
> [-1*Floor(1/t[2]): t in Factorization(2*3^3*5^2*7*11*13)];
[ -1, 0, 0, -1, -1, -1 ]
> &*[-1*Floor(1/t[2]): t in Factorization(2*3^3*5^2*7*11*13)];
0> NaiveDecoding([Modexp(m,d,n): m in mersenne]);
Mersenne numbers!
> /* In base 2 notation, 2^n-1 is a string of n 1's.
> For example, 2^5-1 = 2^4 + 2^3 + 2^2 + 2^1 + 2^0; so in base 2
> it is written as 11111.
> Despite this, my answer was wrong. A Mersenne number is a number of the
> form 2^p - 1 WHERE p IS PRIME. A base 2 repunit is any number of the
> form 2^n - 1.
> Now a repunit, in base 2 or in any other base, can only be prime if the
> number of digits is prime. So the only base 2 repunits that could
> conceivably be prime are the Mersenne numbers. And the base 2 repunits
> that are prime are exactly the Mersenne primes.
> */
> NaiveDecoding([Modexp(m,d,n): m in base3repunits]);
r:=1;
N:=0;
i:=1;
while N lt 7 do
  r:=3*r+1;
  i:=i+1;
  if IsPrime(i) then
    if IsProbablyPrime(r) then
      N:=N+1;
      "The base 3 repunit whose base 3 representation has",i,"digits is prime. Here it 
is:";
      r;
    end if;
  end if;
end while;

> r:=1;
> N:=0;
> i:=1;
> while N lt 7 do
while>   r:=3*r+1;
while>   i:=i+1;
while>   if IsPrime(i) then
while|if>     if IsProbablyPrime(r) then
while|if|if>       N:=N+1;
while|if|if>       "The base 3 repunit whose base 3 representation has",i,"digits is pri\
me. Here it
while|if|if> is:";
while|if|if>       r;
while|if|if>     end if;
while|if>   end if;
while> end while;
The base 3 repunit whose base 3 representation has 3 digits is prime. Here it
is:
13
The base 3 repunit whose base 3 representation has 7 digits is prime. Here it
is:
1093
The base 3 repunit whose base 3 representation has 13 digits is prime. Here it
is:
797161
The base 3 repunit whose base 3 representation has 71 digits is prime. Here it
is:
3754733257489862401973357979128773
The base 3 repunit whose base 3 representation has 103 digits is prime. Here it
is:
6957596529882152968992225251835887181478451547013
The base 3 repunit whose base 3 representation has 541 digits is prime. Here it
is:
66308439547181843673169185149975450335546563790474079717256778420962351634164023840051028\
85034630084417368521192210000003011029079481693984080146458143761582514901416066165280887\
59064995410612765793960501606910585490086339893058591064091241255832207903808201
The base 3 repunit whose base 3 representation has 1091 digits is prime. Here it
is:
17308478920290520561214081551955057184506477986582966242130879304040855943078281750655654\
65877316796954389481170318472442241905892654388323782720732502237578480086224100629849484\
07004907976010144446727825839493126338747635054243392729455736743049576570977803426383945\
19597261348848875851233994488843015578069982923778120456915440994423370412559032927349367\
94717885690746152667525928929590798694320443870675144426955467776443389120197038417789600\
1619872662271102048498620697981918654970696571273721026402375168960030732173
> NaiveDecoding([Modexp(m,d,n): m in amicable]);
n:=1;
f:=0;
while f lt 50 do
  n:=n+1;
  m:=SumOfDivisors(n)-n;
  if m ne n then
    if SumOfDivisors(m) eq n+m then
      f:=f+1;
      "amicable pair found:",m,n;
    end if;
  end if;
end while;
> n:=1;
> f:=0;
> while f lt 50 do
while>   n:=n+1;
while>   m:=SumOfDivisors(n)-n;
while>   if m ne n then
while|if>     if SumOfDivisors(m) eq n+m then
while|if|if>       f:=f+1;
while|if|if>       "amicable pair found:",m,n;
while|if|if>     end if;
while|if>   end if;
while> end while;
amicable pair found: 284 220
amicable pair found: 220 284
amicable pair found: 1210 1184
amicable pair found: 1184 1210
amicable pair found: 2924 2620
amicable pair found: 2620 2924
amicable pair found: 5564 5020
amicable pair found: 5020 5564
amicable pair found: 6368 6232
amicable pair found: 6232 6368
amicable pair found: 10856 10744
amicable pair found: 10744 10856
amicable pair found: 14595 12285
amicable pair found: 12285 14595
amicable pair found: 18416 17296
amicable pair found: 17296 18416
amicable pair found: 76084 63020
amicable pair found: 66992 66928
amicable pair found: 66928 66992
amicable pair found: 71145 67095
amicable pair found: 87633 69615
amicable pair found: 67095 71145
amicable pair found: 63020 76084
amicable pair found: 88730 79750
amicable pair found: 69615 87633
amicable pair found: 79750 88730
amicable pair found: 124155 100485
amicable pair found: 139815 122265
amicable pair found: 123152 122368
amicable pair found: 122368 123152
amicable pair found: 100485 124155
amicable pair found: 122265 139815
amicable pair found: 153176 141664
amicable pair found: 168730 142310
amicable pair found: 141664 153176
amicable pair found: 142310 168730
amicable pair found: 176336 171856
amicable pair found: 180848 176272
amicable pair found: 171856 176336
amicable pair found: 176272 180848
amicable pair found: 203432 185368
amicable pair found: 202444 196724
amicable pair found: 196724 202444
amicable pair found: 185368 203432
amicable pair found: 365084 280540
amicable pair found: 389924 308620
amicable pair found: 430402 319550
amicable pair found: 399592 356408
amicable pair found: 280540 365084
amicable pair found: 308620 389924
> /*
> I had to do "while f lt 50" rather than "while f lt 25" because
> my method was going to find each amicable pair twice.
> */
> UnsetLogFile();