# Probability that a linear Diophantine equation has a solution

Given a linear Diophantine equation $(\mathcal D)$ with two variable $x$ and $y$:

$$(\mathcal D):\quad ax+by=c$$

I want to calculate the probability $P$ that $(\mathcal D)$ has a solution.

We know that $(\mathcal D)$ has a solution if, and only if, $d:=\gcd(a,b)$ divides $c$.

So if $a$ and $b$ are coprime, $(\mathcal D)$ has a solution, so using a result on coprime integers, we have

$$p\geqslant \mathbb P(d=1)=\frac 1{\zeta(2)}=\frac 6{\pi^2}.$$

Since $d$ must divides $c$, for each $d$ there is $\frac 1d$ integers which will work for $c$, so

$$p=\sum_{k=1}^\infty \mathbb P(d=k)\frac 1k.$$

Then we can majorate $\mathbb P(d=k)$ by $\mathbb P(d\geqslant k)=\frac 1k\times \frac 1k$ for $k\geqslant 2$.

So

$$p=\mathbb P(d=1)\times 1+\sum_{k=2}^\infty\mathbb P(d=k)\frac 1k\leqslant \frac 1{\zeta(2)}+\sum_{k=2}^\infty\frac 1{k^2}\frac 1k=\frac 1{\zeta(2)}+\zeta(3)-1.$$

What we get so far with our rough analysis is:

$$\frac 1{\zeta(2)}\leqslant p\leqslant \frac 1{\zeta(2)}+\zeta(3)-1$$

$$0.6079<p<0.8100.$$

I have read somewhere (if I understand the article correctly) that

$$\mathbb P(d=k)=\frac 1{\zeta(2)k^2}.$$

This would mean that

$$p=\sum_{k=1}^\infty \mathbb P(d=k)\frac 1k=\frac 1{\zeta(2)}\sum_{k=1}^\infty\frac 1{k^3}=\fbox{\frac {\zeta(3)}{\zeta(2)}\approx0.7308}.$$

My questions.

• Thanks to comments, I have the new following question: can I formalize this approach to make it rigorous, with a well chosen probabilistic distribution?

• I am really unused to do probabilities, so is the following approach correct?

• Can we generalize the result / the approach to find $p_n$ the probability that a linear Diophantine equation with $n$ variables has a solution?

• If we can't get an exact value for $n$ variables, can we get a numerical approximation of $p_n$?

• What is the limit of $p_n$?

• Just a quick heads-up; you have a typo on line 4. $$d := \operatorname{gcd}(a,d) \mapsto d := \operatorname{gcd}(a,b)$$ – ÍgjøgnumMeg Mar 6 '17 at 14:58
• Cool question: +1! From the Wiki page that you cite: " the probability of k randomly chosen integers being coprime is $\frac{1}{\zeta(k)}$". I don't see it straight away, but I think you can use this to cook up a number in the same way you did for the case $k=2$. Did you try to do so? – Giovanni De Gaetano Mar 6 '17 at 15:00
• The result on coprime integers bases on an assumption that you can generate integers uniformly, doesn't it? – Antoine Mar 6 '17 at 15:01
• @GiovanniDeGaetano Yes, that is what I tried, but I can't quite make it work. – E. Joseph Mar 6 '17 at 15:01
• @Antoine he is using a notion of natural density for integers, as described in Wikipedia. – Giovanni De Gaetano Mar 6 '17 at 15:01

## 1 Answer

The modeling

The text does not mention the adopted signal model. But this model is presented in the article cited. This allows us to take into account the model in which $a,\ b,\ c$ distributed over a uniform discrete law $$a,b,c\in[1,n]\cap\mathcal N,\tag1$$ $$\Pr(a=m) = \Pr(b=m) = \Pr(c=m) = \begin{cases} {1\over n},\quad\text{if}\quad m\in[1,n]\cap\mathcal N,\\ 0,\quad\text{otherwise}. \end{cases}$$

It will be also considered the series of dividers $$p_{0\dots5} = 2,3,4,5,7,9;\quad p_{k} = k\log(k)\text{ for }k>5.\tag2$$ (the table for $k>5$ is there).

The index $k_{max}$ of greatest simple divider for $k\in[1,n]$ must satisfy the inequality $$k_{max}(n)\log(k_{max}(n)) < n,$$ $$k_{max}(n) = [e^{W(n)}],$$ where $W(x)$ is the Lambert W function (see also Wolfram Alpha)

It should be borne in mind that in the model under consideration, prime divisors are computed approximately, and accounting for the powers of prime divisors is bounded by the cases $p_k = 4,\ 9.$ This could reduce the final result.

Will be assumed that $$\Pr(p\mid m) = {1\over p}\Pr(p\le m) = {n-p\over np},\quad\text{if}\quad m\in[1,n]\cap\mathcal N.$$

Calculations for the fixed $\mathbf n$

Let $S(p_k)$ is event "parameters $p,a,b,c$ allow the solution", then the opposite event $\overline{S(p_k)}$ is

$$\overline{S(p_k)} = (p_k\mid a) \wedge (p_k\mid b) \wedge (p_k\nmid c),$$ where the events $(p_k\mid a),\ (p_k\mid b),\ (p_k\nmid c)$ are independent and $$\Pr(p\mid a) = \Pr(p\mid b) = \Pr(p\mid c) = {n-p\over np},\quad \Pr(p\nmid c) = 1-{n-p_k\over np_k},$$ so $$\Pr(\overline S(p_k)) = \left({n-p_k\over np_k}\right)^2\left(1-{n-p_k\over np_k}\right),$$ $$\Pr(S(p_k)) = 1 - \left({n-p_k\over np_k}\right)^2\left(1-{n-p_k\over np_k}\right).\tag3$$

If $k=0\dots5$, then
$$p_k << n\rightarrow \Pr(S(p_k)) = 1-{p_k-1\over p_k^3},$$

The probability that the Diophantine equation has a solution will be asserted as the product of the corresponding probabilities for each of the dividers $p_k$. For example, $$P_{(2\dots9)} = \Pr(S(2)\wedge S(3)\wedge S(4)\wedge S(5)\wedge S(7)\wedge S(9)) = {7\over8} \cdot{25\over27} \cdot{61\over64} \cdot{121\over125} \cdot{337\over343} \cdot{721\over729},$$ $$P_{(2\dots9)}\approx0.726362.$$ Then for $n>>10$ $$P_n \approx \Pr(\bigwedge\limits_{k=0\dots k_{max}(n)} S) = P_{(2\dots9)} \prod\limits_{k=6}^{k_{max}(n)} \Pr(S(k\log(k))),$$ $$P_n\approx 0.726362e^{R(n)},\tag4$$ $$R(n) = \sum\limits_{k=6}^{[e^{W(n)}]} \left.\log\left(1 - \left({n-p_k\over np_k}\right)^2\left(1-{n-p_k\over np_k}\right)\right)\right|_{p_k=k\log(k)} \tag5.$$ For example, $$R(100) \approx -0.017900$$ (see also exp W(100) and R(100)), $$P_{100} \approx 0.726362e^{-0.017900} = 0.713476.$$ For $n>100$ $$P_n\approx 0.713476e^{Q(n)},\tag6$$ $$Q_n = \sum\limits_{k=30}^{[e^{W(n)}]} \lim\limits_{p_k\to k\log(k)}\log\left(1 - \left({n-p_k\over np_k}\right)^2\left(1-{n-p_k\over np_k}\right)\right),$$ $$Q_n \approx -\sum\limits_{k=30}^{[e^{W(n)}]} \lim\limits_{p_k\to k\log(k)}\left({n-p_k\over np_k}\right)^2\left(1-{n-p_k\over np_k}\right),$$ For example, $$Q(100000) \approx -0.001939$$ (see also exp W(100000) and Q(100000)), $$P_{100000} \approx 0.713476e^{-0.001939} = 0.712094.$$

This approach gives $$p\approx 0.712$$

On the other hand, direct multiplication of probabilities for factors $$p\in\{2,3,4,5,7,8,9,11,13,16,17,19,23,25,27,29,31,32,37,41,47\}$$ gives $$\boxed{p<0.641796}.$$

Conclusion

Comparison of the OP estimates with the results of accurate modeling shows that the evaluation of the required probability should be substantially reduced.

In the same time, the stability of the model parameters is achieved due to the identity of the laws of the distribution of the parameters $a,b,c.$

The proposed model can be used for the estimations based on the another distribution laws.