(34S) Prime Factors

[Thomas Klemm]

(05-30-2014 08:24 AM)Thomas Klemm Wrote:  Well, it depends on how NEXTP is implemented.

PRIME?

Code:

/* Test if a number is prime or not using a Miller-Rabin test */
#ifndef TINY_BUILD
static const unsigned char primes[] = {
    2, 3, 5, 7, 11, 13, 17, 19,
    23, 29, 31, 37, 41, 43, 47, 53,
};
#define N_PRIMES    (sizeof(primes) / sizeof(unsigned char))
#define QUICK_CHECK (59*59-1)
#endif

int isPrime(unsigned long long int p) {
#ifndef TINY_BUILD
    int i;
    unsigned long long int s;
#define PRIME_ITERATION 15

    /* Quick check for p <= 2 and evens */
    if (p < 2)  return 0;
    if (p == 2) return 1;
    if ((p&1) == 0) return 0;

    /* We fail for numbers >= 2^63 */
    if ((p & 0x8000000000000000ull) != 0) {
        err(ERR_DOMAIN);
        return 1;
    }

    /* Quick check for divisibility by small primes */
    for (i=1; i<N_PRIMES; i++)
        if (p == primes[i])
            return 1;
        else if ((p % primes[i]) == 0)
            return 0;
    if (p < QUICK_CHECK)
        return 1;

    s = p - 1;
    while ((s&1) == 0)
        s /= 2;

    for(i=0; i<PRIME_ITERATION; i++) {
        unsigned long long int temp = s;
        unsigned long long int mod = expmod(primes[i], temp, p);
        while (temp != p-1 && mod != 1 && mod != p-1) {
            mod = mulmod(mod, mod, p);
            temp += temp;
        }
        if(mod!=p-1 && temp%2==0)
            return 0;
    }
#endif
    return 1;
}

NEXTP

Code:

/* Very minimal routine to return the next prime in sequence
 */
        XLBL"NEXTPRIME"
            FLOOR
            x[<=]1?
                JMP prime_2
            PRIME?
                INC X
            EVEN?
prime_loop::            INC X
            PRIME?
                RTN
            INC X
            JMP prime_loop

prime_2::        CLx
            _INT 2
            RTN