Progress

Below is the table of my progress on this Language Challenge. The numbers below link to my code for the solution, statements of the problems can be found on projecteuler.net. A * indicates a language that was used at some point but has yet to be reused to solve anything, and P&P denotes the problems that have been solved with pencil and paper, no coding necessary.

BrainFuck	1*
Ook!	6*
Whenever	5,2
Cat	3
FALSE	4
GLSL	7
ArnoldC	9
LOLCODE	10
Io	11
Shakespeare	12
Smalltalk	13
Clojure	14
BASIC	16
INTERCAL	17
WIND	18
F#	19
E	20
COBOL	21
SwiftScript	22
Ceylon	23
Erlang	24
Kotlin	26
Befunge	27
Boo	29
K	8,30
ALGOL68	31
Go	12,32
Bash	33
Batch	34
ChucK	35

Whitespace	36
Haskell	2,37
Lua	38
Gosu	39
Rust	40
Ada	41
sML	10,42
Coffeescript	43
Scala	18,44
Rexx	45
Julia	46
x86-64	47
ELM	17,48
OCaml	49
Postscript	7,50
Cobra	51
APL	52
EEL	53
Chapel	54
Elixir	55
Ruby	4,40,56
Racket	7,26,57
WARM	58
C#	26,59
Javascript	9,31,60
Pascal	61
cLisp	13,62
Rebol	63
Tcl	64
Dart	63,65

Python	16,24,39,65,66
Prolog	18,67
Fantom	68
Perl	19,70
Processing	71
J	3,24,35,69,72
Groovy	73
Genie	74
Vala	75
Forth	9,76
Hack	77
R	45,78
CIL	79
Frink	16,80
Dogescript	81
Fortran95	3,15,82
Zimbu	83
Red	84
Idris	85
Squirrel	21,52,55,86
D	22,87
C	48,50,88
PASM	89
JavaBC	90
Java	20,79,91
X10	92
PHP	11,93
C++	14,53,94
P&P	1,5,6,8
P&P	15,25,28,69

Problem 89 - PASM

After solving problem 88, problem 89 is not a particularly difficult problem, it involves doing stuff with Roman Numerals. As I expected the problem to not be too difficult, I looked for a new language that would not be super easy to use to solve this problem. So, I stumbled across and decided to use PASM, the Parrot ASseMbly language for the Parrot Virtual Machine (http://www.parrot.org/). Now, though it is an assembly language, it is an assembly language for a virtual machine, so it a CISC language, and it notably has a number of language one does not find in most assembly languages, most notably a number of string manipulation functions (used below are length, substr, ord. Also readline and a simple print command are nice to have). So, the string manipulation niceties were nice, especially as problem 89 involves reading in and processing a large text file. The main difficulties I had working in PASM was 1) figuring out how to declare constants...even though this ended up being unnecessary, .const and .constant are both documented as existing but don't, so this took a while. And 2) Attempting to use an array. One disadvantage of PASM running on a VM is that, unlike a normal assembly language, memory management is hard, and I could not figure out how to allocate/use an array (I could use a stack, but random access in a stack is...not fun). This difficulty cost a lot of time while writing the program, and also resulted in the GET_CHAR_COUNT method not being too nice.

I think a number of the issues described above were due to my choice to use PASM, while Parrot apparently intends for most users to use (and of course by "use" I mean have their compilers output) PIR (Parrot Intermediate Representation) code which is then converted to PASM. And because this is the intent, good PASM documentation was harder to find than PIR. Luckily, there was a list of opcodes, and a whopping 4 examples that came with Parrot, so those two things combined are what I used to get through most of this.

The last issue was stuff with the input file having Windows newlines while I only tested on Unix newlines to begin with. But oh well, that was fixed relatively easily.

My use of PASM has inspired me to look into soon using a similar language...Java ByteCode (and, if it is sufficiently different, PIR).

Anyway, even if it runs on a VM, it is still basically an assembly language, so my solution below runs very quickly: about 25ms on my machine.

.loadlib 'io_ops'

#ASCII values of important letters
.macro_const M 77
.macro_const D 68
.macro_const C 67
.macro_const L 76
.macro_const X 88
.macro_const V 86
.macro_const I 73

#Impromptu convention:
#arguments passed in I/S 1, returned in I/S 0
.pcc_sub :main main:
    #'The Stack'
    new         P10,'ResizableIntegerArray'
    getstdin    P0
    set         I9, 0
MAINLOOP:
    readline    S1, P0
    unless      S1, FINAL
    #get length of string
    length      I8, S1
    #subtract 2 from length because windows newline
    sub         I8, I8, 2
    inc         I8
    #cutoff one char of newline, as methods expect
    #UNIX newlines
    substr      S1, S1, 0, I8
    dec         I8
    #Now that that is done...
    local_branch P10, GET_STRING_VALUE
    local_branch P10, GET_CHAR_COUNT
    sub         I0, I8, I0
    add         I9, I9, I0 
    branch      MAINLOOP
FINAL:
    say         I9
    end

#Basically a switch statement
#input: I0
#output: I0
GET_CHAR_VALUE:
    eq          I0, .M, M
    eq          I0, .D, D
    eq          I0, .C, C
    eq          I0, .L, L
    eq          I0, .X, X
    eq          I0, .V, V
I:
    set         I0, 1
    local_return P10
V:  
    set         I0, 5
    local_return P10
X:  
    set         I0, 10
    local_return P10
L:  
    set         I0, 50
     local_return P10
C:  
    set         I0, 100
    local_return P10
D:  
    set         I0, 500
    local_return P10
M:  
    set         I0, 1000
    local_return P10

#compute integer corresponding to Roman Numeral
#input: S1
#output: I0
#used: I2, I3, I4, I5, S2
GET_STRING_VALUE:
    length      I2, S1
    dec         I2
    set         I3, 0
    set         I4, 0
    set         I5, 0
LOOP:
    ge          I3, I2, ENDLOOP
    substr      S2, S1, I3, 1 
    ord         I0, S2
    local_branch P10, GET_CHAR_VALUE
    set         I4, I0
    inc         I3
    substr      S2, S1, I3, 1
    ord         I0, S2
    local_branch P10, GET_CHAR_VALUE
    lt          I4, I0, _SUB
    add         I5, I5, I4
    branch      LOOP
_SUB:
    sub         I5, I5, I4
    branch      LOOP    
ENDLOOP:
    set         I0, I5
    local_return P10

#This could be written much more nicely
#If I could figure out how to do arrays
#input: I0
#used: I1, I6, I7
GET_CHAR_COUNT:
    #Move input into I1
    #And compute number of M's
    set         I1, I0
    div         I6, I1, 1000
    mul         I0, I6, 1000
    sub         I1, I1, I0
    set         I7, I6
    #CM
    div         I6, I1, 900
    mul         I0, I6, 900
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #D
    div         I6, I1, 500
    mul         I0, I6, 500
    sub         I1, I1, I0
    add         I7, I7, I6
    #CD
    div         I6, I1, 400
    mul         I0, I6, 400
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #C
    div         I6, I1, 100
    mul         I0, I6, 100
    sub         I1, I1, I0
    add         I7, I7, I6
    #XC
    div         I6, I1, 90
    mul         I0, I6, 90
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #L
    div         I6, I1, 50
    mul         I0, I6, 50
    sub         I1, I1, I0
    add         I7, I7, I6
    #XL
    div         I6, I1, 40
    mul         I0, I6, 40
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #X
    div         I6, I1, 10
    mul         I0, I6, 10
    sub         I1, I1, I0
    add         I7, I7, I6
    #IX
    div         I6, I1, 9
    mul         I0, I6, 9
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #V
    div         I6, I1, 5
    mul         I0, I6, 5
    sub         I1, I1, I0
    add         I7, I7, I6
    #IV
    div         I6, I1, 4
    mul         I0, I6, 4
    sub         I1, I1, I0
    mul         I6, I6, 2
    add         I7, I7, I6
    #I
    add         I0, I7, I1
    local_return P10 

Problem 88 - A struggle completed in C

Problem 88 was quite the problem. I am absolutely certain that the answer I have now is not the best solution, but it is definitely a solution. This problem investigates "Product Sum numbers" - numbers that are both the sum and product of collections of numbers. I got some insight into some of the math going on by reading a paper I found online on the subject: http://www-users.mat.umk.pl/~anow/ps-dvi/si-krl-a.pdf. From what I found in that paper, I wrote up a solution in Python, which gave the correct answer...in 30 minutes. In the search for a faster solution, I realized that I was actually computing all I needed for an answer in the first minute or so of the problem and then used a lot of unnecessary computation to find the answer.

The key to making the problem tractable was to notice that though the problem asks for the minimal Product-Sum numbers for a given k (size of the set for which it is a PS number), the easier domain was to look at numbers and find for what k values they are product sum numbers. When I made this switch, I when from minutes to seconds in the speed of my code.

However, the script that I ran to initially compute the answer was written in python, for ease, but I knew I had to convert to another language eventually for the language challenge. I ended up going with C, and while C works, I ran into...a bit of an issue. In Python I was just using lists and wasn't aware of just how much memory was flying around, but when I wrote a solution in C I found out that the solution actually involves allocating many million-element arrays...so after pushing the size of my arrays and making sure to allocate and free memory well to avoid seg faults, I eventually got to a working solution, which is below. Runs in about 16 seconds in my machine. Much of that time is probably spent in malloc.

#include<stdio.h>
#include<stdlib.h>

int BOUND = 12000;

/** Gets the sums of all factoriations of n, minus the number of terms in the
 *  factorization. This is most of the k-value for which n is a Product-Sum number
 *  k-value = n - factorizationSum*/

/** Admittably not the most efficient of code...this is what I get for
 *  trying to transfer code originally written in python
 */
void getAllFactorizationSums(int n, int includeself,
        int maxlength, int* output, int* outputLength) {

    /**Don't look at exceptionally long factorizations*/
    if (maxlength == 0) {
        output[0]     = -1;
        *outputLength =  1;
        return;
    }
    int*  returnValue = (int*) malloc(sizeof(int) * 15 * BOUND);
    int  rvLength = 0;
    int f;
    for (f = 2; f*f <= n; ++f) {
        if (n % f == 0) {
            getAllFactorizationSums(n / f, 1, maxlength - 1, output, outputLength);
            int j;
            for (j = 0; j < *outputLength; ++j) {
                int fSum = output[j];
                if (fSum != -1) {
                    returnValue[rvLength++] = fSum + f - 1;
                }
            }
        }
    }
    if (includeself) {
        returnValue[rvLength++] = n - 1;
    }
    //copy back to output
    int i;
    for (i = 0; i < rvLength; ++i) {
        output[i] = returnValue[i];
    }
    *outputLength = rvLength;
    free(returnValue);
} 

int min(int x, int y) { return x < y ? x : y; }

int main(int argc, char ** argv) {
    int* productSumsForK = (int*) malloc(sizeof(int) * (BOUND + 1));
    productSumsForK[0] = 0; productSumsForK[1] = 0;
    int i;
    /**Initialize with upper bounds*/
    for (i = 2; i <= BOUND; ++i) {
        productSumsForK[i] = 2 * i;
    }
    /**Allocate lots of space...most of it actually gets used!*/
    int* factorizationSums = (int*) malloc(sizeof(int) * 15 * BOUND);
    int fsLength = 0;
    int N;
    for (N = 4; N <= 2*BOUND; ++N) {
        getAllFactorizationSums(N, 0, 15, factorizationSums, &fsLength);
        for (i = 0; i < fsLength; ++i) {
            int kIndex = N - factorizationSums[i];
            if (kIndex <= BOUND) {
                productSumsForK[kIndex] = min(productSumsForK[kIndex], N);
            }
        }
    }

    //Now, we need to convert to a set to print the answer
    char set[BOUND/4];
    for (i = 0; i < BOUND/4; ++i) { set[i] = 0; }
    for (i = 0; i <= BOUND; ++i) {
        int psNum = productSumsForK[i];
        set[psNum / 8] |= (1 << (psNum % 8));
    }
    int sum = 0;
    for (i = 0; i < 2*BOUND; ++i) {
        if ((set[i / 8] >> (i % 8)) & 1) {
            sum += i;
        }
    }
    printf("%d\n", sum);
    return 0;
}

Problem 50 - Postscript

I have been working out a solution to the next problem, problem 88, and I figured I wanted a pretty capable language for my solution to that one, so I decided to free up such a language. Now, the language i had free to do this freeing up with was Postscript. Prior to this project, my main experiences with Postscript had been writing an interpreter for a subset of the language, and using the language way back when in my initial 24 hour run to solve problem 7. Now, however, I used Postscript to free up C, which means that I used it to solve problem 50. Once I got into the swing of things, the stack-basedness wasn't too bad, the main issues that arose were due to my inability to parse error messages, causing bugs to be hard to remove (and, as would be expected when using an unfamiliar programming paradigm, bugs abounded). This problem gave me a chance to write an implementation of my smarter sieve method from my analysis before...though in my attempt to also generate an independent isPrime method, I messed up at first (an independent isPrime method cannot assume the input number is not divisible by 2 or 3, it actually needs to check that!). But anyway, other than lots of little mistakes like that, things weren't too bad. Postscript isn't the most blazingly fast language though, so even with a number of optimizations that were not present in my C solution to problem 50, the solution below runs in about 50 seconds on my machine, whereas my C solution ran in about 0.5s. But either way, it runs in under a minute, so it is acceptable. Now, to solve problem 88 in C.

%Implementation of "Smartest Sieve" from Sieve Timing :)
%Assumes symbol /primes is defined
/isPrime {
    /x exch def
    /ip true def
    /i 0 def
    /for2 {
        /p primes i get def
        x p mod 0 eq {
            /ip false def
        } {
            p p mul x le {
                /i i 1 add def
                for2
            } if
        } ifelse
    } def
    for2   
    ip     
} def

/primeSieve {
    /N exch def
    /primes N array def
    primes 0 2 put
    primes 1 3 put
    primes 2 5 put
    /primeLength 3 def
    /n 7 def
    /for1 {
        n isPrime {
            primes primeLength n put
            /primeLength primeLength 1 add def
        } if
        /n 2 n mul 3 mod 2 mul n add def
        n N lt {
            for1
        } if
    } def
    for1
} def

%Assumes symbol /primes is defined
/partialSum {
    /sum 0 def
    /high exch def
    /low exch def
    /k low def
    /for3 {
        /sum sum primes k get add def
        /k k 1 add def
        k high lt {
            for3
        } if
    } def
    for3
    sum
} def

/BOUND 1000000 def
%Method defines symbols /primes and /primeLength
%We do not use any return value
BOUND primeSieve
/ans 0 def
/maxChainLength 1 def
/lindex 0 def
/for4 {
    lindex primeLength lt {
        /hindex lindex maxChainLength add def
    hindex primeLength lt {
        /psum lindex hindex partialSum def
        /for5 {
            hindex primeLength lt psum BOUND lt and {
                psum isPrime {
                    /ans psum ans max def
                    /maxChainLength hindex lindex sub def
                } if
                /psum psum primes hindex get add def
                /hindex hindex 1 add def
                for5
            } if
        } def
        for5
        /lindex lindex 1 add def
        for4
    } if
    } if
} def
for4

ans pstack pop

Sieve Timing Continued - A simpler method

After writing up all my sieve timing stuff yesterday, I was thinking of things I could do to add more information to my set of sieve times, and in an attempt to implement some additional optimizations, I realized a method that I probably should have realized a long time ago: instead of computing sub-list of primes based on square roots and such, why not just generate a list of primes in place? Somehow this never came to mind before, but this time it did, and lo and behold, it is about as fast or maybe faster than the recursive method (In the data I have it is much faster at 10^8, but I only ran one trial, so that information is sketchy). Most importantly though, the code is much simpler. The code is below:

/** For some reason this didn't occur to me for a really long time */
void smartestSieve(int bound, int* primes, int* primeslength) {
    primes[0] = 2; primes[1] = 3; primes[2] = 5; int mLength = 3;
    int n;
    for (n = 7; n < bound; n += 2 * ((2*n) % 3)) {
        int isprime = 1;
        int i;
        for (i = 2; i < mLength; ++i) {
            int p = primes[i];
            if (p*p > n) {
                break;
            } else if (n % p == 0) {
                isprime = 0;
                break;
            }
        }
        if (isprime) {
            primes[mLength++] = n;
        }
    }
    *primeslength = mLength;
}

Timing information is below, first, I slotted it into the original graph, removing the two slowest methods:

 And, to emphasize the "part that matters", its comparison with the recursive method:

The information used to generate these graphs is now on the second page of the sheet linked to in the previous post, https://docs.google.com/spreadsheets/d/18-avzFTZDDOPMh7p_FRPOVc9SVfEhxPEoVFGQPTyNWk/edit?usp=sharing

I am now considering looking at the advantages of hand tooled assembly and/or comparing the same methods in Python.

(Something Different) - Sieve Timing

So, most of my posts on this blog are direct solutions to Project Euler problems. However, this time I am writing on something I did which is kind of meta-Project Euler related. A lot of Project Euler problems involve prime numbers, and therefore Prime Sieves often have to be written. For anyone who may not know, a prime sieve is a program that finds primes basically by checking numbers for their divisibility by other numbers: if a number is not divisible by any of the numbers below it (other than 1), then it is prime. However, there are plenty of optimizations to this process that can be done.

For a recent problem (87), I decided to write up a somewhat complicated prime sieve that calls itself recursively to ensure that you only ever "sieve" with prime numbers. That is, to check that a number is prime, instead of checking for divisibility by all numbers less than it, or even all numbers less than the square root, or even all odds less than the square root, you only do divisibility tests with primes. I am pretty happy with my sieve,  which I wrote up first in python quite a while ago. However, after spending the time to rewrite it in D, I wondered whether it was really worth it. Does the improved sieve actually save a significant amount of time over a simpler sieve?

I decided to investigate this question thoroughly. I wrote up a series of sieve methods in C (so that I can feel confident about what code is actually getting executed, unlike in a language like python). The methods are as follows:

1. isPrimeDumbest: checks primality of a number n by dividing by 2..n-1
2. isPrimeSqrt: checks primality of a number n by dividing by 2..sqrt(n)
3. isPrimeSqrt2: checks primality of a number n by dividing by 2 and also all odd numbers 3..sqrt(n).
4. isPrimeSqrt23: checks primality of a number n by dividing by 2 and 3 and then also all numbers in 5..sqrt(n) not divisible by 2 or 3.
5. smartSieve: checks primality of a number n by dividing by all primes in 2..sqrt(n)

Obviously checking only up to the square root will be better than checking all numbers, and obviously checking only odd numbers will give an approximately 2x speed up. The main things I was considering was how much also predividing by 3 helps (it lets you check only 2/7 of numbers, but requires a wee bit extra computation), and whether the overhead of recursion and memory accesses allowed the smart sieve method to be as smart as I wanted it to be.

And in the end, yes, it is as smart as I wanted it to be: to compute all primes below 10^8, in the tests I ran, smartSieve is about 5x faster than isPrimeSqrt, and also more than twice as fast as isPrimeSqrt23. So, in conclusion, smartSieve is worth using, at least if the number of primes concerned is at least around 10^7 (and notably, problem 87, which brought up all of this, required finding all primes less than 5*10^7). Also the tests drove home just how bad it is to not at least shorten things to the square root and truly brute force things...the dumbest method takes about 5 minute at 10^6, while the next slowest method still takes under half a second to compute all primes.

Below are graphs of my data: The horizontal axis is the N, the number below which all primes were computed. It is on a log scale, so I ran tests on all powers of 10 from 10^1 to 10^8 (100 million). The vertical axis is time, measured in seconds. For smaller N, as the code would run too quickly to get any meaningful timing results I ran each sieve 1000000 times for N = 10 to 10^3, then for 10^4 to 10^7, each sieve was run 50 times, except for the dumbest sieve, which I ran only once at 10^5 because its just super slow, it doesn't really matter how much, and for all of these the result I show is the computed average running time of multiple trials. Because 10^8 takes so long and the data is spread out enough, I only ran these tests once instead of computing an average.

Now that I have fully described what the data is, here is the data, which shows exactly what one would expect, except for perhaps the magnitudes. The first graph shows all methods (except no data for the dumbest method above 10^6 because that would take too long to compute), and the second graph shows all but the dumbest method run on small N, so that differences can actually be seen (though the differences are oh-so small on an absolute scale).

In case you want to know the exact actual results: all the data is in a google spreadsheet here: https://docs.google.com/spreadsheets/d/18-avzFTZDDOPMh7p_FRPOVc9SVfEhxPEoVFGQPTyNWk/edit?usp=sharing

Lastly, here is the code I ran to generate the information. In order to get code running as quickly as possible, and because it doesn't hurt, the below was compiled with gcc -O3 (highest amount of optimization). I haven't included the testing code, but its a fairly no-frills python script.

Also, it is notable that I am sure there are other prime sieve methods that may be faster out there in this world, and maybe I should eventually check out others, but this analysis is sufficient for my own purposes for now at least.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

/** Very brute force sieve
 *  Always checks 2..n-1
 */
int isPrimeDumbest(int n) {
    int i;
    for (i = 2; i < n - 1; ++i) {
        if (n % i == 0) {
            return 0;
        }
    }
    return 1;
}

/** Mostly brute force, but checks
 *  2..n-1
 */
int isPrimeSqrt(int n) {
    int i;
    for (i = 2; i*i <= n; ++i) {
        if (n % i == 0) {
            return 0;
        }
    }
    return 1;
}

/** Checks for 2 explicitly, then does a += 2 */
int isPrimeSqrt2(int n) {
    if (n % 2 == 0) { return n == 2; }
    int i;
    for (i = 3; i*i <= n; i+=2) {
        if (n % i == 0) {
            return 0;
        }
    }
    return 1;
}

/** Sieves out both 2 and 3 */
int isPrimeSqrt23(int n) {
    if (n % 2 == 0) { return n == 2; }
    if (n % 3 == 0) { return n == 3; }
    int i;
    for (i = 3; i*i <= n; i += (i % 3 == 1 ? 4 : 2)) {
        if (n % i == 0) {
            return 0;
        }
    }
    return 1;
}

/** Function that takes one of the isPrime Methods above
 *  and sieves:
 *  For consistency, all sieves know 2 is prime
 *  and then check all odds
 */
void mostSieves(int bound, int* primes, int* primeslength, int (*sieveFunc)(int)) {
    primes[0] = 2;
    int c = 1;
    int i;
    for (i = 3; i < bound; i+=2) {
        if (sieveFunc(i)) {
            primes[c++] = i;
        }
    }
    *primeslength = c;
}
/** Recursively determines prime lists
 *  So that mods are only taken for
 *  primes <= sqrt */
void smartSieve(int bound, int* primes, int* primeslength) {
    if (bound < 11) {
        primes[0] = 2; primes[1] = 3; primes[2] = 5; primes[3] = 7;
        *primeslength = 4;
        return;
    }
    int rbound = (int)sqrt(bound * 1.0f);
    smartSieve(rbound, primes, primeslength);
    int partialLength = *primeslength;
    int mLength       = partialLength;
    int current = rbound;
    if (current % 2 == 0) { ++current; }
    while (current < bound) {
        int isprime = 1;
        int i ;
        for (i = 0; i < partialLength; ++i) {
            int p = primes[i];
            if (p*p > current) {
                break;
            } else if (current % p == 0) {
                isprime = 0;
                break;
            }
        }
        if (isprime) {
            primes[mLength++] = current;
        }
        current += 2;
    }
    *primeslength = mLength;
    return; 
}

int main(int argc, char** argv) {
    if (argc < 3) { 
        printf("%s\n", "Not enough arguments");
        return 1;
    }
    int fchoice  = atoi(argv[1]);
    int bound    = atoi(argv[2]);
    int numtimes = atoi(argv[3]);
    /** Allocate an over estimate in the beginning
     *  So that resizing does not effect times
     */
    int* primes = (int*)malloc(bound * sizeof(int));
    // Why not use this place in memory for this...
    int primeslength = 0;
    int i;
    for (i = 0; i < numtimes; ++i) {
    switch(fchoice) {
        case 0:
            mostSieves(bound, primes, &primeslength, &isPrimeDumbest);
            break;
        case 1:
            mostSieves(bound, primes, &primeslength, &isPrimeSqrt);
            break;
        case 2:
            mostSieves(bound, primes, &primeslength, &isPrimeSqrt2);
            break;
        case 3:
            mostSieves(bound, primes, &primeslength, &isPrimeSqrt23);
            break;
        case 4:
            smartSieve(bound, primes, &primeslength);
            break;
        default:
            printf("Not a valid choice\n");
            return 1;
    }
    }
    return 0; 
}

Problem 7 - GLSL

By this point, I have worked in a lot of programming languages. And, for obvious reasons, just about any programming language that I have used in my life for something, I have used in this blog to solve a Project Euler problem. However, this has not been the case for quite all languages. As I do work in the graphics lab at Williams College, naturally I have worked with GLSL, a shader language. Now, as GLSL is intended for shaders and running on the GPU, there are some issues using it for a Project Euler problem.

Perhaps the largest and most obvious issue is how to output an answer in GLSL. Shaders produce images, so clearly in order to output the answer from GLSL, I had to create an image. There are these things called 7-segment displays that display numbers pretty well, without much complicated logic, so I wrote up some GLSL to make a 7-segment display, or, to be more specific, some webGL. webGL is less featured than GLSL, including some annoying features like how in wbGL for loops must have constant bounds. The 7-segment display can be seen on shadertoy: https://www.shadertoy.com/view/4ssSz2.

Further frustrations with webGL include lack of much support for integer arithmetic (almost all of graphics is done with floating point numbers). See the hand-written modulus function and hand power-of-10 calculations for more info on this particular issue...

Then, when I had my code almost all running, I ran an issue with webGL. Remember how I said that for loops must have constant bounds? Well, the reason why is because the loops must get unfolded in the compiler. Now, unfolding a loop 150000 times...that is not the best, indeed, attempting to compile my shader crashed webGL on the machine I was using. So, because that didn't work, I may yet return to webGL to try to solve a problem, keeping in mind all the additional issues.

GLSL has smart loop bounds and doesn't need to do crazy stuff like unfold loops, so my solution below ran and gave the answer. Now, another interesting thing about using GLSL. This code is a pixel shader. What that means is that it is run once per pixel, using information about where that pixel is on screen. So though GPUs are known for being very good at parallelizing things, this problem used parallelism in the dumbest way possible, which is having every pixel compute the answer by itself, and then check if the pixel should be green or black in order to display it in a 7-segment display. Even though about 1000x more work was done than necessary, this runs in a reasonable time (maybe a couple seconds? I didn't set up anything to time it) on my machine. And here is the output, in all of its glory:

#define NUMBER_TO_PRINT esol
#define DISTANCE_BETWEEN_DIGITS 0.05
#define LINE_WIDTH 0.05
#define ON_COLOR  vec4(0.0, 1.0, 0.0, 1.0);
#define OFF_COLOR vec4(0.0, 0.0, 0.0, 1.0);

float digits(float x) {
    if (x < 2.0) { return 1.0; }
    return ceil(log(x) / log(10.0));
}

int imod(int x, int y) {
    return x - y*(x/y);
}

int euler7() {

    int numprimes = 1;
    int ans = 0;
    for (int i = 3; i < 150000; i += 2) {
        bool isprime = (imod(i, 2) != 0);
        for (int j = 3; j*j <= i; j += 2) {
            isprime = isprime && (imod(i, j) > 0);
        }
        numprimes  = numprimes + int(isprime);
        ans       += (int(numprimes == 10001) * i);
        numprimes += int(numprimes == 10001);
    }
    return ans;
}

uniform vec3 iResolution;

void main(void) {

    int esol = euler7();
    vec2 uv = gl_FragCoord.xy / iResolution.xy;
    //Below line switches between glsl and webgl coordinaes
    //uv.y = 1.0 - uv.y;

    float N = digits(float(NUMBER_TO_PRINT));
    float width = (1.0 - (N-1.0)*DISTANCE_BETWEEN_DIGITS) / N;
    float digit = floor(uv.x / (width + DISTANCE_BETWEEN_DIGITS));

    //coordinates within box
    float digitStartx = (width+DISTANCE_BETWEEN_DIGITS)*digit;
    float digitEndx   = digitStartx + width;
    vec2 sub_uv = vec2((uv.x - digitStartx) * 1.0/(digitEndx - digitStartx), (uv.y - .2) * (1.0/.6));
    //Computing a power of 10 was made harder because webgl
    //Could do something less silly in actual glsl
    int p = 1;
    for (int i = 0; i < 10; ++i) {
        p *= (1 + int(i < int(N) - int(digit)) * 9);
    }
    int currentDigit = imod(NUMBER_TO_PRINT, p) / (p/10);

    //Determining the position
    bool tiptop   = sub_uv.y < LINE_WIDTH && sub_uv.y > 0.0;
    bool top      = sub_uv.y < 0.5 - LINE_WIDTH/2.0&& sub_uv.y > 0.0;
    bool left     = sub_uv.x > 0.0                 && sub_uv.x < LINE_WIDTH;
    bool middle_x = sub_uv.x > LINE_WIDTH          && sub_uv.x < 1.0 - LINE_WIDTH;
    bool middle_y = sub_uv.y > 0.5 - LINE_WIDTH/2.0&& sub_uv.y < 0.5 + LINE_WIDTH/2.0;
    bool right    = sub_uv.x > (1.0 - LINE_WIDTH)  && sub_uv.x < 1.0;
    bool bottom   = sub_uv.y > 0.5 + LINE_WIDTH/2.0&& sub_uv.y < 1.0;
    bool botbot   = sub_uv.y > 1.0 - LINE_WIDTH    && sub_uv.y < 1.0;
    

    //shorthand

    int c = currentDigit;
    bool bit =(top && left  && (c != 1 && c != 2 && c != 3 && c != 7)) ||
              (top && right && (c != 5 && c != 6))                     ||
              (tiptop && middle_x && (c != 1 && c != 4))               ||
              (middle_x && middle_y && (c > 1 && c != 7))              ||
              (bottom && left && (c/2*2 == c) && c != 4)               ||
              (bottom && right && (c != 2))                            ||
              (botbot && middle_x && (c != 1 && c != 4 && c != 7));
    gl_FragColor = float(bit && true) * ON_COLOR + float(!(bit && true)) * OFF_COLOR;
}

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);
}

Summertime Problem 22 in SwiftScript

First a note about project euler. The website had its security compromised, and many of the features of the site, including answer checking, are down for now. Therefore, I will likely stick with resolving solved problems, at least for a little while, as this problem is doing. Now, for the problem:

So, I heard about the new language Apple apparently developed to replace Objective-C for their uses, Swift (https://developer.apple.com/swift/), and I thought I might use it, but it seems like all the tools for it are built into Xcode and other Apple things I don't want to have to deal with. However, at the bottom of the page, there is a small note, to compensate for the fact that Apple's Swift is not the first language called Swift: "Looking for the Swift parallel scripting language? Please visit http://swift-lang.org"

So, I visited that website, and found the Swift parallel scripting language. As the name suggests, yes, it is a parallel scripting language. And that basically sums up the issues. One: it is a scripting language. It seems that the developers of Swift only ever intended for it to be used to gather large collections of filenames and commands and run other programs in parallel using it. Two: it is a parallel language, which opens up its own can of worms, such as the fact that all variables can only be written to exactly once.

Problem 22 seemed like a problem that would be doable in SwiftScript, as the largest task of the problem involves reading in a list of names. However, despite having huge amounts of tools for file output, in order to extract a list of strings from a file in SwiftScript, the easiest solution I could find was to treat it as a one-column csv, which then gets parsed as a struct with one string field, which is the name at the top of the list (AARON, which had to have an extra copy added to names.txt), and then getting a normal string array from that.

Another task one would think would be simple: getting the ascii value of a char. SwiftScript has no chars though, only ints and strings, and no way to map from string to int (though there is @toString for the other direction). However, as SwiftScript is built on top of Java, it has a @java directive for calling static methods of Java library classes. So, I used java.lang.Integer.parseInt to tread my names as base-36 integers, in order to extract a mapping a => 1, b => 2, etc.

Sorting the list proved to be incredibly hard in SwiftScript (due to being able to only assign a value to a variable once, and all for loops being executed in parallel, most sorting methods would not work), so to accomplish this task, I used SwiftScripts ability to run other programs and used sort to sort the list.

The last task I had difficulty solving was finding the sum of a list. A recursive funciton (which I use in some places) took far too long to run on a list the size of the number of names this problem must deal with, and you can't just use an accumulator and run through the list, as += does not exist in SwiftScript. However, using a parallel foreach loop to build up an array of running sums does work, as each thread will wait until it has the information to compute the computation.

After completing all of these tasks over a large amount of time, I got to the solution that is below. Runs in 27s on my machine.

type file;
type line {
    string AARON;
}

(int sum) arraysum (int[] arr, int index, int acc) {
    if (index >= @length(arr)) {
        sum = acc;
    } else {
        sum = arraysum(arr, index+1, acc+arr[index]);
    }
}

#Incredibly long winded manner of getting the ascii value
#of a char
(int val) parseInt (string s) {
    int temp = @java("java.lang.Integer", "parseInt", s, 36);
    val = temp - 9;
}

(int val) strvalue (string str) {
    string chars[] = @strsplit(str, "");
    int values[];
    #we start at 1 to chop off strange empty char
    foreach i in [1:@length(chars)-1] {
        values[i-1] = parseInt(chars[i]);
    }
    val = arraysum(values, 0, 0);
}

#Though I might normally find it to be almost like cheating
#to call other programs, this is kind of what SwiftScript
#was made for; and its a program that is pre-installed
#on most UNIX systems.
app (file outfile) sortwords(file infile) {
    sort stdin=@filename(infile) stdout=@filename(outfile);
}

#had to add extra AARON to names.txt to get desired result
#As SwiftScript does not parse strict lists, it requires a 
#category at the top of columns of a csv type file

file original <"names.txt">;
file sorted   <"sorted_names.txt">;
sorted = sortwords(original);
line lines[] = readData("sorted_names.txt");

#Converting the array of lines to strings
string[] names;
foreach i in [0:@length(lines)-1] {
    names[i] =  lines[i].AARON;
}

#Get the name scores
int[] scores;
foreach i in [0:@length(names)-1] {
    scores[i] = (i+1) * strvalue(names[i]);
}

#Workable way to compute the sum of a list
#given the constraints of the language
int[] scoresums;
scoresums[0] = scores[0];
foreach i in [1:@length(names)-1] {
    scoresums[i] = scoresums[i - 1] + scores[i];
}
#Ensures this does not run until
#scoresums is fully written
if (@length(scoresums) == @length(scores)) {
    trace(@toString(scoresums[@length(scoresums)-1]));
}

New Progress Table

Below is the table of my progress on this Language Challenge. The numbers below link to my code for the solution, statements of the problems can be found on projecteuler.net. I changed the arrangement of the table in a way that will hopefully make it more useful: previously, the table was formed in chronological order relative to when I used the language: so the language I used most recently would always be at the end of the table. Now, it is arranged in order Project Euler problems: the languages below are listed in order of the problem they are currently being used to solve. Languages with multiple numbers listed indicate languages that were used to solve problems, and then were swapped out with another language, and then were used to solve another problem. Other changes/details: A * indicates a language that was used at some point but has yet to be reused to solve anything, and P&P denotes the problems that have been solved with pencil and paper, no coding necessary.

BrainFuck	1*
Ook!	6*
Whenever	5,2
Cat	3
FALSE	4
GLSL	7
ArnoldC	9
LOLCODE	10
PHP	11
Shakespeare	12
Smalltalk	13
Clojure	14
BASIC	16
INTERCAL	17
WIND	18
F#	19
E	20
COBOL	21
SwiftScript	22
Ceylon	23
Erlang	24
Kotlin	26
Befunge	27
Boo	29
K	8,30
ALGOL68	31
Go	12,32
Bash	33

Batch	34
ChucK	35
Whitespace	36
Haskell	2,37
Lua	38
Gosu	39
Rust	40
Ada	41
sML	10,42
Coffeescript	43
Scala	18,44
Rexx	45
Julia	46
x86-64	47
ELM	17,48
OCaml	49
Postscript	7,50
Cobra	51
APL	52
C++	14,53
Chapel	54
Elixir	55
Ruby	4,40,56
Racket	7,26,57
WARM	58
C#	26,59
Javascript	9,31,60
Pascal	61

cLisp	13,62
Rebol	63
Tcl	64
Dart	63,65
Python	16,24,39,65,66
Prolog	18,67
Fantom	68
Perl	19,70
Processing	71
J	3,24,35,69,72
Groovy	73
Genie	74
Vala	75
Forth	9,76
Hack	77
R	45,78
Java	20,79
Frink	16,80
Dogescript	81
Fortran95	3,15,82
Zimbu	83
Red	84
Idris	85
Squirrel	21,52,55,86
D	22,87
C	48,50,88
P&P	1,5,6,8
P&P	15,25,28,69