Problem 87 - D

With Project Euler being down, I thought it might be difficult to move on and do more recent problems, but then I realize that between myself and my friends, the next few most recent problems have all been solved by someone. However, on the day when it finally became relevant, Project Euler has re-added the solution checking feature to their website, so there should be no issues related to that moving forward.

Now, for Problem 87, I used D, which I had just recently freed up, having not used it since my initial solution to problem 22. D is a fine language with C-like syntax which compiles directly to native code and lets you do all sorts of stuff that one would want to do in a low or high level language...its a powerful language, which was a nice breath of fresh air having just come from SwiftScript. Problems 87 is  rather simple problem: "how many numbers below 50000000 can be expressed as the sum of a prime square, prime cube, and prime fourth power?": the main issue comes in finding a way to solve it efficiently, which I did by storing all square+cube sums in a set (implemented by hand as an array of bits) and then checking against this set to add in the fourth powers. This was much more efficient than my first solution, which, without going into too many details, took 5min to run.

Other than that, the only issue I had was that I trusted myself too much with my prime sieve. I wrote it from memory based on one I have written in python a number of times, but I made a few small mistakes the first couple of times around that caused it to not quite find all the primes. Of course, this in turn led to some difficulties. However, once I discovered that annoying little error, everything fell into place. Solution runs in 10 seconds on my machine.

import std.stdio;
import std.math;
import std.conv;

/** Writing all the casts is annoying */
int isqrt(int x) {
    return to!int(floor(sqrt(to!float(x))));
}

/** "Basic" Prime Sieve */
int[] sieve(int bound) {
    int[] partial;
    int rbound = isqrt(bound);
    if (bound < 11) {
        return [2, 3, 5, 7];
    } else {
        partial = sieve(rbound);
    }
    int current = rbound;
    if (current % 2 == 0) {
        current += 1;
    }
    int[] newprimes = partial;
    while (current < bound) {
        bool isprime = true;
        int rcurrent = isqrt(current);
        foreach (p; partial) {
            if (p > rcurrent) {
                break;
            }
            if (current % p == 0) {
                isprime = false;
                break;
            }
        }
        if (isprime) {
            newprimes ~= current;
        }
        current += 2;
    }
    return newprimes;
}

int BOUND = 50000000;
/**Set implemented with array of bits
   Not only done for space efficiency, D does not
   allow static arrays with 50000000 elements (at least by default).
   */
uint[50000000/32 + 1] is23sum;

void setIsSum(int elem) {
    if (!isSum(elem)) {
        is23sum[elem / 32] += 1 << (elem % 32);
    }
}

bool isSum(int elem) {
    return ((is23sum[elem / 32] >> (elem % 32)) & 1) > 0;
}
void main() {
    // Get all necessary primes
    int[] primes = sieve(isqrt(BOUND));
    // Precompute arrays
    int[] primes2;
    foreach (p ; primes ) { 
        primes2 ~= pow(p , 2); 
    }
    int[] primes3;
    foreach (p ; primes ) {
        int temp = pow(p, 3);
        if (temp > BOUND) { 
           break;
        } else { 
           primes3 ~= temp;
        }
    }
    int[] primes4;
    foreach (p; primes) {
        int temp = pow(p, 4);
        if (temp > BOUND) {
            break;
        } else {
            primes4 ~= temp;
        }
    }
    is23sum[] = 0;
    // Precompute all numbers that are sums of squares and cubes
    foreach (p2; primes2) {
        foreach (p3; primes3) {
            int s = p2 + p3;
            if (s < BOUND) {
                setIsSum(s);
            }
        }
    }
    int ans = 0;
    // Checking for sum of square, cube, and 4th power
    foreach(i; 27 .. BOUND) {
        foreach (p4; primes4) {
            if ( (i - p4) < 12) {
                break;
            } else if (isSum(i - p4)) {
                ans = ans + 1;
                break;
            }
        }
    }
    writeln(ans);
}