Solovay-Strassen primality test

[Ulf M�ller]

In case anyone cares, here's an implementation of the Solovay-Strassen
primality test.

#include <openssl/bn.h>
#include "bn_prime.h"

/* return Jacobi symbol, or -2 on error */
int BN_jacobi_symbol(BIGNUM *A, BIGNUM *N, BN_CTX *ctx)
        {
        BIGNUM *a, *n, *r, *m4, *m8, *tmp;
        unsigned long e;
        int ok = 0, s = 1;

        a = &(ctx->bn[ctx->tos++]);
        n = &(ctx->bn[ctx->tos++]);
        r = &(ctx->bn[ctx->tos++]);
        m4 = &(ctx->bn[ctx->tos++]);
        m8 = &(ctx->bn[ctx->tos++]);

        if (!BN_set_word(m4, 4)) goto err;
        if (!BN_set_word(m8, 8)) goto err;
        if (!BN_copy(a, A)) goto err;
        if (!BN_copy(n, N)) goto err;

        for (;;)
                {
                if (BN_is_zero(a))
                        {
                        s = 0;
                        break;
                        }
                if (BN_is_one(a))
                        break;
        
                e = 0;
                while (!BN_is_odd(a))
                        {
                        e++;
                        if (!BN_rshift1(a, a)) goto err;
                        }

                if (e % 2)
                        {
                        if (!BN_mod(r, n, m8, ctx)) goto err;
                        if (!BN_is_one(r) && !BN_is_word(r, 7))
                                s = -s;
                        }

                if (!BN_mod(r, n, m4, ctx)) goto err;
                if (BN_is_word(r, 3))
                        {
                        if (!BN_mod(r, a, m4, ctx)) goto err;
                        if (BN_is_word(r, 3))
                                s = -s;
                        }

                if (!BN_is_one(a))
                        if (!BN_mod(n, n, a, ctx)) goto err;

                tmp=n; n=a; a=tmp;
                }
        ok = 1;

 err:
        ctx->tos -= 5;
        return (ok ? s : -2);
        }

static int witness(BIGNUM *a, BIGNUM *n, BN_CTX *ctx, BN_CTX *ctx2)
        {
        int s,ret= -1;
        BIGNUM *r,*n1;

        r = &(ctx->bn[ctx->tos++]);
        n1 = &(ctx->bn[ctx->tos++]);

        if (!BN_copy(n1,n)) goto err;
        if (!BN_rshift1(n1,n1)) goto err; /* n1 = (n-1)/2 */

        if (!BN_mod_exp(r, a, n1, n, ctx2)) goto err;

        if (!BN_add_word(r,1)) goto err;
        if (!BN_is_word(r,2) && BN_cmp(r,n) != 0)
                ret = 1;
        if (ret != 1)
                {
                if ((s = BN_jacobi_symbol(a, n, ctx)) == -2) goto err;
                if ((s == 1 && BN_is_word(r,2)) ||
                    (s == -1 && BN_cmp(r,n) == 0))
                        ret = 0;
                else
                        ret = 1;
                }
 err:
        ctx->tos -= 2;
        return(ret);
        }

int BN_is_prime_solovay_strassen(BIGNUM *a, int checks,
                void (*callback)(int,int,void *),
                BN_CTX *ctx_passed, BN_CTX *ctx2_passed, void *cb_arg,
                int do_trial_division)
        {
        int i,j,ret= -1;
        BIGNUM *check;
        BN_CTX *ctx=NULL,*ctx2=NULL;

        if (checks == BN_prime_checks)
                {
                checks = 50;
                }

        if (!BN_is_odd(a))
                return(0);
        if (do_trial_division)
                {
                for (i = 1; i < NUMPRIMES; i++)
                        if (BN_mod_word(a, primes[i]) == 0) 
                                return 0;
                if (callback != NULL) callback(1,-1,cb_arg);
                }

        if (ctx_passed != NULL)
                ctx=ctx_passed;
        else
                if ((ctx=BN_CTX_new()) == NULL) goto err;
        if (ctx2_passed != NULL)
                ctx2=ctx2_passed;
        else
                if ((ctx2=BN_CTX_new()) == NULL) goto err;

        check= &(ctx->bn[ctx->tos++]);

        for (i=0; i<checks; i++)
                {
                if (!BN_pseudo_rand(check,BN_num_bits(a),0,0)) goto err;
                if (BN_cmp(check, a) >= 0)
                        BN_sub(check, check, a);
                if (BN_is_zero(check)) BN_one(check);
                j=witness(check,a,ctx,ctx2);
                if (j == -1) goto err;
                if (j)
                        {
                        ret=0;
                        goto err;
                        }
                if (callback != NULL) callback(1,i,cb_arg);
                }
        ret=1;
err:
        ctx->tos--;
        if ((ctx_passed == NULL) && (ctx != NULL))
                BN_CTX_free(ctx);
        if ((ctx2_passed == NULL) && (ctx2 != NULL))
                BN_CTX_free(ctx2);
        return(ret);
        }
______________________________________________________________________
OpenSSL Project                                 http://www.openssl.org
Development Mailing List                       [EMAIL PROTECTED]
Automated List Manager                           [EMAIL PROTECTED]