Chapter One:

Prime Numbers

1.1  Recognizing Large Primes

Mod[100!, 101]

Mod[90!, 91]

Mod[2^100, 101]

PowerMod[2, 100, 101]

PowerMod[17, -1, 10001]

Mod[17 %, 10001]

For[n = 3, n < 3000, n = n + 2,
If[PowerMod[2, n - 1, n] == 1 && !PrimeQ[n], Print[n]]]

Do[If[PowerMod[2, n - 1, n] == 1 && !PrimeQ[n], Print[n]],
   {n, 3, 3000, 2}]

FactorInteger[2701]

PrimeQ[1000001]

PrimeQ[10^100 + 1]

?PrimeQ

FactorInteger[2^67 - 1]

FactorInteger[1000001] //Timing

n = 10^30 + 1;
While[!PrimeQ[n], n = n + 2];
n

FirstPrimeAbove[n_] := Block[{k},
		k = n;
		While[!PrimeQ[k], k = k + 1];
		Return[k]]

FirstPrimeAbove[234436]

FirstPrimeAbove[10^30]

Following requires primeQ in Appendix

FirstPrimeAbove[n_] := Block[{k = n},
		While[!primeQ[k], k++];
		k]

FirstPrimeAbove[10^100]

Version 2.0 notes.

In version 2.0 there is a NextPrime[] function provided in the package NumberTheoryFunctions.m. That package contains lots of other useful number theory functions too.

1.2  Public-Key Encryption

p = FirstPrimeAbove[12400123056421444]

q = FirstPrimeAbove[24760095346677660]

n = p q

exponent = 230959;
GCD[exponent, (p-1)(q-1)]

message = 296776666357357111357;

PowerMod[message, exponent, n]

d = PowerMod[exponent, -1, (p - 1)(q - 1)]

PowerMod[%%, d, n]

following cells require that encode and decode from Appendix be loaded

PowerMod[encode["Sell shares now"], exponent, n]

decode[PowerMod[%, d, n]]

1.3  The Distribution of the Primes

Version 2.0 notes:

PrimePi is built in. Also RiemannR, but it is hidden in PrimePi.m and unavailable globally. To make it available, go
into PrimePi.m and add a usage statement for RiemannR
immediately following the usage statement for PrimePi. Then
RiemannR will be available after PrimePi[] has been invoked once.
Both Prime and PrimePi will take arbitrarily large inputs, limited only by your patience.

Prime[10^8]

Table[Prime[n], {n, 100}]

primes = Table[Prime[n], {n, 100}];

Short[primes]

Range[10]

Range[3, 70, 5]

PrimePi from Appendix (due to Ilan Vardi)

PrimePi::usage = "PrimePi[x] returns the number of primes 
less than or equal to x";

PrimePi[x_] := PrimePi[Floor[x]] /; Not[IntegerQ[x]]

PrimePi[n_Integer] := 0 /; n <= 1

PrimePi[4] = 2;

PrimePi[n_Integer] := Block[{t, diff, m},
	t = Round[N[Log[n]]];
	m = Round[N[LogIntegral[n]]];
	diff = n - Prime[m];
	While[ Abs[diff] > t^2 / (1.16)^2,
		m += Round[ diff / t];
		diff = n - Prime[m];] 
	If[diff < 0,
		While[n < Prime[m], m--;],
		While[n >= Prime[m], m++;]; m--];
	Return[m]]  /; n > 1 && n != 4

Attributes[PrimePi] = Listable;

(* Needs PrimePi from Appendix (built into version 2.0 *)
Plot[PrimePi[x], {x, 1, 200}]

Plot[{x^2 Sin[1/x], x^2, -x^2}, {x, -.1, .1}]

Table[{x, PrimePi[x]}, {x, 200}];

Short[%]

Short[Table[{Prime[n], n - 1}, {n, PrimePi[200]}]]

primepoints = Union[%, %%];
Take[primepoints, 10]

ListPlot[primepoints, PlotJoined->True]

{Range[100], PrimePi[Range[100]]};

Transpose[{Range[200], PrimePi[Range[200]]}]  ~Union~
Transpose[
 {Prime[Range[PrimePi[200]]], Range[0, PrimePi[200]-1]}];

Plot[PrimePi[x], {x, 1, 10^6}, PlotStyle->Thickness[.0005]]

Plot[{PrimePi[x], x/Log[x]}, {x, 2, 50000},
  PlotStyle->{{Thickness[.0005]},
			  {Thickness[.0001], Dashing[{.02, .04}]}}]

Plot[{x/(Log[x] - 1.08366) - PrimePi[x],
      LogIntegral[x] - PrimePi[x],
      RiemannR[x] - PrimePi[x]}, {x, 1, 10^7}]

RiemannR  (from Appendix)

RiemannR::usage = "Riemann[x] returns the value of Riemann's
prime-approximating function R(x).";

rlist = Map[1/(#! # N[Zeta[# + 1]]) &, Range[100]];
RiemannR[x_] := 
	Round[1 + (Floor[N[Log[x]] ^ Range[100]]) . rlist]
RiemannR[0] = 0;
Attributes[RiemannR] = {Listable};

(* Needs PrimePi and RiemannR from Appendix *)													


step = 25000 //N;
domain = Range[step, 10^7, step];
options = {PlotJoined->True};

primepilist = Table[PrimePi[x], {x, step, 10^7, step}];
Print["Done primepi"];
legendre = Table[x/(Log[x] - 1.08366), {x, step, 10^7, step}];
Print["Done legendre"];
gauss = LogIntegral[domain];
Print["Done gauss"];
riemann = Table[If[Mod[x, 500000] == 0, Print@x]; RiemannR[x],
	 						{x, step, 10^7, step}];

legendre = Transpose[{domain, legendre - primepilist}];
gauss    = Transpose[{domain,  gauss   - primepilist}];
riemann  = Transpose[{domain, riemann  - primepilist}];

Show[{ListPlot[legendre, options,
		   PlotStyle->{Thickness[1.5/420], GrayLevel[0.5]}],
	 ListPlot[gauss, options, PlotStyle->Thickness[.5/420]],
	 ListPlot[riemann, options, PlotStyle->Thickness[.5/420]]}]

ListPlot[Table[{x, LogIntegral[x] - PrimePi[x]},
			  {x, 0., 2 10^9, 2 10^7}],     PlotJoined->True]

Counting Primes Modulo 4 and 11

PrimePiModn (from Appendix)

PrimePiModn::usage = "PrimePiModn[x, n] divides the first x
primes into the nonzero mod-n residue classes and pairs 
each prime with the mod-n count up to that prime.";

PrimePiModn[x_, n_] := Block[
	{i, temp, modlist = Mod[Prime[Range[x]], n]},
    Table[(temp = Flatten[Position[modlist, i]];
           Transpose[{Prime[temp], Range[Length[temp]]}]),
       {i, n - 1}]]

PrimePiModn[40, 11] // TableForm

list = PrimePiModn[500, 11];
Show[Table[ListPlot[list[[m]], PlotJoined->True], {m, 10}]];

list = PrimePiModn[669, 11];
Show[
	Table[ListPlot[list[[m]], PlotJoined->True,
                              DisplayFunction->Identity],
       {m, 10}], 
  DisplayFunction->$DisplayFunction];

PrimePiMod4 exercise (from Appendix)

ListPlot[Accumulate[Plus, Mod[Prime[Range[3000]], 4] - 2], 
	PlotJoined->True];

1.4  Recursion and Caching

euclid[n_] := euclid[n-1] Prime[n]
euclid[1] = 2;
Table[euclid[n] + 1, {n, 20}]

euclid[n_] := (Print[n]; euclid[n - 1] Prime[n])
euclid[1] = 2;
Table[euclid[n], {n, 8}]

Clear[euclid];
euclid[n_] := (Print[n]; euclid[n] = euclid[n - 1] Prime[n])
euclid[1] = 2;
Table[euclid[n], {n, 8}]

Clear[euclid];
Attributes[euclid] = {Listable};
euclid[n_] := euclid[n] = euclid[n - 1] Prime[n]
euclid[1] = 2;
PrimeQ[euclid[Range[20]] + 1 ]

Select[euclid[Range[20]] + 1, PrimeQ]

Select[Range[40], Function[x, PrimeQ[euclid[x]+1]]]

Select[Range[40], PrimeQ[euclid[#] + 1] &]