Skip to content

Primo Victoria

September 12, 2011

Few weeks back a co-worker of mine recommended I’d check Project Euler out, so I did. It’s based around an ever-increasing number of mathematical/computational problems. When you’ve solved one of them, you get a link to a forum thread where people share and/or discuss their solutions. Older problems have their threads locked, though; to prevent down-bogging and repetition, one would assume.

The problems are quite interesting but rather difficult. I just assumed my Swedish education was to blame, but after solving some of them I found I had employed more elegant solutions than some – and, to be fair, less elegant than some.

Some of the problems can be solved with brute force but I always suspect there to be less computationally expensive ways. Sometimes that pays off, like in problem 16. In case you are linkophobic, the problem was to calculate the sum of the digits of 21000 when written in a non-potency form.

People well versed in Python or otherwise having a casual relationship to boundless integers did one-liner solutions which just brute-forced it by actually calculating the answer, converting to a string and adding the digits together.

Me, I hadn’t been exposed to bignum libraries, although I guess I should have known they existed somewhere. I had previously dabbled in converting 64-bit numbers to strings using 16-bit code, so I had a hazy understanding of how such a thing would be done. However, writing a bignum library – even a skeletal one – in C didn’t seem clever enough, somehow.

Luckily I quickly realised that since we’re always multiplying by two, we’re really just bitshifting the number. From that, it follows that 21000 is a 1 followed by a thousand zeroes when written in binary. So instead of calculating the number, I could just store it like that from the beginning. That way, the only calculations would be the divisions needed to convert it to decimal.

Thus:


%include "string\str.mac"
%include "io\io.mac"
%include "useful.mac"

org 0x100
    mov bx, 10               ; the radix is ten
    xor bp, bp               ; store the sum of all digits in BP
    segmov es, ds            ; make sure DI and SI are operating in the same segment

    cld
while_not_zero:
    xor dx, dx
    mov si, numbers          ; SI should point to the big number
    mov di, numbers          ; dito DI

    mov cx, 63
.each_digit:
    lodsw                    ; DX:AX / BX -> AX mod DX
    div bx                   ; and put the answer back for the next round
    stosw                    ; modulus will never be more than radix, which we divide by
    loop .each_digit         ; making overflow impossible

    add bp, dx               ; add modulus (= decimal digit) to total

    test ax, ax
    jnz while_not_zero

    word2str s, bp, 10       ; print total to screen
    putsline s
    ret

; we'll store 2^1000 in big-endian for simplicity's sake
numbers        dw    256
times    62    dw    0

%include "string\word2str.asm"
%include "io\putsline.asm"

section .bss
    s        resb 128

(The code is macro- and library-augmented, as you may have noticed.)

I felt rather pleased with that one.

Anyway. Onward to the problem I was actually going to talk about.

I haven’t solved a lot of the problems (just seven so far) but I’ve noticed a curious regularity of prime-related problems. I figured I had to tackle primes sooner or later, so I decided to start with one I perceived as being on the easier end of the scale; problem 7.

I started by reading up on primes a bit. Though easily confused by the Wikipedia articles on the subject, I gathered some information on finding primes. There are several sieves for finding primes, but sieves have some logistical difficulties. They would be more suited for when you want to find all primes below a certain number, since you can’t know (I think; maybe wrong on this) beforehand how many primes will be in a given range.

So I would have to devise some sort of checking algorithm. Between starting to think about this problem and actually solving it, I did another – also prime-related. Even before doing that one, I had suspected that all non-primes can be expressed as the product of primes. Turns out that hunch/forgotten-and-half-remembered-knowledge was right.

How did this help then? Read on, not-yet-asleep reader.

The initial mechanism is fairly straight-forward: have a counter that increases for each run of a loop. Check if that counter is a prime. If it is, add to the list. If not, ignore. Either way, increase counter, redo loop until we have the number of primes we wanted.

From the get-go I had the feeling that checking a number against all previous numbers would be a rather horrible, if working, solution. I can’t remember if it was the first refinement I did, but it’s one of the most obvious ones: only check odd numbers. That is, seed the first position in the primes list with 2 and start the counter at 3, then increase the counter by 2 for every loop. Because no even number can ever be a prime.

But it gets better. Since all non-primes are the product of primes, we only have to check the counter against previously known primes, i.e. the list we so far have compiled. This cuts down the number of checks dramatically.

I had also managed to realise that we don’t even have to check against all primes preceding the counter. I figured that since no number can be evenly divisible by anything greater than half of itself, we could put a limit at counter/2 for checking against.

That’s the point at which I ran and achieved the correct answer.

Upon reading the thread of solutions, I was reminded of something I had read but forgotten: that we can limit ourselves to far lower numbers. Namely, the square root of the counter, inclusive (or plus one, as I did below). I’m still a bit fuzzy on the specifics of this, but apparently it’s true. Incorporating this made for dramatic increases in speed.

Another hindsight tweak was to never start checking against the first prime, 2, since we already skipped all even numbers.

So I arrived at this piece of code:


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

int main()
{
    int *primes;
    int found_primes = 1;
    int currnum = 3;
    int i, isprime, lim;
    int numprimes;

    do
    {
        printf("Number of primes to find: ");
        scanf("%d", &numprimes);
    } while (numprimes < 1);

    primes = (int *) malloc(sizeof(int) * numprimes);
    primes[0] = 2;

    while (found_primes < numprimes)
    {
        isprime = 1;
        lim = sqrt(currnum) + 1;

        // if we're testing against primes larger than the square root of currnum plus one, it must be a prime
        // we can start with the second prime, since we're already skipping all even numbers
        for (i = 1; primes[i] < lim; i++)
        {
            // if modulus is non-zero, it must be a non-prime
            if (currnum % primes[i] == 0)
            {
                isprime = 0;
                break;
            }
        }

        if (isprime)
            primes[found_primes++] = currnum;

        currnum += 2;
    }

    printf("%d: %d\n", numprimes, primes[numprimes - 1]);

//  for (i = 0; i < numprimes; i++)
//      printf("%d: %d\n", i + 1, primes[i]);

    free(primes);
    return 0;
}

Pew pew.

The best part is of course that I now have a way of making primes for future problems.

Advertisements

From → General

Leave a Comment

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: