(************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 32098, 960]*) (*NotebookOutlinePosition[ 32933, 986]*) (* CellTagsIndexPosition[ 32889, 982]*) (*WindowFrame->Normal*) Notebook[{ Cell["\<\ This file contains the code that allows one to compute \ Hardy-Littlewood estimates for admissible sets, as illustrated in our paper, \ The Gaussian Zoo (Experimental Mathematics). More info or comments: Stan Wagon (wagon@macalester.edu)\ \>", "Text"], Cell[CellGroupData[{ Cell["Hardy-Littlewood Code", "Section", PageWidth->PaperWidth], Cell["\<\ diamond = {-1, 1, I, -I}; castle = {0, 2, 4, 1 + I, 3 + I, 2*I, 4 + 2*I, 1 + 3*I, 3 + 3*I, 4*I, 2 + \ 4*I, 4 + 4*I}; fifteenexample = {11,13,17,19,23,29,31,37,41,43,47,53,59,61,67}; \ \>", "Input", PageWidth->PaperWidth], Cell[BoxData[{ \(\(\(SetAttributes[GaussianMod, \ Listable];\)\(\n\) \)\), "\n", \(\ \(\(GaussianMod[u_?GaussianIntegerQ, \ z_?GaussianIntegerQ]\ := \ \n\ \ Module[{d, q, \ a\ = \ Re[z], \ b = Im[z], \ x = Re[u], \ y = Im[u]}, \n\ \ {d, \ e}\ = \ ExtendedGCD[a, b]; \n\ \ \ q\ = \ Quotient[y, d]; \n\ \ \ Mod[ x\ - \ a\ q\ e[\([2]\)]\ + \ b\ q\ e[\([1]\)], \n\t\ \((a^2 + b^2)\)/d]\ \ + I\ \((y\ - \ q\ d)\)]\)\(\n\) \)\)}], "Input", PageWidth->PaperWidth], Cell["\<\ ToComplex[e_] := e /. {x_, y_} :> x+I y ToCoords[e_] := e /. x_?NumericQ :> {Re[x], Im[x]} Format[ExponentForm[n_, b_:10]] := Module[{e = N[Log[b, Abs[n]], 20], h, ff, nn = Abs[n]}, h = Round[e]; e = If[h < e + 0.4, h, h - 1]; Sign[n]*If[(ff = IntegerExponent[nn, b]) != 0 && nn/b^ff < Sqrt[nn], nn/b^ff*ToString[b]^ff, SFTF[ToString[b]^e + (nn - b^e)]]] ExponentForm /: N[HoldPattern[ExponentForm[n_, b_:10]]] := N[n] If[$VersionNumber < 4, IntegerExponent[n_, b_Integer] := Module[{m = n, r = 0}, While[Mod[m, b] == 0, m /= b; r++]; r] /; b > 1 && n > 0; IntegerExponent[0, b_Integer] := Infinity]\ \>", "Input",\ PageWidth->PaperWidth], Cell["\<\ FactorFormfor[{}] = 1; FactorFormfor[l_] := SequenceForm @@ (DisplayForm[RowBox[{#, \" \"}]] &) /@ MapThread[If[#2 == 1, If[Head[#1] === Complex, StringJoin[\"(\", ToString[#1], \")\"], ToString[#1]], SuperscriptBox[ToString[#1], #2]] & , Transpose[l]] /. {\"(I)\" -> \"I\", \"(-I)\" -> \"-I\"}; FactorFormfor[x_?Negative] := SequenceForm[\"-\", FactorForm[FactorInteger[-x]]]; HoldPattern[FactorForm[x_]] := FactorFormfor[x]; FactorForm /: N[HoldPattern[FactorForm[l_]], prec_:16] := N[Unfactor[Round[l]], prec]; Unfactor[n_] := Times @@ Apply[Power, Round[n], {1}]\ \>", "Input", PageWidth->PaperWidth], Cell["\<\ HardyLittlewood::usage = \"HardyLittlewood[A, prec:6] returns the \ Hardy-Littlewood constant for the set A of integers or Gaussian integers.\"; ConstellationEstimate::usage = \"ConstellationEstimate[A, prec:6] returns the \ function, of R, say, that gives the estimate for the number of times the \ pattern A occurs in the primes up to R (or the Gaussian primes in the first \ octant of radius R). R can be either numeric or symbolic, and in the latter \ case the result is the symbolic function of R.\"; ZetaGaussian::usage = \"ZetaGaussian[s] gives the version of the Riemann zeta \ function appropriate for the Gaussian integers: the sum of 1/N(z)^s where z \ ranges over Z[i], ignoring associates.\"; PrimeZeta::usage = \"PrimeZeta[s, prec:6] gives the prime zeta function, the \ sum of 1/p^s over prime p. The GaussianIntegers option causes this to compute \ the Gaussian prime zeta function.\"; EulerProduct::usage = \"EulerProduct[k,n] gives the product of \ (p-k)p^(k-1)/(p-1)^k over primes p \[GreaterEqual] n in typeset form. n \ should be larger than k\"; ConstellationEstimateAsymptotic::usage = \"ConstellationEstimateAsymptotic[A] \ is a function that gives the constellation estimate for A in its asymptotic \ form.\";\ \>", "Input", PageWidth->PaperWidth], Cell["\<\ Unprotect[NextPrime]; NextPrime[-3] := -2 NextPrime[-2] := 2 NextPrime[-1] := 2 NextPrime[0] := 2 NextPrime[1] := 2 NextPrime[r_Real] := NextPrime[Floor[r]] NextPrime[n_Integer] := Block[{i = n + 1 + Mod[n, 2]}, If[Abs[n] < 10^20, While[Not[PrimeQ[i]], i += 2], prod = Apply[Times, Prime[Range[Log[ Abs[n] ]]]]; While[True, If[GCD[ i, prod] == 1, If[PrimeQ[i], Break[]]]; i += 2] ]; i ] SFTF[x_] := StandardForm[TraditionalForm[x]] GaussianNorm[a_?GaussianIntegerQ] := Abs[a]^2 GaussianIntegerQ[a_Integer] := True GaussianIntegerQ[Complex[a_Integer, b_Integer]] := True GaussianIntegerQ[a_] := False GaussianPrimeQ[p_?GaussianIntegerQ] := PrimeQ[p, GaussianIntegers -> True] Options[AdmissibleQ] = {GaussianIntegers->False}; AdmissibleQ[c_, opts___] := Module[ {primes, gi, l=Length[c]}, gi = GaussianIntegers /. {opts} /. Options[AdmissibleQ]; If[!FreeQ[c, Complex], gi = True]; If[gi, primes = Select[Flatten[Table[a + b I, {a, 0, Floor[Sqrt[l]]}, {b, 0, Floor[Sqrt[l - a^2]]}]], GaussianPrimeQ]; And @@ (Length[Union[Mod[c, #]]] != GaussianNorm[#] &) /@ primes, And @@ Table[p = Prime[i]; Length[Union[Mod[c, p]]] != p, {i, PrimePi[l]}]]] \ \>", "Input", PageWidth->PaperWidth], Cell[BoxData[{ \(\($MaxExtraPrecision = 300;\)\ \ \ \n\), "\n", \(\(ZetaGaussian[s_]\ := Zeta[s]\ \ 4\^\(-s\)\ \((Zeta[s, 1/4] - Zeta[s, 3/4])\);\)\n\t\t\), "\n", \(\(distancesq[a_, \ b_]\ := \ Abs[b - a]^2;\)\), "\n", \(\(\(maxnorm[set_]\ := Max[Apply[distancesq, DeleteCases[ Distribute[{set, set}, \ List], \ {x_, \ x_}], \ {1}]];\)\(\n\) \)\), "\n", \(\(\(NextNorm[L_]\ := \ Module[{t = L + 1}, \ While[\(! \((\n\t\t\t\t\t\((IntegerQ[tt = Sqrt[t]] && Mod[tt, 4] == 3\ && \ PrimeQ[tt])\) || \n\t\t\((\t Mod[t, 4] == 1\ && \ PrimeQ[t])\))\)\), \ \(t++\)]; \ t];\)\(\n\) \)\), "\n", \(\(\(Gprime[ n_]\ := \ \(FactorInteger[n, \ GaussianIntegers -> True]\)[\([\(-1\), 1]\)];\)\(\n\) \)\), "\n", \(\(\(Options[\ CriticalPrimeInterval]\ = \ {GaussianIntegers -> False};\)\(\n\) \)\), "\n", \(\(CriticalPrimeInterval[A_, \ opts___]\ := \ \((\t imagQ = \(GaussianIntegers\ /. \ {opts}\)\ /. \ Options[\ CriticalPrimeInterval]; \n\t\t\tIf[\(! FreeQ[\ A, \ Complex]\), \ imagQ = True]; \n\t If[imagQ, \ \n\t\tJoin[{1 + I}, \ Select[Range[ Sqrt[maxnorm[ A]]], \((Mod[#, 4] == 3\ && \ \ PrimeQ[#])\)\ &], \n\t Gprime\ /@ \ \tSelect[ Range[maxnorm[ A]], \((Mod[#, 4] == 1\ && \ \ PrimeQ[#])\)\ &]], \n\t\tSelect[ Range[Sqrt[maxnorm[A]]], \ PrimeQ]])\)\t;\)\n\t\t\), "\n", \(\(PrimeZetaGaussian[1, \ ___]\ := \ Infinity;\)\), "\n", \(\(PrimeZeta[1, \ ___]\ := \ Infinity;\)\ \n\), "\n", \(\(\(Options[PrimeZeta]\ = \ {GaussianIntegers -> False};\)\(\n\) \)\), "\n", \(\(\(PrimeZeta[s_, \ prec_Integer : 6, \ opts___Rule]\ := \((\t\timagQ = \(GaussianIntegers\ /. \ {opts}\)\ \ /. \ Options[PrimeZeta]; \n\t\t\tIf[\(! FreeQ[\ A, \ Complex]\), \ imagQ = True]; \n\ Plus @@ \ Table[ MoebiusMu[k]\ *\ Log[\(If[imagQ, \ ZetaGaussian, Zeta]\)[k\ s]]/k, \ \ {k, IndexBoundPrimeZeta[s, \ Ceiling[ Log[10, 1/2\ 10^\((prec + 1)\)\ 2^s]]]}])\);\)\(\n\) \)\), "\n", \(\(IndexBoundHLProductGaussian[prec_, \ len_, \ L_]\ := \ Module[{M = NextNorm[L]}, \n\t\t\tFloor[ ProductLog[\(4\ 10\^\(prec + 1\)\ M\ \((M + 1)\)\ Log[M\/len]\)\/\ \(M - len\)]\/Log[M\/len]]];\)\), "\n", \(\(\(IndexBoundPrimeZeta[s_, \ acc_]\ := Floor[ProductLog[\(\(2\^\(\(1\)\(+\)\(s\)\(\ \)\)\) 10\^\(1 + acc\)\ \ s\ \((1 + s)\)\ Log[2]\)\/\(\((1 - 2\^s)\)\ \((1 - s)\)\)]\/\(s\ Log[2]\)];\)\ \(\n\) \)\), "\n", \(\(IndexBoundHLProductRational[prec_, \ len_, \ L_]\ := \ \n\t If[L == 2, \ \ \n\t\t\t\tMax[ Ceiling[2\ ProductLog[ 10^\((prec/2)\)* Sqrt[len/\((4 - len)\)]*\n\ \ \ \ \ Log[4/len]]/ Log[4/len]\ ], \n\t\t\t\t\t\t\tFloor[\ ProductLog[3\ /2\ \ 10\^\(prec + 1\)\ Log[3/2]]\/Log[3/2]\ ]], \ \n\ \ \ \ \ \ \ \ M = NextPrime[L]; \n\t\t\tCeiling[ 2\ ProductLog[ 1/2\ \ \ 10^\((prec/2)\)*\n\ \ \ \ \ Sqrt[\ len \((M - 1)\)/\((M - 1 - len)\)]* Log[\((M - 1)\)/len]]/ Log[\((M - 1)\)/len]]\ ];\)\)}], "Input", PageWidth->PaperWidth], Cell["\<\ Options[HardyLittlewood] = {GaussianIntegers -> False}; HardyLittlewood[A_List, prec_Integer:6, opts___Rule] := Module[{n = Length[A], primess, imagQ}, imagQ = GaussianIntegers /. {opts} /. Options[HardyLittlewood]; If[!FreeQ[A, Complex], imagQ = True]; primes = CriticalPrimeInterval[A, opts]; max = If[imagQ, Max[GaussianNorm /@ primes], Max[primes]]; If[imagQ, primes = Join[{1 + I}, Select[Range[Sqrt[max]], Mod[#1, 4] == 3 && PrimeQ[#1] & ], Flatten[Gprime /@ Select[Range[max], Mod[#1, 4] == 1 && PrimeQ[#1] & ] /. Complex[a_, b_] :> {a+b*I, b+a*I}]]]; If[ !AdmissibleQ[A, opts], Return[0]]; If[imagQ, primes3mod4 = Select[primes, IntegerQ]; prime2 = If[MemberQ[primes, 1 + I], {1 + I}, {}]; primesimag = Cases[DeleteCases[primes, 1 + I], _Complex]; factor2 = If[prime2 == {1 + I}, 2^(n - 1), 1]; factor3mod4 = Times @@ ((#^2 - Length[Union[Mod[A, #]]])/#^2 * (#^2/(#^2 - 1))^n & ) /@ primes3mod4; factorimag = Times @@ ((gn = GaussianNorm[#1]; (gn - Length[Union[Mod[A, #1]]])/gn*(gn/(gn - 1))^n) & ) /@ primesimag; rat = Times @@ {factor2, factor3mod4, factorimag}; HLFunction[rat, n, NextNorm[max], GaussianIntegers -> True], rat = Times @@ ((# - Length[Union[Mod[A, #]]])/# * (#/(# - 1))^n & ) /@ primes; HLFunction[rat, n, NextPrime[max]]]]\ \>", "Input", PageWidth->PaperWidth], Cell["\<\ Options[HLFunction] = {GaussianIntegers -> False}; Format[HLFunction[rat_, n_Integer, k_, opts___]] := (imagQ = GaussianIntegers /. {opts} /. Options[HLFunction]; DisplayForm[RowBox[Map[MakeBoxes[#]&, {FactorForm[FactorInteger[Numerator[rat]]]/ If[IntegerQ[rat], 1, FactorForm[FactorInteger[Denominator[rat]]]], \ EulerProduct[n, k, GaussianIntegers -> imagQ]}, {1}]]]); \ \>", \ "Input", PageWidth->PaperWidth], Cell["\<\ HLFunction /: N[HLFunction[rat_, n_, k_,opts___Rule], prec_:6] := (N[rat* EulerProduct[n, \ k, opts], prec]) /. -1. -> 1; HLFunction /: Round[HLFunction[rat_, n_, k_,opts___Rule], prec_:6] := Round[(N[rat* \ EulerProduct[n, k, opts], prec])]; (*HLFunction /: (g_?(#=!=Format && NumericFunctionQ[#])&[HLFunction[rat_, n_, \ k_,opts___Rule], prec_:6]) := (g[N[HLFunction[rat, n, k,opts], prec]]) \ *)\ \>", "Input", PageWidth->PaperWidth], Cell["\<\ IntegrateHold /: N[IntegrateHold[f_, {x_, a_, b_}], prec_:6] := Module[{temp = Integrate[f /. x -> ToExpression[x], ToExpression[x]]}, temp /. x -> ToExpression[b]]; Clear[ConstellationEstimate, IntegralFunction]; Options[ConstellationEstimate] = {GaussianIntegers -> False, EvaluateIntegral -> True, \ FormatOutputNicely->True}; fix[expr_] := expr /. xx_Integer /; xx != -1 :> ExponentForm[xx]; Options[IntegralFunction] = {GuassianIntegers -> False}; IntegralFunction[A_, t_, opts___] := IntegralFunction[A, t,opts] =( imagQ = GaussianIntegers /. {opts} /. Options[IntegralFuntion]; If[!FreeQ[A, Complex], imagQ = True]; If[imagQ, Integrate[t/Log[t]^Length[A], t] /. ExpIntegralEi[2*Log[t]] -> LogIntegral[t^2], Integrate[1/Log[t]^Length[A], t]]); ConstellationEstimate[A_, opts___] := Module[{t}, {intQ, imagQ, intformatQ} = {EvaluateIntegral, GaussianIntegers, \ FormatOutputNicely} /. {opts} /. Options[ConstellationEstimate]; functemp = If[intformatQ,Function[z, SFTF[If[!FreeQ[z, _Real],Rationalize, Identity][z] /.eform]], Identity]; $Post = Function[z, zz=If[Head[functemp] =!= Symbol, functemp[z], z];Clear[functemp];zz]; If[!FreeQ[A, Complex], imagQ = True]; If[imagQ, Function[x, Evaluate[HardyLittlewood[A, GaussianIntegers -> True]* (4/Pi^3; \ 2^(Length[A]-2)/Pi^(Length[A]-1))* (If[intQ, IntegralFunction[A,x,GaussianIntegers -> True], IntegrateHold[\"r\"/Log[\"r\"]^Length[A], {\"r\", 0, x}]])]], Function[x, Evaluate[HardyLittlewood[A]* (If[intQ, IntegralFunction[A,x,GaussianIntegers->False], IntegrateHold[1/Log[\"t\"]^Length[A], {\"t\", 0, x}]])]]]];\ \>", \ "Input", PageWidth->PaperWidth], Cell["\<\ Format[IntegrateHold[f_, {x_, a_, b_}]] := DisplayForm[RowBox[{SubsuperscriptBox[\"\[Integral]\", (MakeBoxes[#] &)[a], (MakeBoxes[#] &)[b]], RowBox[{(MakeBoxes[#] &)[f], RowBox[{\"\[DifferentialD]\", ToString[x]}]}]}]]\ \>", "Input", PageWidth->PaperWidth], Cell["\<\ Options[EulerProduct] = {GaussianIntegers -> False}; EulerProduct /: N[EulerProduct[len_, L_, opts___], prec_:6] := (imagQ = GaussianIntegers /. {opts} /. Options[EulerProduct]; If[imagQ, EulerProductNGaussian, EulerProductNRational][Round[len], \ Round[L], Round[prec]]); \ \>", "Input", PageWidth->PaperWidth], Cell["\<\ Format[EulerProduct[n_, k_, opts___]] := (imagQ = GaussianIntegers /. {opts} /. Options[EulerProduct]; base = If[imagQ, \"N(p)\", \"p\"]; DisplayForm[RowBox[{UnderscriptBox[\"\[Product]\", RowBox[{base, \"\[GreaterEqual]\[ThinSpace]\", ToString[k]}]], FractionBox[RowBox[{RowBox[{\"(\", RowBox[{base, \"-\", ToString[n]}], \")\"}], \" \", If[n - 1 > 1, SuperscriptBox[base, ToString[n - 1]], base]}], SuperscriptBox[RowBox[{\"(\", RowBox[{base, \"-\", \"1\"}], \")\"}], ToString[n]]]}]]); \ \>", "Input", PageWidth->PaperWidth], Cell["\<\ EulerProductNRational[len_, L_, prec_] := Module[{term1}, sum1 = 0; cp = Drop[Prime[Range[PrimePi[L]]], -1]; bd = IndexBoundHLProductRational[prec, len, Max[cp]]; Do[smallprimecorrex = -Plus @@ (1/cp^j); factor = (len - len^j)/j; term1 = Re[N[factor*(PrimeZeta[j, If[j == 2, Max[prec, 17], Max[prec, prec + 1 + Ceiling[Log[10, Abs[smallprimecorrex/(term1/factor)]]] ]]] + smallprimecorrex), Max[prec, 17]]]; sum1 += term1; If[traceQ && Mod[j, traceI] == 0, Print[{j, term1, sum1, Abs[term1/sum1]}]], {j, 2, bd}]; N[E^sum1, prec]]; \ \>", "Input", PageWidth->PaperWidth], Cell["\<\ EulerProductNGaussian[len_, L_, prec_] := Module[{term1}, cp = Select[Join[{1 + I}, Select[Range[Sqrt[L]], Mod[#, 4] == 3 && PrimeQ[#] &], Gprime /@ Select[Range[L], Mod[#, 4] == 1 && PrimeQ[#] & ]], GaussianNorm[#] < L & ]; cp1 = GaussianNorm /@ Cases[DeleteCases[cp, 1 + I], _Complex]; cp3 = GaussianNorm /@ Select[cp, IntegerQ]; sum1 = 0; bd = IndexBoundHLProductGaussian[prec, len, Max[GaussianNorm /@ cp]]; If[traceQ, Print[\"The number of terms in the main sum is \", bd]]; Do[If[j > 2 && traceQ && Mod[j, traceI] == 0, Print[{j, term1, sum1, {Precision[term1], Precision[sum1]}, Abs[term1/sum1]}]]; factor = (len - len^j)/j; smallprimecorrex = -(Plus @@ (2/cp1^j) + Plus @@ (1/cp3^j) + If[MemberQ[cp, 1 + I], 1/2^j, 0]); term1 = Re[N[factor*(PrimeZeta[j, If[j == 2, Max[prec, 17], Max[prec, prec + 1 + Ceiling[Log[10, Abs[smallprimecorrex/(term1/factor)]]] ]], GaussianIntegers -> True] + smallprimecorrex), Max[17, prec]]]; sum1 += term1, {j, 2, bd}]; N[E^sum1, prec]] eform = {x_Integer /; Abs[x] > 10009 :> fix[x], x_Real /; N[x - Round[x]] == 0 :> Round[x]}; fixspecialcase[z_] := ToString[4/\[Pi]^3, StandardForm]* \[Pi]^3 /4 z \ \>", \ "Input", PageWidth->PaperWidth], Cell["\<\ Options[ConstellationEstimateAsymptotic] = {GaussianIntegers-> \ False}; ConstellationEstimateAsymptotic[a_, opts___] := Module[{x,n = Length[a], ce}, \ gQ = GaussianIntegers /. {opts} /. Options[ConstellationEstimateAsymptotic]; gQ = gQ || !FreeQ[a, Complex]; ce = ConstellationEstimate[a][x]; If[ce === 0, Return[0 & ]]; If[ !gQ, ce=ce[[1]]; Return[Function[xx, DisplayForm[RowBox[{ MakeBoxes[#]&[ce], StringForm[\"\\!\\((``\\/Log[``]\\^``)\\)\", xx, xx, n], \"+\", \ StringForm[\"O[\\!\\(``\\/Log[``]\\^``\\)]\", xx, xx, n + 1]}]]]] ]; ce = ce[[{1,2,3}]]; Function[xx, DisplayForm[RowBox[{ MakeBoxes[#]&[ce], StringForm[\"\\!\\(``\\^2\\/\\(2 Log[``]\\^``\\)\\)\", xx, xx, n], \"+\", StringForm[\"O[\\!\\(``\\^2\\/Log[``]\\^``\\)]\", xx, xx, n + 1]}]]]]\ \>", "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Examples in \[DoubleStruckCapitalZ]", "Section"], Cell[TextData[{ "The standard usage is to feed in an admissible pattern, such as ", Cell[BoxData[ \(TraditionalForm\`{0, 2, 6}\)]], ", and an upper bound, which can be just the symbol x. The output is \ formatted nicely and the Euler product part is left as an infinite product." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(ConstellationEstimate[{0, 2, \ 6}]\)[x]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[ RowBox[{ RowBox[{"(", FormBox[ TagBox[ RowBox[{ FractionBox[ InterpretationBox[ RowBox[{ TagBox[\(\(3\^2\)\(\ \)\), DisplayForm], "\[InvisibleSpace]", TagBox[\(\(5\^2\)\(\ \)\), DisplayForm]}], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "3", 2], " "}]], DisplayForm[ RowBox[ { SuperscriptBox[ "5", 2], " "}]]], Editable->False], InterpretationBox[ TagBox[\(\(2\^6\)\(\ \)\), DisplayForm], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "2", 6], " "}]]], Editable->False]], TagBox[\(\[Product]\+\(p \ \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) 7\)\(\((p - 3)\)\ p\^2\)\/\((p - 1)\)\ \^3\), DisplayForm]}], DisplayForm], "TraditionalForm"], ")"}], " ", \((\(-\(x\/\(2\ \(log(x)\)\)\)\) - x\/\(2\ \(\(log\^2\)(x)\)\) + \(li(x)\)\/2)\)}], "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell["One may want to suppress some of the fancy formatting.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(ConstellationEstimate[{0, 2, \ 6}, \ FormatOutputNicely \[Rule] False]\)[x]\)], "Input"], Cell[BoxData[ RowBox[{ RowBox[{"(", TagBox[ RowBox[{ FractionBox[ InterpretationBox[ RowBox[{ TagBox[\(\(3\^2\)\(\ \)\), DisplayForm], "\[InvisibleSpace]", TagBox[\(\(5\^2\)\(\ \)\), DisplayForm]}], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "3", 2], " "}]], DisplayForm[ RowBox[ { SuperscriptBox[ "5", 2], " "}]]], Editable->False], InterpretationBox[ TagBox[\(\(2\^6\)\(\ \)\), DisplayForm], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "2", 6], " "}]]], Editable->False]], TagBox[\(\[Product]\+\(p \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) \ 7\)\(\((p - 3)\)\ p\^2\)\/\((p - 1)\)\^3\), DisplayForm]}], DisplayForm], ")"}], " ", \((\(-\(x\/\(2\ Log[x]\^2\)\)\) - x\/\(2\ Log[x]\) + LogIntegral[x]\/2)\)}]], "Output"] }, Open ]], Cell["One may want to suppress evaluation of the integral.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(ConstellationEstimate[{0, 2, \ 6}, \ EvaluateIntegral \[Rule] False]\)[x]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[ RowBox[{ RowBox[{"(", FormBox[ TagBox[ RowBox[{ FractionBox[ InterpretationBox[ RowBox[{ TagBox[\(\(3\^2\)\(\ \)\), DisplayForm], "\[InvisibleSpace]", TagBox[\(\(5\^2\)\(\ \)\), DisplayForm]}], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "3", 2], " "}]], DisplayForm[ RowBox[ { SuperscriptBox[ "5", 2], " "}]]], Editable->False], InterpretationBox[ TagBox[\(\(2\^6\)\(\ \)\), DisplayForm], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "2", 6], " "}]]], Editable->False]], TagBox[\(\[Product]\+\(p \ \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) 7\)\(\((p - 3)\)\ p\^2\)\/\((p - 1)\)\ \^3\), DisplayForm]}], DisplayForm], "TraditionalForm"], ")"}], " ", FormBox[ TagBox[\(\[Integral]\_0\%x\( 1\/Log["t"]\^3\) \[DifferentialD]t\), DisplayForm], "TraditionalForm"]}], "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell["Another function can be used to get the asymptotic form", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(ConstellationEstimateAsymptotic[{0, 2, \ 6}]\)[x]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[ TagBox[ RowBox[{ TagBox[ RowBox[{ FractionBox[ InterpretationBox[ RowBox[{ TagBox[\(\(3\^2\)\(\ \)\), DisplayForm], "\[InvisibleSpace]", TagBox[\(\(5\^2\)\(\ \)\), DisplayForm]}], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "3", 2], " "}]], DisplayForm[ RowBox[ { SuperscriptBox[ "5", 2], " "}]]], Editable->False], InterpretationBox[ TagBox[\(\(2\^6\)\(\ \)\), DisplayForm], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "2", 6], " "}]]], Editable->False]], TagBox[\(\[Product]\+\(p \ \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) 7\)\(\((p - 3)\)\ p\^2\)\/\((p - 1)\)\ \^3\), DisplayForm]}], DisplayForm], InterpretationBox["\<\"\\!\\((\\!\\(TraditionalForm\\`x\\)\\/Log[\ \\!\\(TraditionalForm\\`x\\)]\\^\\!\\(TraditionalForm\\`3\\))\\)\"\>", StringForm[ "\!\((``\/Log[``]\^``)\)", x, x, 3], Editable->False], "+", InterpretationBox["\<\"O[\\!\\(\\!\\(TraditionalForm\\`x\\)\\/Log[\ \\!\\(TraditionalForm\\`x\\)]\\^\\!\\(TraditionalForm\\`4\\)\\)]\"\>", StringForm[ "O[\!\(``\/Log[``]\^``\)]", x, x, 4], Editable->False]}], DisplayForm], "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell["\<\ One can get separately the Euler products that show up in these \ expressions in typeset form.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(EulerProduct[3, 7]\)], "Input"], Cell[BoxData[ TagBox[\(\[Product]\+\(p \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) \ 7\)\(\((p - 3)\)\ p\^2\)\/\((p - 1)\)\^3\), DisplayForm]], "Output"] }, Open ]], Cell["And one can get their numerical values.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(N[EulerProduct[3, 7]]\)], "Input"], Cell[BoxData[ \(0.813013044443062`\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(N[EulerProduct[3, 7], 20]\)], "Input"], Cell[BoxData[ \(0.813012933893467145224574714754785`20\)], "Output"] }, Open ]], Cell["N[] will automatically numericize the infinite product.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(N[\(ConstellationEstimate[{0, 2, \ 6}]\)[x]]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[\(2.8582489843701397`\ \((\(-\(x\/\(2\ \(log \((x)\)\)\)\)\) - x\/\(2\ \(\(log\^2\) \((x)\)\)\) + \(li \((x)\)\)\/2)\)\), "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell["High-precision is possible.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(N[\(ConstellationEstimate[{0, 2, \ 6}]\)[x], \ 20]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[\(2.85824859571922043243014548155794`19.699\ \((\(-\(x\/\(2\ \ \(log \((x)\)\)\)\)\) - x\/\(2\ \(\(log\^2\) \((x)\)\)\) + \(li \((x)\)\)\/2)\)\), "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(N[\(ConstellationEstimate[{0, 2, 6}]\)[10\^10]]\)], "Input"], Cell[BoxData[ TagBox[ FormBox["2.7152836321207583`*^6", "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ "Examples in ", Cell[BoxData[ \(TraditionalForm\`\[DoubleStruckCapitalZ][\[ImaginaryI]]\)]] }], "Section"], Cell["The castle and diamond are pre-defined.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(castle\)], "Input"], Cell[BoxData[ \({0, 2, 4, 1 + \[ImaginaryI], 3 + \[ImaginaryI], 2\ \[ImaginaryI], 4 + 2\ \[ImaginaryI], 1 + 3\ \[ImaginaryI], 3 + 3\ \[ImaginaryI], 4\ \[ImaginaryI], 2 + 4\ \[ImaginaryI], 4 + 4\ \[ImaginaryI]}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(diamond\)], "Input"], Cell[BoxData[ \({\(-1\), 1, \[ImaginaryI], \(-\[ImaginaryI]\)}\)], "Output"] }, Open ]], Cell["Here is the estimate for prime occurrences of the diamond.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(diamondestimate = \(ConstellationEstimate[diamond]\)[x]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[ FractionBox[ RowBox[{"4", " ", RowBox[{"(", FormBox[ TagBox[ RowBox[{ InterpretationBox[ TagBox[\(\(2\^3\)\(\ \)\), DisplayForm], SequenceForm[ DisplayForm[ RowBox[ { SuperscriptBox[ "2", 3], " "}]]], Editable->False], TagBox[\(\[Product]\+\(\(N \((p)\)\) \ \(\(\[GreaterEqual]\)\(\[ThinSpace]\)\) 5\)\(\((N \((p)\) - 4)\)\ \(N \((p)\)\ \)\^3\)\/\((N \((p)\) - 1)\)\^4\), DisplayForm]}], DisplayForm], "TraditionalForm"], ")"}], " ", \((\(-\(\(2\ x\^2\)\/\(3\ \(log(x)\)\)\)\) - x\^2\/\(3\ \(\(log\^2\)(x)\)\) - x\^2\/\(3\ \(\(log\^3\)(x)\)\) + \(4\ \(li(x\^2)\)\)\/3)\)}], \ \(\[Pi]\^3\)], "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(N[diamondestimate]\)], "Input"], Cell[BoxData[ \(0.17797544738231955`\ \((\(-\(\(0.3333333333333333`\ \ x\^2\)\/Log[x]\^3\)\) - \(0.3333333333333333`\ x\^2\)\/Log[x]\^2 - \ \(0.6666666666666666`\ x\^2\)\/Log[x] + 1.3333333333333333`\ LogIntegral[x\^2])\)\)], "Output"] }, Open ]], Cell["\<\ Here is a prediction of the radius at which there should be 1000 \ diamonds.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(FindRoot[ diamondestimate\ \[Equal] \ 1000, \ {x, {100, 10000}}]\)], "Input"], Cell[BoxData[ \({x \[Rule] 7277.460567273188`}\)], "Output"] }, Open ]], Cell["Here is the castle estimate.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(castleestimate = N[\(ConstellationEstimate[castle]\)[x]]\)], "Input"], Cell[BoxData[ TagBox[ FormBox[\(0.19056744720761584`\ \((\(-\(\(4\ x\^2\)\/\(155925\ \(log( x)\)\)\)\) - \(2\ x\^2\)\/\(155925\ \(\(log\^2\)(x)\)\ \) - \(2\ x\^2\)\/\(155925\ \(\(log\^3\)(x)\)\) - x\^2\/\(51975\ \(\(log\^4\)(x)\)\) - \(2\ x\^2\)\/\(51975\ \ \(\(log\^5\)(x)\)\) - x\^2\/\(10395\ \(\(log\^6\)(x)\)\) - x\^2\/\(3465\ \(\(log\^7\)(x)\)\) - x\^2\/\(990\ \(\(log\^8\)(x)\)\) - \(2\ x\^2\)\/\(495\ \ \(\(log\^9\)(x)\)\) - x\^2\/\(55\ \(\(log\^10\)(x)\)\) - x\^2\/\(11\ \(\(log\^11\)(x)\)\) + \(8\ \(li(x\^2)\)\)\/155925)\ \)\), "TraditionalForm"], TraditionalForm, Editable->True]], "Output"] }, Open ]], Cell[TextData[{ "Here is an estimate for one castle. Our paper shows the one we found near \ ", Cell[BoxData[ \(TraditionalForm\`10\^7\)]], "." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(FindRoot[castleestimate \[Equal] 1, \ {x, {10^6, \ 10^9}}]\)], "Input"], Cell[BoxData[ \({x \[Rule] 1.0420046040834072`*^8}\)], "Output"] }, Open ]] }, Closed]] }, FrontEndVersion->"4.1 for Macintosh", ScreenRectangle->{{0, 1024}, {0, 748}}, WindowSize->{690, 659}, WindowMargins->{{50, Automatic}, {Automatic, 28}}, MacintoshSystemPageSetup->"\<\ 00<0001804P000000]P2:?oQon82n@960dL5:0?l0080001804P000000]P2:001 0000I00000400`<300000BL?00400@0000000000000006P801T1T00000000000 00000000000000000000000000000000\>" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[1705, 50, 264, 6, 78, "Text"], Cell[CellGroupData[{ Cell[1994, 60, 65, 1, 56, "Section"], Cell[2062, 63, 236, 9, 132, "Input"], Cell[2301, 74, 563, 12, 171, "Input"], Cell[2867, 88, 693, 19, 357, "Input"], Cell[3563, 109, 673, 20, 312, "Input"], Cell[4239, 131, 1293, 26, 462, "Input"], Cell[5535, 159, 1300, 47, 747, "Input"], Cell[6838, 208, 3619, 75, 1213, "Input"], Cell[10460, 285, 1436, 32, 597, "Input"], Cell[11899, 319, 459, 13, 192, "Input"], Cell[12361, 334, 454, 15, 222, "Input"], Cell[12818, 351, 1752, 47, 822, "Input"], Cell[14573, 400, 289, 7, 102, "Input"], Cell[14865, 409, 330, 9, 162, "Input"], Cell[15198, 420, 593, 13, 207, "Input"], Cell[15794, 435, 642, 18, 282, "Input"], Cell[16439, 455, 1318, 33, 567, "Input"], Cell[17760, 490, 819, 21, 312, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[18616, 516, 54, 0, 36, "Section"], Cell[18673, 518, 303, 6, 46, "Text"], Cell[CellGroupData[{ Cell[19001, 528, 74, 1, 27, "Input"], Cell[19078, 531, 1604, 42, 75, "Output"] }, Open ]], Cell[20697, 576, 70, 0, 30, "Text"], Cell[CellGroupData[{ Cell[20792, 580, 119, 2, 27, "Input"], Cell[20914, 584, 1249, 34, 63, "Output"] }, Open ]], Cell[22178, 621, 68, 0, 30, "Text"], Cell[CellGroupData[{ Cell[22271, 625, 117, 2, 27, "Input"], Cell[22391, 629, 1683, 46, 63, "Output"] }, Open ]], Cell[24089, 678, 71, 0, 30, "Text"], Cell[CellGroupData[{ Cell[24185, 682, 84, 1, 27, "Input"], Cell[24272, 685, 1914, 49, 57, "Output"] }, Open ]], Cell[26201, 737, 118, 3, 30, "Text"], Cell[CellGroupData[{ Cell[26344, 744, 51, 1, 27, "Input"], Cell[26398, 747, 159, 3, 57, "Output"] }, Open ]], Cell[26572, 753, 55, 0, 30, "Text"], Cell[CellGroupData[{ Cell[26652, 757, 54, 1, 27, "Input"], Cell[26709, 760, 52, 1, 26, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[26798, 766, 58, 1, 27, "Input"], Cell[26859, 769, 72, 1, 26, "Output"] }, Open ]], Cell[26946, 773, 71, 0, 30, "Text"], Cell[CellGroupData[{ Cell[27042, 777, 77, 1, 27, "Input"], Cell[27122, 780, 261, 6, 44, "Output"] }, Open ]], Cell[27398, 789, 43, 0, 30, "Text"], Cell[CellGroupData[{ Cell[27466, 793, 83, 1, 27, "Input"], Cell[27552, 796, 285, 7, 44, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[27874, 808, 80, 1, 31, "Input"], Cell[27957, 811, 150, 5, 28, "Output"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[28156, 822, 132, 4, 36, "Section"], Cell[28291, 828, 55, 0, 30, "Text"], Cell[CellGroupData[{ Cell[28371, 832, 39, 1, 27, "Input"], Cell[28413, 835, 248, 4, 26, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[28698, 844, 40, 1, 27, "Input"], Cell[28741, 847, 80, 1, 26, "Output"] }, Open ]], Cell[28836, 851, 74, 0, 30, "Text"], Cell[CellGroupData[{ Cell[28935, 855, 88, 1, 27, "Input"], Cell[29026, 858, 1131, 30, 69, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[30194, 893, 51, 1, 27, "Input"], Cell[30248, 896, 250, 4, 47, "Output"] }, Open ]], Cell[30513, 903, 100, 3, 30, "Text"], Cell[CellGroupData[{ Cell[30638, 910, 104, 2, 27, "Input"], Cell[30745, 914, 64, 1, 26, "Output"] }, Open ]], Cell[30824, 918, 44, 0, 30, "Text"], Cell[CellGroupData[{ Cell[30893, 922, 89, 1, 27, "Input"], Cell[30985, 925, 716, 14, 139, "Output"] }, Open ]], Cell[31716, 942, 167, 6, 30, "Text"], Cell[CellGroupData[{ Cell[31908, 952, 91, 1, 27, "Input"], Cell[32002, 955, 68, 1, 28, "Output"] }, Open ]] }, Closed]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)