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)

Hardy-Littlewood Code

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};

\(\(\(SetAttributes[GaussianMod, \ Listable];\)\(
\)
    \)\)
\(\ \(\(GaussianMod[u_?GaussianIntegerQ, \ 
        z_?GaussianIntegerQ]\  := \ 
\ \ Module[{d, q, \ a\  = \ Re[z], \ 
          b = Im[z], \ x = Re[u], \ y = Im[u]}, 
\ \ {d, \ e}\  = \ 
          ExtendedGCD[a, b]; 
\ \ \ q\  = \ 
          Quotient[y, d]; 
\ \ \ Mod[
            x\  - \ a\ q\ e[\([2]\)]\  + \ 
              b\ q\ e[\([1]\)], 
	\ \((a^2 + b^2)\)/d]\ \  + 
          I\ \((y\  - \ q\ d)\)]\)\(
\)
    \)\)

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]

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}]

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 ≥ 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.";

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]}]]]

\(\($MaxExtraPrecision = 300;\)\ \ \ 
\)
\(\(ZetaGaussian[s_]\  := 
        Zeta[s]\ \ 4\^\(-s\)\ \((Zeta[s, 1/4] - 
              Zeta[s, 3/4])\);\)
		\)
distancesq[a_,  b_]  :=  Abs[b - a]^2;
\(\(\(maxnorm[set_]\  := 
        Max[Apply[distancesq, 
            DeleteCases[
              Distribute[{set, set}, \ List], \ {x_, \ x_}], \ {1}]];\)\(
\)
    \)\)
\(\(\(NextNorm[L_]\  := \ 
        Module[{t = L + 1}, \ 
          While[\(! \((
					\((IntegerQ[tt = Sqrt[t]] && 
                      Mod[tt, 4] == 3\  && \ PrimeQ[tt])\) || 
		\((	
                    Mod[t, 4] == 1\  && \ PrimeQ[t])\))\)\), \ \(t++\)]; \ 
          t];\)\(
\)
    \)\)
\(\(\(Gprime[
          n_]\  := \ \(FactorInteger[n, \ 
            GaussianIntegers -> True]\)[\([\(-1\), 1]\)];\)\(
\)
    \)\)
\(\(\(Options[\ 
          CriticalPrimeInterval]\  = \ {GaussianIntegers -> False};\)\(
\)
    \)\)
\(\(CriticalPrimeInterval[A_, \ opts___]\  := \ \((	
          imagQ = \(GaussianIntegers\  /. \ {opts}\)\  /. \ 
              Options[\ CriticalPrimeInterval]; 
			If[\(! FreeQ[\ A, \ 
                Complex]\), \ imagQ = True]; 
	
          If[imagQ, \ 
		Join[{1 + I}, \ 
              Select[Range[
                  Sqrt[maxnorm[
                      A]]], \((Mod[#, 4] == 3\  && \ \ PrimeQ[#])\)\  &], 
	
              Gprime\  /@ \ 	Select[
                  Range[maxnorm[
                      A]], \((Mod[#, 4] == 
                          1\  && \ \ PrimeQ[#])\)\  &]], 
		Select[
              Range[Sqrt[maxnorm[A]]], \ PrimeQ]])\)	;\)
		\)
PrimeZetaGaussian[1,  ___]  :=  Infinity;
\(\(PrimeZeta[1, \ ___]\  := \ Infinity;\)\ 
\)
\(\(\(Options[PrimeZeta]\  = \ {GaussianIntegers -> False};\)\(
\)
    \)\)
\(\(\(PrimeZeta[s_, \ prec_Integer : 6, \ 
          opts___Rule]\  := \((		imagQ = \(GaussianIntegers\  /. \ {opts}\)\  /. \ Options[PrimeZeta]; 
			If[\(! FreeQ[\ A, \ Complex]\), \ 
            imagQ = True]; 
\ 
          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]]]}])\);\)\(
\)
    \)\)
\(\(IndexBoundHLProductGaussian[prec_, \ len_, \ L_]\  := \ 
        Module[{M = NextNorm[L]}, 
			Floor[
            ProductLog[\(4\ 10\^\(prec + 1\)\ M\ \((M + 1)\)\ Log[M\/len]\)\/\(M - len\)]\/Log[M\/len]]];\)\)
\(\(\(IndexBoundPrimeZeta[s_, \ acc_]\  := 
        Floor[ProductLog[\(\(2\^\(\(1\)\(+\)\(s\)\(\ \)\)\) 10\^\(1 + acc\)\ s\ \((1 + s)\)\ Log[2]\)\/\(\((1 - 2\^s)\)\ \((1 - s)\)\)]\/\(s\ Log[2]\)];\)\(
\)
    \)\)
\(\(IndexBoundHLProductRational[prec_, \ len_, \ L_]\  := \ 
	
        If[L == 2, \ \ 
				Max[
            Ceiling[2\ ProductLog[
                    10^\((prec/2)\)*
                      Sqrt[len/\((4 - len)\)]*
\ \ \ \ \ Log[4/len]]/
                  Log[4/len]\ ], 
							Floor[\ 
              ProductLog[3\ /2\ \ 10\^\(prec + 1\)\ Log[3/2]]\/Log[3/2]\ ]], 
\ \ \ \ \ \ \ \ M = NextPrime[L]; 
			Ceiling[
            2\ ProductLog[
                  1/2\ \ \ 10^\((prec/2)\)*
\ \ \ \ \ Sqrt[\ 
                      len \((M - 1)\)/\((M - 1 - len)\)]*
                    Log[\((M - 1)\)/len]]/
                Log[\((M - 1)\)/len]]\ ];\)\)

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]]]]

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}]]]);

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]]) *)

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}]])]]]];

Format[IntegrateHold[f_, {x_, a_, b_}]] := 
  DisplayForm[RowBox[{SubsuperscriptBox["∫", 
  (MakeBoxes[#] &)[a], 
      (MakeBoxes[#] &)[b]], RowBox[{(MakeBoxes[#] &)[f], 
       RowBox[{"\[DifferentialD]", ToString[x]}]}]}]]

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]]);

Format[EulerProduct[n_, k_, opts___]] := 
   (imagQ = GaussianIntegers /. {opts} /. Options[EulerProduct]; 
  base = If[imagQ, "N(p)", "p"]; 
  DisplayForm[RowBox[{UnderscriptBox["∏", 
    RowBox[{base, "≥ ", 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]]]}]]);

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]];

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/π^3, StandardForm]* π^3 /4 z

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]}]]]]

Examples in ℤ

The standard usage is to feed in an admissible pattern, such as {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.

ConstellationEstimate[{0, 2,  6}][x]

One may want to suppress some of the fancy formatting.

ConstellationEstimate[{0, 2,  6},  
        FormatOutputNicely -> False][x]

One may want to suppress evaluation of the integral.

ConstellationEstimate[{0, 2,  6},  
        EvaluateIntegral -> False][x]

Another function can be used to get the asymptotic form

ConstellationEstimateAsymptotic[{0, 2,  6}][x]

One can get separately the Euler products that show up in these expressions in typeset form.

EulerProduct[3, 7]

And one can get their numerical values.

N[EulerProduct[3, 7]]

N[EulerProduct[3, 7], 20]

N[] will automatically numericize the infinite product.

N[ConstellationEstimate[{0, 2,  6}][x]]

High-precision is possible.

N[ConstellationEstimate[{0, 2,  6}][x],  20]

N[ConstellationEstimate[{0, 2, 6}][(10^10)]]

Examples in ℤ[i]

The castle and diamond are pre-defined.

castle

diamond

Here is the estimate for prime occurrences of the diamond.

diamondestimate = ConstellationEstimate[diamond][x]

N[diamondestimate]

Here is a prediction of the radius at which there should be 1000 diamonds.

FindRoot[
      diamondestimate  ==  1000,  {x, {100, 10000}}]

Here is the castle estimate.

castleestimate = N[ConstellationEstimate[castle][x]]

Here is an estimate for one castle. Our paper shows the one we found near (10^7).

FindRoot[castleestimate == 1,  {x, {10^6,  10^9}}]