I created a code which calculates


where $D$ is the density of $S\cap[a,b]$ (in $A\cap[a,b]$), $[a,b]$ is an interval for $a,b\in\mathbb{R}$, $F_n$ is the Følner Sequence of $A$, and $S\subseteq A$. For more information, click here (replace $G,X,i,g$ with $A,S,n,a$), click (here) and here.

I set $A=\left\{\frac{p}{2^k(2q+1)}:p,q,k\in\mathbb{Z},k\ge 0\right\}$, $F_n=\left\{\frac{p}{2^k(2q+1)}:p,q,k\in\mathbb{Z},k\ge 0,k\le n,|2q+1|\le n\right\}$ and $S=\left\{\frac{j^2+j+1}{k^3+1}:j,k\in\mathbb{Z},k\neq -1\right\}$

In my code $A$ is A[p_,k_,q_], $F_n$ is f[n_,a_,b_], $S$ is S[j_,k_] and $D$ is d.

Clear[A, F, f, p, Ff, S, X, Y, d, j, k];
A[p_, k_, q_] := p/((2^k)*(2*q + 1));
F[p_, n_] := 
  Table[A[p, k, q], {k, 0, Floor[Log[2, n]]}, {q, 0, 
    Floor[(n - 1)/2]}];
f[n_, a_, b_] := 
  p /. Table[
    Solve[a <= A[p, k, q] <= b, p, Integers], {k, 0, 
     Floor[Log[2, n]]}, {q, 0, Floor[(n - 1)/2]}];
Ff[n_, a_, b_] := 
     F[f[n, a, b][[v]][[u]], n][[v]][[u]], {v, 1, 
      Floor[Log[2, n]] + 1}, {u, 1, Floor[(n - 1)/2] + 1}];
S[j_, k_] := (j^2 + j + 1)/(k^3 + 1);
X[n_, a_, b_] := 
      Exists[{j, k}, S[j, k] == # && {j, k} \[Element] Integers]] & /@
     Ff[n, a, b], True];
Y[n_, a_, b_] := Length[Ff[n, a, b]];
d[n_, a_, b_] := N[(X[n, a, b])/Y[n, a, b]];
Ff[4, 1, 2]
Num = X[10, 1, 2]
Den = Y[10, 1, 2]

My result is

{1, 2, 4/3, 5/3, 3/2, 7/6, 11/6, 5/4, 7/4, 13/12, 17/12, 19/12, 23/12}



{0., 0.0653595}

Everything seems correct (Edit: I was wrong it is not correct) except the timing since the calculation must have taken more than 5 minutes.

I want to calculate n at higher numbers as quickly as possible. In fact, I want to calculate Limit[d[n,a,b],n->Infinity].

I was suggested in this post to use methods similar to monte Carlo integration? Are my new methods correct?

Edit: I found out Exists is incorrectly calculating X[n,a,b] for large n

Example: If we set $S=\left\{\frac{m^2}{n^3}:m,n\in\mathbb{Z},n\neq 0\right\}$, I get N[Num/Den] is approximately zero instead of $1$. This makes no sense since $\left\{\frac{m^2}{n^3}:m,n\in\mathbb{Z},n\neq 0\right\}=\left\{\frac{a}{b}:a,b\in\mathbb{Z},b\neq 0\right\}$, so the density should be $1$.

How do we correct Exists?


1 Answer 1


Several remarks:

  1. The timing is correctly calculated. Note that Num and Den are defined using Set (instead of SetDelayed), and so by the time the compiler executes N[Num/Den], these objects have already been calculated. Thus, Timing[N[Num/Den]] only measures the how long it takes to compute an integer division -- essentially zero, as it shows. If you want the full timing, use SetDelayed for Num and Den, or

    Timing[Num = X[10, 1, 2];
           Den = Y[10, 1, 2];
  2. You are getting a wrong density because Mathematica cannot resolve some Diophantine equations (it cannot determine whether they are soluble or not). For example,

    Resolve[Exists[{j, k}, (1 + j + j^2)/(1 + k^3) == 7/6 && {j, k} \[Element] Integers]]

    is returned unevaluated, and Count does not regard it as True. (I don't know whether this is really a false negative, because I don't know whether this equation has solutions; a quick Table search suggests that it does not, but I have no proof).

  3. You can speed up the process by using FindInstance instead of Resolve: if you take e.g. n = 5, a = 1, b = 2, then

    Count[Resolve[Exists[{j, k}, S[j, k] == # && {j, k} \[Element] Integers]] & /@ Ff[n, a, b], True] // AbsoluteTiming
    (* {38.2593, 4} *)


    Length[Cases[FindInstance[S[j, k] == #, {j, k}, Integers] & /@ Ff[n, a, b], {{j -> _, k -> _}}]] // Quiet // AbsoluteTiming
    (* {14.6932, 4} *)

    which is x2.5 faster. If you remove the Quiet you'll see all the cases where Mathematica was unable to resolve the equation; there are several of these. Unfortunately I don't think Mathematica will be able to solve these Diophantine equations unless you do some analytic groundwork first. The black box failed this time, and you'll have to work harder. Maybe ask on math.SE how to tackle this kind of equations.

  • $\begingroup$ Can other software do a better job? $\endgroup$
    – Arbuja
    Nov 25, 2019 at 19:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.