Subversion Repositories planix.SVN

Rev

Rev 33 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include "os.h"
2 
#include <mp.h>
3 
#include "dat.h"
4 
 
5 
// res = b**e
6 
//
7 
// knuth, vol 2, pp 398-400
8 
 
9 
enum {
10 
	Freeb=	0x1,
11 
	Freee=	0x2,
12 
	Freem=	0x4,
13 
};
14 
 
15 
//int expdebug;
16 
 
17 
void
18 
mpexp(mpint *b, mpint *e, mpint *m, mpint *res)
19 
{
20 
	mpint *t[2];
21 
	int tofree;
22 
	mpdigit d, bit;
23 
	int i, j;
24 
 
33 7u83 25 
	assert(m == nil || m->flags & MPnorm);
26 
	assert((e->flags & MPtimesafe) == 0);
27 
	res->flags |= b->flags & MPtimesafe;
28 
 
2 - 29 
	i = mpcmp(e,mpzero);
30 
	if(i==0){
31 
		mpassign(mpone, res);
32 
		return;
33 
	}
34 
	if(i<0)
35 
		sysfatal("mpexp: negative exponent");
36 
 
37 
	t[0] = mpcopy(b);
38 
	t[1] = res;
39 
 
40 
	tofree = 0;
41 
	if(res == b){
42 
		b = mpcopy(b);
43 
		tofree |= Freeb;
44 
	}
45 
	if(res == e){
46 
		e = mpcopy(e);
47 
		tofree |= Freee;
48 
	}
49 
	if(res == m){
50 
		m = mpcopy(m);
51 
		tofree |= Freem;
52 
	}
53 
 
54 
	// skip first bit
55 
	i = e->top-1;
56 
	d = e->p[i];
57 
	for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
58 
		;
59 
	bit >>= 1;
60 
 
61 
	j = 0;
62 
	for(;;){
63 
		for(; bit != 0; bit >>= 1){
33 7u83 64 
			if(m != nil)
65 
				mpmodmul(t[j], t[j], m, t[j^1]);
2 - 66 
			else
33 7u83 67 
				mpmul(t[j], t[j], t[j^1]);
68 
			if(bit & d) {
69 
				if(m != nil)
70 
					mpmodmul(t[j^1], b, m, t[j]);
71 
				else
72 
					mpmul(t[j^1], b, t[j]);
73 
			} else
2 - 74 
				j ^= 1;
75 
		}
76 
		if(--i < 0)
77 
			break;
78 
		bit = mpdighi;
79 
		d = e->p[i];
80 
	}
81 
	if(t[j] == res){
82 
		mpfree(t[j^1]);
83 
	} else {
84 
		mpassign(t[j], res);
85 
		mpfree(t[j]);
86 
	}
87 
 
88 
	if(tofree){
89 
		if(tofree & Freeb)
90 
			mpfree(b);
91 
		if(tofree & Freee)
92 
			mpfree(e);
93 
		if(tofree & Freem)
94 
			mpfree(m);
95 
	}
96 
}