Big numbers not working after 18 bits of size.

[trusktr]

I'm trying to implement my own modpow function, and the SchemeNumber.maxIntegerDigits seems to be the limiting factor. Removing the check altogether for the expt function where the "first argument is not finite" error gets raised on line 2629 allows me to continue.

But I'm not getting the results that I expect. If the numbers are smaller, it works, but if the number are bigger but still smaller than 2^53, then it doesn't (but I have plenty of memory).

Here's my function:

function modpow(x, y, n) {
    // I'm passing in strings or schemeNumbers, no natives.
    x = SchemeNumber(x);
    y = SchemeNumber(y);
    n = SchemeNumber(n);
    var result = 1;
    var binaryRepresentationOfExponent = y.toString(2);

    for (var i=0; i<binaryRepresentationOfExponent.length; i++) {
        if (binaryRepresentationOfExponent.charAt(i) == '1') {
            result = fn.mod( fn["*"]( fn['*'](result, result), x), n);
        }
        else {
            result = fn.mod( fn['*'](result, result), n);
        }
    }

    log(" -- "+x.toString(10)+"^"+y.toString(10)+" mod "+n.toString(10)+" RESULT: " + result.toString(10));
    return result;
}

For example, check out the following outputs of calling my function in a loop, one with small numbers (20 bits or less), one with slightly bigger numbers (25 bits or less) and one with numbers of 30 bits or less (not, all of these numbers are schemeNumbers and they are smaller than the native limit of 2^53):

20 bits:

These ones all work fine:

 -- 139730^394726 mod 789453 RESULT: 501337
 -- 726593^479350 mod 958701 RESULT: 57265
 -- 1026907^396933 mod 793867 RESULT: 792128
 -- 891615^442377 mod 884755 RESULT: 120355
 -- 112542^396142 mod 792285 RESULT: 471156
 -- 538465^514881 mod 1029763 RESULT: 847136
 -- 299308^317502 mod 635005 RESULT: 508104
 -- 131778^394712 mod 789425 RESULT: 208716
 -- 992979^509016 mod 1018033 RESULT: 497588
 -- 987591^417173 mod 834347 RESULT: 73920
 -- 668295^395936 mod 791873 RESULT: 789448
 -- 985849^319283 mod 638567 RESULT: 496576
 -- 403846^378031 mod 756063 RESULT: 238681
 -- 165749^316859 mod 633719 RESULT: 125256
 -- 564165^316539 mod 633079 RESULT: 2156
 -- 107509^479499 mod 958999 RESULT: 506072
 -- 697732^495804 mod 991609 RESULT: 912263
 -- 438769^460206 mod 920413 RESULT: 617465
 -- 505279^284142 mod 568285 RESULT: 209211
 -- 145994^426064 mod 852129 RESULT: 717886
 -- 965249^314503 mod 629007 RESULT: 202976
 -- 275169^482295 mod 964591 RESULT: 23600
 -- 303617^445877 mod 891755 RESULT: 454762

25 bits:

 -- 12565866^8783384 mod 17566769 RESULT: 0
 -- 5421817^12289774 mod 24579549 RESULT: 0
 -- 16398908^13503319 mod 27006639 RESULT: 23855104
 -- 3835155^8520408 mod 17040817 RESULT: 14943832
 -- 6718190^12492678 mod 24985357 RESULT: 0
 -- 1412017^12781916 mod 25563833 RESULT: 18176812
 -- 2047445^10608139 mod 21216279 RESULT: 7553024
 -- 12987550^12585703 mod 25171407 RESULT: 20447232
 -- 21146329^13189640 mod 26379281 RESULT: 0
 -- 8139549^12765618 mod 25531237 RESULT: 21980557
 -- 12431709^14369855 mod 28739711 RESULT: 524288
 -- 1152460^16719699 mod 33439399 RESULT: 19005440
 -- 30056251^11849436 mod 23698873 RESULT: 6104173
 -- 32738224^8965879 mod 17931759 RESULT: 0
 -- 20909850^9929454 mod 19858909 RESULT: 12262975
 -- 13969680^9897159 mod 19794319 RESULT: 15728640
 -- 6306037^13763788 mod 27527577 RESULT: 8555211
 -- 12495844^9082092 mod 18164185 RESULT: 11661601
 -- 18876918^11324268 mod 22648537 RESULT: 5523566
 -- 30256478^10196847 mod 20393695 RESULT: 12845056
 -- 6631520^12079300 mod 24158601 RESULT: 23659506

30 bits:

 -- 10270590^356681627 mod 713363255 RESULT: 0
 -- 498866520^472848363 mod 945696727 RESULT: 0
 -- 144969034^278959505 mod 557919011 RESULT: 0
 -- 311114789^456869444 mod 913738889 RESULT: 0
 -- 207729674^366380073 mod 732760147 RESULT: 0
 -- 24903887^509070672 mod 1018141345 RESULT: 0
 -- 534190065^375754658 mod 751509317 RESULT: 0
 -- 305637^317606645 mod 635213291 RESULT: 419430400
 -- 59483847^473175270 mod 946350541 RESULT: 0
 -- 281579956^313280274 mod 626560549 RESULT: 0
 -- 264150281^297390241 mod 594780483 RESULT: 0
 -- 36209975^276812543 mod 553625087 RESULT: 0
 -- 409668651^393514895 mod 787029791 RESULT: 0
 -- 406661344^283252541 mod 566505083 RESULT: 0
 -- 450666631^398508318 mod 797016637 RESULT: 0
 -- 288335035^320166021 mod 640332043 RESULT: 0
 -- 524629248^518357032 mod 1036714065 RESULT: 0
 -- 155879269^460687891 mod 921375783 RESULT: 0
 -- 56080000^421577698 mod 843155397 RESULT: 0

As you can see, the bigger the numbers get, then the results start to be 0. At around 500 bits, the results start to be NaN instead of 0.

WHat I did was comment out the following (line 2628):

            //if (!SN_isFinite(x))
                //raise("&assertion", "div/mod first argument is not finite", x);

I also tried setting SchemeNumber.maxIntegerDigits = 1e1000000000000; without uncommenting those two lines, but it'd still complain about argument one not being finite.

Might I be doing something wrong in my code? Do I need to do something to lift size restrictions?

[trusktr]

Oops, I just realized that the output is not even correct. I'm not getting expected values. Let me see here...

[trusktr]

Ok, so I've done some further investigation, and I think there might be something not working as it should. I couldn't get my program (a random prime generator, see below) working in JavaScript, so I rewrote it in C, fixed a few programmatic errors I made, then came back to the JavaScript version and fixed it up. The js version works now, but only with numbers up to 18 bits long. At 18 bits, the prime is found almost instantaneously, but starting at 19 bits (much less than JavaScript's 53 bits), the program seems to run indefinitely.

Try it out. Below are both the C and JavaScript versions. The C version works with any number of bits (given you have the memory available). Note that when I say "bits" when referring to the JavaScript version that I simply mean theoretical bits that represent a big number, not actual bits like when referring to the native C version.

To run the C version, pass it a single argument at the command line: the number of bits the prime will have.

To run the JavaScript version, you'd have paste biginteger.js and schemeNumber.js at the top of the code, then npm install bigint-node in the same directory as the file (or install it globally with npm install -g bigint-node). I'm temporarily using bigint-node to simply generate random big integers, then pass their string representations to SchemeNumber(). After all that, you can execute it with Node.js at the command line and it will prompt you for the number of bits (it won't work in a browser): node prime.js.

Try various sizes, and you'll get random primes of the specified size (with a probability of at most 1/(2^128) of generating a non-prime number). At some specific size (19 in my case) it locks up. Note that I've commented out 2629 of schemeNumber.js.

In the meantime, I'll keep investigating the cause of the problem.

C

#include <openssl/bn.h>
#include <stdbool.h>


// BEGIN helper functions.
BIGNUM *mymul(BIGNUM *a, BIGNUM *b) {
    BIGNUM *result = BN_new();
    BN_CTX *ctx = BN_CTX_new();

    BN_mul(result, a, b, ctx);

    return result;
}
BIGNUM *mydiv(BIGNUM *a, BIGNUM *b) {
    BIGNUM *result = BN_new();
    BIGNUM *remainder = BN_new();
    BN_CTX *ctx = BN_CTX_new();

    BN_div(result, remainder, a, b, ctx);

    return result;
}
BIGNUM *myadd(BIGNUM *a, BIGNUM *b) {
    BIGNUM *result = BN_new();

    BN_add(result, a, b);

    return result;
}
BIGNUM *mysub(BIGNUM *a, BIGNUM *b) {
    BIGNUM *result = BN_new();

    BN_sub(result, a, b);

    return result;
}
BIGNUM *mymod(BIGNUM *a, BIGNUM *b) {
    BIGNUM *result = BN_new();
    BN_CTX *ctx = BN_CTX_new();

    BN_mod(result, a, b, ctx);

    return result;
}
BIGNUM *myrshift(BIGNUM *a, int n) {
    BIGNUM *result = BN_new();
    BN_CTX *ctx = BN_CTX_new();

    BN_rshift(result, a, n);

    return result;
}
bool myequal(BIGNUM *a, BIGNUM *b) {
    int result;
    result = BN_cmp(a, b);

    if (result == 0)
        return true;
    else
        return false;
}
// END helper functions.

BIGNUM * modpow2(BIGNUM *x, BIGNUM *y, BIGNUM *n) {
    x = mymod(x, n);

    BIGNUM *one = BN_new(); BN_one(one);
    BIGNUM *two = BN_new(); BN_add(two, one, one);

    if ( myequal(y, one) ) return x;

    if ( !BN_is_odd(y) )
        return mymod( modpow2( mymul(x, x), mydiv( y, two ), n), n);

    else
        return mymod( mymul(x, modpow2( mymul(x, x), mydiv( mysub(y, one), two ), n) ), n);
}

bool isPrime(BIGNUM *number /*p*/, long sizeOfNumberInBits) {

    if (!BN_is_odd(number)) return false;

    bool carmichael = true;
    BIGNUM *x = BN_new();
    BIGNUM *y = BN_new();

    // for i=1 to 128
    int i;
    for (i = 1; i <= 128; i++) {

        // let x E Zp*
        BIGNUM *zero = BN_new(); BN_zero(zero);
        do {
            BN_rand(x, (int)sizeOfNumberInBits, -1, 0); // random integer up to sizeOfNumberInBits bits long
        } while( myequal(x, zero) || BN_cmp(x, number) > 0 || myequal(x, number) );

        // let y = x^((p-1)/2) mod p
        BIGNUM *one = BN_new(); BN_one(one);
        BIGNUM *two = BN_new(); BN_add(two, one, one);
        y = modpow2(x, mydiv(mysub(number, one), two), number);

        // if y != +-1 return false

        if ( !myequal(y, one) && !myequal(y, mysub(number, one) ) ) {
            return false;
        }

        // if y == -1 AKA if y == p-1
        if ( myequal(y, mysub(number, one)) ) carmichael = false;

    // endfor
    }

    return !carmichael;
}

int main(int argc, char* args[]) {
    long numberOfBits = strtol(args[1], NULL, 0);

    printf(" -- number of bits: %i", (int)numberOfBits);

    BIGNUM *number = BN_new();
    do {
        BN_rand(number, (int)numberOfBits, 0, 1);
    } while (!isPrime(number, numberOfBits));

    printf("\n -- Random prime: %s", BN_bn2dec(number));
    return 0;
}

/* Sample output:
 *   ┌─[6:26:57pm/starlancer/trusktr/.../csus_csc152_cryptography/hwk5]
 *   └─╼ genprime() { date && ~/school/csus_csc152_cryptography/hwk5/prime-generator $1; echo; date }
 *
 *   ┌─[6:27:00pm/starlancer/trusktr/.../csus_csc152_cryptography/hwk5]
 *   └─╼ genprime 1024
 *   Wed Dec  4 18:27:01 PST 2013
 *   -- number of bits: 1024
 *   -- Random prime: 141084237503186237301497483293869589783995261992058764171201648326
 *   34526693457266069335781966954751879788733662406265295835505969138750591939516861273
 *   62572339732737730005056561528437974140450298442914896605426555201041598593174235024
 *   83437152132662069886285634899820709143310848051446769098902187551242933816643
 *   Wed Dec  4 18:27:04 PST 2013
 */

JavaScript

SchemeNumber.maxIntegerDigits = "1e1000000000000";

var BigInt = require("bigint-node"); // a big int library
var fn = SchemeNumber.fn; // a big number library
var Readline = require("readline"); // to read input from a stream (e.g. stdin) and put data into a stream (e.g. stdout)

log(" -= Random Prime Generator =- ");

var prompt = Readline.createInterface({
    input: process.stdin,
    output: process.stdout
});

var numberOfBits = 0;
prompt.question(" -- Pick the size of the prime in bits: ", function(answer) {
    numberOfBits = parseInt(answer);
    prompt.close();

    var number = 0;

    do {
        number = BigInt.Random(numberOfBits-1).toStr(10); // random numberOfBits-bit number from bigint-node library.
        number = fn["+"](number, fn.expt("2", ""+(numberOfBits-1)) ); // make the most significant bit equal to 1 while converting to javascript-bignum library format.
        if (number.toString(2).length < numberOfBits) {
            log(" -- WTF: " + number.toString(2));
        }
    } while (!fn['odd?'](number) || !isPrime(number, numberOfBits));

    //number = modpow("8347987234234", "38847792792433", "384197717433");

    log(" -- Your prime number is: " + fn['number->string'](number));
});

function isPrime(number /*p*/, sizeOfNumberInBits) {
    number = SchemeNumber(number);
    if (fn["even?"](number)) return false;

    var carmichael = true;
    var x; var y;

    // for i=1 to 128
    for (var i = 1; i <= 128; i++) {

        // let x E Zp*
        do {
            x = SchemeNumber(BigInt.Random(sizeOfNumberInBits).toStr(10)); // random integer up to sizeOfNumberInBits bits long
        } while( fn["="](x, 0) || fn[">="](x, number) );

            //log("\n -- number: " + fn["number->string"](number));

        // let y = x^((p-1)/2) mod p
        y = modpow(x, fn['/']( fn['-'](number, 1), 2 ), number);

            //log(" -- "+x.toString(10)+" ^ (("+number.toString(10)+"-1)/2) = "+y.toString(10));
            //log(" -- number: "+number.toString(10) );
            //log(" -- -1: "+fn["-"](number, 1) );
            //log(" -- "+y.toString(10)+"==1: "+fn["="](y, 1));
            //log(" -- "+y.toString(10)+"==-1: "+fn["="](y, fn["-"](number, 1) ) );

        // if y != +-1 return false
        if ( !fn['='](y, 1) && !fn['='](y, fn["-"](number, 1)) ) {
            //log(" -- y==1: "+fn["="](y, 1));
            //log(" -- y==-1: "+fn["="](y, fn["-"](number, 1) ) );
            return false;
        }

        if ( fn["="](y, fn["-"](number, 1)) ) carmichael = false;

    // endfor
    }

    return !carmichael;
}

// custom exponentiation functions for this assignment
function pow(x, y) {
    x = new SchemeNumber(x);
    y = new SchemeNumber(y);

    if ( fn['='](y, 1) ) return x;
    if ( fn['even?'](y) ) return pow(fn['*'](x, x), fn["/"](y, 2));
    else return fn["*"]( x, pow( fn["*"](x,x), fn['/']( fn['-'](y,1), 2) ) );
}

function modpow(x, y, n) {
    x = SchemeNumber(x);
    y = SchemeNumber(y);
    n = SchemeNumber(n);
    var result = 1;
    var binaryRepresentationOfExponent = y.toString(2);

    for (var i=0; i<binaryRepresentationOfExponent.length; i++) {
        if (binaryRepresentationOfExponent.charAt(i) == '1') {
            // result = x * result^2 mod n
            result = fn.mod( fn["*"]( fn['*'](result, result), x), n);
        }
        else {
            // result = result^2 mod n
            result = fn.mod( fn['*'](result, result), n);
        }
    }

    return result;
}

function log(value) {
    console.log(value);
}