imod.c

Arithmetic for integers of almost unlimited size for C and C++. Developed and copyrighted by Ra

/* Integer Version 2.1, RD, 19.7.92 	imod.c			*/

/* Replaced several
	if (carry || mvecgt(sum->val, mv, k))
   by
	if (carry || !mvecgt(mv, sum->val, k))
   (or similar)
   in MasMplM, MplasM, MasMmiM, MmiasM.
   Thanks to [e-mail]
   RD, 18.11.93 						*/
/* Completed dMod, RD, 18.11.93					*/

/* Modular arithmetic with Montgomery's method. Based on the large
   integer arithmetic libI. RD, 16.7.93.			*/

#include <imod.h>
#include <idigit.h>
#include <imem.h>
#include <stdlib.h>
#include <timing.h>


typedef DigitType *mvec;
typedef DigitType *mdoubvec;
typedef void *pointer;

/*********************   special names   **************/

#define mvecadd DigitVecAdd
#define mvecsub DigitVecSub
#define	mvecsr1 DigitVecSr1

/***************   special memory management   **********************/

static mvec 
newmvec(mod)
    ModulusType *mod;
{
    DigitType *u;
    int i;
    if (mod->memsingfree)
    {
	u = (DigitType *) (mod->memsingfree);
	mod->memsingfree = *((pointer *) (mod->memsingfree));
	return u;
    }
    else
    {
	i = mod->length * sizeof(DigitType);
	u = (DigitType *) Imalloc((i > sizeof(void *) ? i : sizeof(void *)));
	if (!u)
	    Merror("newmvec: memory full\n");
	return u;
    }
}				/* newmvec */

static void 
delmvec(u, mod)
    mvec u;
    ModulusType *mod;
{
    pointer *v;
    v = (pointer *) u;
    *v = mod->memsingfree;
    mod->memsingfree = (void *) u;
}				/* delmvec */

static mdoubvec 
newmdoubvec(mod)
    ModulusType *mod;
{
    DigitType *u;
    int i;
    if (mod->memdoubfree)
    {
	u = (DigitType *) (mod->memdoubfree);
	mod->memdoubfree = *((pointer *) (mod->memdoubfree));
	return u;
    }
    else
    {
	i = 2 * mod->length * sizeof(DigitType);
	u = (DigitType *) Imalloc((i > sizeof(void *) ? i : sizeof(void *)));
	if (!u)
	    Merror("newmdoubvec: memory full\n");
	return u;
    }
}				/* newmdoubvec */

static void 
delmdoubvec(u, mod)
    mdoubvec u;
    ModulusType *mod;
{
    pointer *v;
    v = (pointer *) u;
    *v = mod->memdoubfree;
    mod->memdoubfree = (void *) u;
}				/* delmdoubvec */

/***************   special vector functions   **********************/

static void 
mvecas(a, b, k)
    mvec a, b;
    int k;
/* a[k]=b[k]; */
{
    for (; k > 0; k--)
	*a++ = *b++;
}				/* mvecas */

static void 
mvecas0(a, k)
    mvec a;
    int k;
/* a[k]=0; */
{
    for (; k > 0; k--)
	*a++ = 0;
}				/* mvecas0 */

static void 
mvecas1(a, k)
    mvec a;
    int k;
/* a[k]=1; */
{
    *a++ = 1;
    k--;
    for (; k > 0; k--)
	*a++ = 0;
}				/* mvecas1 */

BOOLEAN 
mveceq0(a, k)
    mvec a;
    int k;
/* return a[k]==0; */
{
    for (; k > 0; k--)
	if (*a++)
	    return FALSE;
    return TRUE;
}				/* mveceq0 */

BOOLEAN 
mveceq1(a, k)
    mvec a;
    int k;
/* return a[k]==1; */
{
    if (*a++ != 1)
	return FALSE;
    k--;
    for (; k > 0; k--)
	if (*a++)
	    return FALSE;
    return TRUE;
}				/* mveceq1 */

static BOOLEAN 
mvecgt(a, b, k)
    mvec a, b;
    int k;
/* return a[k]>b[k]; lexikographisch */
{
    for (a += k, b += k; k > 0; k--)
	if (*--a > *--b)
	    return TRUE;
	else if (*a < *b)
	    return FALSE;
    return FALSE;
}				/* mvecgt */


static void 
mveccorsub(diff, b, c, mod, k)
    mvec diff, b, c, mod;
    int k;
/* diff[k]=b[k]-c[k]+mod[k]; */
/*
{	DigitVecAdd(diff, b, mod, k);
	DigitVecSub(diff, diff, c, k);
}
*/
{
    DigitType accu, ac = 0, sc = 0, bb, cc, mm, tmp;
    int i;
    for (i = 0; i < k; i++)
    {
	bb = b[i];
	cc = c[i];
	mm = mod[i];
	/* now : bb + ac + mm - sc - cc */
	accu = bb + ac;
	ac = (accu < bb);
	accu += mm;
	ac += (accu < mm);
	tmp = accu - sc;
	sc = tmp > accu;
	accu = tmp - cc;
	sc += accu > tmp;
	diff[i] = accu;
    }
}				/* mveccorsub */

static void 
mvecmul(prod, b, c, k)
    mdoubvec prod;
    mvec b;
    mvec c;
    int k;
/* prod[2k]=b[k]*c[k]; */
{
    int i;

    prod[k] = DigitVecMult(prod, b, c[0], k);
    for (i = 1; i < k; i++)
    {
	prod[k + i] = DigitVecMultAdd(&prod[i], b, c[i], k);
    }
}				/* mvecmul */

static void 
omvecred(t, n, nprime, k)
    mdoubvec t;
    mvec n;
    DigitType nprime;
    int k;
/* Aendere t so modulo n ab, dass seine niedere Haelfte Null ist. */
{
    DigitType *a, m, c = 0, tmp, aa, c1;
    int j;


    for (j = 0; j < k; j++)
    {
	a = &t[j];
	m = *a * nprime;
	tmp = DigitVecMultAdd(a, n, m, k);
	aa = a[k];
	tmp += aa;
	c1 = tmp < aa;
	tmp += c;
	c = c1 + (tmp < c);
	a[k] = tmp;
    }
    if (c || !mvecgt(n, &t[k], k))
	mvecsub(&t[k], &t[k], n, k);
}				/* omvecred */


static void 
mvaddsr1(r, a, mod, k)
    mvec r, a, mod;
    int k;
/* r[k] = (a[k]+mod[k]) >> 1; */
/* r == a moeglich */
{
    DigitType carry;
    carry = DigitVecAdd(r, a, mod, k);
    DigitVecSr1(r, k);
    if (carry)
	r[k - 1] = r[k - 1] | (1UL << (BitsPerDigit - 1));
}				/* mvaddsr1 */

/*********** static function for random Digits ********/

#define NO_RANDOM_BITS 31

static DigitType 
Prandom()
 /* return  random DigitType */
{
    unsigned int x;
    static BOOLEAN init = FALSE;
    int i;
    DigitType ran;

    if (!init)
    {
	init = TRUE;
	x = timeseed();
	srandom(x);
    }

    ran = random();
    i = NO_RANDOM_BITS;
    while (i < BitsPerDigit)
    {
	ran = (ran << NO_RANDOM_BITS) | random();
	i += NO_RANDOM_BITS;
    }
    return ran;
}				/* Prandom */

/*********** static function for inversion ********/

#define MEVEN(A) (!((A)&1))

static BOOLEAN 
mvecinv(inv, val, mod)
    mvec inv, val;
    ModulusType *mod;
{
    int k;
    BOOLEAN ok;
    mvec swap, a, b, x, u, b0;

    k = mod->length;
    if (mveceq0(val, k))
	return FALSE;
    a = newmvec(mod);
    b = newmvec(mod);
    x = newmvec(mod);
    u = newmvec(mod);
    mvecas(a, val, k);
    b0 = mod->vec;
    mvecas(b, b0, k);
    mvecas1(u, k);		/* a == u*a0 + v*b0; */
    mvecas0(x, k);		/* b == x*a0 + y*b0; */
    while (MEVEN(*a))
    {
	mvecsr1(a, k);
	if (MEVEN(*u))
	    mvecsr1(u, k);
	else
	    mvaddsr1(u, u, b0, k);
    }
    while (TRUE)
    {
	if (mvecgt(b, a, k))
	{
	    swap = b;
	    b = a;
	    a = swap;
	    swap = x;
	    x = u;
	    u = swap;
	}
	mvecsub(a, a, b, k);
	if (!mveceq0(a, k))
	{
	    if (mvecgt(x, u, k))
		mveccorsub(u, u, x, b0, k);
	    else
		mvecsub(u, u, x, k);
	    while (MEVEN(*a))
	    {
		mvecsr1(a, k);
		if (MEVEN(*u))
		    mvecsr1(u, k);
		else
		    mvaddsr1(u, u, b0, k);
	    }
	}
	else
	    break;
    }				/* while(TRUE) */
    ok = mveceq1(b, k);
    mvecas(inv, x, k);
    delmvec(u, mod);
    delmvec(x, mod);
    delmvec(b, mod);
    delmvec(a, mod);
    return ok;
}				/* mvecinv */



/***********************************************************
** EXTERN_FUNCTION(void Merror, (const char *));
** 	Error message
*/
void 
Merror(s)
    const char *s;
{
    fprintf(stderr, "M: %s\n", s);
#ifdef unix
    abort();
#else
    exit(-1);
#endif
}				/* Merror */

/***********************************************************
** EXTERN_FUNCTION(void cMod, (ModulusType *, const Integer *));
**	Creator Modulus
*/
void 
cMod(mod, a)
    ModulusType *mod;
    const Integer *a;
{
    Integer r, rsquare, rcube, u, v, d;
    int i;
    if (Ile0(a))
	Merror("cMod: <=0");
    if (Ieq1(a))
	Merror("cMod: 1");
    if (Ieven(a))
	Merror("cMod: even");
    mod->memsingfree = NULL;
    mod->memdoubfree = NULL;
    mod->length = a->length;
    mod->vec = newmvec(mod);
    mvecas(mod->vec, a->vec, mod->length);
    cIasI(&mod->ModIval, a);