> load "tut11data.txt";
Loading "L:\win\Magma\MATH2068\tut11data.txt"
> 1/(100*Log(10));
0.00434294481903251827651128918916
> prob := 1/(100*Log(10));
> /*
> Now the probability that a randomly chosen number less than 10^100 is NOT prime
> is approximately 1 - prob.
> So if we randomly choose k numbers less than 10^100 then the probability
> that none of them are prime is approximately (1 - prob)^k.
> How large does k have to be to make (1 - prob)^k less than 0.5?
> We need Log((1-prob)^k) < Log(0.5).
> That is, we need k*Log(1-prob) < Log(0.5).
> Dividing this inequality through by the NEGATIVE number Log(1-prob)
> we see that we need k > Log(0.5)/Log(1-prob)
> */
> 1-prob;
0.995657055180967481723488710810
> Log(0.5)/Log(1-prob);
159.256211525973648934737844198
> /*
> OK, choose 160 random numbers of 100 digits or less and you have
> about a 50% chance that at least one of them is prime.
> */
> p,b:=choosesafeprime(15);
> p,
> b;
25847 5
> p;
25847
> b;
5
> a:=Random(p-1);
> a;
5743
> prdl(a,b,p);
<5, 1, 0> <125, 3, 0>
<25, 2, 0> <3125, 5, 0>
<125, 3, 0> <15710, 12, 0>
<625, 4, 0> <16507, 48, 0>
<3125, 5, 0> <9875, 97, 0>
<15625, 6, 0> <12587, 194, 1>
<15710, 12, 0> <23394, 776, 4>
<16944, 24, 0> <20853, 776, 6>
<16507, 48, 0> <11042, 1552, 14>
<1975, 96, 0> <1478, 3105, 28>
<9875, 97, 0> <11103, 3107, 28>
<20741, 194, 0> <25216, 12428, 112>
<12587, 194, 1> <21317, 12428, 114>
<16306, 388, 2> <1574, 24856, 230>
<23394, 776, 4> <13503, 24858, 230>
<24883, 776, 5> <5508, 23871, 460>
<20853, 776, 6> <8465, 23873, 460>
<9628, 776, 7> <1749, 21902, 920>
<11042, 1552, 14> <19599, 17960, 1840>
<5465, 3104, 28> <8027, 17960, 1842>
<1478, 3105, 28> <7338, 10076, 3684>
<7390, 3106, 28> <18493, 20154, 7368>
<11103, 3107, 28> <17250, 20154, 7370>
<12266, 6214, 56> <6606, 20154, 7372>
....( 200 lines deleted ) ....
<14450, 1546, 15078> <11315, 25586, 14924>
<10434, 3092, 4310> <14177, 24806, 8004>
<792, 6184, 8620> <5285, 23767, 16008>
<3960, 6185, 8620> <2890, 23769, 16008>
<19800, 6186, 8620> <10434, 21694, 6170>
<10447, 6186, 8621> <3960, 17543, 12340>
<13775, 12372, 17242> <10447, 17544, 12341>
<7798, 24744, 8638> <7798, 18484, 23518>
The discrete log of 5743 modulo 25847 relative to the base 5 is
23571
> /*
> In the notation used in lectures, we compute a sequence of numbers
> x_1, x_2, x_3, ..., and also
> n_1, n_2, n_3, ... and m_1, m_2, m_3, ...
> 
> In each line of the output above, two triples are printed.
> The first is <x_i,n_i,m_i>, the second is <x_j,n_j,m_j>, where j equals 2*i.
> This is done for i = 1, then i = 2, then i = 3, and so on, until
> we find a case where x_i = x_j. 
> When this occurs we can use the values of n_i, m_i, n_j and m_j to work
> out how to express a as a power of b. 
> 
> See the lectures for details.
> 
> In the last line above we have the triples <7798, 24744, 8638> and <7798, 18484, 23518>
> This means that 7798 is congruent mod p to b^24744 * a^8638 and that also
> 7798 is congruent mod p to b^18484 * a^23518.
> */
> Modexp(b,24744,p)*Modexp(a,8638,p) mod p;
7798
> Modexp(b,18484,p)*Modexp(a,23518,p) mod p;
7798
> /* 
> So if a is congruent to b^i (mod p) then
> b^(24744+i*8638) is congruent mod p to b^(18484+i*23518).
> So 24744+i*8638 is congruent mod (p-1) to 18484+i*23518.
> The solution of this is that i must be congruent to 10648 mod (p-1)/2.
> In this example it can be shown by separate argument that  i has to be 
> odd for a to be congruent to b^i mod p. So in fact
> i = 10648 + (p-1)/2 = 10648 + 12923 =23571.
> So in fact Modexp(b,23571,p) should equal a.
> */
> Modexp(5,23571,25847);
5743
> p,b:=choosesafeprime(16);
> a:=Random(p-1);
> p;
3167
> b;
5
> a;
530
> prdl(a,b,p);
<5, 1, 0> <125, 3, 0>
<25, 2, 0> <3125, 5, 0>
<125, 3, 0> <2442, 5, 2>
<625, 4, 0> <1435, 5, 4>
<3125, 5, 0> <208, 11, 8>
<3076, 5, 1> <2033, 13, 8>
<2442, 5, 2> <770, 27, 16>
<2124, 5, 3> <248, 29, 16>
<1435, 5, 4> <1605, 60, 32>
<675, 10, 8> <1684, 240, 128>
<208, 11, 8> <3011, 960, 512>
<1040, 12, 8> <1379, 960, 514>
<2033, 13, 8> <2096, 674, 2056>
<154, 26, 16> <2935, 1349, 946>
<770, 27, 16> <2765, 1350, 947>
<683, 28, 16> <752, 1350, 949>
<248, 29, 16> <2965, 1352, 949>
<1240, 30, 16> <3090, 1353, 950>
<1605, 60, 32> <1805, 1354, 951>
<1254, 120, 64> <339, 2708, 1903>
<1684, 240, 128> <556, 2252, 640>
<1391, 480, 256> <745, 2253, 641>
<3011, 960, 512> <2790, 2255, 641>
<2829, 960, 513> <2013, 2255, 643>
<1379, 960, 514> <848, 2688, 2572>
<1441, 1920, 1028> <1708, 2212, 1978>
<2096, 674, 2056> <2285, 1259, 790>
<587, 1348, 946> <370, 2518, 1582>
<2935, 1349, 946> <2140, 1872, 3164>
<553, 1349, 947> <2070, 1873, 3165>
<2765, 1350, 947> <1473, 580, 3165>
<2296, 1350, 948> <1670, 1161, 3164>
<752, 1350, 949> <1204, 1478, 3158>
<593, 1351, 949> <1282, 2956, 3151>
<2965, 1352, 949> <205, 2746, 3137>
<618, 1352, 950> <1958, 2748, 3137>
<3090, 1353, 950> <334, 1494, 3050>
<361, 1353, 951> <1940, 2990, 2934>
<1805, 1354, 951> <2297, 2462, 2238>
<2349, 2708, 1902> <3018, 1758, 1312>
<339, 2708, 1903> <1025, 1759, 1313>
<1695, 2709, 1903> <1694, 354, 2626>
<556, 2252, 640> <1670, 709, 2086>
<2780, 2253, 640> <1204, 2836, 2012>
<745, 2253, 641> <1282, 2506, 859>
<558, 2254, 641> <205, 1846, 1719>
<2790, 2255, 641> <1958, 1848, 1719>
<2878, 2255, 642> <334, 1060, 544>
<2013, 2255, 643> <1940, 2122, 1088>
<1576, 1344, 1286> <2297, 2156, 1186>
<848, 2688, 2572> <3018, 1146, 2374>
<1073, 2689, 2572> <1025, 1147, 2375>
<1708, 2212, 1978> <1694, 2296, 1584>
<457, 1258, 790> <1670, 1427, 2>
<2285, 1259, 790> <1204, 2542, 8>
<1256, 1259, 791> <1282, 1918, 17>
<370, 2518, 1582> <205, 670, 35>
<1850, 2519, 1582> <1958, 672, 35>
<2140, 1872, 3164> <334, 2688, 140>
<414, 1872, 3165> <1940, 2212, 280>
<2070, 1873, 3165> <2297, 2516, 1120>
<3116, 580, 3164> <3018, 1866, 2242>
<1473, 580, 3165> <1025, 1867, 2243>
<334, 1160, 3164> <1694, 570, 1320>
<1670, 1161, 3164> <1670, 1141, 2640>
The discrete log of 530 modulo 3167 relative to the base 5 is
725
> Modexp(5,725,3167);
530
> p,b:=choosesafeprime(16);
> a:=Random(p-1);
> p;
3803
> b;
2
> a;
1681
> prdl(a,b,p);
<2, 1, 0> <8, 3, 0>
<4, 2, 0> <32, 5, 0>
<8, 3, 0> <128, 7, 0>
<16, 4, 0> <512, 9, 0>
<32, 5, 0> <2048, 11, 0>
<64, 6, 0> <3735, 22, 1>
<128, 7, 0> <2433, 22, 3>
<256, 8, 0> <19, 88, 12>
<512, 9, 0> <76, 90, 12>
<1024, 10, 0> <304, 92, 12>
<2048, 11, 0> <1216, 94, 12>
<3398, 22, 0> <959, 190, 24>
<3735, 22, 1> <1223, 382, 48>
<3585, 22, 2> <797, 766, 96>
<2433, 22, 3> <432, 1534, 192>
<2021, 44, 6> <1728, 1536, 192>
<19, 88, 12> <1258, 3073, 384>
<38, 89, 12> <2064, 2346, 768>
<76, 90, 12> <1472, 891, 1536>
<152, 91, 12> <2624, 1782, 3073>
<304, 92, 12> <295, 1782, 3075>
<608, 93, 12> <1180, 1784, 3075>
<1216, 94, 12> <2008, 3570, 2348>
<2432, 95, 12> <1768, 3339, 894>
<959, 190, 24> <119, 2876, 1789>
<1918, 191, 24> <476, 2878, 1789>
<1223, 382, 48> <1904, 2880, 1789>
<2446, 383, 48> <1914, 1959, 3578>
<797, 766, 96> <2214, 117, 3354>
<1594, 767, 96> <809, 234, 2907>
<432, 1534, 192> <1460, 470, 2012>
<864, 1535, 192> <1293, 1880, 444>
<1728, 1536, 192> <3737, 3718, 1776>
<629, 3072, 384> <2697, 3718, 1778>
<1258, 3073, 384> <962, 3719, 1779>
<2516, 3074, 384> <1457, 3638, 3558>
<2064, 2346, 768> <1550, 3475, 3314>
<736, 890, 1536> <2847, 3148, 2827>
<1472, 891, 1536> <786, 2494, 1854>
<2877, 1782, 3072> <3037, 1188, 3708>
<2624, 1782, 3073> <3697, 2376, 3616>
<3267, 1782, 3074> <1110, 2377, 3617>
<295, 1782, 3075> <3515, 954, 3432>
<590, 1783, 3075> <14, 954, 3434>
<1180, 1784, 3075> <56, 956, 3434>
<2360, 1785, 3075> <224, 958, 3434>
<2008, 3570, 2348> <896, 960, 3434>
<884, 3338, 894> <1532, 1922, 3066>
<1768, 3339, 894> <1146, 43, 2330>
<3561, 2876, 1788> <1321, 88, 858>
<119, 2876, 1789> <295, 176, 1717>
<238, 2877, 1789> <1180, 178, 1717>
<476, 2878, 1789> <2008, 358, 3434>
<952, 2879, 1789> <1768, 717, 3066>
<1904, 2880, 1789> <119, 1434, 2331>
<957, 1958, 3578> <476, 1436, 2331>
<1914, 1959, 3578> <1904, 1438, 2331>
<1107, 116, 3354> <1914, 2877, 860>
<2214, 117, 3354> <2214, 1953, 1720>
The discrete log of 1681 modulo 3803 relative to the base 2 is
164
> Modexp(2,164,3803);
1681
> p,b:=choosesafeprime(20);
> a:=Random(p-1);
> p;
471959
> b;
11
> a;
387374
> prdl(a,b,p);
<11, 1, 0> <1331, 3, 0>
<121, 2, 0> <161051, 5, 0>
<1331, 3, 0> <368978, 10, 1>
<14641, 4, 0> <284948, 20, 4>
<161051, 5, 0> <91333, 41, 8>
<445797, 10, 0> <196236, 43, 8>
<368978, 10, 1> <187099, 87, 16>
<172581, 10, 2> <318, 174, 33>
<284948, 20, 4> <38478, 176, 33>
<8303, 40, 8> <115933, 177, 34>
<91333, 41, 8> <468390, 178, 35>
<60745, 42, 8> <242503, 356, 72>
<196236, 43, 8> <388305, 1424, 288>
<17009, 86, 16> <47291, 2848, 578>
<187099, 87, 16> <58703, 2850, 578>
<364812, 174, 32> <50379, 5702, 1156>
<318, 174, 33> <432351, 5704, 1156>
<3498, 175, 33> <438399, 11408, 2314>
<38478, 176, 33> <262000, 22816, 4630>
<423258, 177, 33> <53671, 45632, 9261>
... (1081 lines deleted) ....
<180580, 260062, 138201> <15081, 289890, 401455>
<73213, 48166, 276402> <366550, 107824, 330952>
<333384, 48167, 276402> <330680, 215648, 189948>
<264610, 48167, 276403> <433726, 215649, 189949>
<30737, 96334, 80848> <355648, 215650, 189950>
<338107, 96335, 80848> <73213, 431300, 379902>
<46969, 96335, 80849> <264610, 431301, 379903>
<44700, 96336, 80849> <338107, 390645, 287848>
<19741, 96337, 80849> <44700, 390646, 287849>
<217151, 96338, 80849> <217151, 390648, 287849>
The discrete log of 387374 modulo 471959 relative to the base 11 is
278032
> Modexp(11,278032,471959);
387374
> Modexp(3,2^(2^1),2^(2^1)+1);
1
> Modexp(3,2^(2^2),2^(2^2)+1);
1
> Modexp(3,2^(2^3),2^(2^3)+1);
1
> Modexp(3,2^(2^4),2^(2^4)+1);
1
> /*
> In each case we have done a command of the form Modexp(3,n-1,n),
> that is, we have computed 3^(n-1) mod n. If n is a prime number
> then Fermat's Little Theorem guarantees that the answer will be 1.
> So if the numbers n that we have used are prime then the above results
> are not at all surprising!
> */
> 2^(2^1)+1;
5
> /* That is prime! */
> 2^(2^2)+1;
17
> /* So is that! */
> 2^(2^3)+1;
257
> /* That is too! */
> 2^(2^4)+1;
65537
> /* ... could be ... */
> IsPrime(65537);
true
> /* It is! */
> 2^(2^5)+1;
4294967297
> IsPrime(2^(2^5)+1);
false
> /* since this value of n is not prime it is highly unlikely
> that 3^(n-1) will be 1 modulo n.
> I can't guess what it's value will be, but I can confidently
> guess that it will not be 1.
> */
> Modexp(3,2^(2^5),2^(2^5)+1);
3029026160
> for n:=3 to 999 by 2 do
for> if Modexp(2,n-1,n) eq 1 then
for|if> if not IsPrime(n) then
for|if|if> print n;
for|if|if> end if;
for|if> end if;
for> end for;
341
561
645
> /*
> Since (2^170)^2-1 is divisible by 341, (2^170 - 1)*(2^170 + 1) is divisible
> by 341. If 341 were prime it would follow that one or other of these factors
> is divisible by 341.
> */
> Modexp(2,170,341);
1
> /*
> OK, 2^170 - 1 is divisible by 341. If 341 were prime it would follow that
> either 2^85 - 1 or 2^85 + 1 is divisible by 341.
> */
> Modexp(2,85,341);
32
> /*
> So neither 2^85 - 1 nor 2^85 + 1 is congruent to 0 mod 341.
> So 341 is not prime. (Of course, it is 31*11.)
> */
> Modexp(2,560,561);
1
> /* (As we knew) */
> Modexp(2,280,561);
1
> Modexp(2,140,561);
67
> /* So 561 is not prime. */
> Modexp(2,644,645);
1
> /* (As we knew) */
> Modexp(2,322,645);
259
> /* So 645 is not prime */
> m:=25326001;
> Factorization(m-1);
[ <2, 4>, <3, 3>, <5, 3>, <7, 1>, <67, 1> ]
> /*
> The highest power of 2 that divides m-1 is 2^4, i.e. 16.
> */
> Modexp(2,m-1,m);
1
> /*
> This is consistent with the possibility that m is prime.
> We say that m is a pseudoprime to the base 2.
> 
> If m were actually prime, we could deduce from the fact
> that 2^(m-1) is congruent to 1 mod m that 2^((m-1)/2) is congruent to
> either 1 or to -1 mod m. Let us check this.
> */
> Modexp(2,(m-1)/2,m);

>> Modexp(2,(m-1)/2,m);
         ^
Runtime error in 'Modexp': Bad argument types
Argument types given: RngIntElt, FldRatElt, RngIntElt
> /*
> Magma doesn't like us using (m-1)/2 in the Modexp command
> because it thinks of (m-1)/2 as a rational number rather
> than an integer. But we can use (m-1) div 2 instead.
> */
> (m-1)/2;
12663000
> (m-1) div 2;
12663000
> Modexp(2,(m-1) div 2,m);
1
> /* This is still consistent with m being prime.
> Note that (m-1) div 2 is still even. If m really were
> prime then 2^((m-1)/2) congruent to 1 would imply 2^((m-1)/4) congruent
> to 1 or to -1.
> */
> Modexp(2,(m-1) div 4,m);
1
> /*
> This is still consistent with the possibility that m is prime.
> And (m-1) div 4 is still even; so we can divide by 2 again.
> */
> Modexp(2,(m-1) div 8,m);
1
> /*
> Still consistent with m being prime. And (m-1) div 8 is still even.
> */
> Modexp(2,(m-1) div 8,m);
1
> Modexp(2,(m-1) div 16,m);
25326000
> /*
> i.e. 2^((m-1)/16) is congruent to -1 mod m. This is still consistent
> with m being prime, and there are no further tests we can do just
> using powers of 2.
> So we say that m is a strong pseudoprime to the base 2.
> */
> Modexp(3,m-1,m);
1
> /* m is a pseudoprime to the base 3 */
> Modexp(3,(m-1) div 2,m);
1
> /* Consistent with m being prime, but we can test it further using
> 3^((m-1)/4).
> */
> Modexp(3,(m-1) div 4,m);
1
> /* Still consistent with m being prime, but we can go on ... */
> Modexp(3,(m-1) div 8,m);
1
> /* Still consistent with m being prime */
> Modexp(3,(m-1) div 16,m);
25326000
> /* Still consistent with m being prime, and we have run out of
> tests we can do using powers of 3.
> m is a strong pseudoprime to the base 3.
> */
> Modexp(5,m-1,m);
1
> /* m is a pseudoprime to the base 5. */
> Modexp(5,(m-1) div 2,m);
1
> Modexp(5,(m-1) div 4,m);
1
> Modexp(5,(m-1) div 8,m);
1
> Modexp(5,(m-1) div 16,m);
1
> /* This is a s far as we can go because (m-1) div 16 is odd.
> m is a strong pseudoprime to the base 5. 
> */
> Modexp(7,m-1,m);
5872860
> /*
> Not consistent with Fermat's Little Theorem. 
> m is not prime, because it is not even a pseudoprime to the base 7.
> */
> for n:=3 to 9999 by 2 do
for> if not IsPrime(n) then
for|if> if Modexp(2,n-1,n) eq 1 then
for|if|if> if Modexp(3,n-1,n) eq 1 then
for|if|if|if> n;
for|if|if|if> end if;
for|if|if> end if;
for|if> end if;
for> end for;
1105
1729
2465
2701
2821
6601
8911
> /*
> These numbers are all pseudoprimes to the base 2, but we will
> show that they are not strong pseudoprimes to the base 2.
> */
> Modexp(2,1104 div 2,1105);
1
> 1104 div 2;
552
> 1104 div 4;
276
> Modexp(2,276,1105);
781
> /* Since 781 is not 1 or -1 mod 1105, we conclude that 1105 is not
> a strong pseudoprime to the base 2.
> */
> Modexp(2,1728,1729);
1
> Modexp(2,864,1729);
1
> Modexp(2,432,1729);
1
> Modexp(2,216,1729);
1
> Modexp(2,108,1729);
1
> Modexp(2,54,1729);
1065
> /* This is not 1 or -1. So 1729 is not a strong pseudoprime to the base 2. */
> Modexp(2,2464.2465);

>> Modexp(2,2464.2465);
         ^
Runtime error in 'Modexp': Bad argument types
Argument types given: RngIntElt, FldReElt
> Modexp(2,2464,2465);
1
> Modexp(2,1232,2465);
1
> Modexp(2,616,2465);
1
> Modexp(2,308,2465);
1886
> /* This is not 1 or -1. So 2465 is not a strong pseudoprime to the base 2. */
> Modexp(2,2700,2701);
1
> Modexp(2,1350,2701);
147
> /* This is not 1 or -1. So 2701 is not a strong pseudoprime to the base 2. */
> Modexp(2,2820,2821);
1
> Modexp(2,1410,2821);
1520
> /* This is not 1 or -1. So 2821 is not a strong pseudoprime to the base 2. */
> Modexp(2,6600,6601);
1
> Modexp(2,3300,6601);
1
> Modexp(2,1650,6601);
4509
> /* This is not 1 or -1. So 6601 is not a strong pseudoprime to the base 2. */
> Modexp(2,8910,8911);
1
> Modexp(2,4455,8911);
6364
> /* This is not 1 or -1. So 8911 is not a strong pseudoprime to the base 2. */
> 8910 div 2;
4455
> Modexp(3,2^(2^5),2^(2^5)+1);
3029026160
> Modexp(3,2^(2^6),2^(2^6)+1);
8752249535465629170
> Modexp(3,2^(2^7),2^(2^7)+1);
47511664169441434718291075092691853899
> Modexp(3,2^(2^8),2^(2^8)+1);
113080593127052224644745291961064595403241347689552251078258028018246279223993
> Modexp(3,2^(2^9),2^(2^9)+1);
13387457852131866017809974335626508736765841341908171621341620739066502578793457441078230\
804865246011339933833061458906559278633032869468345609327807927612
> Modexp(3,2^(2^10),2^(2^10)+1);
14680952897024403602482100551468347348324469982504648594298872798304927867143532626816508\
55751818323747269129462800188929193889626999934059088187841777884701782257649364274500395\
82955979907044211973217849267621360435307454207055855850222255439616231847463892347771894\
224104545724927599853426578286430466818141
> /*
> By Fermat's Little Theorem, none of these Fermat numbers are 
> prime. In each case, 3 is a witness to it.
> */
> UnsetLogFile();