Let two numbers $x$ and $y$ be selected from the set of first $n$ natural numbers with replacement(i.e. the two numbers can be same).The question is to find out the probability that $x^2-y^2$ is divisible by $k\in \mathbb{N}$

For $k=2$

Any number can be expressed as $2p,2p+1$.Now $x^2-y^2=(x-y)(x+y)$ If both numbers are of form $2p+1$ then (x-y) would be divisible by $2$ .if two numbers are of different forms then it will not be divisible by $2$.So the probability in this case is $a^2+(1-a)^2$ where $a$ is probability that number chosen is divisible by $2$ which is $\frac{\lfloor \frac{n}{2} \rfloor}{n}$.However this gets complicated with $k=3$ onwards because numbers in different forms may be divisible.In other words if there a generalisation or way to solve for some large $k$.Thanks.

  • 1
    $\begingroup$ $k=3$ is actually as easy as $k=2$, since any square is either $3p$ or $3p+1$. Looking at the difference of the two squares, we see that it is divisible by $3$ iff $x$ and $y$ are either both divisible by $3$ or neither of them are. $\endgroup$
    – Arthur
    Commented Apr 10, 2017 at 7:15
  • $\begingroup$ @navinstudent You might be interested in the complete analytic solution of your problem provided in my solution. $\endgroup$ Commented Apr 13, 2017 at 17:12

2 Answers 2


A generalization expressed by a set

A good way to generalize this is to use modular arithmetic, or essentially look at the remainder of $\frac{x}{k}$ and $\frac{y}{k}$. As you pointed out in your example for $k=2$, the numbers can only be expressed as $$2*p,2*p+1$$

This can be further generalized into a set of unique expressions that define every number or every $x$ and $y$ for a value $k>1$ and $p≥0$$$S={kp,kp+1...kp+(k-2),kp+(k-1)}$$

Using modular arithmetic, we know that $S\equiv \{0,1,2...k-2,k-1\}\pmod k$

Since $x$ and $y$ are any two terms from the set $S$ for any $p≥0$, we know that $x$ and $y$ are really just equal to either any value from $\{0,1,2...k-2,k-1\}$ since we only need to look at the remainder when divided by $k$ to determine divisibility.

Probability using mod

Now we know $x^2-y^2 = (x-y)(x+y)$ meaning either $(x-y)$ or $(x+y)$ have to be divisible by $k$ in order to fulfill the requirement that $x^2-y^2$ is divisible by $k$.

It should be noted that the probability of choosing a random number being expressed by any expression in set $S$ is uniform and is equal to $\frac{1}{k}$. For example, the probability of choosing a number is expressible by $3p+1$ is $\frac{1}{3}$. This can be proved by showing how every number expressed by $kp$ can always be paired with $kp+1,kp+2,...kp+(k-2),kp+(k-1)$, thus showing how there are an equal number of numbers of each type.

Onto the question....

$(x-y)$ is divisible by $k$ if $(x-y)\equiv 0\pmod k$. This simplifies to $(x$ mod $k)-(y \text{ mod}\ k) = 0$. Thus we now need to find the probability that $x$ and $y$ are both expressible by the same expression in the set S.

Now this comes down to a simple problem of asking what is the probability of choosing two of the same elements from the set $S$ without replacement, which is $\frac{k}{k} *\frac{1}{k} =$ $\frac{k}{k^2}$

$(x+y)$ is divisible by $k$ if $(x+y)\equiv 0\pmod k$. This simplifies to $(x$ mod $k)+(y \text{ mod}\ k) = 0$. Thus we now need to find the probability that $x$ and $y$ are chosen such that the $x+y =$ a multiple of $k$.

Knowing that $x$ and $y$ must come from the set $S$, we can see there is a specific pairing of elements from $S$ that causes the addition of the pair of elements to be equal to $2*kp+k$, a multiple of $k$. The only times that the pair of elements, which will be the expression that expresses $x$ and $y$, add up to $2*kp$ is when one is expressed by $kp+m$ and the other being $kp+(k-m)$ for any whole number $m$. This is because $(kp+m)+(kp+(k-m)) = 2*kp$. For example $3p+1$ can be paired up with $3p+2$ since their sum is $6p+3$, a multiple of 3.

Now we are now looking for the probability the two chosen elements from $S$ are pairs of each other (If both chosen elements are $kp$, it will be a multiple of k still) The total number of pairs that can be chosen is the total number of elements in set $S$ squared, which is $k^2$. The total number of pairs that fulfill the divisibility by $k$ is equal to $\lfloor\frac{k}{2}\rfloor + 1$. The ceiling function is applied to correct errors that occur when $k$ is odd as you obviously can't have a fractional number of pairs. The extra $+1$ is for including the case when both elements are the same.

Thus the probability that $(x+y)$ is divisible by $k$ is $\frac{\lfloor\frac{k}{2}\rfloor + 1}{k^2}$.

However now we encounter a problem of overlap, which are the cases when the two elements chosen from $S$ make $(x-y)$ and $(x+y)$ divisible by k. To get rid of the over count we need to subtract the number of times this happens for a given set $S$. We know we need to subtract once for the case when both elements are $kp$. We also need to subtract another one depending if $k$ is even or not. This is because when $k$ is even, we have the case where the element $kp+\frac{k}{2}$ can pair with itself to fulfill everything. Thus the final answer is $$\frac{k}{k^2}+\frac{\lfloor\frac{k}{2}\rfloor + 1}{k^2} - \frac{(k+1 \text{ mod }2)+1}{k^2}$$

The $\frac{(k+1 \text{ mod }2)+1}{k^2}$ tells us to subtract one more and subtract another if $k$ is even, hence we shifted the modular part by 1 since $( \text{ even number mod }2 = 0)$

  • $\begingroup$ Thanks for your answer.I think you have approached the problem by taking that the two numbers are chosen from the infinite set of natural numbers but how will the answer depend on n if the two numbers are chosen from the set of first n natural numbers with replacement. Before that would you please explain more on how you applied the ceiling function to get the probability that x+y is divisible by k.Thanks. $\endgroup$
    – Navin
    Commented Apr 10, 2017 at 7:47
  • $\begingroup$ The number of pairings that work is always equal to k divided by 2. However, there is a special case when k is odd. We can just check a few explicit exmaples to prove whether or not to use ceiling or floor. If k is 3, we can pair 0,0 and 1,2. If k is 4, we pair 0,0 1,3 2,2. So in fact I wrote it wrong, it should be the floor function as 3/2 floored is 1 and adding the 0,0 pair is 2. For k=4, the floor of k/2 is 2, and adding the 0,0 pair is 3. Sorry for the mistake, will change. I will also work on your other part of the question $\endgroup$
    – Stone
    Commented Apr 10, 2017 at 8:03
  • 2
    $\begingroup$ Suppose $k = p \cdot q$. $x^2 - y^2$ is also divisible by $k$ if $(x - y)$ is divisible by $p$ and $(x + y)$ is divisible by $q$. So it is not necessary that $( x- y)$ is divisible by $k$ or $(x + y)$ is divisible by $k$. $\endgroup$
    – PSPACEhard
    Commented Apr 10, 2017 at 8:10
  • $\begingroup$ If you can find an example where such a case appears and either x-y doesn't be equal 0 or x+y doesn't equal a multiple of k, then your point would be valid. However I haven't been able to find such an example $\endgroup$
    – Stone
    Commented Apr 10, 2017 at 8:43
  • $\begingroup$ Suppose $n = 8, k = 8$ and $x = 3$, $y = 1$. Would this be a valid case? $\endgroup$
    – PSPACEhard
    Commented Apr 10, 2017 at 12:42

EDIT 12.04.17

After some days of studying this interesting problem, I provide now a complete analytic solution of the problem by retrieving a related problem in OEIS.

I give exact values of the first few probabilities, and give an exact general formula for the probababily if the divisor $n$ is prime.

A Monte-Carlo-Simulation is also presented.

Notice that my findings came up in opposite temporal order.

Analytic solution to the problem

Looking up the sequence of the first few terms of

$$s(n) = p(n) n^2 \tag{1}$$

which are

$${1, 2, 5, 8, 9, 10, 13, 24, 21, 18, 21, 40, 25, 26, 45}$$

in the online-encyclopedia of integer sequences https://oeis.org/ gives us the entry A062803 "Number of solutions to $x^2 == y^2 mod(n)$", created initially by Ahmed Fares (ahmedfares(AT)my-deja.com), Jul 19 2001. A formula was devised by Vladeta Jovovic on Sep 22, 2003: "Multiplicative with a(2^e)=e*2^e and a(p^e)=((p-1)*e+p)*p^(e-1) for an odd prime p."

Let us explore this a little bit closer. Our problem is transformed into the question of the number $s(n)$ of solutions to the congruence

$$x^2-y^2 == 0, \; mod(n)\tag{1}$$

Our probabilities are then given by

$$p(n) = s(n)/n^2\tag{2}$$

Suppose the number $n$ has the representation

$$n = \prod_{i=1}^{k} p_i^{a_i}\tag{3}$$

where $p_i$ is the $i$-th prime number appearing in $n$ in ascending order, $a_i$ ist multiplicity (exponent) and $k$ is the number of different prime factors in $n$. Notice that in number theory $k$ is traditionally called $\omega(n)$, the number of distinct prime factors. It is implemented in Mathematica as PrimeNu[n].

Then the statement of being multiplicative means that we can apply the formula to each of the prime power factors separately, and multiply the result together.

This gives for $n$ even

$$s(n_{even}) = (a_1 2^{a_1}) \prod_{i=2}^{k} ((p_i-1)a_i+p_i) p_i^{a_i-1}\tag{4.1}$$

and for $n$ odd

$$s(n_{odd}) = \prod_{i=1}^{k} ((p_i-1)a_i+p_i) p_i^{a_i-1}\tag{4.2}$$

It is easy to see that for an odd prime $n$, $s(n) = 2n-1$, as claimed earlier.

The formula simplifies if $n$ is square-free. Then all non-vanishing $a_i$ are equal to $1$, and we find

$$s(n_{even}) = 2 \prod_{i=2}^{k} (2 p_i-1)\tag{5.1}$$ $$s(n_{odd}) = \prod_{i=1}^{k} (2 p_i-1)\tag{5.2}$$

For those who are interested in encoding this formula, here is an example in Mathematica

s[n_] := Module[{fi, pi, ai, pout},   
  fi = FactorInteger[n];  
  pi = #[[1]] & /@ fi;  
  ai = #[[2]] & /@ fi;  
  pout = If[OddQ[n],   
    Product[((pi[[i]] - 1) ai[[i]] + pi[[i]]) pi[[i]]^(ai[[i]] - 1), {i, 1, Length[pi]}],    
    ai[[1]] 2^ai[[1]] Product[((pi[[i]] - 1) ai[[i]] + pi[[i]])
       pi[[i]]^(ai[[i]] - 1), {i, 2, Length[pi]}]]]  

The prime factor decomposition of $n$ is done by the function FactorInteger[]. From this we extract the $p_i$ and $a_i$, and then apply the formula of Jovovic.

Exact values of the probabilities

We make the (natural) assumption that all possible remainders of a randomly chosen number $x$ modulo $n$ have equal probability.

Then we can calculate the exact values of the probabilities with the following piece of code (written here in Mathematica)

h[n_] := 1/n^2 Count[Flatten[Table[Mod[x^2 - y^2, n], {x, 0, n - 1}, {y, 0, n - 1}]], 0]


For a given divisor $n$ the expression $z = x^2-y^2$ needs to be considered only for $x$ and $y$ bewteen $0$ and $n-1$. The Table[] lists all elements $x^2-y^2 mod(n)$, and Flatten[] puts them in one array. Now Count[.,0] counts the zeroes in this array. Dividing this by $n^2$ gives the probability.

The result for $n = 1..30$ in the format $\{n,p(n)\}$ are

$$h(n)_{tab} = \left( \begin{array}{ccccc} \{1,1\} & \left\{2,\frac{1}{2}\right\} & \left\{3,\frac{5}{9}\right\} & \left\{4,\frac{1}{2}\right\} & \left\{5,\frac{9}{25}\right\} \\ \left\{6,\frac{5}{18}\right\} & \left\{7,\frac{13}{49}\right\} & \left\{8,\frac{3}{8}\right\} & \left\{9,\frac{7}{27}\right\} & \left\{10,\frac{9}{50}\right\} \\ \left\{11,\frac{21}{121}\right\} & \left\{12,\frac{5}{18}\right\} & \left\{13,\frac{25}{169}\right\} & \left\{14,\frac{13}{98}\right\} & \left\{15,\frac{1}{5}\right\} \\ \left\{16,\frac{1}{4}\right\} & \left\{17,\frac{33}{289}\right\} & \left\{18,\frac{7}{54}\right\} & \left\{19,\frac{37}{361}\right\} & \left\{20,\frac{9}{50}\right\} \\ \left\{21,\frac{65}{441}\right\} & \left\{22,\frac{21}{242}\right\} & \left\{23,\frac{45}{529}\right\} & \left\{24,\frac{5}{24}\right\} & \left\{25,\frac{13}{125}\right\} \\ \left\{26,\frac{25}{338}\right\} & \left\{27,\frac{1}{9}\right\} & \left\{28,\frac{13}{98}\right\} & \left\{29,\frac{57}{841}\right\} & \left\{30,\frac{1}{10}\right\} \\ \end{array} \right)$$

The simulation results (see below) are in reasonable agreement with these results.

If $n$ is an odd prime number the probability is given by

$$p(n)=\frac{2 n-1}{n^2}$$

If $n$ is composite I have not found the exact formula. The problem here seems to be related to quadratic residues.


EDIT 11.04.17

We distinguish two possible basic sets of integers from which to select for the divisibility test:

  1. Set with repetition

We create a set $m$ consisting of all numbers $z = x^2-y^2$ of integers where $1<= x <= n_{max}, 1<= y <= n_{max}$.

  1. Set without repetition

The set $m_0$ is obtained by removing from $m$ all duplicates.

Then, for a given divisor $n$ we estimate the probability of divisibility by the ratio of the number of elements of the set for which $\frac{z}{n}$ is integer relative to all elements of the set.

The resulting probabilities for $m_{max} = 10^3$ and divisors in the range $n = 1..30$ are

Case 1 (with repetition)

List of results in the format $(n, p(n))$

$$((1, 1.0), (2, 0.5000005), (3, 0.55533378), (4, 0.5000005), (5, \ 0.35968128), (6, 0.27766739), (7, 0.26530612), (8, 0.37475112), (9, \ 0.25933407), (10, 0.17984114), (11, 0.17355372), (12, 0.27766739), \ (13, 0.14792899), (14, 0.13265356), (15, 0.19973533), (16, \ 0.25000075), (17, 0.11417454), (18, 0.12966754), (19, 0.10246197), \ (20, 0.17984114), (21, 0.14736213), (22, 0.08677736), (23, \ 0.08502487), (24, 0.20808462), (25, 0.10419251), (26, 0.073965), (27, \ 0.11103881), (28, 0.13265356), (29, 0.06774345), (30, 0.09986916))$$

The graph

enter image description here

Case 2 (no repetition)

List of results in the format $(n, p(n))$

$$((1, 0.5125414), (2, 0.19847984), (3, 0.22136804), (4, 0.19847984), \ (5, 0.14290505), (6, 0.08393305), (7, 0.10653981), (8, 0.11713062), \ (9, 0.08960969), (10, 0.05436621), (11, 0.07129035), (12, \ 0.08393305), (13, 0.06139016), (14, 0.04070156), (15, 0.05862769), \ (16, 0.06700892), (17, 0.04836223), (18, 0.03377941), (19, \ 0.04369956), (20, 0.05436621), (21, 0.04403089), (22, 0.02740017), \ (23, 0.03676344), (24, 0.04804486), (25, 0.03682731), (26, \ 0.0236357), (27, 0.03531034), (28, 0.04070156), (29, 0.02972352), \ (30, 0.02204489))$$

The graph

enter image description here

Notice the remarkable peridocity with a period of $4$ in the "artificial" case 2. I'm sure there is a simple explanation for this, and that some readers can give it.

  • $\begingroup$ However, your answer for probability when dividing by 2 is different from mine. I got 1/4 + 1/4 = 1/2 whereas you got 0.2. The two quarters represent the probability of choosing x and y such that they are congruent in terms of modulo 2. $\endgroup$
    – Stone
    Commented Apr 11, 2017 at 0:13
  • $\begingroup$ @Stone you are right. This difference is due to my perhaps artificial removal of double occurences in the set m. I have added the "pure case" to my answer $\endgroup$ Commented Apr 11, 2017 at 8:41

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .