Calculate number of pairs a,b where a² - b² = n

Question

In one site there was one question to make one program where you will be given a positive integer n and we have to calculate the number of such pairs (a, b), where n = a² - b² and both a and b are positive integers.

for example 15 = 4² - 1² = 8² - 7².

Currently this question is closed and in my case I am still not able to do it properly. Still I won't be trying to find solution here but I just want to know that where I am doing wrong with my code.

At first I thought of simply increment a, b by 1 and check every equation. But time limit is 4 seconds and 1 ≤ n ≤ 10¹³ so it was taking a very long time for larger values of n.

Since a² - b² = (a+b)(a-b) then I took one example 15 and solved for a = n and b = a-1 and observed that I got some special series:

(2n+1): n=1,2,3 ...   (for a = 2 and b = 1)

4n: n=2,3,4 ...       (for a = 3 and b = 1/b = 2)

6n+3: n = 2,3,4 ...   (for a = 4 and b = 1/b = 2/b = 3)

8n: n = 3,4,5 ...

Later I used this concept and made one program. It is working fine with 6-digit numbers(approx).

My program:

n = input()
count = 0
for i in range(1,n):
    first = (n-((2*i)-1))
    second = ((2*i)+2*(i-1))
    if first % second == 0 and first/second >= i:
        count += 1
    elif n%(4*i) == 0 and n/(4*i) >= (i+1):
        count +=1
print count

But for some larger values it is giving some error(16/30). If I make n in range to 100000 then I get few error as compared to previous one (27/30). But still it is taking time for larger input. So does my method is right or do I have to think of something else. If this is wrong then why ?

Explanation

Let us say I want to find for no. 15 then according to formula (a+b)(a-b) = n total no. of possibilities :

(2+1)(2-1) = 3 # BLOCK 1

---

(3+1)(3-1) = 8 #BLOCK 2

(3+2)(3-2) = 5

---

(4+1)(4-1) = 15

(4+2)(4-2) = 12 #BLOCK 3

(4+3)(4-3) = 7

---

And so on ... till (8+7)(8-7) = 15

So if I go on like this then the last line of each block (in this case it is n = 3,5,7) we found that it is changing as (2n+1) where n = 1,2,3 .... (here we can think n as if from which block we are starting the count)

Now for every last second line (in this case it is 8,12) the equation will be 4n where n = 2,3,4 ....(since first block has only 1 line so there can't be any last second line so we will count from second block)

Repeat this process again and again .... We get:

2n+1, 4n, 6n+3, 8n, 10n+5, 12n, 14n+7, 16n ... and so on

If we simply try to find these equation every time for different lines then we find out that we can use this concept in program. If any no. 'n' is divisible by this equation it means it has occurred on those lines of block. Means each time it is divisible it will increment the count.

Answer 1

You were correct to try a^2 - b^2 = (a+b)(a-b) venue: take your n and break it into p and q so that p > q and p * q = n.

Usually there are several ways of doing it and almost each way defines a pair you need. a would be (p + q)/2 and b would be (p - q)/2. Thus, you only need pairs where either both p and q are even, or they are both odd.

You don't really need calculating a and b, or p and q, just number of pairs.

Thus, break your n into prime divisors and count distinct ways to divide those divisors into 2 different sets so that either both sets contain at least one 2, or neither do.

Edit:

Ok, some more math:

Case with n=1 is trivial: no solutions. For all other cases we can assume that n has nonempty set of prime divisors.

Let's start with case when all prime divisors are distinct and 2 is not present.

Each divisor can be placed either in set 1 or in set 2. Since sets will be obviously different, we're guaranteed p > q and since we're working with n's divisors, we can be sure that p * q = n.

Divisors' placement can be described bijectively by a string of bits with the length as number of divisors. Hence, we'll have 2^n placements, 2^(<number of divisors>-1) pairs of p and q ( -1 because swapping sets gives same p and q).

Now what if some prime divisor (except for 2, which is still not present) is present k>1 times?

For all other divisors, situation is exactly as above, but our special divisor will provide *(k+1), instead of *2 to number of pairs. Well, actually, *2 of other divisors is just a special case with k = 1.

So, if 2 is not present, we go through n's prime factorisation and for each divisor and its respective k number of occurences, we multiply our number of placements by k+1. If result is even - there were odd ks, meaning that in no placement our 2 sets were equal. We still counted every pair of p and q twice (once as set1, set2 and once as set2, set1) so we divide k+1 product by 2 and that's our answer. If result is not even - we have even quantity of every divisor, meaning that n is a square. That means we counted twice every pair of p and q but we've also counted pair of sqrt(n) once - which we shouldn't have! We would have to substract 1 before dividing by 2.

If 2 is present, things don't get much complicated either. If there's only one 2 we return 0 - because there's no a and b that work. If there's more than one - we say that one goes to each set every time so we simply drop 2 twos and stop concerning ourselves with them, repeating steps above.

And now for some code:

from collections import Counter
from functools import reduce
from operator import mul

def product(stuff): return reduce(mul, stuff)

def square_pairs_count(n):
    """Return number of pairs (a,b) so that a*a - b*b = n

    (Algo explanation goes here)
    """
    if n % 2 == 0 and n % 4 != 0:
        return 0 # (a+b) is odd and (a-b) is even (or vice versa)? No way!
    factorisation = Counter(prime_factorisation(n))
    if factorisation[2] > 1:
        # We use 2 '2's to guarantee that both (a+b) and (a-b) are even
        factorisation[2] -= 2
    factor_count_product = product(k + 1 for k in factorisation.values())

    if factor_count_product % 2 == 1:
        # n is a square and we counted root-root pair as (a+b) and (a-b), undo that
        factor_count_product -= 1

    return factor_count_product // 2 # we've counted every pair twice, so //2

print(square_pairs_count(155235236))

Of course, this code could be called comprehensible only in context of all math above. In the wild, most of the code would be a docstring explaining what happens.

Number of divisors grows pretty slowly with growth of n. It all boils down to how fast you produce them prime factors.

---

Your idea with "blocks":

You can notice that "block" number c contains a = c + 1 and all bs that would give n > 0 (but not necessarily integer).

If you shift block numbering so that (2+1)(2-1) is actually block 2 then you will have block number equal to a used in this block with block 1 being empty. (because no b and n would satisfy our needs).

You could iterate through each pair like:

for a in range(2, max_a(n)):
    for b in range(1, a):
        # test our (a-b)(a+b)

However, that means iterating a lot.

You could notice that each block can produce no more than 1 solution, so we could iterate through blocks instead. However, we'd have to check if a block really produces a solution:

for a in range(2, max_a(n)):
    b_squared = a*a - n
    # check if b_squared is actually a square

The trick will be in determining max_a(n).

Surely, if 2*a - 1 > n, then a^2 - (a-1)^2 > n and thus for any b we're interested in, a^2 - b^2 > n.That makes our max_a something like (n + 1)//2, or n//2 + 1 for sake of corner cases.

That's too much! we'll have to check around 10^13 blocks, and that's a lot. Well, checking couple hundred thousand primes is not a simple thing either, but it's way faster.

So, unless you can lower the boundary for a, the code would still take lots of time.

---

Modulo sieving

If you are interested in generative solutions, I'd recommend picking some number N, rewriting equation modulo N and looking for moduli N for a and b that are possible with given n. If, for example, N is around hundred buth there's only 10 possible modulo pairs that work - you would have to iterate only through 10 (a,b) pairs out of 10 000.

For example, if N = 3 and n mod 3 = 2 then a mod 3 = 0 and either b mod 3 = 1 or b mod 3 = 2. Of each 9 consecutive possible options for (a,b) only 2 remain. Of each 3 consecutive possible options for a only 1 remains.

For example, if N = 5 and n mod 5 = 1 then only a^2 mod 5 = 1, b^2 mod 5 = 0 and a^2 mod 5 = 0, b^2 mod 5 = 4. Of 25 possible options for (a,b) pair this leaves only 4. Of 5 possible options for a, it leaves only 3.

Moreover, if you have calculated allowed modulo pairs for some M and some N you can then calculate allowed modulo pairs for M*N and instead of iterating possible pairs with steps of size N (or M), go with M*N steps.

Combining 2 examples above, you'll get that you only need to check 3 out of 15 options for a.

Using these 2 tricks (modulo sieving and combining) you can iterate with ever-increasing steps without missing anything. However, problem does not ask for generating those pairs, only for counting their number, so that's a bit excessive.

If you find that option interesting, take a closer look at sieves and especially wheels part.

Answer 2

You were correct to look at: \$a^2-b^2=(a-b)(a+b)\$, but you only ended up using it for specific cases. I came up with a way of creating all \$(a, b)\$ pairs that a 7th grader might be able to discover.

Please forgive me for any spelling mistakes or faults in my answer, as this is my first ever answer.

Main

The first step is to factor this binomial: \$a^2-b^2=(a-b)(a+b)=n\$. Since \$a-b\$ (\$a>b\$) and \$a+b\$ are both natural numbers, \$a-b\$ and \$a+b\$ must be factors of \$n\$.

I'm gonna define \$a-b\$ to be \$k\$ and from that arises that \$a+b=\frac{n}{k}\$.

We can now derive a formula for \$b=a-k\$, since we just defined \$k=a-b\$.

Plugging that into \$a+b=\frac{n}{k}\$ yields the following result: $$2a-k=\frac{n}{k}$$, and rearranging we get: $$2a=\frac{n}{k}+k$$ From this we can determine that \$\frac{n}{k}+k\$ must be even, because \$a\$ is defined to be a natural number and any natural number multiplied by \$2\$ is even.

Note the case of having an odd number of factors. In that case we'll possibly encounter that \$k=\frac{n}{k}\$. We'll prove that this \$k\$ value will not yield a valid \$(a, b)\$ pair. Recall that we stumbled upon the formula \$a-b=k\$ and defined \$a+b=\frac{n}{k}\$. In this special case \$a+b=k\$. These statements cannot both be true, since that would imply that \$b=0\$, which is a contradiction to the fact that \$a\$ and \$b\$ are natural numbers. Therefore the special case \$k=\frac{n}{k}\$ will not yield a valid \$(a, b)\$ pair.

We just found that for a solution to exist there must a factor \$k\$ of \$n\$, such that: \$\frac{n}{k}+k\$ is even and \$\frac{n}{k}\neq k\$.

Let's prove that there is always a valid \$(a, b)\$ pair if the above statement is true. If \$\frac{n}{k}+k\$ is even, then \$a\$ will be a natural number, since \$\frac{n}{k}+k\$ is divisible by 2. Now because \$a\$ and \$k\$ are always natural numbers and \$a>k\$, \$b\$ must also be a natural number \$(b=a-k)\$. And now we have proven that all \$k\$ values, that satisfy the statement \$\frac{n}{k}+k\$ is even, except the special case of \$k=\frac{n}{k}\$, will always yield a valid \$(a, b)\$ pair.

Since \$k\$ and \$\frac{n}{k}\$ are factors of \$n\$, these will be the only possible \$(a, b)\$ pairs for \$n\$.

Wrapping up

I came up with a simple way of finding all valid \$(a,b)\$ pairs for a given \$n\$ value.

Firstly write down all the factors of \$n\$ in an ascending order. Then create pairs of factors and 'opposite factors' (a factor that multiplies with the other factor to give \$n\$). If the sum of the pair's factor and 'opposite factor' is even, then the factor (\$k\$), when substituted into \$2a=\frac{n}{k}+k\$, will yield a valid \$(a, b)\$ pair. If there are no such factor pairs, then there are no valid \$(a, b)\$ pairs.

For example, lets take \$n=36\$.

Firstly write down all its factors in an ascending order like so: $$1, 2, 3, 4, 6, 9, 12, 18, 36$$ Then find the factor pairs that satisfy the required equalities \$(2 \$ and \$ 18)\$ and plug the factor in to \$2a=\frac{n}{k}+k\$ to find \$a\$ (\$10\$) and then into \$b=a-k\$ to find \$b\$ (\$8\$). Plugging \$a\$ and \$b\$ we get the true statement: $$10^2-8^2=36$$

Code

To count the number of valid \$(a, b)\$ pairs, we generate all of the factors of \$n\$ and then check if it will qualify. If it does, then increase the count by \$1\$. Do this until you run out of factors.

This code is very unoptimized, but it gets the job done.

This will most likely not work for your specifications though.

    from math import sqrt, floor

    def differences(n):
        divisors = get_divisors(n)
        count = 0

        l = len(divisors)
        for i in range(0, floor(l/2)):
            if (divisors[i] + divisors[l-i-1]) % 2 == 0:
            count += 1

            a = int((n+divisors[i]*divisors[i])/(2*divisors[i]))
            b = a-divisors[i]
            print(("[" + str(divisors[i]) + "]").ljust(12) + ("y=" + str(a) + "^2-" + str(b) + "^2").ljust(35) + "| Result: " + str(a*a-b*b))

        print()
        return count

    def get_divisors(n):
        divisors = []
        endDiv = []

        # List to store half of the divisors
        for i in range(1, int(sqrt(n) + 1)):
            if (n % i == 0):
            # Check if divisors are equal
            if (n / i == i):
                divisors.append(i)
            else:
                # Otherwise print both
                divisors.append(i)
                endDiv.append(int(n / i))

        endDiv.reverse() 
        return divisors + endDiv

    N = 36
    print("N: " + str(N) + " | Count: " + str(differences(N)))