> PrimesUpTo(1000);
[ 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 ]
> #PrimesUpTo(1000);
168
> #PrimesUpTo(2000);
303
> #PrimesUpTo(3000);
430
> #PrimesUpTo(4000);
550
> #PrimesUpTo(5000);
669
> #PrimesUpTo(6000);
783
> #PrimesUpTo(7000);
900
> #PrimesUpTo(8000);
1007
> #PrimesUpTo(9000);
1117
> #PrimesUpTo(10000);
1229
> #PrimesUpTo(10)*Log(10)/10;
0.921034037197618273607196581873
> #PrimesUpTo(100)*Log(100)/100;
1.15129254649702284200899572734
> #PrimesUpTo(1000)*Log(1000)/1000;
1.16050288686899902474506769316
> #PrimesUpTo(10000)*Log(10000)/10000;
1.13195083171587285826324459912
> #PrimesUpTo(100000)*Log(100000)/100000;
1.10431981059994431005502870167
> #PrimesUpTo(1000000)*Log(1000000)/1000000;
1.08448994777907958862426575926
> #PrimesUpTo(10000000)*Log(10000000)/10000000;
1.07117478896182292064729492007
> #PrimesUpTo(100000000)*Log(100000000)/100000000;
1.06129923175648075811311015652
> /* If they say that this sequence has limit equal to 1, then I believe
> them. But we have used a large value of n, and still Pi(n)Log(n)/n is
> not all that close to 1. */
> safeprimes:= function(N)
function> SP:={};
function> p:=3;
function> while p le N do
function|while> if IsPrime((p-1) div 2) then
function|while|if> Include(~SP,p);
function|while|if> end if;
function|while> p:=NextPrime(p);
function|while> end while;
function> return SP;
function> end function;
> for n:=1000 to 10000 by 1000 do
for> "When n is",n,"there are",#safeprimes(n),"safe primes less than n";
for> end for;
When n is 1000 there are 25 safe primes less than n
When n is 2000 there are 37 safe primes less than n
When n is 3000 there are 50 safe primes less than n
When n is 4000 there are 60 safe primes less than n
When n is 5000 there are 72 safe primes less than n
When n is 6000 there are 83 safe primes less than n
When n is 7000 there are 91 safe primes less than n
When n is 8000 there are 101 safe primes less than n
When n is 9000 there are 110 safe primes less than n
When n is 10000 there are 115 safe primes less than n
> safeprimes(1000);
{ 5, 7, 11, 23, 47, 59, 83, 107, 167, 179, 227, 263, 347, 359, 383, 467, 479, 503, 563, 
587, 719, 839, 863, 887, 983 }
> #PrimesUpTo(10^8);
5761455
> #PrimesUpTo(10^7);
664579
> #safeprimes(10^7);
30657
> for n:=1000 to 10000 by 1000 do
for> "When n is",n,"the approximation to the number of safe primes is",(n/Log(n))/Log(n/\
2);
for> end for;
When n is 1000 the approximation to the number of safe primes is 
23.2942809922057436653517233978
When n is 2000 the approximation to the number of safe primes is 
38.0914840235445736153903243034
When n is 3000 the approximation to the number of safe primes is 
51.2362184944720591708406411600
When n is 4000 the approximation to the number of safe primes is 
63.4494996450652928943178527997
When n is 5000 the approximation to the number of safe primes is 
75.0312376655233046942013989963
When n is 6000 the approximation to the number of safe primes is 
86.1431398178331847883782887007
When n is 7000 the approximation to the number of safe primes is 
96.8851449663931676720957532646
When n is 8000 the approximation to the number of safe primes is 
107.324556818455943910036673338
When n is 9000 the approximation to the number of safe primes is 
117.509477394852063683529749672
When n is 10000 the approximation to the number of safe primes is 
127.475822181930986440403606364
> C:=RealField()!(&*[p*(p-2)/(p-1)^2 : p in PrimesUpTo(10000) | p gt 2 ]);
> C;
0.660168296505503361995760236662
> n:=10^300;
> prob:=C/(Log(n/2))^2;
> prob;
1.38628583627580210489813405592E-6
> /* Notice the E-6 on the end of this, meaning "times 10 to the power -6".
> The probability that a randomly chosen number less than 10^300 is a safe
> prime is about 0.000001386..., which is rather close to zero.
> The probablity that a randomly chosen number less than 10^300 is NOT a safe
> prime is 1 minus the above --- about 0.999998614... (rather close to 1).
> So if k numbers less than 10^300 are chosen at random the probability
> that none of them are safe primes is about (0.999998614)^k.
> How large does k have to be to make this less than 0.5?  */
> 1-prob;
0.999998613714163724197895101866
> /* (1-prob)^k is less than 0.5 if Log((1-prob)^k) is less than Log(0.5)
> Now Log((1-prob)^k) equals k*Log(1-prob), and for this to be less than
> Log(0.5) we need k greater than Log(0.5)/Log(1-prob). (Note that Log(1-prob)
> and Log(0.5) are both negative numbers, since 1-prob and 0.5 are both 
> in the interval from 0 to 1. Remember that if an inequality like x < y is
> divided through by a negative number r then the inequlity is reversed: (x/r) > (y/r)
> */
> Log(0.5)/Log(1-prob);
500002.728132809881233820059365
> /* OK, I would have to make about 500000 random choices to have about
> 50% chance of getting a safe prime. */
> F:=FiniteField(97);
> a:=F!50;
> b:=F!77;
> a;
50
> b;
77
> /* a and b are 50 and 77 -- except that they are going to be treated
> as residues mod 97 rather than integers. Arithmetic involving a and b will 
> be residue arithmetic mod 97  */
> Parent(a);
Finite field of size 97
> /* This is saying that a belongs to the set of residues mod 97, which
> mathematicians refer to as the finite field with 97 elements. */
> Parent(50);
Integer Ring
> /* The ordinary number 50 belongs to the set of Integers (which is an
> example of a mathematical system called a "ring")  */
> a*b;
67
> 50*77;
3850
> 50*77 mod 97;
67
> (50+77) mod 97;
30
> a+b;
30
> a^111;
33
> 50^111;
38518598887744717061119558851698546370762032964307763904798775911331176757812500000000000\
00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000\
00000000000
> 50^111 mod 97;
33
> /* But 50^111 mod 97 is much less efficient than Modexp(50,111,97), which
> gets magma to use residue arithmetic rather than calculating the
> large number 50^111.  */
> Modexp(50,111,97);
33
> /* Of course, a^111 is exactly the same as Modexp(50,111,97) */
> b^(-1);
63
> InverseMod(77,97);
63
> P<x>:=PolynomialRing(F);
> f5:=(x-1)*(x-2)*(x-3)*(x-4)/((5-1)*(5-2)*(5-3)*(5-4));
> f4:=(x-1)*(x-2)*(x-3)*(x-5)/((4-1)*(4-2)*(4-3)*(4-5));
> f3:=(x-1)*(x-2)*(x-4)*(x-5)/((3-1)*(3-2)*(3-4)*(3-5));
> f2:=(x-1)*(x-3)*(x-4)*(x-5)/((2-1)*(2-3)*(2-4)*(2-5));
> f1:=(x-2)*(x-3)*(x-4)*(x-5)/((1-2)*(1-3)*(1-4)*(1-5));
> Evaluate(f5,1);
0
> Evaluate(f5,2);
0
> Evaluate(f5,3);
0
> Evaluate(f5,4);
0
> Evaluate(f5,5);
1
> [ [ Evaluate(fi,xj) : xj in [1..5] ] : fi in [f1,f2,f3,f4,f5] ];
[
    [ 1, 0, 0, 0, 0 ],
    [ 0, 1, 0, 0, 0 ],
    [ 0, 0, 1, 0, 0 ],
    [ 0, 0, 0, 1, 0 ],
    [ 0, 0, 0, 0, 1 ]
]
> a1:=73; a2:=50; a3:=36; a4:=82; a5:=17;
> f:=a1*f1+a2*f2+a3*f3+a4*f4+a5*f5;
> f;
15*x^4 + 4*x^3 + 42*x^2 + 83*x + 26
> /* This polynomial should take the value a1 at x = 1, a2 at x = 2, and so on */
> Evaluate(f,1);
73
> Evaluate(f,2);
50
> Evaluate(f,3);
36
> Evaluate(f,4);
82
> Evaluate(f,5);
17
> g:=33+Random(F)*x+Random(F)*x^2+Random(F)*x^3+Random(F)*x^4;
> g;
44*x^4 + 47*x^3 + 49*x^2 + 90*x + 33
> /* (I have chosen a random quartic polynomial with constant term 33) */
> 
> /* A polynomial (over a field) is completely determined by its
> values at d+1 points, if d is the degree of the polynomial. So
> for a quartic polynomial like g above, we can determine g once
> we know its values at 5 points.
> */
> a1:=Evaluate(g,1);
> a2:=Evaluate(g,2);
> a3:=Evaluate(g,3);
> a4:=Evaluate(g,4);
> a5:=Evaluate(g,5);
> a1;
69
> a2;
34
> a3;
48
> a4;
26
> a5;
66
> /* Putting these values into the Lagrange Interpolation Formula
> will give us the formual for g. */
> a1*f1+a2*f2+a3*f3+a4*f4+a5*f5;
44*x^4 + 47*x^3 + 49*x^2 + 90*x + 33
> /* Sure enough, that is the same thing we had above. */
> p:=447555834974539;
> F:=FiniteField(447555834974539);
> P<x>:=PolynomialRing(F);
> f4:=(x-1)*(x-3)/((4-1)*(4-3));
> f3:=(x-1)*(x-4)/((3-1)*(3-4));
> f1:=(x-3)*(x-4)/((1-3)*(1-4));
> f:=196231291191342*f1+195412581909834*f3+163633523397347*f4;
> Evaluate(f,0);
165270941960363
> /* Or, rather than work out f as I have done above, I could have
> just worked out f(0) directly: */
> m4:=F!((-1)*(-3)/((4-1)*(4-3)));
> m3:=F!((-1)*(-4)/((3-1)*(3-4)));
> m1:=F!((-3)*(-4)/((1-3)*(1-4)));
> 196231291191342*m1+195412581909834*m3+163633523397347*m4;
165270941960363
> UnsetLogFile();