Project Euler #127 abc-hits

Question

Problem (Rephrased from here):

• The radical of \$n\$, \$rad(n)\$, is the product of distinct prime factors of \$n\$. For example, \$504 = 2^3 × 3^2 × 7\$, so \$rad(504) = 2 × 3 × 7 = 42\$.

We shall define the triplet of positive integers \$(a, b, c)\$ to be an abc-hit if:

• \$GCD(a, b) = GCD(a, c) = GCD(b, c) = 1\$
• \$a < b\$
• \$a + b = c\$
• \$rad(abc) < c\$

If \$(a,b,c)\$ is an abc-hit such that \$c<120000\$, Find sum of all \$c\$.

I have this code which runs in ~25.727 seconds, which I would like to write in a better way.

Firstly I noted these:

• If GCD(a,b)=1 then GCD(a,a+b)=1 and GCD(b,a+b)=1
• rad(a*b*c)=rad(a)*rad(b)*rad(c) if GCD(a,b,c)=1

Code:

private static double[] rad;
private static int[] primes;

public static void main(String[] args) {
    long t1 = System.currentTimeMillis();
    System.out.println(solve());
    long t2 = System.currentTimeMillis();
    System.out.println((double) (t2 - t1) / 1000.0);
}

private static long solve() {
    final int lim = 1000;
    long sumc = 0;
    sievePrimesArray(lim);
    sieveRadArray(lim);
    for (int a = 1; a < lim; a++) {
        // The value of `b` is minimum `a+1`, so `notcoprime` is an
        // array which tells if `i` and `a` have some common prime factor.

        boolean[] notcoprime = new boolean[lim - a - 1];
        for (int i = 0; i < primes.length && primes[i] <= a; i++) {
            if (a % primes[i] == 0) {
                // We need `primes[i]*j>=a+1`

                for (int j = (a + 1) / primes[i] + (((a + 1) / primes[i] < a + 1) ? 1 : 0);
                     primes[i] * j < lim; j++) {
                    notcoprime[primes[i] * j - (a + 1)] = true;
                }
            }
        }
        for (int b = a + 1; b + a < lim; b++) {
            int c = a + b;
            if (!notcoprime[b - (a + 1)]) {

            // Here to avoid overflow I used `double` and since `BigInteger` was slow.
            // The maximum value of `rad` rad reached is `119999` and
            // `rad[a-1]*rad[b-1]<=119999*119999=14399760001` which is still in `double`

                if (rad[a-1] * rad[b-1] < (double) c / rad[c-1]) {
                    sumc += c;
                }
            }
        }
    }
    return sumc;
}

Rather than finding all prime factors of a number, I multiply rad[x*p] by p where p is prime, once.

public static void sieveRadArray(int lim) {
    rad = new double[lim];
    Arrays.fill(rad, 1);
    for (int x : primes) {
        for (int j = 1; x * j < lim; j++) {
            rad[j * x - 1] *= x;
        }
    }
}

private static void sievePrimesArray(int lim) {
    boolean[] isPrime = new boolean[lim + 1];
    Arrays.fill(isPrime, true);
    for (int i = 2; i * i <= lim; i++) {
        if (isPrime[i]) {
            for (int j = i; i * j <= lim; j++) {
                isPrime[i * j] = false;
            }
        }
    }
    ArrayList<Integer> prime_list = new ArrayList<Integer>();
    prime_list.add(2);
    for (int j = 3; j < lim; j++) {
        if (isPrime[j]) {
            prime_list.add(j);
        }
    }
    primes = prime_list.stream().mapToInt(i -> i).toArray();
}

Answer 1

• Kudos for doing math groundwork, particularly for noticing multiplicative property of \$rad\$.

• I see no reason to resort to double. 14399760001 is well within the long range.

• Computing notcoprime costs a lot. Replacing it with

  for (int a = 1; a < lim; a++) {
      for (int b = a + 1; b + a < lim; b++) {
          if(gcd(a, b) == 1) {
              ....
  

  reduces the execution time by the factor of 1.5 even for the most naive gcd implementation.

• On my system your code as is runs for about 0.06 sec. The reported 25+ sec must have some other explanation. What system do you have?