![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to lmi.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to lmi.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/lmi.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #include "inline.h"
5:
6: typedef struct oGEN_LM {
7: int id;
8: N dense;
9: N sparse;
10: int bit;
11: unsigned int lower;
12: } *GEN_LM;
13:
14: void pwrlm0(N,N,N *);
15: void gen_simpn(N,N*);
16: void gen_simpn_force(N,N*);
17:
18: int lm_lazy;
19:
20: static GEN_LM current_mod_lm;
21:
22: void random_lm(r)
23: LM *r;
24: {
25: N n,t;
26:
27: if ( !current_mod_lm )
28: error("random_lm : current_mod_lm is not set");
29: randomn(current_mod_lm->bit+1,&n);
30: gen_simpn_force(n,&t);
31: MKLM(t,*r);
32: }
33:
34: void ntosparsen(p,bits)
35: N p;
36: N *bits;
37: {
38: int l,i,j,nz;
39: N r;
40: unsigned int *pb,*rb;
41: if ( !p )
42: *bits = 0;
43: else {
44: l = n_bits(p);
45: pb = BD(p);
46: for ( i = l-1, nz = 0; i >= 0; i-- )
47: if ( pb[i>>5]&(1<<i&31) ) nz++;
48: *bits = r = NALLOC(nz); PL(r) = nz;
49: rb = BD(r);
50: for ( i = l-1, j = 0; i >= 0; i-- )
51: if ( pb[i>>5]&(1<<i&31) )
52: rb[j++] = i;
53: }
54: }
55:
56: void setmod_lm(p)
57: N p;
58: {
59: int i;
60:
61: if ( !current_mod_lm )
62: current_mod_lm = (GEN_LM)MALLOC(sizeof(struct oGEN_LM));
63: bzero((char *)current_mod_lm,sizeof(struct oGEN_LM));
64: current_mod_lm->dense = p;
65: ntosparsen(p,¤t_mod_lm->sparse);
66: current_mod_lm->bit = n_bits(p);
67: current_mod_lm->id = UP2_DENSE; /* use ordinary rep. by default */
68: if ( !(current_mod_lm->bit%32) ) {
69: for ( i = PL(p)-1; i >= 1; i-- )
70: if ( BD(p)[i] != 0xffffffff )
71: break;
72: if ( !i ) {
73: current_mod_lm->lower = (unsigned int)((((UL)1)<<32)-((UL)BD(p)[0]));
74: current_mod_lm->id = UP2_SPARSE; /* XXX */
75: }
76: }
77: }
78:
79: void getmod_lm(p)
80: N *p;
81: {
82: if ( !current_mod_lm )
83: *p = 0;
84: else
85: *p = current_mod_lm->dense;
86: }
87:
88: void simplm(n,r)
89: LM n;
90: LM *r;
91: {
92: N rem;
93:
94: if ( !n )
95: *r = 0;
96: else if ( NID(n) != N_LM )
97: *r = n;
98: else {
99: gen_simpn(n->body,&rem);
100: MKLM(rem,*r);
101: }
102: }
103:
104: void qtolm(q,l)
105: Q q;
106: LM *l;
107: {
108: N rn;
109: LM a,b,c;
110:
111: if ( !q || (OID(q)==O_N && ((NID(q) == N_LM) || (NID(q) == N_GFPN))) ) { /* XXX */
112: *l = (LM)q;
113: } else if ( OID(q) == O_N && NID(q) == N_Q ) {
114: gen_simpn(NM(q),&rn); MKLM(rn,a);
115: if ( SGN(q) < 0 ) {
116: chsgnlm(a,&c); a = c;
117: }
118: if ( DN(q) ) {
119: gen_simpn_force(DN(q),&rn); MKLM(rn,b);
120: if ( !b )
121: error("qtolm : invalid denominator");
122: divlm(a,b,l);
123: } else
124: *l = a;
125: } else
126: error("qtolm : invalid argument");
127: }
128:
129: #define NZLM(a) ((a)&&(NID(a)==N_LM))
130:
131: void addlm(a,b,c)
132: LM a,b;
133: LM *c;
134: {
135: int s;
136: N t,t1;
137: LM z;
138:
139: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
140: if ( !a )
141: *c = b;
142: else if ( !b )
143: *c = a;
144: else {
145: addn(a->body,b->body,&t);
146: gen_simpn(t,&t1);
147: MKLM(t1,*c);
148: }
149: }
150:
151: void sublm(a,b,c)
152: LM a,b;
153: LM *c;
154: {
155: int s;
156: N t,t1,mod;
157: LM z;
158:
159: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
160: if ( !b )
161: *c = a;
162: else if ( !a )
163: chsgnlm(b,c);
164: else {
165: s = subn(a->body,b->body,&t);
166: if ( !s )
167: *c = 0;
168: else {
169: MKLM(t,z);
170: if ( s < 0 )
171: chsgnlm(z,c);
172: else
173: *c = z;
174: }
175: }
176: }
177:
178: void mullm(a,b,c)
179: LM a,b;
180: LM *c;
181: {
182: N t,q,r;
183: LM z;
184:
185: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
186: if ( !a || !b )
187: *c = 0;
188: else {
189: kmuln(a->body,b->body,&t);
190: gen_simpn(t,&r);
191: MKLM(r,*c);
192: }
193: }
194:
195: void divlm(a,b,c)
196: LM a,b;
197: LM *c;
198: {
199: LM r;
200: Q t,m,i;
201: LM z;
202:
203: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
204: if ( !b )
205: error("divlm: division by 0");
206: else if ( !a )
207: *c = 0;
208: else {
209: NTOQ(b->body,1,t); NTOQ(current_mod_lm->dense,1,m);
210: invl(t,m,&i);
211: MKLM(NM(i),r);
212: mullm(a,r,c);
213: }
214: }
215:
216: void chsgnlm(a,c)
217: LM a,*c;
218: {
219: LM t;
220: N s,u;
221:
222: if ( !a )
223: *c = a;
224: else {
225: qtolm((Q)a,&t); a = t;
226: gen_simpn_force(a->body,&s);
227: subn(current_mod_lm->dense,s,&u);
228: MKLM(u,*c);
229: }
230: }
231:
232: void pwrlm(a,b,c)
233: LM a;
234: Q b;
235: LM *c;
236: {
237: LM t;
238: N s;
239:
240: if ( !b ) {
241: MKLM(ONEN,*c);
242: } else if ( !a )
243: *c = 0;
244: else {
245: qtolm((Q)a,&t); a = t;
246: pwrlm0(a->body,NM(b),&s);
247: MKLM(s,*c);
248: }
249: }
250:
251: void pwrlm0(a,n,c)
252: N a,n;
253: N *c;
254: {
255: N n1,t,t1,t2,r1;
256: int r;
257:
258: if ( !n )
259: *c = ONEN;
260: else if ( UNIN(n) )
261: *c = a;
262: else {
263: r = divin(n,2,&n1);
264: pwrlm0(a,n1,&t);
265: kmuln(t,t,&t1); gen_simpn(t1,&r1);
266: if ( r ) {
267: kmuln(r1,a,&t2); gen_simpn(t2,&r1);
268: }
269: *c = r1;
270: }
271: }
272:
273: int cmplm(a,b)
274: LM a,b;
275: {
276: LM z;
277:
278: qtolm((Q)a,&z); a = z; qtolm((Q)b,&z); b = z;
279: if ( !a )
280: if ( !b )
281: return 0;
282: else
283: return -1;
284: else if ( !b )
285: return 1;
286: else
287: return cmpn(a->body,b->body);
288: }
289:
290: void remn_special(N,N,int,unsigned int ,N *);
291:
292: void gen_simpn(a,b)
293: N a;
294: N *b;
295: {
296: if ( !current_mod_lm )
297: error("gen_simpn: current_mod_lm is not set");
298: if ( lm_lazy )
299: *b = a;
300: else if ( current_mod_lm->id == UP2_SPARSE )
301: remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
302: else
303: remn(a,current_mod_lm->dense,b);
304: }
305:
306: void gen_simpn_force(a,b)
307: N a;
308: N *b;
309: {
310: if ( !current_mod_lm )
311: error("gen_simpn_force: current_mod_lm is not set");
312: else if ( current_mod_lm->id == UP2_SPARSE )
313: remn_special(a,current_mod_lm->dense,current_mod_lm->bit,current_mod_lm->lower,b);
314: else
315: remn(a,current_mod_lm->dense,b);
316: }
317: