> n:=CRT([573,2],[1001,999]);
> n;
214787
> n mod 1001;
573
> n mod 999;
2
> /* It checks! */
> m:= 425257806103696493;
> Factorization(m);
[ <463262383, 1>, <917963171, 1> ]
> p:=Factorization(m)[1,1];
> q:=Factorization(m)[2,1];
> p;
463262383
> q;
917963171
> n:=361971017277558996;
> n mod p;
16
> n mod q;
36
> /* So the square roots of n mod p are 4 and -4, and the square roots
> of n mod q are 6 and -6. */
> sqrt1:=CRT([4,6],[p,q]);
> sqrt2:=CRT([-4,6],[p,q]);
> sqrt3:=CRT([4,-6],[p,q]);
> sqrt4:=CRT([-4,-6],[p,q]);
> sqrt1;
36197101727755902
> sqrt2;
180985508638779486
> sqrt3;
244272297464917007
> sqrt4;
389060704375940591
> Modexp(sqrt1,2,m);
361971017277558996
> Modexp(sqrt2,2,m);
361971017277558996
> Modexp(sqrt3,2,m);
361971017277558996
> Modexp(sqrt4,2,m);
361971017277558996
> /* We have indeed found four square roots of n */
> k:=139642041051408538;
> sp1:=Modexp(k, (p+1) div 4, p);
> sq1:=Modexp(k, (q+1) div 4, q);
> sp1;
171119576
> sq1;
691671103
> squareroot1:=CRT([sp1,sq1],[p,q]);
> squareroot2:=CRT([-sp1,sq1],[p,q]);
> squareroot3:=CRT([sp1,-sq1],[p,q]);
> squareroot4:=CRT([-sp1,-sq1],[p,q]);
> Modexp(squareroot1,2,m);
139642041051408538
> Modexp(squareroot2,2,m);
139642041051408538
> Modexp(squareroot3,2,m);
139642041051408538
> Modexp(squareroot4,2,m);
139642041051408538
> /* It worked, we got four square roots of k. */
> load "tut12data.txt";
Loading "L:\win\Magma\MATH2068\tut12data.txt"
> p1;
98193403447007846364293549496534840626586107516549603723217892945966367816514457509916143\
00548372406884587748387062551285546165836862826501114696264357987998850335889756854554162\
9801435915684340758467
> p2;
21113364761575055130453152887554511875646730393620606418588730835077058605450797046261043\
11278484682177015020265630401938641762620145861114474213244984415671567916416192451305701\
8918326407315610702719
> p1 mod 4;
3
> p2 mod 4;
3
> message;
[ 289670920632665415146023305042226552693218914161246512772055900330540760736835879927322495465050285551\
45080936642242897437461748000789241785040438500804447624424607886412731896010957111015170263369055557176\
30439849103681399901538447053099828878983282788320636760613854941327093012128378169565575563619073761614\
3229239160530648919140398259123589507989007988305925916516517517965760763988935073796, 
17478219172562037232986477227997605794911353656937344929804479154703815613797749018323153379869743653153\
527816674906227026382346347520369453959166651459144310776107529758045800444871202397617234859713316 ]
> message[1];
28967092063266541514602330504222655269321891416124651277205590033054076073683587992732249546505028555145\
08093664224289743746174800078924178504043850080444762442460788641273189601095711101517026336905555717630\
43984910368139990153844705309982887898328278832063676061385494132709301212837816956557556361907376161432\
29239160530648919140398259123589507989007988305925916516517517965760763988935073796
> message[2];
17478219172562037232986477227997605794911353656937344929804479154703815613797749018323153379869743653153\
527816674906227026382346347520369453959166651459144310776107529758045800444871202397617234859713316
> m11:=Modexp(message[1], (p1+1) div 4, p1);
> m12:=Modexp(message[1], (p2+1) div 4, p2);
> m1:=CRT([m11,m12],[p1,p2]);
> NaiveDecoding([m1]);
ERROR! -52 is not a valid ascii code number.

     You have not calculated the encoded form of the plaintext correctly.

> /* OK , I must have chosen the wrong square root. */
> m1:=CRT([m11,-m12],[p1,p2]);
> NaiveDecoding([m1]);
ERROR! 315 is not a valid ascii code number.

     You have not calculated the encoded form of the plaintext correctly.

> /* Still wrong! There are still two more possibilities ... */
> m1:=CRT([-m11,-m12],[p1,p2]);
> NaiveDecoding([m1]);
ERROR! 389 is not a valid ascii code number.

     You have not calculated the encoded form of the plaintext correctly.
> /* The last one had better be right! */
> m1:=CRT([-m11,m12],[p1,p2]);
> NaiveDecoding([m1]);
Fashion is a form of ugliness so intolerable that we have to alter
> /* Sounds like something that Oscar Wilde might have said ... */
> m21:=Modexp(message[2], (p1+1) div 4, p1);
> m22:=Modexp(message[2], (p2+1) div 4, p2);
> m2:=CRT([m21,m22],[p1,p2]);
> NaiveDecoding([m2]);
ERROR! 687 is not a valid ascii code number.

     You have not calculated the encoded form of the plaintext correctly.

> /* OK , I'll guess again. */
> m2:=CRT([m21,-m22],[p1,p2]);
> NaiveDecoding([m2]);
 it every six months.
Oscar Wilde.
> /* Bingo! */
> NaiveDecoding([m1,m2]);
Fashion is a form of ugliness so intolerable that we have to alter it every six months.
Oscar Wilde.
> Modexp(500,32768,65537);
65536
> /* Note that if p = 65537 then (p-1)/2 = 32768. So what the above calculation tells 
> is that 500 raised to the power (p-1)/2 is congruent to 65536 mod p.
> That is, 500^((p-1)/2) is congruent to -1 mod p (since 65536 is -1 mod 65537).
> Recall that a is a nonsquare mod p if a^((p-1)/2) is -1 mod p. So we have shown
> that 500 is a nonsquare.
> 
> By Fermat's Little Theorem we know that the order of 500 mod 65537 has to
> be a divisor of 65536. Since 65536 is a power of 2 its only divisors are smaller
> powers of 2. So the only divisor of 65536 that is not also a divisor of 32768
> is 65536 itself. Now we know that the order of 500 mod 65537 is not a divisor
> of 32768 since 500^32768 is not congruent to 1 mod p -- it is congruent to -1
> mod p. So the order of 500 mod p must be 65536, which means that 500 is a primitive
> root mod p. */
> 
> /* Since 500 is a primitive root there must be an i in the set {0,1,2,...,65535}
> such that 23 = 500^i (using residue arithmetic mod 65537). Our aim is to find i. */
> 
> Modexp(23,32768,65537);
65536
> /* So (working in residue arithemetic, and remembering that 65536 = -1) we have
> -1 = 23^32768 = (500^i)^32768 = (500^32768)^i = (-1)^i.
> This shows us that i is odd. We put i = 2*j + 1. 
> Now 23 = = 500^i = 500^(2*j+1) gives 23 * 500^(-1) = 500^(2*j). 
> Let us calculate 23 * 500^(-1):   */
> 23*Modexp(500,-1,65537) mod 65537;
15860
> 
> /* OK, 15860 = 500^(2*j). So 15860^16384 = (500^(2*j))^16384 = (500^32768)^j = (-1)^j. */
> Modexp(15860,16384,65537);
65536
> /* So -1 = 15860^16384 = (-1)^j. So j is odd.
> Let j = 2*k + 1. 
> Remember that we are trying to find i, and that i = 2*j + 1. So now
> we have i = 2*(2*k+1)+1 = 4*k + 3.
> And 23 = 500^i = 500^(4*k+3); so 23 * 500^(-3) = 500^(4*k)  */
> 23*Modexp(500,-3,65537) mod 65537;
57206
> /* OK, 57206 = 500^(4*k). So 57206^8192 = (500^(4*k))^8192 = (500^32768)^k = (-1)^k. */
> Modexp(57206,8192,65537);
1
> /* So 1 = 57206^8192 = (-1)^k. So k is even.
> Let k = 2*m. Thus i = 4*k +3 = 8*m +3.
> And 23 * 500^(-3) = 500^(8*m).
> We already know that 23 * 500^(-3) is 57206. So 57206 = 500^(8*m), and
> 57206^4096 = (500^(8*m))^4096 = (500^32768)^m = (-1)^m   */
> Modexp(57206,4096,65537);
1
> /* So 1 = 57206^4096 = (-1)^m. So m is even.
> Let m = 2*n. Thus i = 8*m +3 = 16*n + 3.
> And 23 * 500^(-3) = 500^(16*n).
> So 57206 = 500^(16*n), and
> 57206^2048 = (500^(16*n))^2048 = (500^32768)^n = (-1)^n   */
> Modexp(57206,2048,65537);
1
> /* So 1 = 57206^2048 = (-1)^n. So n is even.
> Let n = 2*l. Thus i = 32*l +3.
> And 23 * 500^(-3) = 500^(32*l).
> So 57206 = 500^(32*l), and
> 57206^1024 = (500^(32*l))^1024 = (500^32768)^l = (-1)^l   */
> Modexp(57206,1024,65537);
1
> /* So 1 = 57206^1024 = (-1)^l. So l is even.
> Let l = 2*q. Thus i = 64*q +3.
> And 23 * 500^(-3) = 500^(64*q).
> So 57206 = 500^(64*q), and
> 57206^512 = (500^(64*q))^512 = (500^32768)^q = (-1)^q   */
> Modexp(57206,512,65537);
65536
> /* So -1 = 57206^512 = (-1)^q. So q is odd.
> Let q = 2*r + 1. Thus i = 64*(2*r+1) +3 = 128*r + 67.
> And 23 = 500^(128*r+67), giving 23 * 500^(-67) = 500^(128*r). */
> 23*Modexp(500,-67,65537) mod 65537;
49554
> /* So 49554 = 500^(128*r) and
> 49554^256 = (500^(128*r))^256 = (500^32768)^r = (-1)^r   */
> Modexp(49554,256,65537);
65536
> /* So -1 = 49554^256 = (-1)^r. So r is odd.
> Let r = 2*s + 1. Thus i = 128*(2*s+1) + 67 = 256*s + 195.
> And 23 = 500^(256*s+195), giving 23 * 500^(-195) = 500^(256*s). */
> 23*Modexp(500,-195,65537) mod 65537;
9589
> /* So 9589 = 500^(256*s) and
> 9589^128 = (500^(256*s))^128 = (500^32768)^s = (-1)^s   */
> Modexp(9589,128,65537);
1
> /* So 1 = 9589^128 = (-1)^s. So s is even.
> Let s = 2*t. Thus i = 512*t + 195.
> And 23 * 500^(-195) = 500^(512*t). 
> So 9589 = 500^(512*t) and
> 9589^64 = (500^(512*t))^64 = (500^32768)^t = (-1)^t   */
> Modexp(9589,64,65537);
65536
> /* So -1 = 9589^64 = (-1)^t. So t is odd.
> Let t = 2*u + 1. Thus i = 512*(2*u+1) + 195 = 1024*u + 707.
> And 23 = 500^(1024*u+707), giving 23 * 500^(-707) = 500^(1024*u). */
> 23*Modexp(500,-707,65537) mod 65537;
16
> /* So 16 = 500^(1024*u) and
> 16^32 = (500^(1024*u))^32 =(500^32768)^u = (-1)^u   */
> Modexp(16,32,65537);
1
> /* So 1 = 16^32 = (-1)^u. So u is even.
> Let u = 2*v. Thus i = 2048*v + 707.
> And 23 * 500^(-707) = 500^(2048*v).
> So 16 = 500^(2048*v) and
> 16^16 = (500^(2048*v))^16 = (500^32768)^v = (-1)^v  */
> Modexp(16,16,65537);
1
> /* So 1 = 16^16 = (-1)^v. So v is even.
> Let v = 2*w. Thus i = 4096*w + 707.
> And 23 * 500^(-707) = 500^(4096*w).
> So 16 = 500^(4096*w) and
> 16^8 = (500^(4096*w))^8 = (500^32768)^w = (-1)^w  */
> Modexp(16,8,65537);
1
> /* So 1 = 16^8 = (-1)^w. So w is even.
> Let w = 2*x. Thus i = 8192*x + 707.
> And 23 * 500^(-707) = 500^(8192*x).
> So 16 = 500^(8192*x) and
> 16^4 = (500^(8192*x))^4 = (500^32768)^x = (-1)^x  */
> Modexp(16,4,65537);
65536
> /* So -1 = 16^4 = (-1)^x. So x is odd.
> Let x = 2*y + 1. Thus i = 8192*(2*y+1) + 707 = 16384*y + 8899.
> And 23 = 500^(16384*y+8899), giving 23 * 500^(-8899) = 500^(16384*y). */
> 23*Modexp(500,-8899,65537) mod 65537;
65281
> /* So 65281 = 500^(16384*y) and
> 65281^2 = (500^(16384*y)^2 = (500^32768)^y = (-1)^y  */
> Modexp(65281,2,65537);
65536
> /* So -1 = 65281^2 = (-1)^y. So y is odd.
> Let y = 2*z + 1. Thus i = 16384*(2*z+1) + 8899 = 32768*z + 25283.
> And 23 = 500^(32768*z+25283), giving 23 * 500^(-25283) = 500^(32768*z). 
> We had better evaluate 23 * 500^(-25283). */
> 23*Modexp(500,-25283,65537) mod 65537;
1
> /* This tells us that we are finished. 23 * 500^(-25283) = 1
> says that 23 = 500^25283. So the number i we are looking for
> (the solution of 23 = 500^i) is i = 25283. */
> 
> /* The above seems a bit tedious by hand, but actually it is
> computationally efficient. We computed the residue of i modulo 2,
> then computed the residue of i modulo 4, then modulo 8, then
> modulo 16, and so on up to 65536. Since the modulus doubles at
> every step and 65536 = 2^16, it takes 16 steps to reach the end.
> Now 16 is the log (to the base 2) of 65536. The number of steps
> is directly proportional to the logarithm of the number we are
> dealing with.
> If we had to use 65536 steps, that would have been very bad.
> If the number of steps were the square root of 65536, that also would
> have been very bad. But the number of steps is proportional to
> the log, which is proportional to the number of bits. 
> 65536 is a 16 bit number.
> Our algorithm is thus a polynomial time algorithm, which
> is computationally good. */
> 
> UnsetLogFile();