X
X/*
X * Mulitipliziere INT array der Laenge l mit einer INT (d = n * m)
X * return neue Laenge
X */
Xint n_mult( n, m, d, l)
Xregister INT *n;
Xregister INT m;
XINT *d;
X{
X	int i;
X	register LONG mul;
X	
X	for (i=l,mul=0; i; i--) {
X		mul += (LONG)m * (LONG)*n++;
X		*d++ = TOINT(mul);
X		mul /= (MAXINT +1);
X	}
X	
X	if (mul) {		/* carry  ? */
X		l++;
X		*d = mul;
X	}
X	
X	return( l );
X}
X
X/*
X * Mulitipliziere eine NUMBER mit einer INT (d = n * m)
X */
Xvoid a_imult( n, m, d )
XNUMBER *n;
XINT m;
XNUMBER *d;
X{
X	if (m == 0)
X		d->n_len=0;
X	else if (m == 1)
X		a_assign( d, n );
X	else
X		d->n_len = n_mult( n->n_part, m, d->n_part, n->n_len );
X}
X  
X/*
X * Multipliziere zwei NUMBER (d = m1 * m2) 
X */
Xvoid a_mult( m1, m2, d )
XNUMBER *m1,*m2,*d;
X{
X	static VLONG id[ MAXLEN ];		/* Zwischenspeicher	*/
X	register VLONG *vp;			/* Pointer darin	*/
X	register LONG tp1;			/* Zwischenspeicher fuer m1 */
X	register INT *p2;
X	INT *p1;
X	int l1,l2,ld,lc,l,i,j;
X	
X	l1 = m1->n_len;
X	l2 = m2->n_len;
X	l = l1 + l2;
X	if (l >= MAXLEN)
X		abort();
X
X	for (i=l, vp=id; i--;)
X		*vp++ = 0;
X	
X			/* ohne Uebertrag in Zwischenspeicher multiplizieren */
X	for ( p1 = m1->n_part, i=0; i < l1 ; i++, p1++ ) {
X
X		tp1 = (LONG)*p1;		
X		vp = &id[i];
X		for ( p2 = m2->n_part, j = l2; j--;)
X			*vp++ += (VLONG)(tp1 * (LONG)*p2++);
X	}
X
X			/* jetzt alle Uebertraege beruecksichtigen	*/
X	ld = 0;
X	for (lc=0, vp=id, p1=d->n_part; lc++ < l;) {
X		register VLONG tmp;
X
X		tmp = *vp++;
X		if ( *p1++ = TOINT(tmp) )
X			ld = lc;
X		*vp += tmp / ((VLONG)MAXINT +1);
X	}
X	
X	d->n_len = ld;
X}
X
X
X/*
X * Dividiere Zwei NUMBER mit Rest (q= d1 / z2[0] Rest r)
X * z2[i] = z2[0] * 2^i,  i=0..MAXBIT-1
X * r = 0 : kein Rest
X * q = 0 : kein Quotient
X */
Xvoid n_div( d1, z2, q, r )
XNUMBER *d1,*z2,*q,*r;
X{
X	static	NUMBER dummy_rest;  /* Dummy Variable, falls r = 0 */
X	static	NUMBER dummy_quot;  /* Dummy Variable, falla q = 0 */
X	INT *i1,*i1e,*i3;
X	int l2,ld,l,lq;
X#if MAXINT != 1
X	INT z;
X	int pw,l2t;
X#endif
X
X	if (!z2->n_len)
X		abort();
X		
X	if (!r)
X		r = &dummy_rest;
X	if (!q)
X		q = &dummy_quot;
X	
X	a_assign( r, d1 );	/* Kopie von d1 in den Rest		*/
X	
X	l2= z2->n_len;		/* Laenge von z2[0]			*/
X	l = r->n_len - l2;	/* Laenge des noch ''rechts'' liegenden
X					Stuecks von d1			*/
X	lq= l +1;		/* Laenge des Quotienten		*/
X	i3= q->n_part + l;
X	i1= r->n_part + l;
X	ld = l2;		/* aktuelle Laenge des ''Vergleichsstuecks''
X					von d1				*/
X	i1e= i1 + (ld-1);
X	
X	for (; l >= 0; ld++, i1--, i1e--, l--, i3--) {
X		*i3 = 0;
X
X		if (ld == l2 && ! *i1e) {
X			ld--;
X			continue;
X		}
X		
X		if ( ld > l2 || (ld == l2 && n_cmp( i1, z2->n_part, l2) >= 0 ) ) {
X#if MAXINT != 1
X				/* nach 2er-Potenzen zerlegen	*/
X			for (pw=MAXBIT-1, z=(INT)HIGHBIT; pw >= 0; pw--, z /= 2) {
X				if ( ld > (l2t= z2[pw].n_len)
X					|| (ld == l2t
X					    && n_cmp( i1, z2[pw].n_part, ld) >= 0) ) {
X					ld = n_sub( i1, z2[pw].n_part, i1, ld, l2t );
X					(*i3) += z;
X				}
X			}
X#else
X				/* bei MAXINT == 1 alles viel einfacher	*/
X			ld = n_sub( i1, z2->n_part, i1, ld, l2 );
X			(*i3) ++;
X#endif
X		}
X	}
X	
X			/* Korrektur, falls l von Anfang an Negativ war */
X	l ++;
X	lq -= l;
X	ld += l;
X	
X	if (lq>0 && !q->n_part[lq -1])	/* evtl. Laenge korrigieren	*/
X		lq--;
X	
X	q->n_len = lq;
X	r->n_len = ld -1;
X}
X
X/*
X * Dividiere Zwei NUMBER mit Rest (q= d1 / z2[0] Rest r)
X * z2[i] = z2[0] * 2^i,  i=0..MAXBIT-1
X * r = 0 : kein Rest
X * q = 0 : kein Quotient
X */
Xvoid a_div( d1, d2, q, r )
XNUMBER *d1,*d2,*q,*r;
X{
X#if MAXINT != 1
X	NUMBER z2[MAXBIT];
X	INT z;
X	int i;
X	
X	a_assign( &z2[0], d2 );
X	for (i=1,z=2; i < MAXBIT; i++, z *= 2)
X		a_imult( d2, z, &z2[i] );
X	
X	d2 = z2;
X#endif
X
X	n_div( d1, d2, q, r );
X}
X
X/*
X * Dividiere eine NUMBER durch 2
X */
Xvoid a_div2( n )
XNUMBER *n;
X{
X#if MAXBIT == LOWBITS
X	register INT *p;
X	int i;
X
X#if MAXINT != 1
X	register INT h;
X	register int c;
X
X	c=0;
X	i= n->n_len;
X	p= &n->n_part[i-1];
X	
X	for (; i--;) {
X		if (c) {
X			c = (h= *p) & 1;
X			h /= 2;
X			h |= HIGHBIT;
X		}
X		else {
X			c = (h= *p) & 1;
X			h /= 2;
X		}
X		
X		*p-- = h;
X	}
X	
X	if ( (i= n->n_len) && n->n_part[i-1] == 0 )
X		n->n_len = i-1;
X
X#else  /* MAXBIT != 1 */
X	p = n->n_part;
X	i = n->n_len;
X	
X	if (i) {
X		n->n_len = i-1;
X		for (; --i ; p++)
X			p[0] = p[1];
X	}
X#endif /* MAXBIT != 1 */
X#else  /* MAXBIT == LOWBITS */
X	a_div( n, &a_two, n, NUM0P );
X#endif /* MAXBIT == LOWBITS */
X}
X
X
X/*
X *	MODULO-FUNKTIONEN
X */
X
Xstatic NUMBER mod_z2[ MAXBIT ];
X
X/*
X * Init
X */
Xvoid m_init( n, o )
XNUMBER *n,*o;
X{
X	INT z;
X	int i;
X	
X	if (o)
X		a_assign( o, &mod_z2[0] );
X	
X	if (! a_cmp( n, &mod_z2[0] ) )
X		return;
X	
X	for (i=0,z=1; i < MAXBIT; i++, z *= 2)
X		a_imult( n, z, &mod_z2[i] );
X}
X
Xvoid m_add( s1, s2, d )
XNUMBER *s1, *s2, *d;
X{
X	a_add( s1, s2, d );
X	if (a_cmp( d, mod_z2 ) >= 0)
X		a_sub( d, mod_z2, d );
X}
X
Xvoid m_mult( m1, m2, d )
XNUMBER *m1,*m2,*d;
X{
X	a_mult( m1, m2, d );
X	n_div( d, mod_z2, NUM0P, d );
X}
X
X/*
X * Restklassen Exponent
X */
Xvoid m_exp( x, n, z )
XNUMBER *x,*n,*z;
X{
X	NUMBER xt,nt;
X	
X	a_assign( &nt, n );
X	a_assign( &xt, x );
X	a_assign( z, &a_one );
X	
X	while (nt.n_len) {
X		while ( ! (nt.n_part[0] & 1) ) {
X			m_mult( &xt, &xt, &xt );
X			a_div2( &nt );
X		}
X		m_mult( &xt, z, z );
X		a_sub( &nt, &a_one, &nt );
X	}
X}
X
X/*
X * GGT
X */
Xvoid a_ggt( a, b, f )
XNUMBER *a,*b,*f;
X{
X	NUMBER t[2];
X	int at,bt, tmp;
X	
X	a_assign( &t[0], a ); at= 0;
X	a_assign( &t[1], b ); bt= 1;
X	
X	if ( a_cmp( &t[at], &t[bt] ) < 0 ) {
X		tmp= at; at= bt; bt= tmp;
X	}
X				/* euklidischer Algorithmus		*/
X	while ( t[bt].n_len ) {
X		a_div( &t[at], &t[bt], NUM0P, &t[at] );
X		tmp= at; at= bt; bt= tmp;
X	}
X	
X	a_assign( f, &t[at] );
X}
X	
X/*
X * die untersten b bits der Dualdarstellung von n
X * die bits muessen in ein int passen
X */
Xint n_bits(n,b)
XNUMBER *n;
X{
X	INT *p;
X	int l;
X	unsigned r;
X	int m = (1<<b) -1;
X
X	if ( n->n_len == 0)
X		return(0);
X	
X	if (LOWBITS >= b)
X		return( n->n_part[0] & m );
X
X#if LOWBITS != 0
X	l = (b-1) / LOWBITS;
X#else
X	l = n->n_len -1;
X#endif
X	for (p= &n->n_part[l],r=0; l-- >= 0 && b > 0; b-= LOWBITS, p--) {
X		r *= (unsigned)(MAXINT+1);
X		r += (unsigned)*p;
X	}
X
X	return( r & m );
X}
X
X/*
X * Anzahl der bits von n bei Dualdarstellung
X */
Xint n_bitlen( n )
XNUMBER *n;
X{
X	NUMBER b;
X	int i;
X	
X	a_assign( &b, &a_one );
X	
X	for (i=0; a_cmp( &b, n ) <= 0; a_mult( &b, &a_two, &b ), i++)
X		;
X	
X	return(i);
X}
SHAR_EOF
if test 11614 -ne "`wc -c < 'arith.c'`"
then
	echo shar: "error transmitting 'arith.c'" '(should have been 11614 characters)'
fi
fi
echo shar: "extracting 'arith.h'" '(1464 characters)'
if test -f 'arith.h'
then
	echo shar: "will not over-write existing file 'arith.h'"
else
sed 's/^X//' << \SHAR_EOF > 'arith.h'
X/*******************************************************************************
X*									       *
X*	Copyright (c) Martin Nicolay,  22. Nov. 1988			       *
X*									       *
X*	Wenn diese (oder sinngemaess uebersetzte) Copyright-Angabe enthalten   *
X*	bleibt, darf diese Source fuer jeden nichtkomerziellen Zweck weiter    *
X*	verwendet werden.						       *
X*									       *
X*	martin@trillian.megalon.de					       *
X*									       *
X*******************************************************************************/
X
X#ifndef	_arith_h_
X#define	_arith_h_
X
X#ifndef	_conf_h_
X#include	"conf.h"
X#endif
X
Xextern NUMBER a_one,a_two;
X
X/*
X * Prototypes
X */
X
Xvoid	a_add		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	a_assign	P(( NUMBER*, NUMBER* ));
Xint	a_cmp		P(( NUMBER*, NUMBER* ));
Xvoid	a_div		P(( NUMBER*, NUMBER*, NUMBER*, NUMBER* ));
Xvoid	a_div2		P(( NUMBER* ));
Xvoid	a_ggt		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	a_imult		P(( NUMBER*, INT, NUMBER* ));
Xvoid	a_mult		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	a_sub		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	m_init		P(( NUMBER*, NUMBER* ));
Xvoid	m_add		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	m_mult		P(( NUMBER*, NUMBER*, NUMBER* ));
Xvoid	m_exp		P(( NUMBER*, NUMBER*, NUMBER* ));
Xint	n_bits		P(( NUMBER*, int));
Xvoid	n_div		P(( NUMBER*, NUMBER*, NUMBER*, NUMBER* ));
Xint	n_cmp		P(( INT*, INT*, int ));
Xint	n_mult		P(( INT*, INT, INT*, int ));
Xint	n_sub		P(( INT*, INT*, INT*, int, int ));
Xint	n_bitlen	P(( NUMBER* ));
X
X#endif
SHAR_EOF
if test 1464 -ne "`wc -c < 'arith.h'`"
then
	echo shar: "error transmitting 'arith.h'" '(should have been 1464 characters)'
fi
fi
echo shar: "extracting 'conf.h'" '(2050 characters)'
if test -f 'conf.h'
then
	echo shar: "will not over-write existing file 'conf.h'"
else
sed 's/^X//' << \SHAR_EOF > 'conf.h'
X/*******************************************************************************
X*									       *
X*	Copyright (c) Martin Nicolay,  22. Nov. 1988			       *
X*									       *
X*	Wenn diese (oder sinngemaess uebersetzte) Copyright-Angabe enthalten   *
X*	bleibt, darf diese Source fuer jeden nichtkomerziellen Zweck weiter    *
X*	verwendet werden.						       *
X*									       *
X*	martin@trillian.megalon.de					       *
X*									       *
X*******************************************************************************/
X
X#ifndef	_conf_h_
X#define	_conf_h_
X
Xtypedef	unsigned char INT;		/* muss MAXINT fassen		*/
Xtypedef	unsigned int LONG;		/* muss 2*MAXINT+1 fassen	*/
Xtypedef	unsigned long VLONG;		/* muse MAXLEN * (MAXINT^2 +1) fassen*/
X
X#if	defined( M_XENIX )
X#define	P(x)	x			/* Funktions Prototypen an	*/
X#else
X#define	P(x)	()			/* Funktions Prototypen aus	*/
X#endif
X
X/*
X *	(MAXINT+1)-adic Zahlen
X */
X
X/*
X *	MAXINT		Maximale Zahl pro Elemenmt (muss int sein)
X *	MAXBIT		Maximales Bit von MAXINT
X *	LOWBITS		Anzahl der consekutiven low Bits von MAXINT
X *	HIGHBIT		Hoechsten Bit von MAXINT
X *	TOINT		muss (INT)( (x) % MAXINT) ergeben
X *	MAXLEN		Laenge der INT Array in jeder NUMBER
X */
X
X#define MAXINT		0xFF
X
X#if MAXINT == 99
X#define	MAXBIT		7
X#define	LOWBITS 	2
X#endif
X#if MAXINT == 9
X#define	MAXBIT		4
X#define	LOWBITS 	1
X#endif
X#if MAXINT == 1
X#define MAXBIT		1
X#endif
X#if MAXINT == 0xFF
X#define MAXBIT		8
X#define	TOINT(x)	((INT)(x))
X#endif
X
X#ifndef	MAXBIT
X#include	"<< ERROR: MAXBIT must be defined >>"
X#endif
X#ifndef	LOWBITS
X#if MAXINT == (1 << MAXBIT) - 1
X#define	LOWBITS		MAXBIT
X#else
X#include	"<< ERROR: LOWBITS must be defined >>"
X#endif
X#endif
X
X#define	MAXLEN		(300*8/(MAXBIT + 1))
X#define	STRLEN		(MAXLEN*MAXBIT/4)
X#define	HIGHBIT		(1 << (MAXBIT-1) )
X#ifndef	TOINT
X#if LOWBITS == MAXBIT
X#define	TOINT(x)	((INT)((x)&MAXINT))
X#else
X#define	TOINT(x)	((INT)((x)%(MAXINT+1)))
X#endif
X#endif
X
Xtypedef struct {
X	int	n_len;			/* Hoechster benutzter Index	*/
X	INT	n_part[MAXLEN];
X} NUMBER;
X
X#define	NUM0P	((NUMBER *)0)		/* Abkuerzung			*/
X
X#endif
SHAR_EOF
if test 2050 -ne "`wc -c < 'conf.h'`"
then
	echo shar: "error transmitting 'conf.h'" '(should have been 2050 characters)'
fi
fi
echo shar: "extracting 'genprim.c'" '(1416 characters)'
if test -f 'genprim.c'
then
	echo shar: "will not over-write existing file 'genprim.c'"
else
sed 's/^X//' << \SHAR_EOF > 'genprim.c'
X/*******************************************************************************
X*									       *
X*	Copyright (c) Martin Nicolay,  22. Nov. 1988			       *
X*									       *
X*	Wenn diese (oder sinngemaess uebersetzte) Copyright-Angabe enthalten   *
X*	bleibt, darf diese Source fuer jeden nichtkomerziellen Zweck weiter    *
X*	verwendet werden.						       *
X*									       *
X*	martin@trillian.megalon.de					       *
X*									       *