Code for Electronic Supplement

PrimePi Overload

Off[General::spell1];
Unprotect[PrimePi];
PrimePi[10.^23|10^23] = 1925320391606803968923;
PrimePi[10.^22|10^22] = 201467286689315906290;
PrimePi[10.^21|10^21] = 21127269486018731928;
PrimePi[10.^20|10^20] = 2220819602560918840;
PrimePi[10.^19|10^19] = 234057667276344607;
PrimePi[10.^18|10^18] = 24739954287740860;
PrimePi[10.^17|10^17] = 2623557157654233;
PrimePi[10.^16|10^16] = 279238341033925;
PrimePi[10.^15|10^15] = 29844570422669;
PrimePi[10.^14|10^14] = 3204941750802;
PrimePi[10.^13|10^13] = 346065536839;
PrimePi[10.^12|10^12] = 37607912018;
PrimePi[10.^11|10^11] = 4118054813;
PrimePi[10.^10|10^10] = 455052511;
Protect[PrimePi];

FactorForm

FactorForm::usage = "FactorForm[f] puts the list f, which should be a list of the form returned by FactorInteger, into traditional mathematical notation.";

FactorForm[l_] :=  SequenceForm@@
 (  (DisplayForm[RowBox[{#1, " "}]] & ) /@ 
 (MapThread[
  If[#2 == 1,
  If[Head[#1]===Complex,
    StringJoin["(", ToString[#1], ")"], 
    ToString@#1], SuperscriptBox[ToString@#1, #2]]&, 
   Transpose[l]]))

Prime Numbers

Code for Chapter Opening Art

.  Basic Number Theory Functions

Prime[1000]

(10^(Range[5]))

Prime[(10^(Range[9]))]

PrimePi[(10^6)]

Prime[78499]

FactorInteger[11737654214175]

Transpose[%]

((%[[1]])^(%[[2]]))

Apply[Times,%]

FactorInteger[100!]//FactorForm

Divisors[123456789]

GCD[76895,16302982080]

LCM[165,150]

Mod[1238719479147974,100]

Divisible[100!+1,101]

WilsonQ[n_] := Divisible[(n-1)!+1,n]
w=Select[Range[2, 200], WilsonQ]

Table[Mod[w[[i]]!,i], {i, Length[w]}]

Select[Range[2,200],Divisible[(#-1)!+1,#]&]

Table[Mod[(n-1)!+1,n],{n,2,200}]

PowerMod[9,1000000,7]

PowerMod[9999999999,(10^100),700]

PowerMod[57,-1,100]

57 93

PowerMod[444,(1/2),1000]

x/.Solve[{(x^2)==444,Modulus==1000}]

Union[PowerMod[Range[1,540],540,541]]

PseudoprimeQ[n_]:=!PrimeQ[n]&&PowerMod[2,n-1,n]==1
Select[Range[10000],PseudoprimeQ]

PowerMod[3,EulerPhi[1016],1016]

EulerPhi[267]

good[n_]:=EulerPhi[n]==176
Select[Range[1000],good]

Reduce[n>0&&EulerPhi[n]==400000,n,Integers]

Short[n/.{ToRules[Reduce[EulerPhi[n]==400000&&n>0,n,Integers]]}]

Union[Table[Length[Reduce[n>0&&EulerPhi[n]==k,n,Integers]],{k,500}]]

.  Where the Primes Are

Plot[PrimePi[x],{x,1,200},PlotStyle->{Thickness[0.004],Black},Exclusions->Range[200]]

Plot[PrimePi[x],{x,1,(10^6)},PlotStyle->{Thickness[0.004],Black}]

Plot[{PrimePi[x],(x/(Log[x]))},{x,2,50000},PlotStyle->{{Black},{Black,Thickness[0.006],Dashed}}]

RiemannRData=(1./(Range[200]! Range[200] Zeta[Range[200]+1]));
RiemannR[x_?NumberQ]:=1+((Log[x])^(Range[200])) . RiemannRData;

x=(10^23);Grid[{
	{"Legendre",Round[(x/(Log[x]-1.08366))]},
	{"Chebyshev", Round[(x/(Log[x]-1))]},
          {"Gauss",Round[LogIntegral[x]]},
	{"Riemann",Round[RiemannR[x]]},
	{"π((10^23))", PrimePi[x]}},Dividers->All]

Plot[{LogIntegral[x]-PrimePi[x],RiemannR[x]-PrimePi[x],(x/(Log[x]-1))-PrimePi[x]},{x,2,(10^6)},MaxRecursion->3,Frame->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}}]

ListLinePlot[
Table[{i,Log[10,Abs[#[(10^i)]-PrimePi[(10^i)]]]},{i,2,23}]&/@{((#)/(Log[#]-1.08366))&,((#)/(Log[#]-1))&,LogIntegral,RiemannR},Frame->True,
PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue},{Thick,Cyan}},
FrameTicks->{{Table[{i,(("10")^i)},{i,0,18,3}],None},{Table[{i,(("10")^i)},{i,3,23,4}],None}}, Axes->None, PlotRange->{-1,18.1},Epilog ->{Text["Gauss", {18,8.2},{0,0}, {1, .4}], Text["Riemann", {18,5.5},{0,0}, {1, .4}],Text["Legendre", {18,14.5}, {0,0}, {1,.65}], Text["Chebyshev", {18.3,12.5},{0,0}, {1, .65}]}]

{PrimePi[(10^9)], Round[RiemannR[(10^9)]]}

ListLinePlot[{
Table[{i,Log[10,LogIntegral[(10^i)]-PrimePi[(10^i)]]},{i,2,23}],
Table[{i,Log[10, Log[(10^i)](Sqrt[(10^i)])]},{i,2,23}]},Frame->True,PlotStyle->{{Thickness[0.004],Black}},FrameTicks->{{Table[{i,(("10")^i)},{i,0,9,3}],None},{Table[{i,(("10")^i)},{i,3,23,4}],None}},Axes->None,PlotRange->{-1,12.5},Epilog->{Text["Error using LogIntegral",{17,6.1},{0,0},{1,.44}],Text["Riemann Hypothesis bound",{14,9.3},{0,0},{1,0.51}]}]

.  The Prime Number Race

Prime[Range[50]]

Mod[Prime[Range[2,50]],4,-1]

Accumulate[{a,b,c}]

Accumulate[Mod[Prime[Range[2,50]],4,-1]]

race=ListLinePlot[Accumulate[Mod[Prime[Range[2,3000]],4, -1]],Frame->True,Axes->False,GridLines->{{},{0}}]

Show[race,PlotRange->{{2920,2980},{-8,1}},
	FrameTicks->{{Automatic,None},{{2940,2946,2960}, None}}]

Prime[2947]

ListLinePlot[Accumulate[Mod[Prime[Range[2,3000]],4,-1]],Frame->True,Axes->False,GridLines->{{},{0}}, FrameTicks->{{{-30,-15,0,5},{}},{Append[Table[{n-1,Prime[n]},{n,500,2500,500}],{2945,Prime[2946]}],{}}}, AspectRatio->0.3, PlotRange->{-35, 5}]

PrimePiMod[x_?NumberQ,n_Integer,a_Integer]:=Count[Prime[Range[PrimePi[x]]],_?(Mod[#,n]==a&)];
{PrimePi[100], PrimePiMod[100,4,1],PrimePiMod[100,4,3]}

Plot[Evaluate[Table[PrimePiMod[x,11,i],{i,10}]],{x,1,200},PlotStyle->Thick]

AttachPosns[m_]:=Transpose[{m,Range[Length[m]]}]
PrimePiMod[x_,n_]:=AttachPosns/@Table[Select[Prime[Range[PrimePi[x]]],Mod[#,n]==i&],{i,n-1}]

PrimePiMod[100,5]

ListLinePlot[PrimePiMod[5000,11],PlotStyle->Thickness[0.004],Frame->True]

.  Euclid and Fibonacci

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

PrimeProduct[n_]:=(AppendTo[nCount,n];PrimeProduct[n-1] Prime[n])
PrimeProduct[1]=2;

nCount={};
Table[PrimeProduct[n]+1,{n,8}];
nCount

Clear[PrimeProduct]
PrimeProduct[n_]:=PrimeProduct[n]=PrimeProduct[n-1] Prime[n]
PrimeProduct[1]=2;

PrimeProduct[5]

?PrimeProduct

Fib[0|1]=1;
Fib[n_]:=Fib[n-1]+Fib[n-2]
Timing[Fib[26]]

Clear[Fib];
Fib[0|1]=1;
Fib[n_]:=Fib[n]=Fib[n-1]+Fib[n-2]
Timing[Fib[26]]

Attributes[PrimeProduct]=Listable;

PrimeQ[PrimeProduct[Range[20]]+1]

Select[Range[100],PrimeQ[PrimeProduct[#]+1]&]

.  Strong Pseudoprimes

OddPart[n_]:=With[{p=IntegerExponent[n,2]},{(n/(2^p)),p}];
Map[OddPart,{100,101,102}]

spspSequence[b_,n_]:=Module[{m,s},{m,s}=OddPart[n-1];NestList[Mod[((#)^2),n,-1]&,Mod[PowerMod[b,m,n],n,-1],s]];

spspSequence[2,7919]

spspSequence[2,Prime[10000]]

spspSequence[2,341]

StrongPseudoprimeQ[b_][n_]:=Module[{bSeq=spspSequence[b,n]},(bSeq[[1]]==1||MemberQ[bSeq,-1])&&!PrimeQ[n]]/;n>1&&OddQ[n]
StrongPseudoprimeQ[_][1]=False;

Select[Range[10000],StrongPseudoprimeQ[2]]

StrongPseudoprimeQ[3][1373653]

FactorInteger[1373653]

Table[StrongPseudoprimeQ[Prime[i]][3215031751],{i,1,10}]

FactorInteger[3215031751]

n=18215745452589259639 4337082250616490391 867416450123298079;
Select[Range[101],!StrongPseudoprimeQ[#][n]&]

Do[If[!StrongPseudoprimeQ[i][n],Print[i];Break[]],{i,200}]