Code for Electronic Supplement

FrobeniusNumberNijenhuisDijkstra

FrobeniusNumberNijenhuisDijkstra[A_List] :=Module[{a,b,j,k,i,l,Y,X,w,v,ww},Clear[h,ptr];
      Clear[h, ptr, P];
      {a, b} = {First[A], Rest[A]};
      If[GCD @@ A ≠ 1, Return[∞]];
      NDup[loc_] :=(For[j=Quotient[i=loc,2],w<h[j],{i,j}={j,Quotient[j,2]},
            ptr[Mod[h[i]=h[j],a]]=i];h[i]=w;ptr[wwm]=i);
      h[l=1]=0;   ptr[_]=0;h[0]=0;
      While[l>0, (X = h[1];
          ptr[Mod[X,a]] = -1;
          Y = h[l--];
          If[l>0,For[j=1;i=2,If[i<l&&h[i+1]<h[i],i++];i≤l&&Y>h[i],
              j=i;i*=2,ptr[Mod[h[j]=h[i],a]]=j];h[j]=Y;ptr[Mod[Y,a]]=j]);
        ww = X+b;
        wwmod = Mod[ww,a];
        Do[w=ww[[p]]; wwm=wwmod[[p]]; v=ptr[wwm];
          Which[v==0, NDup[++l], v>0 && w<h[v], NDup[v]], {p, Length[b]}]];
      X-a];

FrobeniusNumberNijenhuisDijkstraDecreasing[A_] := Module[{a,b,j,i,l,Y,X,w,v,ww}, Clear[h,ptr,P];
      {a, b} = {First[A], Rest[A]};
   If[GCD @@ A != 1, Return[∞]];
      NDup[loc_] := (For[j=Quotient[i=loc,2],w<h[j],{i,j}={j,Quotient[j,2]},
            ptr[Mod[h[i]=h[j],a]]=i];h[i]=w;ptr[wwm]=i);
      h[l=1] = 0; ptr[_] = 0; h[0] = 0;
      P[0] = Length[b];
      While[l > 0,(X = h[1]; ptr[Mod[X,a]] = -1; Y = h[l--];
          If[l>0,For[j=1;i=2,If[i<l&&h[i+1]<h[i],i++];i≤l&&Y>h[i],
              j=i;i*=2,ptr[Mod[h[j]=h[i],a]]=j];h[j]=Y;
            ptr[Mod[Y,a]]=j]);
        ww = X+b; (*the new weights*)
        wwmod = Mod[ww,a]; (*the new vertices*)
        Do[w = ww[[p]]; wwm = wwmod[[p]]; v = ptr[wwm];
     Which[v==0, NDup[++l]; P[wwm] = p,
           v > 0 && w < h[v], NDup[v]; P[wwm]=p],
        {p, P[Mod[X,a]]}]];
      X-a];

ReachableGaussianPrimes

DiskBackground::usage = "DiskBackground is an option to ReachableGaussianPrimes that sets the color of the background disk.";

PrintTheCounts::usage = "PrintTheCounts is an option to ReachableGaussianPrimes that causes a count of the primes to be printed, but only if the UseLines option is set to False."; 
ReachableGaussianPrimes::usage = "ReachableGaussianPrimes[d, rad] displays the Gaussian primes within distance r from the origin that one can reach from 1 + I using steps of size at most d and stepping on the Gaussian primes only. The unreachable primes are shown in a contrasting style.";

ReachablesOnly::usage="ReachablesOnly is an option to ReachableGaussianPrimes that causes the unreachable primes to not be shown. If UseLines is set to False, then ReachableOnly automatically becomes True.";

ReachableStyle::usage = "ReachableStyle is an option to ReachableGaussianPrimes that sets the style of the reachable points."; 

LineStyle::usage = "LineStyle is an option to ReachableGaussianPrimes that sets the style of the lines in the network."; 

UnreachableStyle::usage="UnreachableStyle is an option to ReachableGaussianPrimes that sets the style of the unreachable points.";
UseLines::usage = "UseLines is an option to ReachableGaussianPrimes that shows the edges connecting the reachable primes.";

GaussianIntegerQ[_Integer|Complex[_Integer,_Integer]]:=True;
GaussianIntegerQ[_]:=False;
GaussianPrimeQ[p_?GaussianIntegerQ]:=PrimeQ[p,GaussianIntegers->True];
GaussianPrimesInExpandedFirstOctant[rad_,k_]:=Select[Flatten[Table[{a,b},{a,-Ceiling[k],rad},{b,-Ceiling[k],Min[a+Floor[k (Sqrt[2])],(Sqrt[(rad^2)-(a^2)])]}],1],GaussianPrimeQ[#1[[1]]+#1[[2]] I]&];
vectors[k_]:=vectors[k]=Select[Flatten[Table[{i,j},{i,-Floor[k],Floor[k]},{j,-Floor[k],Floor[k]}],1],#1≠{0,0}&&((#1[[1]])^2)+((#1[[2]])^2)≤(k^2)&&Mod[#1[[1]]-#1[[2]],2]==0&];
neighbors[pt_,k_]:=Select[(pt+#1&)/@vectors[k],GaussianPrimeQ[#1[[1]]+#1[[2]] I]&];
neighbors[{1,1},k_]:=Select[Join[({2,1}+#1&)/@vectors[k],{{2,1},{1,2},{-1,1},{1,-1}}],(((#1[[1]]-1))^2)+(((#1[[2]]-1))^2)≤k&&GaussianPrimeQ[#1[[1]]+#1[[2]] I]&];

symmetrize[{a_Integer,b_Integer}]:={{a,b},{-a,b},{a,-b},{-a,-b},{b,a},{-b,a},{b,-a},{-b,-a}};
symmetrize[{p_,q_}]:=Transpose[{symmetrize[p],symmetrize[q]}];
symmetrize[pts_]:=Flatten[symmetrize/@pts,1]/;Depth[pts]==3;
symmetrize[{}]={};
graphedges[Gp_,d_]:=(Sort[{Gp,#1}]&)/@neighbors[Gp,d];
fix[pairs_]:=Scan[Which[mark[#1[[1]]],mark[#1[[2]]]=True,mark[#1[[2]]],mark[#1[[1]]]=True]&,pairs];

Options[ReachableGaussianPrimes]={ReachablesOnly->False,UseLines->True,ReachableStyle->{PointSize[0.01],RGBColor[0.5,0,0]},PrintTheCounts->False,UnreachableStyle->{PointSize[0.006],RGBColor[0,0,0.7]},DiskBackground->RGBColor[1,1,0.6], LineStyle->{Thickness[0.002], RGBColor[0.5,0,0]}};

(* it is important that mark be not Moduled *) 
ReachableGaussianPrimes[d_,rad_,opts___]:=Module[{unreachablePrimes,roadNetwork,unmarkedPairs,roads},{bg,reachsty,unreachsty,lineQ,roQ,prQ,lsty}={DiskBackground,ReachableStyle,UnreachableStyle,UseLines,ReachablesOnly,PrintTheCounts,LineStyle}/.{opts}/.Options[ReachableGaussianPrimes];If[!lineQ,roQ=True];
If[!ListQ[unreachsty],unreachsty={unreachsty}];
If[!ListQ[reachsty],reachsty={reachsty}];
If[!ListQ[lsty],lsty={lsty}];

  roads = Union @@((graphedges[#, d]&) /@  GaussianPrimesInExpandedFirstOctant[rad, d]);
    
  Clear[mark];   mark[_] := False;   mark[{1, 1}] = True; 
  unmarkedPairs = roads; 
  FixedPoint[(fix[unmarkedPairs]; 
    Length[unmarkedPairs =  Select[unmarkedPairs,!mark[#[[1]]] || !mark[#[[2]]]&]])&,Length[unmarkedPairs]]; 
 firstOctant = Select[roads, mark[#[[1]]]&]; 
 If[!lineQ, firstOctant = Union[Cases[firstOctant,
     {x_Integer, y_} /; 0 ≤ y ≤ x, Infinity]]]; 
 roadNetwork = Union @@(symmetrize /@ firstOctant); 
 unreachablePrimes =   If[roQ, {}, Union[symmetrize[Union[Select[Flatten[roads, 1], 
     (!mark[#]) && Norm[#]≤ rad & ]]]]]; 
If[!lineQ && prQ, Print[StringForm[ "There are `` primes in the first octant and `` in the entire component.",Length[firstOctant], Length[roadNetwork]]]]; 
 Graphics[{{EdgeForm[Thickness[0.002]],bg, Disk[{0, 0}, rad]}, 
  Append[unreachsty, Point @ unreachablePrimes], 

  If[lineQ,Append[lsty,  Line[roadNetwork]],Append[reachsty,  Point/@ roadNetwork]]}, 
     Sequence @@ FilterRules[{opts}, First /@ Options[Graphics]], 
      PlotRange -> All]]

Benford Law Demo

GetData[property_]:=DeleteCases[N[Table[CountryData[country,property],{country,CountryData[]}]],x_/;x==0||!NumericQ[x]];
GetData["ExchangeRate"]=Union[Table[CountryData[c, "ExchangeRate"], {c, CountryData[]}]];
Benford[d_]:=Log10[((d+1)/d)];
BenfordVector=N[Benford/@Range[9]];
BenfordΔ[proportions_]:=((Norm[proportions-BenfordVector])/(Norm[{0,0,0,0,0,0,0,0,1}-BenfordVector]));
BenfordBars=Graphics[
{LightPurple,EdgeForm[{GrayLevel[0.2], Red,Thickness[0.002]}],Table[Rectangle[{i-0.3, 0},{i+0.3,BenfordVector[[i]]}],{i,9}]}];
FirstDigit[x_]:=RealDigits[N[x]][[1,1]];
DigitFrequencies[d_] := With[{dd=Sort[Tally[FirstDigit/@d]]},
Last /@ Union[dd,Table[{mm,0},{mm,Complement[Range[9],First/@dd]}]]];
DigitProportions[d__] :=With[{ddd= DigitFrequencies[d]},N[ ddd/ Total[ ddd]]];
BenfordDeviance[data_] := BenfordΔ[DigitProportions[data]];
NumericProperties=Select[CountryData["Properties"], NumericQ[CountryData["USA", #]] &];
units[property_] := CountryData["USA",property,"Units"];

spread[data_]:=With[{d=Log10[Sort[Abs[data]]], n=Length[data]},
middle=d[[Round[0.75 n]]]-d[[Round[0.25 n]]];
d1=DeleteCases[d,h_/;h<d[[Round[0.25 n]]]-1.5 middle||h>d[[Round[0.75 n]]]+1.5 middle];
d1[[-1]]-d1[[1]]];

SpreadAndDeviation1=Table[With[{data= GetData[p]},
{spread[data],BenfordDeviance[data], p}],{p,NumericProperties}];

Options[BenfordDemo] = {ImageSizes->{500,200,245,200}, FontSizeLarge->12, FontSizeSmall->10};

BenfordDemo[opts___]:=({sz,fsl,fss}={ImageSizes, FontSizeLarge, FontSizeSmall}/. {opts}/. Options[BenfordDemo];
Manipulate[
currentdata=Abs[GetData[property1]];
Pane[Column[{Show[Graphics[{PointSize[0.02],Red,Point[Most[Select[SpreadAndDeviation1,#1[[3]]==property1&,1][[1]]]]}],imageScatter,Frame->True,Axes->False,FrameLabel->{"Logarithmic spread","Deviation from Benford's law"},PlotRange->{{0,8},{0,0.7}},FrameTicks->{Range[0,9],Range[0,1,.1],{},{}},LabelStyle->fsl,FrameTicksStyle->fss,AspectRatio->0.3,ImagePadding->{{40,10}, {40,5}},ImageSize->sz[[1]],
PlotLabel->Style["Place the mouse over a point to see the three graphs corresponding to that point.",fss]],
GraphicsRow[DeleteCases[{
If[trueQ,mm1=Min[currentdata];mm2=Max[currentdata];
Histogram[currentdata,Automatic, "Probability",ChartElementFunction->"FadingRectangle",ChartStyle->Blue,PlotRange->All,Axes->False,Frame->True,PlotLabel->"Raw data",
 AspectRatio-> 0.7,ImageSize->sz[[2]],
FrameTicks->
{{Automatic,None},{ N@FindDivisions[{mm1, mm2},3]/. (z_/; z==Round[z]&&z<10001):> Round[z], None}}, FrameLabel->None, LabelStyle->fsl, FrameTicksStyle->fss],{}],

If[BenfordQ,Show[ BenfordBars,
Graphics[{Black,Red,PointSize[0.05],
EdgeForm[{Black, Thickness[0.0003]}], FaceForm[],Rectangle[#-{.4,.005},#+{.4,.005}]& /@ Transpose[{Range[9],DigitProportions[currentdata]}]}, PlotRange->All],
 PlotRange->{{0.5,9.5}, {0, .4}},AspectRatio->.7,
FrameTicks-> {{{}, {{Log10@2,"(Subscript[log, 10])(2)"},{Log10@1.5,"(Subscript[log, 10])((3/2))"}}},{Range[9],{}}},Axes->False, Frame->True,
PlotLabel-> Style[If[!MemberQ[{None, "People"},units[property1]], StringForm["``\n(in ``)", property1,units[property1]], property1],fsl],LabelStyle->fsl, FrameTicksStyle->fss,ImageSize->sz[[3]], AspectRatio-> 0.5, FrameLabel->{"Digits","Probability",None,None}],{}],

If[logQ,Show[Histogram[Log10[Abs[currentdata]],Automatic, "ProbabilityDensity",
ChartElementFunction->"FadingRectangle",ChartStyle->Orange], PlotRange->{{-5,14},All}, Axes->False, Frame->True,
FrameTicks->
{{#,(10^(ToString[#]))}&/@ Range[-2,14,5],Automatic,{},{}},
  PlotLabel-> "Distribution of (Subscript[log, 10])(data)",LabelStyle->fsl, FrameTicksStyle->fss,ImageSize->sz[[4]], AspectRatio-> 0.7,
 FrameLabel->{units[property1];None, None }], {}]}, {}], Spacings->-30 ]}, Alignment->Center], ImageSize-> {1.05 Total[Rest[sz]],0.6 Total[Rest[sz]]},ImageSize->{520,340}],

{{property1,"Area", "property"},NumericProperties},  
{{trueQ,True, "true distribution"}, {True, False}, None},
{{BenfordQ,True, "Benford comparison"}, {True, False}, None},
{{logQ,True, "(Subscript[log, 10]) distribution"}, {True, False}, None},
Evaluate[Sequence@@ FilterRules[{opts},Options[Manipulate]]],
ControlPlacement->Top,SaveDefinitions->True,TrackedSymbols:>{ property1  }, Initialization->(property1="Area";imageScatter=ListPlot[Table[Tooltip[{d[[1]],d[[2]]},With[{property0=d[[3]]},Dynamic[property1=property0]]],{d,SpreadAndDeviation1}],Frame->True,Axes->False,PlotStyle->PointSize[0.01],FrameLabel->{"Logarithmic spread","Benford law deviation"},PlotRange->{{0,8},{0,0.7}},GridLines->{{},{}},ImageSize->400,FrameTicks->{Range[0,9],Range[0,1,.1],{},{}},LabelStyle->fsl,FrameTicksStyle->10,AspectRatio->0.3])])

Miscellany

Code for Chapter Opening Art

ReachableGaussianPrimes[(Sqrt[8]),100,ImageSize->400]

.  An Educational Integral

int1=(∫_0^1)(x^20) (E^x)\[DifferentialD]x;
N[int1]

J[n_]:=E-n J[n-1];
J[0]=E-1;
N[J[20]]

NIntegrate[(x^20) (E^x),{x,0,1}]

Expand[(∫_0^1)(x^20) (E^x)\[DifferentialD]x]

N[{2432902008176640000,895014631192902121 E},30]

N[(∫_0^1)(x^20) (E^x)\[DifferentialD]x,30]

J[20]

Clear[J];
J[n_]:=((E-J[n+1])/(n+1));

J[100]=1234.5678;
J[20]

.  Making the Alternating Harmonic Series Disappear

HarmonicRearrangement[T_, maxTerm_:60] := Module[{ans={}, sum=0., nextpos=1, nextneg=2},
While[Min[nextpos,nextneg]<maxTerm,
AppendTo[ans,0];While[sum≤T,sum+=(1/nextpos);ans[[-1]]++;nextpos+=2];
AppendTo[ans,0];While[sum≥T,sum-=(1/nextneg);ans[[-1]]++;nextneg+=2]];
ans];

HarmonicRearrangement[2]

TypesetSeries[code_] := ({odd,even}=Transpose[Partition[code,2]];PrependTo[odd,1];
PrependTo[even,1];DisplayForm[Row[Append[Rest[
Flatten[Map[FractionBox["1", ToString[#]]&,Riffle[Range[#[[1]], #[[2]]-2,2]& /@ Partition[2 Accumulate[ odd]-1,2,1],
Range[#[[1]], #[[2]]-2,2]& /@ Partition[2 Accumulate[ even],2,1]],{2}]/. FractionBox["1", s_] :>{If[OddQ[ToExpression[s]]," + "," - "],FractionBox["1", s] },2]]/. FractionBox["1", "1"]->"1", " + ···"]]])

e=TypesetSeries[HarmonicRearrangement[2,15]]

TypesetToRational[e_]:=ToExpression[
StringJoin[Drop[e[[1,1]],-3]/.FractionBox[a_,b_]:>Sequence[ToString[a],"/",ToString[b]]]]

TypesetToRational[e]
N[%]

rearrπ=HarmonicRearrangement[π,40]

asymptoticRatio[target_]:=(E^(2 target))/4;

asymptoticRatio[π]
N[%]

asymptoticRatio/@{0,Log[2],Log[3],1,π}

rearr0=HarmonicRearrangement[0]

TypesetSeries[HarmonicRearrangement[0,10]]==0

.  Bulletproof Prime Numbers

primeProd=Times@@Prime[Range[500]];
CertifiedRandomPrime[d_]:=Module[{p,r,q=CertifiedRandomPrime[Ceiling[(d/2)]+1]},
r=RandomInteger[{Ceiling[((10^(d-1))/q)],Floor[((10^d)/q)]}];While[p=r q+1;GCD[p,primeProd]≠1||PowerMod[2,p-1,p]≠1||GCD[p,PowerMod[2,r,p]-1]≠1,r++];
p];
CertifiedRandomPrime[d_/;d<8]:=RandomPrime[(10^({d-1,d}))]

CertifiedRandomPrime[100]//Timing

(p=RandomPrime[{(10^100),(10^101)}])//Timing

Needs["PrimalityProving`"]
ProvablePrimeQ[p]//Timing

.  Gaussian Moats

PrimeQ[{2,3}, GaussianIntegers->True]

FactorInteger[2,GaussianIntegers->True]

FactorInteger[1030507+2040608I]

Nonresidue[p_?OddQ]:=Select[Range[2,100],JacobiSymbol[#,p]==-1&,1][[1]];
Sum2Squares[p_/;Mod[p,4]==1] := NestWhile[{#[[2]],Mod@@#}&,
{p,PowerMod[Nonresidue[p],(p-1)/4,p]},((#[[1]])^2)>p&];
Sum2Squares[73]

p=0;While[Mod[p,4]≠1,p=RandomPrime[(10^50)]];p
ans=Sum2Squares[p]
((Norm[ans])^2)

n=40;
gPrimes=Select[Tuples[Range[-n,n],2],PrimeQ[#[[1]]+#[[2]] I,GaussianIntegers->True]&];
Graphics[{RGBColor[0.4,0,0],Rectangle/@(gPrimes/.z_?NumericQ:>z-(1/2))},Frame->True,FrameTicks->{Range[-n,n,10],({#,# I}&)/@Range[-n,n,10],None,None}]

GaussianPrimesInExpandedFirstOctant[rad_,k_]:=Select[Flatten[Table[{a,b},{a,-Ceiling[k],rad},{b,-Ceiling[k],Min[a+Floor[k (Sqrt[2])],(Sqrt[(rad^2)-(a^2)])]}],1],PrimeQ[#[[1]]+#[[2]] I,GaussianIntegers->True]&];

Graphics[{Line[{{15,0},{0,0},{11,11}}],PointSize[0.02],Point[GaussianPrimesInExpandedFirstOctant[15,2]]},Frame->True,PlotRange->{{-3,15.5},{-3,11.5}},FrameTicks->{Range[-2,14,2],({#1,#1 I}&)/@Range[-2,10,2],{},{}}]

neighbors[{1,1},2]={{2,1},{1,2},{-1,1},{1,-1}};
neighbors[pt_,2]:=Select[(pt+#&)/@{{1,1},{1,-1},{-1,1},{-1,-1},{2,0},{0,2},{0,-2},{-2,0}},PrimeQ[#[[1]]+#[[2]] I,GaussianIntegers->True]&]

primes=GaussianPrimesInExpandedFirstOctant[50,2];
edges[p_,d_]:=Table[Sort[{p,q}],{q,neighbors[p,d]}];
roads=Union@@Table[edges[p,2],{p,primes}];
Length/@{primes,roads}

Graphics[{Thickness[0.004],{RGBColor[0.5,0,0],Line[roads]},{Gray,Line[{{-3,-3},{40,40}}]}},Frame->True,FrameTicks->{Range[0,50,10],({#,# I}&)/@Range[0,50,10],None,None},GridLines->{{{0,{Gray}}},{{0,{Gray}}}}]

Needs["Combinatorica`"];

verts2=Union @@ roads;
edgesInGraph=List/@(roads/.{a_Integer,b_}:>Position[verts2,{a,b},1][[1,1]]);
ShowGraph[g=Graph[edgesInGraph,List/@verts2],VertexStyle->{RGBColor[0,0,0.7],PointSize[0.008]},EdgeStyle->{{RGBColor[0.5,0,0],Thickness[0.003]}}]

root=Position[verts2,{1,1}][[1,1]]

reachable=Select[ConnectedComponents[g],MemberQ[#,17]&,root][[1]]

ShowGraph[InduceSubgraph[g, reachable], VertexStyle->{RGBColor[0,0,0.7], PointSize[0.008]},  EdgeStyle->{{RGBColor[0.5,0,0],Thickness[0.003]}}]

symmetrize[{a_Integer,b_}]:={{a,b},{-a,b},{a,-b},{-a,-b},{b,a},{-b,a},{b,-a},{-b,-a}};
symmetrize[{p_List,q_List}]:=Transpose[{symmetrize[p],symmetrize[q]}];
symmetrize[pts_List]:=Flatten[symmetrize/@pts,1]/;Depth[pts]==3;

reachablePrimes=verts2[[reachable]];
firstOctant=Select[roads,MemberQ[reachablePrimes,#[[1]]]&];
roadNetwork=Union@@symmetrize/@firstOctant;
unreachablePrimes=Union[symmetrize[Complement[primes,reachablePrimes]]];

Graphics[{{EdgeForm[Thickness[0.002]],RGBColor[1,1,0.6],Disk[{0,0},50]},
{RGBColor[0,0,0.7],PointSize[0.005],Point[unreachablePrimes]},
Thickness[0.002],RGBColor[0.5,0,0],Line[roadNetwork]}]

ReachableGaussianPrimes[(Sqrt[8]),100]

.  Frobenius Number by Graphs

cols={Red,Blue};
A={6,9,20};a=A[[1]];B=Rest[A];
edgeRules=Flatten[Table[{i->Mod[i+j,a,1],j},{j,B},{i,a}],1];
pts=Table[{Cos[2 π t/a],Sin[2 π t/a]},{t,0,a-1}];
GraphPlot[edgeRules,DirectedEdges->True,VertexCoordinateRules->Table[i->pts[[i]],{i,a}],VertexRenderingFunction->({Circle[#1,0.07],Text[#2/.a->0,#1]}&),EdgeRenderingFunction->({cols[[Position[B,#3,1][[1,1]]]],Thickness[0.006],Arrow[#1,0.15]}&),BaseStyle->Bold,ImageSize->300]

<<Combinatorica`

FrobeniusGraph[A_]:=Module[{B=Rest[A],a=A[[1]]},verts=Table[{{Cos[j 2 π/a],Sin[j 2 π/a]}},{j,0,a-1}];edges=Flatten[Table[{{i,Mod[i-1+j,a]+1},EdgeWeight->j,EdgeColor->If[j==A[[2]],Red,Blue]},{j,B},{i,a}],1];Graph[edges,verts,EdgeDirection->True]];
g=FrobeniusGraph[A];
ShowGraph[g,VertexStyle->White,VertexLabel->RotateLeft[Range[0,a-1]],VertexLabelPosition->Center,EdgeColor->Red]/.Arrow[e_]:>Arrow[e,0.04]

{parents, wts}=Dijkstra[g,6]

tree=DeleteCases[Transpose[{parents,Range[6]}],{i_,i_}]

edgeRules1=Flatten[Table[{i->Mod[i+j,a,1],j},{j,B},{i,a}],1];
GraphPlot[edgeRules1,DirectedEdges->True,VertexRenderingFunction->({White,EdgeForm[Black],Disk[#1,0.07],Black,Text[#2/. 6->0,#1]}&),VertexCoordinateRules->Table[i->pts[[i]],{i,a}],EdgeRenderingFunction->({If[MemberQ[tree,#2],{cols[[Position[B,#3,1][[1,1]]]],Thickness[0.008],Arrow[#1,0.2]},{Thickness[0.0015],GrayLevel[0.2],Arrowheads[0.025],Arrow[#1,0.2]}]}&),BaseStyle->Bold]

pts1={{2,2}, {0.2,1}, {2,0}, {0,2}, {1.8,1}, {0,0}};
GraphPlot[edgeRules1,DirectedEdges->True,VertexRenderingFunction-> (Text[#2/. 6->0,#1]&),VertexCoordinateRules->Table[i->pts1[[i]],{i,a}],EdgeRenderingFunction->({If[MemberQ[tree,#2],{cols[[Position[B,#3,1][[1,1]]]],Thickness[0.008],Arrow[#1,0.2]},{Thickness[0.005],GrayLevel[0.2],Arrowheads[0.04],Arrow[#1,0.2]}]}&)]

Aprime=Prime[{100,1000,10000, 100000}]
FrobeniusNumberNijenhuisDijkstra[Aprime]//Timing

FrobeniusNumber[Aprime]//Timing

FrobeniusNumberBF[A_]:=Module[{L,Q,h,t,u,w,a=A[[1]]},
Q=Array[0&,a];
L=Append[Table[∞,{a-1}],0]; (* initialize queue and labels *)
h=t=a;                                                      (* place a on the queue *) 
While[t<∞,If[t==h,t=∞];
{v,Q[[h]],h}={h,0,Q[[h]]};            (* pull vertex off head of queue *) 
Do[{u,w}={Mod[b+v,a],b+L[[v]]}; (* form new discovered vertex and weight *)
    If[w<L[[u]],L[[u]]=w;                   (* relax *) 
       If[Q[[u]]==0,If[t<∞,t=Q[[t]]=u,t=h=u]] ],    (* add u to queue *) 
{b,Rest[A]}]];
Max[L]-a];

FrobeniusNumberBF[{6,9,20}]

FrobeniusNumberBF[Aprime]//Timing

a=(10^6);
A=Prepend[RandomInteger[{a+1,10 a},14],a];
FrobeniusNumber[A]//Timing

a=(10^20);
A=Prepend[RandomInteger[{a+1,10 a},6],a];
FrobeniusNumber[A]//Timing

soln=FrobeniusSolve[A, (10^25),1]

soln.A

.  Benford's Law of First Digits

GetData[property_]:=DeleteCases[N[Table[CountryData[country,property],{country,CountryData[]}]],x_/;x==0||!NumericQ[x]];

Benford[d_]:=Log10[((d+1)/d)];
BenfordVector=N[Benford/@Range[9]]

FirstDigit[x_]:=RealDigits[N[x]][[1,1]];
DigitFrequencies[d_] := (dd=Sort[Tally[FirstDigit/@d]];
Last /@ Union[dd,Table[{mm,0},{mm,Complement[Range[9],First/@dd]}]]);
DigitProportions[d__] :=With[{ddd= DigitFrequencies[d]},N[ ddd/ Total[ ddd]]];

DigitProportions[GetData["Area"]]

BenfordBars=Graphics[
{LightPurple,EdgeForm[{GrayLevel[0.2], Red,Thickness[0.002]}],Table[Rectangle[{i-0.3, 0},{i+0.3,BenfordVector[[i]]}],{i,9}]}];
Show[BenfordBars,Graphics[{
EdgeForm[{Thickness[0.0002]}],FaceForm[],Rectangle[#-{0.4,.0035},#+{0.4,.0035}]& /@ Transpose[{Range[9],DigitProportions[GetData["Area"]]}]}], Frame->True,
FrameTicks->{{Table[{Log10[(i+1)/i],StringForm[If[i==1,"(Subscript[log, 10])(``)","(Subscript[log, 10])(``/``)"], i+1,i]},{i, Range[4]}], Range[0.05,0.3,0.05]},{Range[9],{}}},AspectRatio-> 0.6,PlotLabel->"Country areas and Benford's Law"]

BenfordΔ[proportions_]:=((Norm[proportions-BenfordVector])/(Norm[{0,0,0,0,0,0,0,0,1}-BenfordVector]));

BenfordΔ[DigitProportions[GetData["Area"]]]

cityPops=DeleteCases[Table[CityData[c, "Population"], {c, CityData[]}], 0];

Length[cityPops]

cityProportions=DigitProportions[cityPops]

Show[BenfordBars,Graphics[{
EdgeForm[{ Thickness[0.0003]}], FaceForm[],Rectangle[#-{.4,.0035},#+{.4,.0035}]& /@ Transpose[{Range[9],cityProportions}]}], Frame->True, Axes->False, 
FrameTicks->{{Table[{Log10[(i+1)/i],StringForm[If[i==1,"(Subscript[log, 10])(``)","(Subscript[log, 10])(``/``)"], i+1,i]},{i, Range[4]}], Range[0.05,0.3,0.05]},{Range[9],{}}},AspectRatio-> 0.6,PlotLabel->"163 369 city populations and Benford's Law"]

BenfordΔ[cityProportions]

g[x_]:=(E^(-(((x-3.5))^4.)))/(2 Gamma[5/4])
f[x_]:= (PDF[NormalDistribution[5,1.4]][x] +  g[x])/2.
Show[Plot[f[x],{x,0,9},PlotStyle->Thickness[0.008]],Table[Plot[f[x],{x,i,i+Log10[2]},PlotStyle->Thick,Filling->0,FillingStyle->Red],{i,0,9}],AspectRatio->0.4,Frame->True,Axes->False,PlotRangePadding->0.02,GridLines->{{},{0}},FrameTicks->{{Automatic,{}},Reverse[{Range[-2,9,1],({#1,(10^(ToString[#1]))/.(10^("0"))->1/.(10^("1"))->10}&)/@Range[9]}]}]

Plot[((f[Log10[x]])/(x Log[10])),{x,100,10000},PlotRange->All,PlotStyle->Thick,Frame->True,Axes->False]

((NIntegrate[f[x],{x,3,3+Log10[2]}])/(NIntegrate[f[x],{x,3,4}]))

Table[((NIntegrate[f[x],{x,i,i+Log10[2]}])/(NIntegrate[f[x],{x,i,i+1}])),{i,0,9}]

(∑_(i=0)^9)NIntegrate[f[x],{x,i,i+Log10[2]}]

g[x_]:= PDF[NormalDistribution[5,14/10]][x];

(∫_(-∞)^(-35))g[x]\[DifferentialD]x+(∫_45^∞)g[x]\[DifferentialD]x
N[%]

oneTerm=(∫_i^(i+δ))g[x]\[DifferentialD]x

deviation=((∑_(i=-35)^45)oneTerm)-δ;

N[N[deviation/.δ->Log10[2],200]]

$MaxExtraPrecision=200;
errors=N[Table[{δ,N[(∑_(i=-35)^45)oneTerm-δ,20]},{δ,(1/40),1-(1/40),(1/40)}]];

ListLinePlot[errors, Frame->True, Axes->False, PlotStyle->Thick, GridLines->{{Log10[2]},{0}}, FrameTicks->{{Automatic, {}},{{0,{ Log10[2],"(Subscript[log, 10])2"},.5, .8,1},{}}}, PlotRange->{{0,1}, All}]

f1[x_]:= PDF[NormalDistribution[5,0.3]][x];
Show[Plot[f1[x],{x,3,6},PlotStyle->Thickness[0.008], PlotRange->All],Table[Plot[f1[x],{x,i,i+Log10[2]},PlotStyle->Thick,Filling->0,FillingStyle->Red],{i,3,6}],
AspectRatio->0.4,Frame->True,Axes->False,PlotRangePadding->0.08,GridLines->{{},{0}},FrameTicks->{{Automatic,{}},Reverse[{Range[-2,9,1],({#1,(10^(ToString[#1]))/.(10^("0"))->1/.(10^("1"))->10}&)/@Range[9]}]}]

(∑_(i=3)^8)NIntegrate[f1[x],{x,i,i+Log10[2]}]

cityPopsLog=N[Log10[cityPops]];
Max[cityPopsLog]

histData=Transpose[{Range[0,7,.5],((2. Last/@Sort[Tally[cityPopsLog,Floor[2 #1]==Floor[2 #2]&]])/(Length[cityPops]))}]

ListLinePlot[histData, Frame->True, Axes->False, PlotStyle->Thick]

cityPopsLogPL = Interpolation[histData, InterpolationOrder->1];

(∫_0^7.)cityPopsLogPL[x]\[DifferentialD]x

Show[
Table[Plot[cityPopsLogPL[x],{x,i,i+Log10[2]},PlotStyle->Thick, Filling->0, FillingStyle->Red],{i,0,7}],
Plot[cityPopsLogPL[x],{x, 0, 7}, PlotStyle->{Blue,Thickness[0.008]}],
 Frame->True,FrameLabel->{"163369 city populations","Probability"},
FrameTicks->{Join[({#1,(10^(ToString[#1]))/.(10^("0"))->1/.(10^("1"))->10}&)/@Range[0,7],{}],Range[0, 0.4, 0.1],{{2+Log10[2],200},{3+Log10[2],2000},{4+Log10[2],2·(10^("4"))},{5+Log10[2],2·(10^("5"))}},{}}, GridLines->{{2,2.301, 3,4,5,3.301,4.301, 5.301},{0}}, PlotRange->All, PlotRangePadding -> 0.05,Axes->None,AspectRatio-> 0.5, LabelStyle->12,FrameTicksStyle->10]

Table[(((∫_i^(i+Log10[2]))cityPopsLogPL[x]\[DifferentialD]x)/((∫_i^(i+1))cityPopsLogPL[x]\[DifferentialD]x)),{i,0,6}]

(∑_(i=0)^6)(∫_i^(i+Log10[2]))cityPopsLogPL[x]\[DifferentialD]x

GetData["ExchangeRate"]=Union[Table[CountryData[c, "ExchangeRate"], {c, CountryData[]}]];

CountryData["India", "ElectricalGridSocketImages"]

NumericProperties=Select[CountryData["Properties"], NumericQ[CountryData["USA", #]] &]

spread[data_]:=(d=Log10[Sort[Abs[data]]];
n=Length[d];
middle=d[[Round[0.75 n]]]-d[[Round[0.25 n]]];
d1=DeleteCases[d,h_/;h<d[[Round[0.25 n]]]-1.5 middle||h>d[[Round[0.75 n]]]+1.5 middle];
d1[[-1]]-d1[[1]]);

SpreadAndDeviation=Table[With[{data= GetData[p]},
{spread[data],BenfordΔ[DigitProportions[data]], p}],{p,NumericProperties}];

ListPlot[Table[Tooltip[{d[[1]],d[[2]]},d[[3]]],{d,SpreadAndDeviation}],Frame->True,Axes->False,FrameLabel->{"Logarithmic spread","Deviation from Benford Law"},PlotRange->All]

BenfordDemo[]

properties = {"IndependenceYear","MedianAge","InternetHosts","Area", "ManufacturingValueAdded"};
propertyLabels={"Independence year","Median age","Internet hosts","Area", "Manufacturing value added"};
yticks={Range[0,20,5],Range[0,10,1],Range[0,1,0.1],Range[0,1,0.1],Range[0,1,0.1]};

GraphicsGrid[Table[propData=GetData[p];pLab=propertyLabels[[Position[properties, p][[1,1]]]];
ytic=yticks[[Position[properties, p][[1,1]]]];
{Show[Graphics[{Thickness[0.012],RGBColor[1,0.2,0.2],Line@Table[{i,BenfordVector[[i]]}, {i,9}]},PlotRange->{{0.8,9.2}, {-.02, .45}}],

ListPlot[DigitProportions[propData], PlotRange->All,PlotStyle->{Black,PointSize[0.04]}],
 PlotRange->{{0.8,9.2}, {-.02, .45}},FrameTicks-> {{{.3,.1},{}},{Range[9],{}}},Axes->False, Frame->True, PlotLabel-> Style[StringForm["``: Δ = `` ",pLab,
NumberForm[BenfordΔ[DigitProportions[propData]],2]],11], AspectRatio->0.62],

Show[Histogram[Log10[Abs[propData]],Automatic, "ProbabilityDensity",
ChartElementFunction->"FadingRectangle",ChartStyle->Orange], PlotRange->{{-5,14},All}, 
Axes->False, Frame->True,FrameTicks->
{{{},ytic},{{#,(10^(ToString[#]))/. 10^"0" -> "1"}&/@ {-3,0,3,6,9,12},{}}}, PlotLabel-> Style[StringForm["``: spread = `` ",pLab,
NumberForm[spread[propData],2]],11]]},{p, properties}],ImageSize->500]