![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ddM.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / asm|
Return to ddM.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / asm |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/asm/ddM.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
2: #include "ca.h"
3: #include "base.h"
4: #include "inline.h"
5:
6: void ksquareummain(int,UM,UM);
7: void kmulummain(int,UM,UM,UM);
8: void c_copyum(UM,int,int *);
9: void copyum(UM,UM);
10: void extractum(UM,int,int,UM);
11: void ksquareum(int,UM,UM);
12: void kmulum(int,UM,UM,UM);
13:
14: /*
15: * mod is declared as 'int', because several xxxum functions contains signed
16: * integer addition/subtraction. So mod should be less than 2^31.
17: */
18:
19: void mulum(mod,p1,p2,pr)
20: int mod;
21: UM p1,p2,pr;
22: {
23: int *pc1,*pcr;
24: int *c1,*c2,*cr;
25: unsigned int mul;
26: int i,j,d1,d2;
27:
28: if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
29: DEG(pr) = -1;
30: return;
31: }
32: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
33: bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
34: for ( i = 0; i <= d2; i++, cr++ )
35: if ( mul = *c2++ )
36: for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ ) {
37: DMAR(*pc1,mul,*pcr,mod,*pcr)
38: }
39: DEG(pr) = d1 + d2;
40: }
41:
42: void mulsum(mod,p,n,pr)
43: int mod,n;
44: UM p,pr;
45: {
46: int *sp,*dp;
47: int i;
48:
49: for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
50: i >= 0; i--, dp--, sp-- ) {
51: DMAR(*sp,n,0,mod,*dp)
52: }
53: }
54:
55: int divum(mod,p1,p2,pq)
56: int mod;
57: UM p1,p2,pq;
58: {
59: int *pc1,*pct;
60: int *c1,*c2,*ct;
61: unsigned int inv,hd,tmp;
62: int i,j, d1,d2,dd;
63:
64: if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
65: DEG(pq) = -1;
66: return d1;
67: }
68: c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
69: if ( ( hd = c2[d2] ) != 1 ) {
70: inv = invm(hd,mod);
71: for ( pc1 = c2 + d2; pc1 >= c2; pc1-- ) {
72: DMAR(*pc1,inv,0,mod,*pc1)
73: }
74: } else
75: inv = 1;
76: for ( i = dd, ct = c1+d1; i >= 0; i-- )
77: if ( tmp = *ct-- ) {
78: tmp = mod - tmp;
79: for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- ) {
80: DMAR(*pc1,tmp,*pct,mod,*pct)
81: }
82: }
83: if ( inv != 1 ) {
84: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ ) {
85: DMAR(*pc1,inv,0,mod,*pc1)
86: }
87: for ( pc1 = c2, pct = c2+d2, inv = hd; pc1 <= pct; pc1++ ) {
88: DMAR(*pc1,inv,0,mod,*pc1)
89: }
90: }
91: for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
92: for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
93: *pct-- = *pc1--;
94: return i;
95: }
96:
97: void diffum(mod,f,fd)
98: int mod;
99: UM f,fd;
100: {
101: int *dp,*sp;
102: int i;
103: UL ltmp;
104:
105: for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
106: i >= 1; i--, dp--, sp-- ) {
107: DMAR(*sp,i,0,mod,*dp)
108: }
109: degum(fd,DEG(f) - 1);
110: }
111:
112: unsigned int pwrm(mod,a,n)
113: int mod,a;
114: int n;
115: {
116: unsigned int s,t;
117:
118: if ( !n )
119: return 1;
120: else if ( n == 1 )
121: return a;
122: else {
123: t = pwrm(mod,a,n/2);
124: DMAR(t,t,0,mod,s)
125: if ( n % 2 ) {
126: DMAR(s,a,0,mod,t)
127: return t;
128: } else
129: return s;
130: }
131: }
132:
133: unsigned int invm(s,mod)
134: unsigned int s;
135: int mod;
136: {
137: unsigned int r,a2,q;
138: unsigned int f1,f2,a1;
139:
140: for ( f1 = s, f2 = mod, a1 = 1, a2 = 0; ; ) {
141: q = f1/f2; r = f1 - f2*q; f1 = f2;
142: if ( !(f2 = r) )
143: break;
144: DMAR(a2,q,0,mod,r)
145: /* r = ( a1 - r + mod ) % mod; */
146: if ( a1 >= r )
147: r = a1 - r;
148: else {
149: r = (mod - r) + a1;
150: }
151: a1 = a2; a2 = r;
152: }
153: /* return( ( a2 >= 0 ? a2 : a2 + mod ) ); */
154: return a2;
155: }
156:
157: unsigned int rem(n,m)
158: N n;
159: unsigned int m;
160: {
161: unsigned int *x;
162: unsigned int t,r;
163: int i;
164:
165: if ( !n )
166: return 0;
167: for ( i = PL(n)-1, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
168: #if defined(sparc)
169: r = dsar(m,r,*x);
170: #else
171: DSAB(m,r,*x,t,r)
172: #endif
173: }
174: return r;
175: }
176:
177: #ifndef sparc
178: void addpadic(mod,n,n1,n2)
179: int mod;
180: int n;
181: unsigned int *n1,*n2;
182: {
183: unsigned int carry,tmp;
184: int i;
185:
186: for ( i = 0, carry = 0; i < n; i++ ) {
187: tmp = *n1++ + *n2 + carry;
188: DQR(tmp,mod,carry,*n2++)
189: /*
190: carry = tmp / mod;
191: *n2++ = tmp - ( carry * mod );
192: */
193: }
194: }
195: #endif
196:
197: void mulpadic(mod,n,n1,n2,nr)
198: int mod;
199: int n;
200: unsigned int *n1;
201: unsigned int *n2,*nr;
202: {
203: unsigned int *pn1,*pnr;
204: unsigned int carry,mul;
205: int i,j;
206:
207: bzero((char *)nr,(int)(n*sizeof(int)));
208: for ( j = 0; j < n; j++, n2++, nr++ )
209: if ( mul = *n2 )
210: for ( i = j, carry = 0, pn1 = n1, pnr = nr;
211: i < n; i++, pn1++, pnr++ ) {
212: carry += *pnr;
213: DMAB(mod,*pn1,mul,carry,carry,*pnr)
214: }
215: }
216:
217: extern up_kara_mag;
218:
219: void kmulum(mod,n1,n2,nr)
220: UM n1,n2,nr;
221: {
222: UM n,t,s,m,carry;
223: int d,d1,d2,len,i,l;
224: unsigned int *r,*r0;
225:
226: if ( !n1 || !n2 ) {
227: nr->d = -1; return;
228: }
229: d1 = DEG(n1)+1; d2 = DEG(n2)+1;
230: if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {
231: mulum(mod,n1,n2,nr); return;
232: }
233: if ( d1 < d2 ) {
234: n = n1; n1 = n2; n2 = n;
235: d = d1; d1 = d2; d2 = d;
236: }
237: if ( d2 > (d1+1)/2 ) {
238: kmulummain(mod,n1,n2,nr); return;
239: }
240: d = (d1/d2)+((d1%d2)!=0);
241: len = (d+1)*d2;
242: r0 = (unsigned int *)ALLOCA(len*sizeof(int));
243: bzero((char *)r0,len*sizeof(int));
244: m = W_UMALLOC(d2+1);
245: carry = W_UMALLOC(d2+1);
246: t = W_UMALLOC(d1+d2+1);
247: s = W_UMALLOC(d1+d2+1);
248: for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
249: extractum(n1,i*d2,d2,m);
250: if ( m ) {
251: kmulum(mod,m,n2,t);
252: addum(mod,t,carry,s);
253: c_copyum(s,d2,r);
254: extractum(s,d2,d2,carry);
255: } else {
256: c_copyum(carry,d2,r);
257: carry = 0;
258: }
259: }
260: c_copyum(carry,d2,r);
261: for ( l = len - 1; !r0[l]; l-- );
262: l++;
263: DEG(nr) = l-1;
264: bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
265: }
266:
267: void ksquareum(mod,n1,nr)
268: int mod;
269: UM n1,nr;
270: {
271: int d1;
272:
273: if ( !n1 ) {
274: nr->d = -1; return;
275: }
276: d1 = DEG(n1)+1;
277: if ( (d1 < up_kara_mag) ) {
278: pwrum(mod,n1,2,nr); return;
279: }
280: ksquareummain(mod,n1,nr);
281: }
282:
283: void extractum(n,index,len,nr)
284: UM n;
285: int index,len;
286: UM nr;
287: {
288: int *m;
289: int l;
290:
291: if ( !n ) {
292: nr->d = -1; return;
293: }
294: m = COEF(n)+index;
295: if ( (l = (DEG(n)+1)-index) >= len ) {
296: for ( l = len - 1; (l >= 0) && !m[l]; l-- );
297: l++;
298: }
299: if ( l <= 0 ) {
300: nr->d = -1; return;
301: } else {
302: DEG(nr) = l-1;
303: bcopy((char *)m,(char *)COEF(nr),l*sizeof(Q));
304: }
305: }
306:
307: void copyum(n1,n2)
308: UM n1,n2;
309: {
310: n2->d = n1->d;
311: bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(int));
312: }
313:
314: void c_copyum(n,len,p)
315: UM n;
316: int len;
317: int *p;
318: {
319: if ( n )
320: bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(int));
321: }
322:
323: void kmulummain(mod,n1,n2,nr)
324: int mod;
325: UM n1,n2,nr;
326: {
327: int d1,d2,h,len;
328: UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
329:
330: d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
331: n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
332: n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);
333: lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);
334: mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);
335: mid = W_UMALLOC(d1+d2+1);
336: s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);
337: extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
338: extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
339: kmulum(mod,n1hi,n2hi,hi); kmulum(mod,n1lo,n2lo,lo);
340: len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
341: bzero((char *)COEF(t1),len*sizeof(int));
342: if ( lo )
343: bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
344: if ( hi )
345: bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
346:
347: addum(mod,hi,lo,mid1);
348: subum(mod,n1hi,n1lo,s1); subum(mod,n2lo,n2hi,s2); kmulum(mod,s1,s2,mid2);
349: addum(mod,mid1,mid2,mid);
350: if ( mid ) {
351: len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
352: bzero((char *)COEF(t2),len*sizeof(int));
353: bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
354: addum(mod,t1,t2,nr);
355: } else
356: copyum(t1,nr);
357: }
358:
359: void ksquareummain(mod,n1,nr)
360: int mod;
361: UM n1,nr;
362: {
363: int d1,h,len;
364: UM n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
365:
366: d1 = DEG(n1)+1; h = (d1+1)/2;
367: n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
368: lo = W_UMALLOC(2*d1+1); hi = W_UMALLOC(2*d1+1);
369: mid1 = W_UMALLOC(2*d1+1); mid2 = W_UMALLOC(2*d1+1);
370: mid = W_UMALLOC(2*d1+1);
371: s1 = W_UMALLOC(2*d1+1);
372: extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
373: ksquareum(mod,n1hi,hi); ksquareum(mod,n1lo,lo);
374: len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
375: bzero((char *)COEF(t1),len*sizeof(int));
376: if ( lo )
377: bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
378: if ( hi )
379: bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
380:
381: addum(mod,hi,lo,mid1);
382: subum(mod,n1hi,n1lo,s1); ksquareum(mod,s1,mid2);
383: subum(mod,mid1,mid2,mid);
384: if ( mid ) {
385: len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
386: bzero((char *)COEF(t2),len*sizeof(int));
387: bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
388: addum(mod,t1,t2,nr);
389: } else
390: copyum(t1,nr);
391: }