This is a slight modification of Euclid’s factoring algorithm:

Given a positive integer a, we can find all pairs of positive integers b and c such that (a-1)*(a+1)=b²*c.

If c is a square number, then essentially you lose sight of the root of the equation and you need to do a GCD until it can be re-expressed in your form, such that GCD( (a-1) MOD b²*c ) or GCD( (a+1) MOD b²*c )is a prime factor of b²*c. See also Shor’s algorithm.

Solving for a and b, or for a and c, is diophantine and NP-hard as it reduces to factoring. This might have some relation as to why when c is a square, there is no remainder 1 solution. Remeber, the Shor algorithm solves for the period r, of the function A**x MOD N =1 such that gcd(a**(r/2)+/-1) MOD N is a factor of N. Same idea, essentially.

Other related concepts might have to do with the separation of twin primes, but I don’t know enough about that to say much.

PS Also, this kind of function you present is known as an elliptical curve to those with a pen paper, two points and a string.