shithub: femtolisp

ref: 76a7c66de3702ff0486900f7d4fa626407b953c3
dir: /3rd/mp/mpexp.c/

#include "platform.h"
#include "mp.h"

// res = b**e
//
// knuth, vol 2, pp 398-400

enum {
	Freeb=	0x1,
	Freee=	0x2,
	Freem=	0x4,
};

//int expdebug;

void
mpexp(mpint *b, mpint *e, mpint *m, mpint *res)
{
	mpint *t[2];
	int tofree;
	mpdigit d, bit;
	int i, j;

	assert(m == nil || m->flags & MPnorm);
	assert((e->flags & MPtimesafe) == 0);
	res->flags |= b->flags & MPtimesafe;

	i = mpcmp(e,mpzero);
	if(i==0){
		mpassign(mpone, res);
		return;
	}
	if(i < 0){
		fprintf(stderr, "mpexp: negative exponent\n");
		exit(2);
	}

	t[0] = mpcopy(b);
	t[1] = res;

	tofree = 0;
	if(res == b){
		b = mpcopy(b);
		tofree |= Freeb;
	}
	if(res == e){
		e = mpcopy(e);
		tofree |= Freee;
	}
	if(res == m){
		m = mpcopy(m);
		tofree |= Freem;
	}

	// skip first bit
	i = e->top-1;
	d = e->p[i];
	for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
		;
	bit >>= 1;

	j = 0;
	for(;;){
		for(; bit != 0; bit >>= 1){
			if(m != nil)
				mpmodmul(t[j], t[j], m, t[j^1]);
			else
				mpmul(t[j], t[j], t[j^1]);
			if(bit & d) {
				if(m != nil)
					mpmodmul(t[j^1], b, m, t[j]);
				else
					mpmul(t[j^1], b, t[j]);
			} else
				j ^= 1;
		}
		if(--i < 0)
			break;
		bit = mpdighi;
		d = e->p[i];
	}
	if(t[j] == res){
		mpfree(t[j^1]);
	} else {
		mpassign(t[j], res);
		mpfree(t[j]);
	}

	if(tofree){
		if(tofree & Freeb)
			mpfree(b);
		if(tofree & Freee)
			mpfree(e);
		if(tofree & Freem)
			mpfree(m);
	}
}