Problem 79 - CIL, who needs high level languages?

After using Parrot Assembly and Java Bytecode, I figured I might as well continue my low level / "intermediate representation" language stint. CIL is the "Common Intermediate Language": it is essentially Java Bytecode for the .Net VM. So, programming in it was not too bad, maybe even a little bit nicer than Java Bytecode at times: one standard "call" method was nice instead of invokestatic, invokevirtual, etc. And most of the language is very similar to Java Bytecode - both are stack based, so they use very similar syntax. And the problem I chose to solve was pretty similar pro grammatically to problem 90 that I solved in Java Bytecode: both just involved lots of set operations, which are easy to express as ors, ands, etc (also, my code contains an xor, so clearly I am doing fancy stuff...). Also, CIL was nice as, though Java Bytecode exists only for Java, CIL is supposed to be able to support more than just C#, so there was no need in this case to create a class to wrap the main method in.

There was, however, one very important way in which programming in CIL was far less fun than Java Bytecode. I could not find as much good documentation, and the error messages are horrible. At least using mono (basically the open source .Net VM, while Microsoft also has one) there were essentially two errors one could get: "Unrecoverable Syntax Error" if an error happens at compile time, and "Invalid IL code" if the error happens at runtime. So, even though my code was basically correct very early (a couple of or's had to be changed to and's once I had things running...) Forgetting a single branch statement to go tho the top of a loop and a single pop instruction to clear the stack caused very hard to debug run time errors that took lots of "remove code until the error goes away" style debugging, which was not fun.

So basically, CIL is a fine stack-based assembly languages, as long as you don't cause any errors.
Also, it is worth noting that despite the title of this post...high level languages are actually quite nice. Variable names are really nice, even if I have gotten pretty good at remembering which register holds what variable.

Below code runs in 50ms on my machine.

.assembly extern mscorlib {}
.assembly e79 {}

.namespace E79 
{
    .method public static void main() cil managed
    {
        .entrypoint
        .maxstack 100
        .locals([0] int32[],
                [1] int32,
                [2] int32,
                [3] int32,
                [4] int32,
                [5] int32,
                [6] int32)
        ldc.i4 10
        newarr int32
        stloc.0
        ldc.i4.0
        stloc.3
        ldc.i4.0
        stloc.1
LOOP_INIT:
        ldloc.1
        ldloc.0
        ldlen
        bge OVERLOOP_INIT
        ldloc.0
        ldloc.1
        ldc.i4.0
        stelem.i4
        ldloc.1
        ldc.i4.1
        add
        stloc.1
        br LOOP_INIT
OVERLOOP_INIT:

LOOP_MAIN:
        call string [mscorlib]System.Console::ReadLine()
        dup
        brnull OVERLOOP_MAIN
        call int32 [mscorlib]System.Int32::Parse(string)
        stloc.2
        ldloc.2
        ldc.i4 100
        div
        stloc 4
        ldloc.2
        ldc.i4 10
        div
        ldc.i4 10
        rem
        stloc 5
        ldloc.2
        ldc.i4 10
        rem
        stloc 6
        ldloc.0
        ldloc 5
        ldloc.0
        ldloc 5
        ldelem.i4
        ldc.i4.1
        ldloc 4
        shl
        or
        stelem.i4
        ldloc.0
        ldloc 6
        ldloc.0
        ldloc 6
        ldelem.i4
        ldc.i4.1
        ldloc 5
        shl
        or
        ldc.i4.1
        ldloc 4
        shl
        or
        stelem.i4
        ldloc.3
        ldc.i4.1
        ldloc 4
        shl
        or
        ldc.i4.1
        ldloc 5
        shl
        or
        ldc.i4.1
        ldloc 6
        shl
        or
        stloc.3
        br LOOP_MAIN
OVERLOOP_MAIN:
        pop
        ldc.i4.0
        stloc 4
        ldc.i4.0
        stloc.1
LOOP_FINAL:
        ldloc.1
        ldc.i4 10
        bge OVERLOOP_FINAL
        ldloc.0
        ldloc.1
        ldelem.i4
        ldloc.3
        ldloc.1
        shr
        ldc.i4.1
        and
        ldc.i4.1
        xor
        or
        brzero IFPART
ELSEPART:
        ldloc.1
        ldc.i4.1
        add
        stloc.1
        br LOOP_FINAL
IFPART:
        ldloc 4
        ldc.i4 10 
        mul
        ldloc.1
        add
        stloc 4
        ldc.i4.1
        ldloc.1
        shl
        not
        ldloc.3
        and
        stloc.3
        ldc.i4.0
        stloc 5
LOOP_CLEAR:
        ldloc 5
        ldc.i4 10
        bge OVERLOOP_CLEAR
        ldloc.0
        ldloc 5
        ldloc.0
        ldloc 5
        ldelem.i4
        ldc.i4.1
        ldloc.1
        shl
        not
        and
        stelem.i4
        ldloc 5
        ldc.i4.1
        add
        stloc 5
        br LOOP_CLEAR
OVERLOOP_CLEAR:
        ldc.i4.0
        stloc.1
        br LOOP_FINAL
OVERLOOP_FINAL:
        ldloc.s 4
        call void [mscorlib]System.Console::WriteLine(int32)
        ret
    }
}

Problem 90 - Java Bytecode!

So, after solving problem 89 in Parrot Assembly, I claimed that I was interesting in trying out Java Bytecode. So, I did indeed use it. Turns out Java Bytecode is not exactly the world's most user friendly language to write in. Its not a standard assembly language, but is rather stack based, with some commands taking arguments directly, but most only directing with the stack, which makes things certainly different. Also, thank god I didn't have to invoke too many methods that weren't static in my own class...I wrote a helper function that is a wrapper for System.out.println, it was so frustrating to call it for debugging purposes!

Anyway, problem 90 was a rather simple problem, being capable of being expressed in terms of bitwise operations. Just lots of ors and ands and shifts, treating integers as sets, and iterating over all pairs of possible sets of {0...9}. This is only 1024 * 1024 =~ 1 million sets, and as the operations per pair of sets is exceedingly simple (and most can be thrown out very quickly for not having enough elements to represent a die), this was an entirely reasonable amount of work. Almost all the difficulty with this problem came from the problem, and not the language (well, the only problem related issue was me forgetting that, despite 6 and 9 being identified, a die could still have both, but I only forgot that exceedingly briefly, thanks to the example on the project euler page).

Anyway, below is my code...I would apologize for the lack of comments, as (I know) it makes it hard to read, but at least the ByteCode assembler I used didn't seem to have any syntax for allowing comments...so though unfortunate, it is what it is. Code runs in about 150ms, though it is notable that much of that is probably just starting up the JVM. Maybe my next language will be wonderfully low level again, or maybe I will go up to more sophisticated languages soon. Anyway, here is the beautiful code:

.class public e90
.super java/lang/Object

.method public static printInt(I)V
    .limit stack  100
    .limit locals 100
    getstatic     java/lang/System out Ljava/io/PrintStream;
    iload_0
    invokevirtual java/io/PrintStream println (I)V
    return
.end method

.method public static sixBitsSet(I)I
    .limit stack  100
    .limit locals 100
    iconst_1
    istore_1
    iconst_0
    istore_2
BITSETLOOP:
    iload_1
    iload_0
    if_icmpgt   BITSETLOOP_END
    iload_1
    iload_0
    iand
    ifeq        NOADD
    iinc        2 1
NOADD:
    iload_1
    iconst_1
    ishl
    istore_1
    goto        BITSETLOOP
BITSETLOOP_END:
    bipush      6
    iload_2
    if_icmpeq   RETURN_TRUE
    iconst_0
    ireturn
RETURN_TRUE:
    iconst_1
    ireturn
.end method
    
.method public static setContains(II)I
    .limit stack  100
    .limit locals 100
    iconst_1
    iload_1
    ishl
    iload_0
    iand
    ifgt        RETURN_SUCCESS
    iconst_0
    ireturn
RETURN_SUCCESS:
    iconst_1
    ireturn
.end method
          
.method public static check(III)I
    .limit stack  100
    .limit locals 100
    iload_0
    invokestatic e90/sixBitsSet(I)I
    iload_1
    invokestatic e90/sixBitsSet(I)I
    iadd
    iconst_2
    if_icmpne   RETURN_FAILURE
    iload_0
    sipush      512
    iand
    ifeq        SKIP_ADDSIX1
    iload_0
    sipush      64
    ior
    istore_0
SKIP_ADDSIX1:
    iload_1
    sipush      512
    iand
    ifeq        SKIP_ADDSIX2
    iload_1
    sipush      64
    ior
    istore_1
SKIP_ADDSIX2:
    iload_2
    bipush      10
    irem
    istore_3
    iload_3
    bipush      9
    if_icmpne   SKIP_SIXIFY
    bipush      6
    istore_3
SKIP_SIXIFY:
    iload_0
    iload_3
    invokestatic e90/setContains(II)I
    istore      4
    iload_1
    iload_3
    invokestatic e90/setContains(II)I
    istore      5
    iload_2
    bipush      10
    idiv
    istore_3
    iload_3
    bipush      9
    if_icmpne   SKIP_SIXIFY2
    bipush      6
    istore      3
SKIP_SIXIFY2:
    iload_0
    iload_3
    invokestatic e90/setContains(II)I
    istore      6
    iload_1
    iload_3
    invokestatic e90/setContains(II)I
    istore      7
    iload       4
    iload       7
    iand
    iload       5
    iload       6
    iand
    ior
    ireturn
RETURN_FAILURE:
    iconst_0
    ireturn        
.end method   
              
.method public static main : ([Ljava/lang/String;)V
    .limit stack  100
    .limit locals 100
    
    iconst_0
    istore      6
    bipush      9
    newarray    int
    astore      4        
    iconst_0    
    istore_0 
INITLOOP:
    iload_0
    bipush      9
    if_icmpge    INITLOOP_END
    aload       4
    iload_0
    iinc        0 1
    iload_0
    iload_0
    imul
    iastore
    goto        INITLOOP
INITLOOP_END:
    iconst_0
    istore_0
OUTERLOOP:      
    iload_0     
    sipush      1024
    if_icmpge   OUTERLOOP_END
    iload_0     
    istore_1    
INNERLOOP:      
    iload_1     
    sipush      1024
    if_icmpge   INNERLOOP_END
    iconst_0
    istore      5
INNERMOSTLOOP:
    iload       5
    bipush      9
    if_icmpge   INNERMOSTLOOP_END_SUCCESS
    iload_0
    iload_1
    aload       4
    iload       5
    iaload
    invokestatic e90/check(III)I
    ifeq        INNERMOSTLOOP_END_FAILURE
    iinc        5 1
    goto        INNERMOSTLOOP
INNERMOSTLOOP_END_SUCCESS:
    iinc        6 1
INNERMOSTLOOP_END_FAILURE:
    iinc        1 1
    goto        INNERLOOP
INNERLOOP_END:    
    iinc        0 1
    goto        OUTERLOOP
OUTERLOOP_END:
    iload       6
    invokestatic e90/printInt(I)V  
    return       
.end method       

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