Counting number of integer solutions of z³ = x²y² for given z

Question

I have been working on this problem for some time now. The problem asks to find two numbers x and y such that the product of the squares of x and y equal the cube of k. My approach was to solve for x, given that an input of 1 would give a number, lim, that when squared is equal to k cubed. The code first tests if the square root of the function arg is an integer or not. It not then 0 is returned. If not, then the lim is found. The code then iterates through a range of numbers, i, starting at 1 up to the square root of lim, testing to see if the quotient of lim and i is an integer or not. Since the problem asks for the number of divisors, I found it best to use reduce and count the number of successful x and y pairs, rather than iterate through a list to find its size.

def c(k: Long) = {
  val sqrtk = math.sqrt(k)

  sqrtk.isWhole match {
    case false => 0
    case true =>
      val lim = k * sqrtk

      (1 to math.sqrt(lim).toInt).foldLeft(0)((acc, i) =>
        if ((lim / i).isWhole) acc + 2 else acc
      )
  }
}

I am pretty sure that the problem lies in the use of the two square roots and the division being performed, but I have so far been unable to make this code perform any better. Initially I was unable to even pass the basic tests and now I can pass about 48 of the random tests, and am hopefully close to passing all the tests. I have tried numerous other solutions, including making my own square root function, and trying to find all divisors with prime numbers, etc., but none have been as fast as the code I have posted above. Like the title says, this code is written in Scala.

Thanks in advance for any assistance that can be offered.

NOTE:

• The combinations (x^2, y^2) must be whole/integer values.
• The fixed tests for the problem have 33 k values, the largest of them is 10000000095. These tests are followed by so far up to 48 random k values. The max time allowed to finish all k values is 16 seconds, and times out before finishing all the tests.

Answer 1

Here is a straightforward approach:

• z * z * z = x * x * y * y means that z * sqrt(z) = x * y
• If sqrt(z) is not an integer, it is irrational, therefore there are 0 solutions
• Otherwise, let q = sqrt(z) be the integer square root.
• Factorize q = f1^e1 * ... fN^eN with prime factors f_i and positive integer exponents e_i.
• Then x * y = f1^(3 * e1) * ... * fN^(3 * eN)
• There are thus (3 * e1 + 1) * ... * (3 * eN + 1) essentially different ways to partition the factors between x and y

So, assuming some kind of factorization algorithm factorize, that would be:

def c(k: Long): Long = {
  val q = math.sqrt(k).toLong
  if (q * q == k) {
    val fs = factorize(q)
    fs.map(e => 3 * e._2 + 1).product
  } else {
    0
  }
}

---

Full code with an overly simplistic factorize:

def c(k: Long): Long = {
  val q = math.sqrt(k).toLong
  if (q * q == k) {
    val fs = factorize(q)
    fs.map(e => 3 * e._2 + 1).product
  } else {
    0
  }
}

@main def main(): Unit = {
  println(c(1))
  println(c(4))
  println(c(4096576))
  println(c(2019))
}

def factorize(i: Long): List[(Long, Int)] = {
  var odd = i
  var twoExp = 0
  while (odd % 2 == 0) {
    odd /= 2
    twoExp += 1
  }

  def factorizeRec(n: Long, divisor: Long, limit: Long): List[(Long, Int)] = {
    require(n >= 3)
    var unfactoredPart = n
    var e = 0
    while (unfactoredPart % divisor == 0) {
      unfactoredPart /= divisor
      e += 1
    }
    if (unfactoredPart == 1) {
      if (e > 0) {
        List((divisor, e))
      } else {
        Nil
      }
    } else {
      if (e > 0) {
        val newLimit = (math.sqrt(unfactoredPart).toInt + 1)
        (divisor, e) :: factorizeRec(unfactoredPart, divisor + 2, newLimit)
      } else {
        factorizeRec(unfactoredPart, divisor + 2, limit)
      }
    }
  }

  val sqrtLim = (math.sqrt(i).toInt + 1)
  val unevenFactors = 
    if (odd >= 3) then factorizeRec(odd, 3, sqrtLim) 
    else Nil

  if (twoExp > 0) {
    (2, twoExp) :: unevenFactors
  } else {
    unevenFactors
  }
} 

which prints

1
4
160
0

as expected. I didn't test whether it's fast enough, didn't bother to register.

If that's still too slow, steal Pollard's Rho somewhere.

Answer 2

As I mentioned in comment, any time you are comparing sqrt() to a whole number, you are probably off track, even if the logic is correct.

There are some clues in the math of this problem. First, everything is multiplied, there are no additions anywhere. What is the implication? The prime factors of the LHS must be the factors of the RHS. Further, we know that the count of each of the factors must be identical.

So then you get concerned about the prime factors of the limiting value of 10B, which is a large number of primes. And we can observe that LHS is a cubic, so we know any prime must occur 3 times in the factorization of z^3. The cube root of 10B is still large, but could be iterated. But there is a nuanced relationship here between the powers of x, y, z. z occurs as a cubic and both x and y as squares, so it is impossible for the occurrences of factor a to occur just 3 times. In fact, it must occur 6 times (LCM of 2,3) if it occurs at all.

So, we need to consider only the prime factors up to 10B^(1/6) which is ~ 68, or about 20 factors.

This works, is written in python, but should be pretty easy to map to scala. It tests 100K random trials in less ~550ms

Code:

#prime factorize

# solving z^3 = y^2 * x^2

from random import randint


def test(k):
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, ]
    if k==1 : return 1, {}
    # factor k^3
    k = k**3
    factors = {}  # a dictionary mapping of prime:occurrences ...
    for p in primes:
        while k % (p**6) == 0:
            factors[p] = factors.get(p,0) + 1
            k /= p**6
        if k == 1: break
    if k != 1: return 0, {}

    # now loop through the factors & counts to construct possible pairs.
    # for each occurrence of each factor, we can put 0, 1, 2, ..., count*3 occurrences in "x"
    # so if there is 1 occurrence, we can put 0, 1, 2, 3 multiples of the factor in the first digit
    # or, more generally, 3 * count + 1
    if len(factors) == 0:
        return 0, {}
    tot = 1
    for v in factors.values():
        tot *= 3*v+1
    return tot, factors


for k in (1, 4, 16, 2019, 4096576, 10_000_000_000):
        print(k, test(k))

# test over randomized list
vals = [randint(1,10_000_000_000) for _ in range(100_000)]
for val in vals:
    result = test(val)
    if result[0]:
        print(val, result)

Output:

1 (1, {})
4 (4, {2: 1})
16 (7, {2: 2})
2019 (0, {})
4096576 (160, {2: 3, 11: 1, 23: 1})
10000000000 (0, {})
[Finished in 563ms]

Answer 3

SOLUTION: The insight I needed to solve this was the relationship between z³ and the factors of x²y². In this case, the factors for the square root of k, combined with the expression (3 * e1 + 1) * ... * (3 * eN + 1), gives the solution to the problem.

def c(k: Long) = {
  val sqrtk = math.sqrt(k)

  @annotation.tailrec
  def factorList(x: Long, n: Long = 2, ls: List[Long] = Nil): List[Long] = {
    n * n > x match {
      case false if x % n == 0 => factorList(x / n, n, n :: ls)
      case false               => factorList(x, n + 1, ls)
      case true                => x :: ls
    }
  }

  sqrtk.isWhole match {
    case false => 0
    case true =>
      factorList(sqrtk.toLong)
        .groupBy(identity)
        .map { case (_, v) => 3 * v.size + 1 }
        .product
  }
}