5
\$\begingroup\$

This is the second version of the task. The original task had a defect that the given range of integers was too small. This was pointed out by @harold that other methods couldn't defeat the way of hard-coding all the possible prime divisors. Now you have a larger range of integers to test from 0 to 2000000. Also, the result of your program must be sane for all possible input.

Multiplication now stores the overflowed bits in the r0 register. You can name labels with minimum restriction, and the added dbg instruction dumps the value of all registers.


Write a program in assembly that will find out whether the input is a prime number. The program will run through mjvm which has a simple set of instructions and a determined number of clock cycles for each instruction. The input will be from 0 to 2000000 inclusive, and the objective is to write the most optimized program spending the least number of clock cycles to determine the primality of all inputs.

This is an example program doing a primality test in mjvm assembly.

.     ind  r1 .    ; r1 = input ; originally meant input as decimal
.     lti  r1 2    ; r0 = r1 < 2 ; test less than immediate
.     cjmp false . ; if (r0) goto false
.     eqi  r1 2    ; r0 = r1 == 2
.     cjmp true .  ; if (r0) goto true
.     mv   r2 r1   ; r2 = r1
.     andi r2 1    ; r2 &= 1
.     eqi  r2 0    ; r0 = r2 == 0
.     cjmp false . ; if (r0) goto false
.     mvi  r3 3    ; r3 = 3
loop  mv   r2 r3   ; r2 = r3
.     mul  r2 r2   ; r2 *= r2
.     gt   r2 r1   ; r0 = r2 > r1
.     cjmp true .  ; if (r0) goto true
.     mv   r2 r1   ; r2 = r1
.     div  r2 r3   ; r2 /= r3, r0 = r2 % r3
.     lnot r0 .    ; r0 = !r0 ; logical not
.     cjmp false . ; if (r0) goto false
.     addi r3 2    ; r3 += 2
.     jmp  loop .  ; goto loop
true  mvi  r1 0xffffffff ; r1 = -1
.     outd r1 .    ; print r1 as decimal
.     end  .  .    ; end program
false mvi  r1 0
.     outd r1 .
.     end  .  .

For input from 0 to 2000000, this program consumes total 4043390047 cycles. While some obvious optimizations have been done, there is still a lot of room for improvement.

Now if you're interested, this is the whole list of instructions. Each instruction takes a fixed number of arguments, which can either be a register, an immediate (constant), or a label (address). The type of the argument is also fixed; r means register, i means immediate, and L means label. An instruction with an i suffix takes an immediate as the third argument while the same instruction without an i takes a register as the third argument. Also, I'll explain later, but r0 is a special register.

end - end program
dbg - dumps the value of all registers and pause the program
ind  r - read the input in decimal digits and store the value in `r`
outd r - output the value in `r` as decimal digits
not  r - bitwise not; set `r` to `~r` (all bits reversed)
lnot r - logical not; if `r` is 0, set `r` to ~0 (all bits on); otherwise, set `r` to 0
jmp  L - jump to `L`
cjmp L - if `r0` is non-zero, jump to L

mv(i)  r1 r2(i) - copy `r2(i)` to `r1`
cmv(i) r1 r2(i) - if `r0` is non-zero, copy `r2(i)` to `r1`

eq(i)  r1 r2(i) - if `r1` equals `r2(i)`, set `r0` to ~0; otherwise set `r0` to 0

neq(i), lt(i), le(i), gt(i), ge(i) - analogus to eq(i)

and(i) r1 r2(i) - set `r1` to `r1 & r2(i)`; bitwise `and`

or(i), xor(i) - analogus to and(i)

shl(i), shr(i) - shift left and shift right; analogus to and(i);

add(i), sub(i) - analogus to and(i)

mul(i) r1 r2(i) - from the result of `r1 * r2(i)`, the high bits are stored in `r0`
                  and the low bits are stored in `r1`;
                  if the first argument is `r0`, the low bits are stored in `r0`;
                  the high bits are 0 when there is no overflow

div(i) r1 r2(i) - set `r0` to `r1 % r2(i)` and set `r1` to `r1 / r2(i)`;
                  if the first argument is `r0`, `r0` is `r1 / r2(i)`

As you can see, the result of comparison instructions and the remainder part of division is stored in r0. The value is also used to perform a conditional jump or move.

These instructions consume 2 clock cycles.

cjmp, cmv(i), mul(i)

These instructions consume 32 clock cycles.

div(i)

These instructions consume 10 clock cycles.

ind, outd

dbg does not consume any clock cycle.

All other instructions, including end consume 1 clock cycle.

When the program starts, a cycle is consumed per each instruction of the whole program. That means, for example, if your program has 15 instructions in total, initially 15 cycles will be counted.

There are 16 registers from r0 to rf each holding a 32-bit unsigned integer. The number after r is hexadecimal. The initial value of all registers is 0. There is no memory other than the registers.

The result of an operation after an overflow is guaranteed to be result % pow(2, 32). This also applies to an immediate value outside the range of unsigned 32-bits. Division by zero is undefined, but most likely you'll get a hardware interrupt.

Labels can be any sequence of characters except whitespaces and .. The maximum length of each label is 15 characters, and there can be total 256 labels in a single program.

An immediate operand can be hexadecimal starting with 0x or 0X, or octal starting with 0.

You start the program by reading an input with ind. The input will be automatically given by the VM. You end the program by a sequence of outd and end. The VM will determine whether your final output is correct.

If the input is not prime, the output shall be 0. If the input is prime, the output shall be ~0. 0 and 1 are not prime numbers. All the comparison instructions result as ~0 for a true value, and you can often use this for efficient bitwise logic.

In the assembly code file, each line should start with 4 instructions divided by one or more whitespaces, and . is used as a placeholder. Everything else in the line is ignored and thus can be used as comments.

The current implementation of the VM may have bugs, you can report any faulty code in the comments. An alternative implementation is also welcomed. I will post here if there is any.

This is the code of the VM to measure your program's correctness and efficiency. You can compile with any C compiler, and the program will interpret your assembly file given as the first argument.

#include <stdint.h>
#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>

static uint32_t prime(uint32_t n) {
    if (n < 2) return 0;
    if (n == 2) return -1;
    if (!(n % 2)) return 0;
    for (int d = 3; d * d <= n; d += 2) {
        if (!(n % d)) return 0;
    }
    return -1;
}

enum {
    NOP, DBG, END,
    IND, OUTD,
    NOT, LNOT,
    JMP, CJMP,
    MV, CMV,
    EQ, NEQ, LT, LE, GT, GE,
    AND, OR, XOR, SHL, SHR,
    ADD, SUB, MUL, DIV,
    MVI, CMVI,
    EQI, NEQI, LTI, LEI, GTI, GEI,
    ANDI, ORI, XORI, SHLI, SHRI,
    ADDI, SUBI, MULI, DIVI,
};

_Noreturn static void verr_(const char *m, va_list a) {
    vprintf(m, a);
    putchar('\n');
    exit(1);
}

_Noreturn static void err(const char *m, ...) {
    va_list a;
    va_start(a, m);
    verr_(m, a);
}

static void chk(int ok, const char *m, ...) {
    if (!ok) {
        va_list a;
        va_start(a, m);
        verr_(m, a);
    }
}

#define LBL_DIC_CAP 0x100

typedef struct {
    struct {
        char id[0x10];
        int i;
    } _[LBL_DIC_CAP];
    int n;
} lblDic_t;

static void lblDic_init(lblDic_t *_) {
    _->n = 0;
}

static void lblDic_add(lblDic_t *_, const char *id, int i) {
    chk(_->n < LBL_DIC_CAP, "too many labels");
    strcpy(_->_[_->n].id, id);
    _->_[_->n++].i = i;
}

static int cmp_(const void *a, const void *b) {
    return strcmp(a, b);
}

static void lblDic_sort(lblDic_t *_) {
    qsort(_, _->n, sizeof(*_->_), cmp_);
}

static int lblDic_search(lblDic_t *_, const char *id) {
    int o = &_->_->i - (int *)_;
    int *p = bsearch(id, _, _->n, sizeof(*_->_), cmp_);
    if (!p) return -1;
    return p[o];
}

static int ins(const char *s) {
    if (!strcmp(s, "dbg")) return DBG;
    if (!strcmp(s, "end")) return END;
    if (!strcmp(s, "ind")) return IND;
    if (!strcmp(s, "outd")) return OUTD;
    if (!strcmp(s, "not")) return NOT;
    if (!strcmp(s, "lnot")) return LNOT;
    if (!strcmp(s, "jmp")) return JMP;
    if (!strcmp(s, "cjmp")) return CJMP;
    if (!strcmp(s, "mv")) return MV;
    if (!strcmp(s, "cmv")) return CMV;
    if (!strcmp(s, "eq")) return EQ;
    if (!strcmp(s, "neq")) return NEQ;
    if (!strcmp(s, "lt")) return LT;
    if (!strcmp(s, "le")) return LE;
    if (!strcmp(s, "gt")) return GT;
    if (!strcmp(s, "ge")) return GE;
    if (!strcmp(s, "and")) return AND;
    if (!strcmp(s, "or")) return OR;
    if (!strcmp(s, "xor")) return XOR;
    if (!strcmp(s, "shl")) return SHL;
    if (!strcmp(s, "shr")) return SHR;
    if (!strcmp(s, "add")) return ADD;
    if (!strcmp(s, "sub")) return SUB;
    if (!strcmp(s, "mul")) return MUL;
    if (!strcmp(s, "div")) return DIV;
    if (!strcmp(s, "mvi")) return MVI;
    if (!strcmp(s, "cmvi")) return CMVI;
    if (!strcmp(s, "eqi")) return EQI;
    if (!strcmp(s, "neqi")) return NEQI;
    if (!strcmp(s, "lti")) return LTI;
    if (!strcmp(s, "lei")) return LEI;
    if (!strcmp(s, "gti")) return GTI;
    if (!strcmp(s, "gei")) return GEI;
    if (!strcmp(s, "andi")) return ANDI;
    if (!strcmp(s, "ori")) return ORI;
    if (!strcmp(s, "xori")) return XORI;
    if (!strcmp(s, "shli")) return SHLI;
    if (!strcmp(s, "shri")) return SHRI;
    if (!strcmp(s, "addi")) return ADDI;
    if (!strcmp(s, "subi")) return SUBI;
    if (!strcmp(s, "muli")) return MULI;
    if (!strcmp(s, "divi")) return DIVI;
    err("invalid instruction: %s", s);
}

static int reg(const char *s) {
    chk(*s == 'r', "invalid register: %s", s);
    if ('0' <= s[1] && s[1] <= '9') return s[1] - '0';
    if ('a' <= s[1] && s[1] <= 'f') return s[1] - 'a' + 10;
    err("invalid register: %s", s);
}

int main(int argc, char **argv) {
    chk(argc == 2, "must have 1 argument");
    int sz;
    uint32_t p[0x10000][3];
    uint32_t r[0x10];
    lblDic_t l;
    lblDic_init(&l);
    FILE *f = fopen(argv[1], "r");
    chk((int)f, "failed to open file: %s", argv[1]);
    for (int i = 0;; ++i) {
        char s[0x100];
        int m = fscanf(f, "%s", s);
        if (m < 0) break;
        chk(strlen(s) < 0x10, "%s is too long", s);
        if (*s != '.') {
            if (lblDic_search(&l, s) < 0) {
                lblDic_add(&l, s, i);
                lblDic_sort(&l);
            }
        }
        int c;
        while ((c = getc(f)) != '\n' && m > 0);
    }
    rewind(f);
    char s[4][0x10];
    for (int i = 0;; ++i) {
        int m = fscanf(f, "%s %s %s %s", *s, s[1], s[2], s[3]);
        if (m < 0) {
            sz = i;
            break;
        }
        chk(m == 4, "parse error at line %d", i + 1);
        *p[i] = ins(s[1]);
        if (*p[i] <= END) {
            p[i][1] = NOP;
        } else if (*p[i] == JMP || *p[i] == CJMP) {
            p[i][1] = lblDic_search(&l, s[2]);
            chk(p[i][1] != -1u, "unknown label: %s", s[2]);
        } else {
            p[i][1] = reg(s[2]);
        }
        if (*p[i] <= CJMP) {
            p[i][2] = NOP;
        } else if (*p[i] >= MVI) {
            p[i][2] = strtoul(s[3], NULL, 0);
        } else {
            p[i][2] = reg(s[3]);
        }
        while ((m = getc(f)) != '\n' && m > 0);
    }
    chk(!fclose(f), "sth's very twisted with your system");
    for (int i = 0; i < l.n; ++i) {
        printf("%16s%6d\n", l._[i].id, l._[i].i);
    }
    for (int i = 0; i < sz; ++i) {
        printf("%6d%4u%4u  %u\n", i, *p[i], p[i][1], p[i][2]);
    }
    long long c = 0;
    for (int n = 0; n <= 2000000; ++n) {
        memset(r, 0, sizeof(r));
        uint32_t o = 1;
        c += sz;
        for (uint32_t *i = *p; i < p[sz]; i += 3) {
        start:
            switch(*i) {
                case DBG: {
                    for (int i = 0; i < 010; ++i) {
                        printf("r%x%11u  r%x%11u\n", i, r[i], i + 010, r[i + 010]);
                    }
                    fputs("press ENTER to continue", stdout);
                    getchar();
                    break;
                }
                case END: ++c; goto end;
                case IND: c += 10; r[i[1]] = n; break;
                case OUTD: c += 10; o = r[i[1]]; break;
                case NOT: ++c; r[i[1]] = ~r[i[1]]; break;
                case LNOT: ++c; r[i[1]] = r[i[1]] ? 0 : -1; break;
                case JMP: ++c; i = p[i[1]]; goto start;
                case CJMP: c += 2; if (*r) {i = p[i[1]]; goto start;} break;
                case MV: ++c; r[i[1]] = r[i[2]]; break;
                case CMV: c += 2; if (*r) r[i[1]] = r[i[2]]; break;
                case EQ: ++c; *r = r[i[1]] == r[i[2]] ? -1 : 0; break;
                case NEQ: ++c; *r = r[i[1]] != r[i[2]] ? -1 : 0; break;
                case LT: ++c; *r = r[i[1]] < r[i[2]] ? -1 : 0; break;
                case LE: ++c; *r = r[i[1]] <= r[i[2]] ? -1 : 0; break;
                case GT: ++c; *r = r[i[1]] > r[i[2]] ? -1 : 0; break;
                case GE: ++c; *r = r[i[1]] >= r[i[2]] ? -1 : 0; break;
                case AND: ++c; r[i[1]] &= r[i[2]]; break;
                case OR: ++c; r[i[1]] |= r[i[2]]; break;
                case XOR: ++c; r[i[1]] ^= r[i[2]]; break;
                case SHL: ++c; r[i[1]] <<= r[i[2]]; break;
                case SHR: ++c; r[i[1]] >>= r[i[2]]; break;
                case ADD: ++c; r[i[1]] += r[i[2]]; break;
                case SUB: ++c; r[i[1]] -= r[i[2]]; break;
                case MUL: {
                    c += 2;
                    uint64_t p = (uint64_t)r[i[1]] * r[i[2]];
                    *r = p >> 0x20;
                    r[i[1]] = p;
                    break;
                }
                case DIV: {
                    c += 32;
                    uint32_t rm = r[i[1]] % r[i[2]];
                    uint32_t q = r[i[1]] / r[i[2]];
                    *r = rm;
                    r[i[1]] = q;
                    break;
                }
                case MVI: ++c; r[i[1]] = i[2]; break;
                case CMVI: c += 2; if (*r) r[i[1]] = i[2]; break;
                case EQI: ++c; *r = r[i[1]] == i[2] ? -1 : 0; break;
                case NEQI: ++c; *r = r[i[1]] != i[2] ? -1 : 0; break;
                case LTI: ++c; *r = r[i[1]] < i[2] ? -1 : 0; break;
                case LEI: ++c; *r = r[i[1]] <= i[2] ? -1 : 0; break;
                case GTI: ++c; *r = r[i[1]] > i[2] ? -1 : 0; break;
                case GEI: ++c; *r = r[i[1]] >= i[2] ? -1 : 0; break;
                case ANDI: ++c; r[i[1]] &= i[2]; break;
                case ORI: ++c; r[i[1]] |= i[2]; break;
                case XORI: ++c; r[i[1]] ^= i[2]; break;
                case SHLI: ++c; r[i[1]] <<= i[2]; break;
                case SHRI: ++c; r[i[1]] >>= i[2]; break;
                case ADDI: ++c; r[i[1]] += i[2]; break;
                case SUBI: ++c; r[i[1]] -= i[2]; break;
                case MULI: {
                    c += 2;
                    uint64_t p = (uint64_t)r[i[1]] * i[2];
                    *r = p >> 0x20;
                    r[i[1]] = p;
                    break;
                }
                case DIVI: {
                    c += 32;
                    uint32_t rm = r[i[1]] % i[2];
                    uint32_t q = r[i[1]] / i[2];
                    *r = rm;
                    r[i[1]] = q;
                    break;
                }
            }
        }
    end:
        chk(o == prime(n), "wrong result for %d", n);
    }
    printf("total cycles: %lld\n", c);
    return 0;
}
\$\endgroup\$
11
  • \$\begingroup\$ @harold Do you think keeping neg will make some useful optimizations available? \$\endgroup\$
    – xiver77
    Jan 1 at 1:13
  • \$\begingroup\$ Maybe, but nothing particularly interesting. Reverse-subtract would be more useful in its place, but I'm not actually planning to use it for this task.. \$\endgroup\$
    – harold
    Jan 1 at 1:16
  • \$\begingroup\$ Related \$\endgroup\$
    – Luis Mendo
    Jan 1 at 1:25
  • \$\begingroup\$ Are we allowed to assume the input is between 2 and 10000, and return the wrong answer for anything else? \$\endgroup\$ Jan 1 at 3:27
  • 1
    \$\begingroup\$ @graffe I'm actually devising a task, this time in a real machine with real programming languages. I'm currently testing for adequate ranges. I'll put on sandbox as soon as things are arranged. \$\endgroup\$
    – xiver77
    Jan 6 at 9:54
5
\$\begingroup\$

Stuff at the bottom is for the earlier version of the question.

Here's a simple technique to optimize the basic trial division from the example code, at 2736719982 cycles. It hopefully won't be the fastest (it's a bit boring), but at least it's a new target to beat. Techniques used were:

  • Rewritten loop logic to have the branch at the bottom, and no jmp.

  • Change the trial division base from { 2 } to { 2, 3 }. There are an explicit div-by-3 test, and a special case for small numbers, before the loop to enable that. The alternating increment pattern of 2, 4, 2, 4 .. is achieved with xori r6 6 which changes a 2 to a 4, and a 4 to a 2. The base { 2, 3, 5 } with its increment pattern of [4, 2, 4, 2, 4, 6, 2, 6] by contrast wouldn't be anywhere near that nice to implement.

  • Reuse an outd/end pair by using a jmp, saving code-size that way was worth over a million cycles.

  • Also, the result of your program must be sane for all possible input.

    I did not verify this, so maybe I cheated, who knows. Actually I tried to test more inputs, but that was too slow and I gave up.

Code:

. ind  r1 .
. lti  r1 2 
. cjmp F .  
. eqi  r1 2 
. cjmp T .  
. mv   r2 r1
. andi r2 1
. eqi  r2 0
. cjmp F .
. lei  r1 8
. cjmp T .
. mv   r2 r1
. muli r2 -1431655765
. lei  r2 1431655765
. cjmp F .
. mvi  r3 5
. mvi  r6 4
L mv   r2 r1
. div  r2 r3
. lnot r0 .
. cjmp F .
. xori r6 6
. add  r3 r6
. mv   r2 r3
. mul  r2 r2
. le   r2 r1
. cjmp L .
T mvi  r1 -1
. jmp X .
F mvi  r1 0
X outd r1 .
. end  . .

For previous version of the question:

I have tried a couple of approaches. So far, the one that is winning is a relatively boring one.. namely an unrolled trial division, at 886783 795197 (using Anders Kaseorg's suggestion) cycles:

. ind  r2 .
. eqi  r2 2
. cjmp T .
. mv   r0 r2
. andi r0 1
. lnot r0 .
. cjmp F .
. lei  r2 68
. cjmp A .
. mv   r1 r2
. muli r2 -1431655765
. lei  r2 1431655765
. cjmp F .
. muli r2 1717986919
. lei  r2 858993459
. cjmp F .
. muli r2 -1840700269
. lei  r2 613566756
. cjmp F .
. muli r2 390451573
. lei  r2 390451572
. cjmp F .
. muli r2 1982292599
. lei  r2 330382099
. cjmp F .
. muli r2 1010580541
. lei  r2 252645135
. cjmp F .
. muli r2 -1356305461
. lei  r2 226050910
. cjmp F .
. muli r2 1493901669
. lei  r2 186737708
. cjmp F .
. muli r2 592409283
. lei  r2 148102320
. cjmp F .
. muli r2 -2078209981
. lei  r2 138547332
. cjmp F .
. muli r2 -1741202957
. lei  r2 116080197
. cjmp F .
. muli r2 -104755299
. lei  r2 104755299
. cjmp F .
. muli r2 -1598127365
. lei  r2 99882960
. cjmp F .
. muli r2 1644881093
. lei  r2 91382282
. cjmp F .
. muli r2 -1215556781
. lei  r2 81037118
. cjmp F .
. muli r2 1019144783
. lei  r2 72796055
. cjmp F .
. muli r2 2112278999
. lei  r2 70409299
. cjmp F .
. lti  r1 4489
. cjmp T .
. muli r2 -769247873
. lei  r2 64103989
. cjmp F .
. lti  r1 5041
. cjmp T .
. muli r2 1935759909
. lei  r2 60492497
. cjmp F .
. lti  r1 5329
. cjmp T .
. muli r2 1882725391
. lei  r2 58835168
. cjmp F .
. lti  r1 6241
. cjmp T .
. muli r2 -1196066841
. lei  r2 54366674
. cjmp F .
. lti  r1 6889
. cjmp T .
. muli r2 1397158037
. lei  r2 51746593
. cjmp F .
. lti  r1 7921
. cjmp T .
. muli r2 579096715
. lei  r2 48258059
. cjmp F .
. lti  r1 9409
. cjmp T .
. muli r2 132834041
. lei  r2 44278013
. cjmp F .
T mvi  r1 -1
. outd r1 .
. end  . .
F mvi  r1 0
. outd r1 .
. end  . .
A mvi  r1 1689570158
. shri r2 1
. shr  r1 r2
. andi r1 1
. mvi  r0 0
. sub  r0 r1
. outd r0 .
. end  . .

First I test whether the input is 2 (output True), then whether it is even (output False), then there is a small bitmap of odd primes, and then there is a bunch of divisibility tests by primes from 3 up through 97, which is enough to cover everything up to 10000.

The divisibility tests are basically of the form if (n * mulinv(d) <= 0xFFFFFFFF / d) output False where mulinv is the multiplicative inverse modulo 232, except that instead of using the actual n each time, the result of the multiplication is carried forward and that multiplication is implicitly "undone" by multiplying by d1 * mulinv(d2) where d1 is the previous divisor and d2 is the current one. Some of them have an extra test first to see if the input is smaller than the square of the divisor, and then output True.

Oddly it seems that the main effect of the small "bitmap of primes"-test for small numbers is that it gets rid of some of those extra "input is smaller than square of divisor"-tests, which is not why I added it, but so be it.

I also tried a deterministic Miller Rabin test with bases 2 and 3, which I was very hopeful for, but it was actually quite bad: over 3 million cycles initially, then I added some of those divisibility tests from above so the main MR-test would run less often, but it turned out that the best thing to do was to add all of them .. and then remove the actual MR-test.

\$\endgroup\$
6
  • \$\begingroup\$ I guess I was wrong to set such a small range from 2 to 10000. Your solution seems to be the ideal one with minor possible optimizations if any. I will think what will work as a better range. Your current answer will be mentioned when I make a change to the task. \$\endgroup\$
    – xiver77
    Jan 1 at 4:32
  • \$\begingroup\$ @xiver77 I'm also inclined to blame the small range, but with a caveat: writing the MR test was only doable because all the numbers were below 2^16 (which makes a mul \ div pair feasible for a modular multiplication, the limited range means the product doesn't overflow). Within the current constraints of the VM, a "full range" MR test (or pretty much any other number theory trick) would be quite difficult. \$\endgroup\$
    – harold
    Jan 1 at 4:36
  • \$\begingroup\$ You can save a few cycles by keeping a running product. Instead of computing n * mulinv(d1) followed by n * mulinv(d2), you can multiply the first result by the constant d1 * mulinv(d2). \$\endgroup\$ Jan 1 at 4:41
  • \$\begingroup\$ @AndersKaseorg good, I've added that trick \$\endgroup\$
    – harold
    Jan 1 at 4:49
  • \$\begingroup\$ As a simpler alternative for Miller-Rabin, one could try a simple Fermat test with base 2. The exceptions (A001567) can be covered by dividing by 10 primes (3, 17, 19, 23, 29, 31, 37, 43, 53, 73), or when combined with base 3 test, just 3 primes (5, 7, 37). \$\endgroup\$
    – Bubbler
    Jan 1 at 8:20
3
\$\begingroup\$

hard-coded divisibility test to its extreme - 2156770452 cycles

During an attempt to optimize @harold's solution, I discovered that every time I add an hard-coded divisibility test, the number of cycles decrease. The reason is actually quite simple. A div instruction consumes 32 cycles, but a hard-coded divisibility test can be achieved by a sequence of mv -> muli -> lei -> cjmp, which is total 6 cycles plus 3 cycles for the increased code size. div can only be optimized out when the value of the operand is specified.

The result is a program with 917 instructions, and the few instructions at the end are just there to produce correct output outside the given range of the task. Of course I did not write this by hand, I generated it with a Python script.


Out of curiosity, I wanted to check if the hard-coded division method scales well at least to the range of 32-bit unsigned integers. I found this paper that describes an efficient algorithm (FJ32_256) doing a primality test in the range of 32-bit unsigned integers. I decided to use this algorithm for comparison, then generated a C program (pastebin) with the same Python script tuned to generate C instead of mjvm assembly (pastebin), for a larger range of possible input from \$0\$ to \$2^{32}\$

The test was done through these 5 ranges

[0, 20000000)
[100000000, 120000000)
[200000000, 220000000)
[300000000, 320000000)
[400000000, 420000000)

This is the result where the numbers are the counted (real) clock cycles from my 2.4GHz CPU.

range 0 - hard-coded 1966207882
range 0 - FJ32_256   1896766512

range 1 - hard-coded 1877253992
range 1 - FJ32_256   2118426310

range 2 - hard-coded 1874180510
range 2 - FJ32_256   2196823877

range 3 - hard-coded 5235061704
range 3 - FJ32_256   2244187362

range 4 - hard-coded 5173315658
range 4 - FJ32_256   2266020651

From this non-rigorous test, we can see that hard-coded divisibility test actually performs well up to 220 million, and then the performance drops due to the increased number of divisibility tests for larger numbers. FJ32_256 on the other hand has a very stable performance throughout the whole range.


This is the actual generated assembly code.

. ind r1 .
. lti r1 2
. cjmp F .
. eqi r1 2
. cjmp T .
. mv r0 r1
. andi r0 1
. lnot r0 .
. cjmp F .
. lti r1 9
. cjmp T .
. mv r2 r1
. muli r2 0xaaaaaaab
. lei r2 0x55555555
. cjmp F .
. mv r2 r1
. muli r2 0xcccccccd
. lei r2 0x33333333
. cjmp F .
. mv r2 r1
. muli r2 0xb6db6db7
. lei r2 0x24924924
. cjmp F .
. lti r1 121
. cjmp T .
. mv r2 r1
. muli r2 0xba2e8ba3
. lei r2 0x1745d174
. cjmp F .
. mv r2 r1
. muli r2 0xc4ec4ec5
. lei r2 0x13b13b13
. cjmp F .
. mv r2 r1
. muli r2 0xf0f0f0f1
. lei r2 0xf0f0f0f
. cjmp F .
. mv r2 r1
. muli r2 0x286bca1b
. lei r2 0xd79435e
. cjmp F .
. mv r2 r1
. muli r2 0xe9bd37a7
. lei r2 0xb21642c
. cjmp F .
. mv r2 r1
. muli r2 0x4f72c235
. lei r2 0x8d3dcb0
. cjmp F .
. mv r2 r1
. muli r2 0xbdef7bdf
. lei r2 0x8421084
. cjmp F .
. mv r2 r1
. muli r2 0x914c1bad
. lei r2 0x6eb3e45
. cjmp F .
. mv r2 r1
. muli r2 0xc18f9c19
. lei r2 0x63e7063
. cjmp F .
. mv r2 r1
. muli r2 0x2fa0be83
. lei r2 0x5f417d0
. cjmp F .
. mv r2 r1
. muli r2 0x677d46cf
. lei r2 0x572620a
. cjmp F .
. mv r2 r1
. muli r2 0x8c13521d
. lei r2 0x4d4873e
. cjmp F .
. mv r2 r1
. muli r2 0xa08ad8f3
. lei r2 0x456c797
. cjmp F .
. mv r2 r1
. muli r2 0xc10c9715
. lei r2 0x4325c53
. cjmp F .
. mv r2 r1
. muli r2 0x7a44c6b
. lei r2 0x3d22635
. cjmp F .
. mv r2 r1
. muli r2 0xe327a977
. lei r2 0x39b0ad1
. cjmp F .
. mv r2 r1
. muli r2 0xc7e3f1f9
. lei r2 0x381c0e0
. cjmp F .
. mv r2 r1
. muli r2 0x613716af
. lei r2 0x33d91d2
. cjmp F .
. mv r2 r1
. muli r2 0x2b2e43db
. lei r2 0x3159721
. cjmp F .
. mv r2 r1
. muli r2 0xfa3f47e9
. lei r2 0x2e05c0b
. cjmp F .
. mv r2 r1
. muli r2 0x5f02a3a1
. lei r2 0x2a3a0fd
. cjmp F .
. mv r2 r1
. muli r2 0x7c32b16d
. lei r2 0x288df0c
. cjmp F .
. mv r2 r1
. muli r2 0xd3431b57
. lei r2 0x27c4597
. cjmp F .
. mv r2 r1
. muli r2 0x8d28ac43
. lei r2 0x2647c69
. cjmp F .
. mv r2 r1
. muli r2 0xda6c0965
. lei r2 0x2593f69
. cjmp F .
. mv r2 r1
. muli r2 0xfdbc091
. lei r2 0x243f6f0
. cjmp F .
. lti r1 16129
. cjmp T .
. mv r2 r1
. muli r2 0xefdfbf7f
. lei r2 0x2040810
. cjmp F .
. mv r2 r1
. muli r2 0xc9484e2b
. lei r2 0x1f44659
. cjmp F .
. mv r2 r1
. muli r2 0x77975b9
. lei r2 0x1de5d6e
. cjmp F .
. mv r2 r1
. muli r2 0x70586723
. lei r2 0x1d77b65
. cjmp F .
. mv r2 r1
. muli r2 0x8ce2cabd
. lei r2 0x1b7d6c3
. cjmp F .
. mv r2 r1
. muli r2 0xbf937f27
. lei r2 0x1b20364
. cjmp F .
. mv r2 r1
. muli r2 0x2c0685b5
. lei r2 0x1a16d3f
. cjmp F .
. mv r2 r1
. muli r2 0x451ab30b
. lei r2 0x1920fb4
. cjmp F .
. mv r2 r1
. muli r2 0xdb35a717
. lei r2 0x1886e5f
. cjmp F .
. mv r2 r1
. muli r2 0xd516325
. lei r2 0x17ad220
. cjmp F .
. mv r2 r1
. muli r2 0xd962ae7b
. lei r2 0x16e1f76
. cjmp F .
. mv r2 r1
. muli r2 0x10f8ed9d
. lei r2 0x16a13cd
. cjmp F .
. mv r2 r1
. muli r2 0xee936f3f
. lei r2 0x1571ed3
. cjmp F .
. mv r2 r1
. muli r2 0x90948f41
. lei r2 0x1539094
. cjmp F .
. mv r2 r1
. muli r2 0x3d137e0d
. lei r2 0x14cab88
. cjmp F .
. mv r2 r1
. muli r2 0xef46c0f7
. lei r2 0x149539e
. cjmp F .
. mv r2 r1
. muli r2 0x6e68575b
. lei r2 0x13698df
. cjmp F .
. mv r2 r1
. muli r2 0xdb43bb1f
. lei r2 0x125e227
. cjmp F .
. mv r2 r1
. muli r2 0x9ba144cb
. lei r2 0x120b470
. cjmp F .
. mv r2 r1
. muli r2 0x478bbced
. lei r2 0x11e2ef3
. cjmp F .
. mv r2 r1
. muli r2 0x1fdcd759
. lei r2 0x1194538
. cjmp F .
. mv r2 r1
. muli r2 0x437b2e0f
. lei r2 0x112358e
. cjmp F .
. mv r2 r1
. muli r2 0x10fef011
. lei r2 0x10fef01
. cjmp F .
. mv r2 r1
. muli r2 0x9a020a33
. lei r2 0x105197f
. cjmp F .
. mv r2 r1
. muli r2 0xff00ff01
. lei r2 0xff00ff
. cjmp F .
. mv r2 r1
. muli r2 0x70e99cb7
. lei r2 0xf92fb2
. cjmp F .
. mv r2 r1
. muli r2 0x6205b5c5
. lei r2 0xf3a0d5
. cjmp F .
. mv r2 r1
. muli r2 0xa27acdef
. lei r2 0xf1d48b
. cjmp F .
. mv r2 r1
. muli r2 0x25e4463d
. lei r2 0xec9791
. cjmp F .
. mv r2 r1
. muli r2 0x749cb29
. lei r2 0xe93965
. cjmp F .
. mv r2 r1
. muli r2 0xc9b97113
. lei r2 0xe79372
. cjmp F .
. mv r2 r1
. muli r2 0x84ce32ad
. lei r2 0xdfac1f
. cjmp F .
. mv r2 r1
. muli r2 0xc74be1fb
. lei r2 0xd578e9
. cjmp F .
. mv r2 r1
. muli r2 0xa7198487
. lei r2 0xd2ba08
. cjmp F .
. mv r2 r1
. muli r2 0x39409d09
. lei r2 0xd16154
. cjmp F .
. mv r2 r1
. muli r2 0x6f71de15
. lei r2 0xcebcf8
. cjmp F .
. mv r2 r1
. muli r2 0xbfce8063
. lei r2 0xc5fe74
. cjmp F .
. mv r2 r1
. muli r2 0xf61fe7b1
. lei r2 0xc27806
. cjmp F .
. mv r2 r1
. muli r2 0x70e046d3
. lei r2 0xbcdd53
. cjmp F .
. mv r2 r1
. muli r2 0xf1545af5
. lei r2 0xbbc840
. cjmp F .
. mv r2 r1
. muli r2 0x9a7862a1
. lei r2 0xb9a786
. cjmp F .
. mv r2 r1
. muli r2 0x2a128a57
. lei r2 0xb68d31
. cjmp F .
. mv r2 r1
. muli r2 0xb7747d8f
. lei r2 0xb2927c
. cjmp F .
. mv r2 r1
. muli r2 0xbb5e06dd
. lei r2 0xafb321
. cjmp F .
. mv r2 r1
. muli r2 0x12e9b5b3
. lei r2 0xaceb0f
. cjmp F .
. mv r2 r1
. muli r2 0xec9dbe7f
. lei r2 0xab1cbd
. cjmp F .
. mv r2 r1
. muli r2 0xec41cf4d
. lei r2 0xa87917
. cjmp F .
. mv r2 r1
. muli r2 0xaec02945
. lei r2 0xa513fd
. cjmp F .
. mv r2 r1
. muli r2 0x8382df71
. lei r2 0xa36e71
. cjmp F .
. mv r2 r1
. muli r2 0x84b1c2a9
. lei r2 0xa03c16
. cjmp F .
. mv r2 r1
. muli r2 0x75eb3a0b
. lei r2 0x9c6916
. cjmp F .
. mv r2 r1
. muli r2 0xfa86fe2d
. lei r2 0x9baade
. cjmp F .
. mv r2 r1
. muli r2 0x3f8df54f
. lei r2 0x980e41
. cjmp F .
. mv r2 r1
. muli r2 0x975a751
. lei r2 0x975a75
. cjmp F .
. mv r2 r1
. muli r2 0xc3efac07
. lei r2 0x9548e4
. cjmp F .
. mv r2 r1
. muli r2 0xa8299b73
. lei r2 0x93efd1
. cjmp F .
. mv r2 r1
. muli r2 0x9ba70e41
. lei r2 0x91f5bc
. cjmp F .
. mv r2 r1
. muli r2 0x23d9e879
. lei r2 0x8f67a1
. cjmp F .
. mv r2 r1
. muli r2 0xc494d305
. lei r2 0x8e2917
. cjmp F .
. mv r2 r1
. muli r2 0xab67652f
. lei r2 0x8d8be3
. cjmp F .
. mv r2 r1
. muli r2 0xfb10fe5b
. lei r2 0x8c5584
. cjmp F .
. mv r2 r1
. muli r2 0xbf54fa1f
. lei r2 0x88d180
. cjmp F .
. mv r2 r1
. muli r2 0xb98f81d7
. lei r2 0x869222
. cjmp F .
. mv r2 r1
. muli r2 0xe90f1ec3
. lei r2 0x85797b
. cjmp F .
. mv r2 r1
. muli r2 0xbed87f3b
. lei r2 0x8355ac
. cjmp F .
. mv r2 r1
. muli r2 0x16e70fc7
. lei r2 0x824a4e
. cjmp F .
. mv r2 r1
. muli r2 0x9dece355
. lei r2 0x80c121
. cjmp F .
. mv r2 r1
. muli r2 0x73f62c39
. lei r2 0x7dc9f3
. cjmp F .
. mv r2 r1
. muli r2 0xad46f9a3
. lei r2 0x7d4ece
. cjmp F .
. mv r2 r1
. muli r2 0x24e8d035
. lei r2 0x79237d
. cjmp F .
. mv r2 r1
. muli r2 0x2319bd8b
. lei r2 0x77cf53
. cjmp F .
. mv r2 r1
. muli r2 0xc7ed9da5
. lei r2 0x75a8ac
. cjmp F .
. mv r2 r1
. muli r2 0xfea2c8fb
. lei r2 0x7467ac
. cjmp F .
. mv r2 r1
. muli r2 0xce0f4c09
. lei r2 0x732d70
. cjmp F .
. mv r2 r1
. muli r2 0x544986f3
. lei r2 0x72c62a
. cjmp F .
. mv r2 r1
. muli r2 0x55a10dc1
. lei r2 0x7194a1
. cjmp F .
. mv r2 r1
. muli r2 0x85e33763
. lei r2 0x6fa549
. cjmp F .
. mv r2 r1
. muli r2 0xd84886b1
. lei r2 0x6e8419
. cjmp F .
. mv r2 r1
. muli r2 0x31260967
. lei r2 0x6d68b5
. cjmp F .
. mv r2 r1
. muli r2 0xd1ff25e9
. lei r2 0x6d0b80
. cjmp F .
. mv r2 r1
. muli r2 0x5b84d99f
. lei r2 0x6bf790
. cjmp F .
. mv r2 r1
. muli r2 0x1335df6d
. lei r2 0x6ae907
. cjmp F .
. mv r2 r1
. muli r2 0x75d5add9
. lei r2 0x6a3799
. cjmp F .
. mv r2 r1
. muli r2 0x3c619a43
. lei r2 0x69dfbd
. cjmp F .
. mv r2 r1
. muli r2 0x4767747
. lei r2 0x67dc4c
. cjmp F .
. mv r2 r1
. muli r2 0x663d81
. lei r2 0x663d80
. cjmp F .
. mv r2 r1
. muli r2 0x671ddc2b
. lei r2 0x65ec17
. cjmp F .
. mv r2 r1
. muli r2 0xc1e12337
. lei r2 0x654ac8
. cjmp F .
. mv r2 r1
. muli r2 0x9cd09045
. lei r2 0x645c85
. cjmp F .
. mv r2 r1
. muli r2 0x91496b9b
. lei r2 0x637299
. cjmp F .
. mv r2 r1
. muli r2 0xc7d7b8bd
. lei r2 0x632591
. cjmp F .
. mv r2 r1
. muli r2 0x9f006161
. lei r2 0x6160ff
. cjmp F .
. mv r2 r1
. muli r2 0x5e28152d
. lei r2 0x60cdb5
. cjmp F .
. mv r2 r1
. muli r2 0xbfe803
. lei r2 0x5ff401
. cjmp F .
. mv r2 r1
. muli r2 0x9e907c7b
. lei r2 0x5ed79e
. cjmp F .
. mv r2 r1
. muli r2 0x76528895
. lei r2 0x5d7d42
. cjmp F .
. mv r2 r1
. muli r2 0x1ce2c0d
. lei r2 0x5c6f35
. cjmp F .
. mv r2 r1
. muli r2 0xbed7c42f
. lei r2 0x5b2618
. cjmp F .
. mv r2 r1
. muli r2 0xd4b010e7
. lei r2 0x5a2553
. cjmp F .
. mv r2 r1
. muli r2 0x1ebbe575
. lei r2 0x59686c
. cjmp F .
. mv r2 r1
. muli r2 0xb47b52cb
. lei r2 0x58ae97
. cjmp F .
. mv r2 r1
. muli r2 0x64f3f0d7
. lei r2 0x58345f
. cjmp F .
. mv r2 r1
. muli r2 0x316d6c0f
. lei r2 0x5743d5
. cjmp F .
. mv r2 r1
. muli r2 0x91c1195d
. lei r2 0x5692c4
. cjmp F .
. mv r2 r1
. muli r2 0xa27b1f49
. lei r2 0x561e46
. cjmp F .
. mv r2 r1
. muli r2 0xe508fd01
. lei r2 0x5538ed
. cjmp F .
. mv r2 r1
. muli r2 0x133551cd
. lei r2 0x54c807
. cjmp F .
. mv r2 r1
. muli r2 0x2d8a3f1b
. lei r2 0x5345ef
. cjmp F .
. mv r2 r1
. muli r2 0xc34ad735
. lei r2 0x523a75
. cjmp F .
. mv r2 r1
. muli r2 0xa714919
. lei r2 0x510237
. cjmp F .
. mv r2 r1
. muli r2 0x24eea383
. lei r2 0x50cf12
. cjmp F .
. mv r2 r1
. muli r2 0x42ba771d
. lei r2 0x4fd319
. cjmp F .
. mv r2 r1
. muli r2 0x7772287
. lei r2 0x4fa170
. cjmp F .
. mv r2 r1
. muli r2 0x5e69ddf3
. lei r2 0x4f3ed6
. cjmp F .
. mv r2 r1
. muli r2 0x3b4a6c15
. lei r2 0x4f0de5
. cjmp F .
. mv r2 r1
. muli r2 0xc606b677
. lei r2 0x4e1cae
. cjmp F .
. mv r2 r1
. muli r2 0x46d3e1fd
. lei r2 0x4cd47b
. cjmp F .
. mv r2 r1
. muli r2 0x484a14e9
. lei r2 0x4c78ae
. cjmp F .
. mv r2 r1
. muli r2 0x1ce874d3
. lei r2 0x4c4b19
. cjmp F .
. mv r2 r1
. muli r2 0x473189f
. lei r2 0x4bf093
. cjmp F .
. mv r2 r1
. muli r2 0x372b7e65
. lei r2 0x4aba3c
. cjmp F .
. mv r2 r1
. muli r2 0x4f9e5d91
. lei r2 0x4a6360
. cjmp F .
. mv r2 r1
. muli r2 0x446bd9bb
. lei r2 0x4a383e
. cjmp F .
. mv r2 r1
. muli r2 0xe777c647
. lei r2 0x49e28f
. cjmp F .
. mv r2 r1
. muli r2 0xf61f0c23
. lei r2 0x48417b
. cjmp F .
. mv r2 r1
. muli r2 0xa5cbbb6f
. lei r2 0x47f043
. cjmp F .
. mv r2 r1
. muli r2 0x69daac27
. lei r2 0x474ff2
. cjmp F .
. mv r2 r1
. muli r2 0x637aa061
. lei r2 0x468b6f
. cjmp F .
. mv r2 r1
. muli r2 0x1fb15099
. lei r2 0x45f13f
. cjmp F .
. mv r2 r1
. muli r2 0x712c5825
. lei r2 0x45a522
. cjmp F .
. mv r2 r1
. muli r2 0xff30637b
. lei r2 0x45342c
. cjmp F .
. mv r2 r1
. muli r2 0x1131289
. lei r2 0x44c4a2
. cjmp F .
. mv r2 r1
. muli r2 0xf5acdf7
. lei r2 0x43c5c2
. cjmp F .
. mv r2 r1
. muli r2 0x4d3f89e3
. lei r2 0x437e49
. cjmp F .
. mv r2 r1
. muli r2 0xd2253531
. lei r2 0x43142d
. cjmp F .
. mv r2 r1
. muli r2 0x7bf69fe7
. lei r2 0x42ab5c
. cjmp F .
. mv r2 r1
. muli r2 0xcfb1781f
. lei r2 0x422195
. cjmp F .
. mv r2 r1
. muli r2 0x318e81ed
. lei r2 0x41bbb2
. cjmp F .
. mv r2 r1
. muli r2 0x9f148d11
. lei r2 0x40f391
. cjmp F .
. mv r2 r1
. muli r2 0x2c7a505d
. lei r2 0x40b1e9
. cjmp F .
. mv r2 r1
. muli r2 0x28728f33
. lei r2 0x405064
. cjmp F .
. mv r2 r1
. muli r2 0xe5ec7155
. lei r2 0x403024
. cjmp F .
. mv r2 r1
. muli r2 0x9fe829b7
. lei r2 0x3f90c2
. cjmp F .
. mv r2 r1
. muli r2 0x6a50ca39
. lei r2 0x3f7141
. cjmp F .
. mv r2 r1
. muli r2 0xb6d26aef
. lei r2 0x3f1377
. cjmp F .
. mv r2 r1
. muli r2 0xa8251829
. lei r2 0x3e7988
. cjmp F .
. mv r2 r1
. muli r2 0x1b863613
. lei r2 0x3e5b19
. cjmp F .
. mv r2 r1
. muli r2 0x20d077ad
. lei r2 0x3dc4a5
. cjmp F .
. mv r2 r1
. muli r2 0x2e3d2b97
. lei r2 0x3da6e4
. cjmp F .
. mv r2 r1
. muli r2 0x3dc8eba5
. lei r2 0x3d4e4f
. cjmp F .
. mv r2 r1
. muli r2 0x3229ebbf
. lei r2 0x3c4a6b
. cjmp F .
. mv r2 r1
. muli r2 0x7e01686b
. lei r2 0x3c11d5
. cjmp F .
. mv r2 r1
. muli r2 0x2c086e8d
. lei r2 0x3bf5b1
. cjmp F .
. mv r2 r1
. muli r2 0x9f632df9
. lei r2 0x3bbdb9
. cjmp F .
. mv r2 r1
. muli r2 0xdff892af
. lei r2 0x3b6a88
. cjmp F .
. mv r2 r1
. muli r2 0x3f04d8fd
. lei r2 0x3b183c
. cjmp F .
. mv r2 r1
. muli r2 0x9cfdeff5
. lei r2 0x3aabe3
. cjmp F .
. mv r2 r1
. muli r2 0xbda9d4b
. lei r2 0x3a5ba3
. cjmp F .
. mv r2 r1
. muli r2 0x24f5cbd9
. lei r2 0x3a0c3e
. cjmp F .
. mv r2 r1
. muli r2 0x94cbbb7f
. lei r2 0x38f035
. cjmp F .
. mv r2 r1
. muli r2 0xb4f43b81
. lei r2 0x38d6ec
. cjmp F .
. mv r2 r1
. muli r2 0x34d4323
. lei r2 0x3859cf
. cjmp F .
. mv r2 r1
. muli r2 0x74f5b99b
. lei r2 0x37f741
. cjmp F .
. mv r2 r1
. muli r2 0xc68ea1b5
. lei r2 0x377df0
. cjmp F .
. mv r2 r1
. muli r2 0x96c0cf0b
. lei r2 0x373622
. cjmp F .
. mv r2 r1
. muli r2 0xe33edf99
. lei r2 0x36ef0c
. cjmp F .
. mv r2 r1
. muli r2 0xb897c451
. lei r2 0x36915f
. cjmp F .
. mv r2 r1
. muli r2 0x8fb91695
. lei r2 0x36072c
. cjmp F .
. mv r2 r1
. muli r2 0xe5ea8b41
. lei r2 0x35d9b7
. cjmp F .
. mv r2 r1
. muli r2 0x93fd7cf7
. lei r2 0x359615
. cjmp F .
. mv r2 r1
. muli r2 0xb3f8805
. lei r2 0x35531c
. cjmp F .
. mv r2 r1
. muli r2 0xa912822f
. lei r2 0x353cee
. cjmp F .
. mv r2 r1
. muli r2 0x13a9147d
. lei r2 0x34fad3
. cjmp F .
. mv r2 r1
. muli r2 0xc4c3f21
. lei r2 0x347884
. cjmp F .
. mv r2 r1
. muli r2 0x7a6883c3
. lei r2 0x340dd3
. cjmp F .
. mv r2 r1
. muli r2 0x2ab33855
. lei r2 0x3351fd
. cjmp F .
. mv r2 r1
. muli r2 0x82e6faff
. lei r2 0x333d72
. cjmp F .
. mv r2 r1
. muli r2 0x17be8dab
. lei r2 0x33148d
. cjmp F .
. mv r2 r1
. muli r2 0x7cb91939
. lei r2 0x32d7ae
. cjmp F .
. mv r2 r1
. muli r2 0x42a09ea3
. lei r2 0x32c385
. cjmp F .
. mv r2 r1
. muli r2 0x6d8b8bf1
. lei r2 0x328766
. cjmp F .
. mv r2 r1
. muli r2 0xba8a223d
. lei r2 0x325fa1
. cjmp F .
. mv r2 r1
. muli r2 0x9c3182a7
. lei r2 0x324bd6
. cjmp F .
. mv r2 r1
. muli r2 0x474bcd13
. lei r2 0x32246e
. cjmp F .
. mv r2 r1
. muli r2 0xc5311a97
. lei r2 0x31afa5
. cjmp F .
. mv r2 r1
. muli r2 0xc00c6719
. lei r2 0x319c63
. cjmp F .
. mv r2 r1
. muli r2 0x4b0b61cf
. lei r2 0x3162f7
. cjmp F .
. mv r2 r1
. muli r2 0x9b8e63b1
. lei r2 0x30271f
. cjmp F .
. mv r2 r1
. muli r2 0x41eb667
. lei r2 0x2ff104
. cjmp F .
. mv r2 r1
. muli r2 0x896a76f5
. lei r2 0x2fbb62
. cjmp F .
. mv r2 r1
. muli r2 0x4a55a46d
. lei r2 0x2f7499
. cjmp F .
. mv r2 r1
. muli r2 0x58046447
. lei r2 0x2ed84a
. cjmp F .
. mv r2 r1
. muli r2 0x69be3a81
. lei r2 0x2e832d
. cjmp F .
. mvi r3 1423
L mv r2 r3
. mul r2 r2
. gt r2 r1
. cjmp T .
. mv r2 r1
. div r2 r3
. lnot r0 .
. cjmp F .
. addi r3 2
. jmp L .
T not r4 .
F outd r4 .
. end . .

And this is the Python script used to generate the code, which actually steals a bit of optimized output from gcc.

import sympy
import os

p = []
m = []
c = []
M = sympy.sqrt(2000000)
i = 3
while i < M:
    if (sympy.isprime(i)):
        p.append(i)
    i += 2
f = open("prime.c", "w")
f.write("int main() {}\n")
for i in range(len(p)):
    f.write("unsigned d" + str(p[i]) + "(unsigned a) {return !(a % " +
    str(p[i]) + ");}\n")
f.close()
os.system("gcc -S -masm=intel -O3 -oprime.s prime.c")
f = open("prime.s", "r")
s = f.read()
f.close()

def h(d):
    d = int(d)
    if d < 0:
        return str(hex(d + (1 << 32)))
    else:
        return str(hex(d))

def g(s):
    for i in range(len(s)):
        if s[i] == '-':
            s = s[i:]
            break;
        if s[i].isdigit():
            s = s[i:]
            break;
    d = ""
    i = 0
    while i < len(s):
        if s[i] == '-' or s[i].isdigit():
            d += s[i]
        else:
            break
        i += 1
    return h(d), s[i:]

for i in range(len(p)):
    s = s[s.find("d" + str(p[i]) + ":"):]
    s = s[s.find("imul"):]
    n, s = g(s)
    m.append(n)
    n, s = g(s)
    c.append(n)
f = open("prime.s", "w")
f.write(". ind r1 .\n. lti r1 2\n. cjmp F .\n. eqi r1 2\n. cjmp T .\n" +
". mv r0 r1\n. andi r0 1\n. lnot r0 .\n. cjmp F .\n. lti r1 9\n. cjmp T .\n")
l = 9
for i in range(len(p)):
    if p[i] > l:
        l = p[i] * p[i]
        f.write(". lti r1 " + str(l) + "\n. cjmp T .\n")
    f.write(". mv r2 r1\n. muli r2 " + m[i] + "\n. lei r2 " + c[i] +
    "\n. cjmp F .\n")
f.write(". mvi r3 1423\nL mv r2 r3\n. mul r2 r2\n. gt r2 r1\n. cjmp T .\n" +
". mv r2 r1\n. div r2 r3\n. lnot r0 .\n. cjmp F .\n. addi r3 2\n. jmp L .\n" +
"T not r4 .\nF outd r4 .\n. end . .\n")
f.close()
\$\endgroup\$
3
  • \$\begingroup\$ You only need to test primes up to the square root of the input. Once you've tested the small primes and the square of the next prime to be tested is bigger than the input, you can skip to the end. \$\endgroup\$ Jan 7 at 20:19
  • \$\begingroup\$ @LevelRiverSt I actually thought about that a while after I posted this. That's exactly why there is a bit leap in CPU time going from 200 million to 300 million (the tests happen at 9, 121, 16129, 260467321). On real CPU, adding the test every line gained a lot in the higher ranges but was slower in the lower ranges. I didn't test what works better on this VM, but I just guessed maybe the increased code size won't pay off because the given range is only up to 2 million. \$\endgroup\$
    – xiver77
    Jan 8 at 0:54
  • \$\begingroup\$ Sorry I didn't realise you already had some tests in there. I figured you wouldn't want to do the input size test every time. I would have thought the optimum would be somewhere in the range of one input size test per every 5 to 10 trial divisions. I doubt the test at 9 is doing anything for you. I would say get rid of that, and add several more at higher values. \$\endgroup\$ Jan 8 at 1:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.