Annotation of OpenXM_contrib2/asir2000/engine/up_gf2n.c, Revision 1.1.1.1
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/up_gf2n.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include <math.h>
4:
5: extern int debug_up;
6: extern int up_lazy;
7: extern GEN_UP2 current_mod_gf2n;
8:
9: void squarep_gf2n(vl,n1,nr)
10: VL vl;
11: P n1;
12: P *nr;
13: {
14: UP b1,br;
15:
16: if ( !n1 )
17: *nr = 0;
18: else if ( OID(n1) == O_N )
19: mulp(vl,n1,n1,nr);
20: else {
21: ptoup(n1,&b1);
22: squareup_gf2n(b1,&br);
23: uptop(br,nr);
24: }
25: }
26:
27: void squareup_gf2n(n1,nr)
28: UP n1;
29: UP *nr;
30: {
31: UP r;
32: GF2N *c1,*c;
33: int i,d1,d;
34:
35: if ( !n1 )
36: *nr = 0;
37: else if ( !n1->d ) {
38: *nr = r = UPALLOC(0); r->d = 0;
39: squaregf2n((GF2N)n1->c[0],(GF2N *)(&r->c[0]));
40: } else {
41: d1 = n1->d;
42: d = 2*d1;
43: *nr = r = UPALLOC(d); r->d = d;
44: c1 = (GF2N *)n1->c; c = (GF2N *)r->c;
45: bzero((char *)c,(d+1)*sizeof(GF2N *));
46: for ( i = 0; i <= d1; i++ )
47: squaregf2n(c1[i],&c[2*i]);
48: }
49: }
50:
51: /* x^(2^n) mod f */
52:
53: void powermodup_gf2n(f,xp)
54: UP f;
55: UP *xp;
56: {
57: UP x,t,invf;
58: int k,n;
59: GF2N lm;
60: struct oEGT eg_sq,eg_rem,eg_mul,eg_inv,eg0,eg1,eg2;
61:
62: n = degup2(current_mod_gf2n->dense);
63: MKGF2N(ONEUP2,lm);
64: x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;
65:
66: reverseup(f,f->d,&t);
67: invmodup(t,f->d,&invf);
68: for ( k = 0; k < n; k++ ) {
69: squareup_gf2n(x,&t);
70: rembymulup_special(t,f,invf,&x);
71: /* remup(t,f,&x); */
72: }
73: *xp = x;
74: }
75:
76: /* g^d mod f */
77:
78: void generic_powermodup_gf2n(g,f,d,xp)
79: UP g,f;
80: Q d;
81: UP *xp;
82: {
83: N e;
84: UP x,y,t,invf,s;
85: int k;
86: GF2N lm;
87:
88: e = NM(d);
89: MKGF2N(ONEUP2,lm);
90: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
91: remup(g,f,&x);
92: if ( !x ) {
93: *xp = !d ? y : 0;
94: return;
95: } else if ( !x->d ) {
96: pwrup(x,d,xp);
97: return;
98: }
99: reverseup(f,f->d,&t);
100: invmodup(t,f->d,&invf);
101: for ( k = n_bits(e)-1; k >= 0; k-- ) {
102: squareup_gf2n(y,&t);
103: rembymulup_special(t,f,invf,&s);
104: y = s;
105: if ( e->b[k/32] & (1<<(k%32)) ) {
106: mulup(y,x,&t);
107: remup(t,f,&s);
108: y = s;
109: }
110: }
111: *xp = y;
112: }
113:
114: /* g+g^2+...+g^(2^(nd-1)) mod f; where e = deg(mod) */
115:
116: void tracemodup_gf2n(g,f,d,xp)
117: UP g,f;
118: Q d;
119: UP *xp;
120: {
121: UP x,t,s,u,invf;
122: int en,i;
123:
124: en = QTOS(d)*degup2(current_mod_gf2n->dense);
125: remup(g,f,&x);
126: if ( !x ) {
127: *xp = 0;
128: return;
129: }
130: reverseup(f,f->d,&t);
131: invmodup(t,f->d,&invf);
132: for ( i = 1, t = s = x; i < en; i++ ) {
133: squareup_gf2n(t,&u);
134: rembymulup_special(u,f,invf,&t);
135: addup(s,t,&u); s = u;
136: }
137: *xp = s;
138: }
139:
140: void tracemodup_gf2n_slow(g,f,d,xp)
141: UP g,f;
142: Q d;
143: UP *xp;
144: {
145: UP x,t,s,u;
146: int en,i;
147:
148: en = QTOS(d)*degup2(current_mod_gf2n->dense);
149: remup(g,f,&x);
150: if ( !x ) {
151: *xp = 0;
152: return;
153: }
154: for ( i = 1, t = s = x; i < en; i++ ) {
155: squareup_gf2n(t,&u);
156: remup(u,f,&t);
157: addup(s,t,&u); s = u;
158: }
159: *xp = s;
160: }
161:
162: static struct oEGT eg_trace_tab,eg_trace_mul;
163:
164: void tracemodup_gf2n_tab(g,f,d,xp)
165: UP g,f;
166: Q d;
167: UP *xp;
168: {
169: UP x0,x2,t,s,u;
170: int en,i;
171: UP *tab;
172: GF2N one;
173: struct oEGT eg1,eg2;
174:
175: en = QTOS(d)*degup2(current_mod_gf2n->dense);
176: remup(g,f,&t); g = t;
177: if ( !g ) {
178: *xp = 0;
179: return;
180: }
181:
182: MKGF2N(ONEUP2,one);
183: x0 = UPALLOC(0); x0->d = 0; x0->c[0] = (Num)one;
184: x2 = UPALLOC(2); x2->d = 2; x2->c[2] = (Num)one;
185:
186: tab = (UP *)ALLOCA(en*sizeof(UP));
187: tab[0] = x0;
188: remup(x2,f,&tab[1]);
189:
190: for ( i = 2; i < en; i++ ) {
191: mulup(tab[i-1],tab[1],&t); remup(t,f,&tab[i]);
192: }
193:
194: for ( i = 1, t = s = g; i < en; i++ ) {
195: square_rem_tab_up_gf2n(t,tab,&u); t = u;
196: addup(s,t,&u); s = u;
197: }
198: *xp = s;
199: }
200:
201: void square_rem_tab_up_gf2n(f,tab,rp)
202: UP f;
203: UP *tab;
204: UP *rp;
205: {
206: UP s,t,u,n;
207: Num *c;
208: int i,d;
209:
210: n = UPALLOC(0); n->d = 0;
211: if ( !f )
212: *rp = 0;
213: else {
214: d = f->d; c = f->c;
215: up_lazy = 1;
216: for ( i = 0, s = 0; i <= d; i++ ) {
217: squaregf2n((GF2N)c[i],(GF2N *)(&n->c[0]));
218: mulup(tab[i],n,&t); addup(s,t,&u); s = u;
219: }
220: up_lazy = 0;
221: simpup(s,rp);
222: }
223: }
224:
225: void powertabup_gf2n(f,xp,tab)
226: UP f;
227: UP xp;
228: UP *tab;
229: {
230: UP y,t,invf;
231: int i,d;
232: GF2N lm;
233:
234: d = f->d;
235: MKGF2N(ONEUP2,lm);
236: y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
237: tab[0] = y;
238: tab[1] = xp;
239:
240: reverseup(f,f->d,&t);
241: invmodup(t,f->d,&invf);
242:
243: for ( i = 2; i < d; i++ ) {
244: if ( debug_up )
245: fprintf(stderr,".");
246: if ( !(i%2) )
247: squareup_gf2n(tab[i/2],&t);
248: else
249: kmulup(tab[i-1],xp,&t);
250: rembymulup_special(t,f,invf,&tab[i]);
251: /* remup(t,f,&tab[i]); */
252: }
253: }
254:
255: void find_root_gf2n(f,r)
256: UP f;
257: GF2N *r;
258: {
259: UP g,ut,c,t,h,rem;
260: int n;
261: GF2N rn;
262: struct oEGT eg0,eg1,eg_trace;
263:
264: n = degup2(current_mod_gf2n->dense);
265: g = f;
266: while ( g->d > 1 ) {
267: ut = UPALLOC(1); ut->c[0] = 0;
268: randomgf2n(&rn);
269: if ( !rn )
270: continue;
271: ut->c[1] = (Num)rn; ut->d = 1;
272: tracemodup_gf2n_tab(ut,f,ONE,&c);
273: gcdup(c,g,&h);
274: if ( h->d && h->d < g->d ) {
275: if ( 2*h->d > g->d ) {
276: qrup(g,h,&t,&rem); g = t;
277: if ( rem )
278: error("find_root_gf2n : cannot happen");
279: } else
280: g = h;
281: }
282: monicup(g,&t); g = t;
283: printf("deg(g)=%d\n",g->d);
284: }
285: divgf2n((GF2N)g->c[0],(GF2N)g->c[1],r);
286: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>