Warning: Attempt to read property "date" on null in /usr/local/www/websvn.planix.org/blame.php on line 247

Warning: Attempt to read property "msg" on null in /usr/local/www/websvn.planix.org/blame.php on line 247
WebSVN – planix.SVN – Blame – /os/trunk/sys/src/libmp/port/cnfield.c – Rev 26

Subversion Repositories planix.SVN

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
26 7u83 1 
#include "os.h"
2 
#include <mp.h>
3 
#include "dat.h"
4 
 
5 
/*
6 
 * fast reduction for crandall numbers of the form: 2^n - c
7 
 */
8 
 
9 
enum {
10 
	MAXDIG = 1024 / Dbits,
11 
};
12 
 
13 
typedef struct CNfield CNfield;
14 
struct CNfield
15 
{
16 
	Mfield;
17 
 
18 
	mpint	m[1];
19 
 
20 
	int	s;
21 
	mpdigit	c;
22 
};
23 
 
24 
static int
25 
cnreduce(Mfield *m, mpint *a, mpint *r)
26 
{
27 
	mpdigit q[MAXDIG-1], t[MAXDIG], d;
28 
	CNfield *f = (CNfield*)m;
29 
	int qn, tn, k;
30 
 
31 
	k = f->top;
32 
	if((a->top - k) >= MAXDIG)
33 
		return -1;
34 
 
35 
	mpleft(a, f->s, r);
36 
	if(r->top <= k)
37 
		mpbits(r, (k+1)*Dbits);
38 
 
39 
	/* q = hi(r) */
40 
	qn = r->top - k;
41 
	memmove(q, r->p+k, qn*Dbytes);
42 
 
43 
	/* r = lo(r) */
44 
	r->top = k;
45 
	r->sign = 1;
46 
 
47 
	do {
48 
		/* t = q*c */
49 
		tn = qn+1;
50 
		memset(t, 0, tn*Dbytes);
51 
		mpvecdigmuladd(q, qn, f->c, t);
52 
 
53 
		/* q = hi(t) */
54 
		qn = tn - k;
55 
		if(qn <= 0) qn = 0;
56 
		else memmove(q, t+k, qn*Dbytes);
57 
 
58 
		/* r += lo(t) */
59 
		if(tn > k)
60 
			tn = k;
61 
		mpvecadd(r->p, k, t, tn, r->p);
62 
 
63 
		/* if(r >= m) r -= m */
64 
		mpvecsub(r->p, k+1, f->m->p, k, t);
65 
		d = t[k];
66 
		for(tn = 0; tn < k; tn++)
67 
			r->p[tn] = (r->p[tn] & d) | (t[tn] & ~d);
68 
	} while(qn > 0);
69 
 
70 
	if(f->s != 0)
71 
		mpright(r, f->s, r);
72 
	mpnorm(r);
73 
 
74 
	return 0;
75 
}
76 
 
77 
Mfield*
78 
cnfield(mpint *N)
79 
{
80 
	mpint *M, *C;
81 
	CNfield *f;
82 
	mpdigit d;
83 
	int s;
84 
 
85 
	if(N->top <= 2 || N->top >= MAXDIG)
86 
		return nil;
87 
	f = nil;
88 
	d = N->p[N->top-1];
89 
	for(s = 0; (d & (mpdigit)1<<Dbits-1) == 0; s++)
90 
		d <<= 1;
91 
	C = mpnew(0);
92 
	M = mpcopy(N);
93 
	mpleft(N, s, M);
94 
	mpleft(mpone, M->top*Dbits, C);
95 
	mpsub(C, M, C);
96 
	if(C->top != 1)
97 
		goto out;
98 
	f = mallocz(sizeof(CNfield) + M->top*sizeof(mpdigit), 1);
99 
	if(f == nil)
100 
		goto out;
101 
	f->s = s;
102 
	f->c = C->p[0];
103 
	f->m->size = M->top;
104 
	f->m->p = (mpdigit*)&f[1];
105 
	mpassign(M, f->m);
106 
	mpassign(N, f);
107 
	f->reduce = cnreduce;
108 
	f->flags |= MPfield;
109 
out:
110 
	mpfree(M);
111 
	mpfree(C);
112 
 
113 
	return f;
114 
}