Annotation of OpenXM_contrib2/asir2000/engine/igcdhack.c, Revision 1.3
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 ! noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.3 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/igcdhack.c,v 1.2 2000/08/21 08:31:27 noro Exp $
1.2 noro 49: */
1.1 noro 50: #if 0
51: #include "base.h"
52:
53: #define SWAP(a,b,Type) { Type temp=a; a = b; b = temp; }
54: #define SIGNED_VAL(a,s) ((s)>0?(a): -(a))
55: #endif
56:
57: #undef FULLSET
58:
59:
60: static int negate_nbd( b, l )
61: int *b, l;
62: /* It is assumed that l > 0 and b[0] != 0 */
63: {
64: int i, nz;
65:
66: *b = BASE - *b;
67: for ( nz = i = 1 ; i < l; ) {
68: b++, i++;
69: if ( (*b ^= BMASK) != 0 ) nz = i;
70: }
71: return nz;
72: }
73:
74: /****************/
75: /****************/
76:
77: static int abs_U_V_maxrshift( u, lu, v, lv, a )
78: int *u, lu, *v, lv, *a;
79: /*
80: * Compute | U - V | >> (as much as possible) in a[], where U = u[0:lu-1]
81: * and V = v[0:lv-1]. Value returned is the length of the result
82: * multiplied by the sign of (U-V).
83: */
84: {
85: int sign, rsh, t, b, i, i_1, l, *p = a;
86:
87: /* first, determine U <, ==, > V */
88: sign = 1;
89: if ( lu == lv ) {
90: int i;
91:
92: for ( i = lu -1; i >= 0; i-- ) {
93: if ( u[i] == v[i] ) continue;
94: lu = lv = i + 1;
95: if ( u[i] > v[i] ) goto U_V;
96: else goto swap;
97: }
98: return 0;
99: } else if ( lu > lv ) goto U_V;
100: swap:
101: SWAP( lu, lv, int ); SWAP( u, v, int * );
102: sign = -1;
103: U_V: /* U > V, i.e., (lu > lv) || (lu == lv && u[lu-1] > v[lv-1]) */
104: for ( i = 0; i < lv; i++ ) if ( u[i] != v[i] ) break;
105: if ( i != 0 ) u += i, lu -= i, v += i, lv -= i;
106: if ( lv == 0 ) { /* no more pv[], and lu > 0 */
107: i = trailingzeros_nbd( u, &rsh );
108: i = rshift_nbd( u, lu, rsh, i, a );
109: return SIGNED_VAL( i, sign );
110: }
111: /**/
112: if ( (t = u[0] - v[0]) >= 0 ) b = 0;
113: else b = 1, t += BASE;
114: TRAILINGZEROS( t, rsh );
115: /**/
116: if ( rsh != 0 ) {
117: int w, lsh = BSH - rsh;
118:
119: /* subtract v[*] from u[*] */
120: for ( i = 1, l = 0; i < lv; i++, t = w >> rsh ) {
121: if ( (w = u[i] - v[i] - b) >= 0 ) b = 0;
122: else b = 1, w += BASE;
123: if ( (*p++ = t | ((w << lsh) & BMASK)) != 0 ) l = i;
124: }
125: /* subtraction of V completed. Note there may exist non-0 borrow b */
126: if ( i >= lu ) /* lu == lv, and b==0 is implied */ goto chk_nz;
127: /* we still have (lu-i) elms in U */
128: if ( b != 0 ) { /* non-zero borrow */
129: if ( (w = u[i++] -1) != 0 ) t |= ((w << lsh) & BMASK);
130: if ( i >= lu ) { /* lu == lv+1. The M.S. word w may beome 0 */
131: if ( w != 0 ) {
132: *p++ = t;
133: l = lv;
134: t = w >> rsh;
135: } else lu--;
136: goto chk_nz;
137: }
138: *p++ = t;
139: if ( w < 0 ) {
140: while ( (w = u[i++]) == 0 ) *p++ = BMASK;
141: w -= 1;
142: *p++ = (BMASK >> rsh) | ((w << lsh) & BMASK);
143: }
144: t = w >> rsh;
145: }
146: for ( ; i < lu; i++, t = w >> rsh ) *p++ = t | (((w = u[i]) << lsh) & BMASK);
147: l = i -1;
148: chk_nz:
149: if ( t != 0 ) *p = t, l = /*1 + (p - a)*/ lu;
150: goto ret_nz;
151: }
152: /* rsh == 0 : need not be right-shifted */
153: *p++ = t, l = 1;
154: for ( i = 1; i < lv; ) {
155: if ( (t = u[i] - v[i] - b) >= 0 ) b = 0;
156: else b = 1, t += BASE;
157: i++;
158: if ( (*p++ = t) != 0 ) l = i;
159: }
160: if ( i >= lu ) /* lu == lv, and b==0 is implied */ goto ret_nz;
161: if ( b != 0 ) { /* non-zero borrow */
162: t = u[i++] -1;
163: if ( i >= lu ) /* lu == lv+1. The M.S. word w may beome 0 */ goto chk_nz;
164: if ( t < 0 ) {
165: do *p++ = BMASK; while ( (t = u[i++]) == 0 );
166: t -= 1;
167: if ( i >= lu ) {
168: l = lu;
169: if ( t == 0 ) l--;
170: else *p = t;
171: goto ret_nz;
172: }
173: }
174: *p++ = t;
175: }
176: for ( ; i < lu; i++ ) *p++ = u[i];
177: l = lu;
178: /**/
179: ret_nz:
180: return SIGNED_VAL( l, sign );
181: }
182:
183: /****************/
184: /****************/
185:
186: #ifdef CALL
187: static int aV_plus_c(), aV_plus_c_sh();
188: #endif
189: static int abs_U_aV_c(), abs_U_aV_c_sh();
190: static int abs_aVplusc_U(), abs_aVplusc_U_sh();
191:
192: static int abs_U_aV_maxrshift( u, lu, a, v, lv, ans )
193: int *u, lu, a, *v, lv, *ans;
194: {
195: int i, t, c, rsh, l, *p = ans, w, lsh;
196: int hi;
197:
198: if ( a == 1 ) return( (l = abs_U_V_maxrshift( u, lu, v, lv, ans )) >= 0 ? l : -l );
199: /**/
200: for ( l = lu <= lv ? lu : lv, i = c = 0; i < l; ) {
201: DMA( a, v[i], c, hi, t, not_used );
202: c = hi;
203: if ( (t -= u[i++]) != 0 ) goto non0w;
204: }
205: /**/
206: if ( i < lv ) { /* no more u[] elm.'s, and we still have non-zero v[] elm's */
207: do {
208: DMA( a, v[i], c, hi, t, not_used );
209: i++, c = hi;
210: } while ( t == 0 );
211: TRAILINGZEROS( t, rsh );
212: v += i, lv -= i;
213: #ifdef CALL
214: if ( rsh == 0 ) {
215: *p++ = t;
216: return( aV_plus_c( a, v, lv, c, p ) + 1 );
217: } else return aV_plus_c_sh( a, v, lv, c, rsh, t, ans );
218: #else
219: i = 0;
220: if ( rsh == 0 ) {
221: *p++ = t;
222: for ( ; i < lv; i++, c = hi, *p++ = t ) { DMA( a, v[i], c, hi, t, not_used ); }
223: if ( c == 0 ) return( i + 1 );
224: *p = c;
225: return( i + 2 );
226: } else {
227: for ( lsh = BSH - rsh; i < lv; i++, c = hi, t = w >> rsh ) {
228: DMA( a, v[i], c, hi, w, not_used );
229: *p++ = t | ((w << lsh) & BMASK);
230: }
231: if ( c != 0 ) {
232: *p++ = t | ((c << lsh) & BMASK);
233: i++;
234: t = c >> rsh;
235: }
236: if ( t == 0 ) return i;
237: *p = t;
238: return( i + 1 );
239: }
240: #endif
241: }
242: /* no more v[] elm.'s */
243: if ( i >= lu ) {
244: if ( c == 0 ) return 0;
245: c_sh:
246: while ( (c&1) == 0 ) c >>= 1;
247: *p = c;
248: return 1;
249: }
250: t = u[i++] - c;
251: if ( i >= lu ) {
252: if ( t == 0 ) return 0;
253: c = t < 0 ? -t : t;
254: goto c_sh;
255: } else if ( t == 0 ) {
256: while ( (t = u[i++]) == 0 ) ;
257: c = 0;
258: } else if ( t < 0 ) t += BASE, c = 1;
259: else c = 0;
260: /* diff turned out to be > 0 */
261: u += i, lu -= i;
262: TRAILINGZEROS( t, rsh );
263: i = 0;
264: if ( rsh == 0 ) {
265: *p++ = t;
266: if ( c != 0 ) {
267: for ( ; i < lu; *p++ = BMASK ) if ( (t = u[i++]) != 0 ) break;
268: t--;
269: if ( i >= lu && t == 0 ) return i;
270: *p++ = t;
271: }
272: for ( ; i < lu; i++ ) *p++ = u[i];
273: } else { /* rsh != 0 */
274: lsh = BSH - rsh;
275: if ( c != 0 ) { /* there still exists u[] elms.'s because diff is known to be > 0 */
276: if ( (w = u[i++]) == 0 ) {
277: *p++ = t | ((BMASK << lsh) & BMASK);
278: for ( ; (w = u[i++]) == 0; *p++ = BMASK ) ;
279: t = BMASK >> rsh;
280: }
281: w--;
282: *p++ = t | ((w << lsh) & BMASK);
283: t = w >> rsh;
284: }
285: for ( ; i < lu; i++, t = w >> rsh ) *p++ = t | (((w = u[i]) << lsh) & BMASK);
286: if ( t == 0 ) return( i );
287: *p = t;
288: }
289: return( i + 1 );
290: /**/
291: non0w:
292: u += i, lu -= i, v += i, lv -= i;
293: /* | t + BASE*(c + a*v[0:lv-1] - u[0:lu-1] | */
294: if ( lv < lu ) { /* a*V < U if lv<lu-1, and very likely if lv==lu-1 */
295: t = -t; /* We shall compute U - a*V */
296: if ( t < 0 ) t += BASE, c++;
297: } else /* a*V > U if lv>lu, and very likely even if lv==lu */
298: if ( t < 0 ) t += BASE, c--;
299: TRAILINGZEROS( t, rsh );
300: return( rsh == 0 ? (lv < lu ? abs_U_aV_c( t, u, lu, a, v, lv, c, ans ) :
301: abs_aVplusc_U( t, a, v, lv, c, u, lu, ans )) :
302: lv < lu ? abs_U_aV_c_sh( t, u, lu, a, v, lv, c, rsh, ans ) :
303: abs_aVplusc_U_sh( t, a, v, lv, c, u, lu, rsh, ans ) );
304: }
305:
306: /*********/
307:
308: static int aV_plus_c( a, v, lv, c, p )
309: int a, *v, lv, c, *p;
310: /*
311: * Compute (a*V + c) in p[], where V = v[0:lv-1].
312: */
313: {
314: int i, t;
315: int hi;
316:
317: for ( i = 0; i < lv; i++, *p++ = t, c = hi ) { DMA( a, v[i], c, hi, t, not_used); }
318: if ( c == 0 ) return i;
319: *p = c;
320: return( i + 1 );
321: }
322:
323: static int aV_plus_c_sh( a, v, lv, c, rsh, t, p )
324: int a, *v, lv, c, rsh, *p;
325: /*
326: * Compute (t + BASE*((a*V + c) >> rsh)) in p[], where V = v[0:lv-1].
327: */
328: {
329: int i, w, lsh = BSH - rsh;
330: int hi;
331:
332: for ( i = 0; i < lv; i++, c = hi, t = w >> rsh ) {
333: DMA( a, v[i], c, hi, w, not_used );
334: *p++ = t | ((w << lsh) & BMASK);
335: }
336: if ( c != 0 ) {
337: *p++ = t | ((c << lsh) & BMASK);
338: i++;
339: t = c >> rsh;
340: }
341: if ( t == 0 ) return i;
342: *p = t;
343: return( i + 1 );
344: }
345:
346: /*********/
347:
348: static int abs_aVplusc_U( t, a, v, lv, c, u, lu, ans )
349: int t, a, *v, lv, c, *u, lu, *ans;
350: /* compute | t + BASE*(a*V+c - U) | in ans[], where V=v[0:lv-1], U=u[0:lu-1],
351: * and a > 1, -1 <= c < BASE, lv >= lu >=0 && t != 0 are assumed.
352: */
353: {
354: int i, l, b, *p = ans;
355: int hi;
356:
357: *p++ = t, l = 1;
358: if ( c < 0 ) c = 0, b = 1;
359: else b = 0;
360: for ( i = 0; i < lu; *p++ = t ) {
361: DMA( a, v[i], c, hi, t, not_used );
362: t -= u[i++], c = hi;
363: if ( b != 0 ) t--;
364: b = 0;
365: if ( t > 0 ) l = i + 1;
366: else if ( t < 0 ) {
367: t += BASE, l = i + 1;
368: if ( c != 0 ) c--;
369: else b = 1;
370: }
371: }
372: /* at this point, i == lu and we have written (i+1) entries in ans[] */
373: if ( i >= lv ) { /* no more elms in v[] */
374: if ( b != 0 ) {
375: if ( i == 0 ) { ans[i] = BASE - t; return 1; }
376: l = negate_nbd( ans, i ); /* <================ */
377: /* if ( (t = BMASK ^ ans[i]) == 0 ) return l;*/
378: if ( (t ^= BMASK) == 0 ) return l;
379: ans[i] = t;
380: return( i + 1 );
381: }
382: if ( c == 0 ) return l;
383: *p = c;
384: return( i + 2 );
385: }
386: /* diff turned out to be > 0 */
387: if ( b != 0 ) { /* c == 0 */
388: while ( (b = v[i++]) == 0 ) *p++ = BMASK;
389: DM( a, b, hi, t );
390: if ( t != 0 ) t--, c = hi;
391: else t = BMASK, c = hi - 1;
392: *p++ - t;
393: }
394: #ifdef CALL
395: return( aV_plus_c( a, &v[i], lv - i, c, p ) + i +1 );
396: #else
397: for ( ; i < lv; i++, *p++ = t, c = hi ) { DMA( a, v[i], c, hi, t, not_used );}
398: if ( c == 0 ) return( i + 1 );
399: *p = c;
400: return( i + 2 );
401: #endif
402: }
403:
404: static int abs_aVplusc_U_sh( t, a, v, lv, c, u, lu, rsh, ans )
405: int t, a, *v, lv, c, *u, lu, rsh, *ans;
406: {
407: int i, l, w, b, lsh = BSH - rsh, *p = ans;
408: int hi;
409:
410: if ( c < 0 ) c = 0, b = 1;
411: else b = 0;
412: for ( i = 0; i < lu; t = w >> rsh ) {
413: DMA( a, v[i], c, hi, w, not_used );
414: w -= u[i++], c = hi;
415: if ( b != 0 ) w--;
416: b = 0;
417: if ( w < 0 ) {
418: w += BASE;
419: if ( c != 0 ) c--;
420: else b = 1;
421: }
422: if ( (*p++ = t | ((w << lsh) & BMASK)) != 0 ) l = i;
423: }
424: /* at this point, i == lu and we have written i entries in ans[] */
425: if ( i >= lv ) {
426: if ( b != 0 ) { /* diff turned out to be < 0 */
427: if ( i == 0 ) {
428: *p = (BASE >> rsh) - t;
429: return 1;
430: }
431: l = negate_nbd( ans, i ); /* <================ */
432: if ( (t ^= (BMASK >> rsh)) == 0 ) return l;
433: *p = t;
434: return( i + 1 );
435: }
436: if ( c == 0 ) {
437: if ( t == 0 ) return l;
438: *p = t;
439: return( i + 1 );
440: }
441: *p++ = t | ((c << lsh) & BMASK);
442: if ( (t = c >> rsh) == 0 ) return( i + 1 );
443: *p = t;
444: return( i + 2 );
445: }
446: /* lv >= lu+1 = i+1, and diff is > 0 */
447: if ( b != 0 ) { /* c == 0 */
448: if ( (w = v[i++]) == 0 ) {
449: *p++ = t | ((BMASK << lsh) & BMASK);
450: while ( (w = v[i++]) == 0 ) *p++ = BMASK;
451: t = BMASK >> rsh;
452: } else {
453: DM( a, w, c, b );
454: if ( b != 0 ) b--;
455: else b = BMASK, c--;
456: *p++ = t | ((b << lsh) & BMASK);
457: t = b >> rsh;
458: }
459: }
460: /**/
461: #ifdef CALL
462: return( aV_plus_c_sh( a, &v[i], lv - i, c, rsh, t, p ) + i );
463: #else
464: for ( ; i < lv; i++, t = w >> rsh, c = hi ) {
465: DMA( a, v[i], c, hi, w, not_used );
466: *p++ = t | ((w << lsh) & BMASK);
467: }
468: if ( c != 0 ) {
469: *p++ = t | ((c << lsh) & BMASK);
470: i++, t = c >> rsh;
471: }
472: if ( t == 0 ) return i;
473: *p = t;
474: return( i + 1 );
475: #endif
476: }
477:
478: /*********/
479:
480: static int abs_U_aV_c( t, u, lu, a, v, lv, c, ans )
481: int t, *u, lu, a, *v, lv, c, *ans;
482: /* compute | t + BASE*(U - a*V - c) | in ans[], where U = u[0:lu-1], V = v[0:lv-1],
483: * c is <= BASE-1, and lu >= lv+1 && t != 0 assumed.
484: */
485: {
486: int i, w, l, *p = ans;
487: int hi;
488:
489: *p++ = t, l = 1;
490: for ( i = 0; i < lv; *p++ = t ) {
491: DMA( a, v[i], c, hi, w, not_used );
492: t = u[i] - w;
493: i++, c = hi;
494: if ( t != 0 ) l = i + 1;
495: if ( t < 0 ) t += BASE, c++;
496: }
497: /* at this point, i == lv */
498: if ( lu == lv+1 ) {
499: if ( (t = u[i]) == c ) return l;
500: else if ( t > c ) { *p = t - c; return( lu + 1 ); }
501: l = negate_nbd( ans, lu ); /* <================ */
502: if ( (c -= (t + 1)) == 0 ) return l;
503: *p = c;
504: return( lu + 1 );
505: }
506: /* lu >= lv+2 and the diff turned out to be >= 0 */
507: if ( c != 0 ) {
508: /* note c <= BASE. If c == BASE, we must care about the length of the result
509: of the case when lu == lv+2 && u[lu-2:lu-1]=u[i:i+1] == BASE */
510: if ( c >= BASE && i+2 == lu && u[i+1] == 1 ) {
511: /* m.s. word of the diff becomes 0 */
512: if ( (t = u[i] - (c & BMASK)) == 0 ) return l;
513: *p = t;
514: return lu;
515: }
516: for ( ; 1; c = 1, *p++ = t + BASE ) if ( (t = u[i++] - c) >= 0 ) break;
517: if ( t == 0 && i >= lu ) return lu;
518: *p++ = t;
519: }
520: for ( ; i < lu; i++ ) *p++ = u[i];
521: return( lu + 1 );
522: }
523:
524: static int abs_U_aV_c_sh( t, u, lu, a, v, lv, c, rsh, ans )
525: int t, *u, lu, a, *v, lv, c, rsh, *ans;
526: /* r-shift version of the above */
527: {
528: int i, w, l, lsh = BSH - rsh, *p = ans;
529: int hi;
530:
531: for ( l = i = 0; i < lv; t = w >> rsh ) {
532: DMA( a, v[i], c, hi, w, not_used );
533: w = u[i] - w, c = hi;
534: i++;
535: if ( w < 0 ) w += BASE, c++;
536: if ( (*p++ = t | ((w << lsh) & BMASK)) != 0 ) l = i;
537: }
538: /* at this point, i == lv, and we have written lv elm.s in ans[] */
539: if ( lu == lv+1 ) {
540: if ( (w = u[i] - c) == 0 ) {
541: if ( t == 0 ) return l;
542: *p = t;
543: return lu; /* == lv+1 */
544: } else if ( w > 0 ) {
545: l0: *p++ = t | ((w << lsh) & BMASK);
546: if ( (w >>= rsh) == 0 ) return lu; /* == lv+1 */
547: *p = w;
548: return( lu + 1 );
549: }
550: /* diff turned out to be < 0 */
551: w = - w - 1;
552: if ( lv == 0 ) { /* lu == 1 */
553: t = (BASE >> rsh) - t;
554: if ( w != 0 ) goto l0;
555: *p = t;
556: return 1;
557: }
558: l = negate_nbd( ans, lv ); /* <================ */
559: t ^= (BMASK >> rsh);
560: if ( w != 0 ) goto l0;
561: if ( t == 0 ) return l;
562: *p = t;
563: return lu; /* == lv + 1 */
564: }
565: /* now, lu >= lv+2 == i+2 */
566: if ( c != 0 ) {
567: /* note c <= BASE. If c == BASE, we must care about the length of the result
568: of the case when lu == lv+2 && u[lu-2:lu-1]=u[i:i+1] == BASE */
569: if ( c >= BASE && i+2 == lu && t == 0 && u[i] == 0 && u[i+1] == 1 )
570: return l; /* m.s. word of the diff has become 0 */
571: for ( ; 1; t = w >> rsh, c = 1 ) {
572: if ( (w = u[i++] - c) >= 0 ) break;
573: w += BASE;
574: *p++ = t | ((w << lsh) & BMASK);
575: }
576: t |= ((w << lsh) & BMASK);
577: w >>= rsh;
578: if ( i >= lu ) {
579: if ( w != 0 ) {
580: *p++ = t;
581: *p = w;
582: return( lu + 1 );
583: } else if ( t == 0 ) return( i - 1 );
584: else { *p = t; return i; }
585: }
586: *p++ = t;
587: t = w;
588: }
589: for ( ; i < lu; i++, t = w >> rsh ) *p++ = t | (((w = u[i]) << lsh) & BMASK);
590: if ( t == 0 ) return i;
591: *p = t;
592: return( i + 1 );
593: }
594:
595: /****************/
596: /****************/
597: #ifdef FULLSET
598: static int abs_axU_V_maxrshift(), abs_axU_bxV_maxrshift();
599:
600: static int abs_aU_bV_maxrshift( a, u, lu, b, v, lv, ans )
601: int a, *u, lu, b, *v, lv, *ans;
602: /*
603: * Compute | a*U - b*V | >> (as much as possible) in ans[], where U=u[0:lu-1]
604: * and V=v[0:lv-1].
605: */
606: {
607: if ( a == 1 ) return abs_U_aV_maxrshift( u, lu, b, v, lv, ans );
608: else if ( b == 1 ) return abs_U_aV_maxrshift( v, lv, a, u, lu, ans );
609: else if ( a == b ) return abs_axU_V_maxrshift( a, u, lu, v, lv, ans );
610: return abs_axU_bxV_maxrshift( a, u, lu, b, v, lv, ans );
611: }
612:
613: static int abs_axU_V_maxrshift( a, u, lu, v, lv, ans )
614: int a, *u, lu, *v, lv, *ans;
615: /*
616: * Compute a*| U - V | >> (as much as possible) in ans[], where U=u[0:lu-1]
617: * and V=v[0:lv-1], 0 < a < 2^BASE.
618: */
619: {
620: int i, l, b, c, t, rsh, *p = ans, lsh, w, x;
621: int hi;
622:
623: if ( lu > lv ) goto U_V;
624: else if ( lu < lv ) goto swap;
625: else {
626: for ( i = lu; i-- > 0; ) {
627: if ( u[i] == v[i] ) continue;
628: lu = lv = i + 1;
629: if ( u[i] > v[i] ) goto U_V;
630: else goto swap;
631: }
632: return 0;
633: }
634: swap:
635: SWAP( lu, lv, int ); SWAP( u, v, int * );
636: U_V:
637: for ( b = c = i = 0; i < lv; ) {
638: if ( (w = u[i] - v[i] - b) < 0 ) w += BASE, b = 1;
639: else b = 0;
640: DMA( a, w, c, hi, t, not_used );
641: i++, c = hi;
642: if ( t != 0 ) goto non0w;
643: }
644: while ( i < lu ) {
645: if ( (w = u[i] - b) < 0 ) w += BASE, b = 1;
646: else b = 0;
647: DMA( a, w, c, hi, t, not_used );
648: i++, c = hi;
649: if ( t != 0 ) goto non0w;
650: }
651: if ( b != 0 ) c -= a;
652: else if ( c == 0 ) return 0;
653: while ( (c&1) == 0 ) c >>= 1;
654: *p = c;
655: return 1;
656: /**/
657: non0w:
658: u += i, lu -= i, v += i, lv -= i;
659: TRAILINGZEROS( t, rsh );
660: i = 0;
661: if ( rsh == 0 ) {
662: *p++ = t, l = 0;
663: for ( ; i < lv; c = hi ) {
664: if ( (w = u[i] - v[i] - b) < 0 ) w += BASE, b = 1;
665: else b = 0;
666: DMA( a, w, c, hi, t, not_used );
667: i++;
668: if ( (*p++ = t) != 0 ) l = i;
669: }
670: for ( ; i < lu; c = hi ) {
671: if ( (w = u[i] - b) < 0 ) w += BASE, b = 1;
672: else b = 1;
673: DMA( a, w, c, hi, t, not_used );
674: i++;
675: if ( (*p++ = t) != 0 ) l = i;
676: }
677: if ( b != 0 ) c -= a;
678: if ( c == 0 ) return( l + 1 );
679: *p = c;
680: return( i + 2 );
681: } else {
682: for ( lsh = BSH - rsh, l = 0; i < lv; c = hi, t = x >> rsh ) {
683: if ( (w = u[i] - v[i] - b) < 0 ) w += BASE, b = 1;
684: else b = 0;
685: DMA( a, w, c, hi, x, not_used );
686: i++;
687: if ( (*p++ = t | ((x << lsh) & BMASK)) != 0 ) l = i;
688: }
689: for ( ; i < lu; i++, c = hi, t = x >> rsh ) {
690: if ( (w = u[i] - b) < 0 ) w += BASE, b = 1;
691: else b = 0;
692: DMA( a, w, c, hi, x, not_used );
693: i++;
694: if ( (*p++ = t | ((x << lsh) & BMASK)) != 0 ) l = i;
695: }
696: if ( b != 0 ) c -= a;
697: if ( c != 0 ) {
698: *p++ = t | ((c << lsh) & BMASK);
699: i++, t = c >> rsh;
700: } else if ( t == 0 ) return l;
701: if ( t == 0 ) return( i );
702: *p = t;
703: return( i + 1 );
704: }
705: }
706: #endif /* FULLSET */
707:
708: static int abs_axU_bxV_maxrshift( a, u, lu, b, v, lv, ans )
709: int a, *u, lu, b, *v, lv, *ans;
710: /*
711: * Compute | a*U - b*V | >> (as much as possible) in ans[], where U=u[0:lu-1]
712: * and V=v[0:lv-1], 0 < a, b < 2^BASE.
713: */
714: {
715: int i, l, c, d, s, t, w, rsh, *p = ans, lsh;
716: int hi;
717: static int bw_int32();
718:
719: BitWidth( a, c ); BitWidth( u[lu -1], s );
720: BitWidth( b, d ); BitWidth( v[lv -1], t );
721: if ( lu < lv || lu == lv && (c + s) < (d + t) ) {
722: SWAP( a, b, int ); SWAP( lu, lv, int ); SWAP( u, v, int * );
723: }
724: for ( i = c = d = 0; i < lv; ) {
725: DMA( a, u[i], c, hi, t, not_used ); /* 0 <= c <= 2^BSH -2 */
726: c = hi;
727: DMA( b, v[i], d, hi, s, not_used ); /* 0 <= d <= 2^BSH -1 */
728: i++, d = hi;
729: if ( (t -= s) == 0 ) continue;
730: else if ( t < 0 ) {
731: t += BASE;
732: if ( c != 0 ) c--;
733: else d++;
734: }
735: goto non0w;
736: }
737: if ( i >= lu ) { /* no more elm.'s in u[] and v[] */
738: only_c:
739: if ( (c -= d) == 0 ) return 0;
740: else if ( c < 0 ) c = -c;
741: sh_onlyc:
742: /* TRAILINGZEROS( c, rsh ); */
743: while ( (c&1) == 0 ) c >>= 1;
744: *p = c;
745: return 1;
746: }
747: /* there is at least one more elm. in u[] */
748: DMA( a, u[i], c, hi, t, not_used );
749: i++, c = hi;
750: if ( i >= lu ) { /* in this case, diff still can be <= 0 */
751: if ( hi == 0 ) { c = t; goto only_c; }
752: if ( (t -= d) == 0 ) { c = hi; goto sh_onlyc; }
753: else if ( t < 0 ) t += BASE, hi--;
754: TRAILINGZEROS( t, rsh );
755: if ( rsh == 0 ) *p++ = t;
756: else {
757: *p++ = t | ((hi << (BSH - rsh)) & BMASK);
758: hi >>= rsh;
759: }
760: if ( hi == 0 ) return 1;
761: *p = hi;
762: return 2;
763: } /* diff turned out to be > 0 */ else if ( (t -= d) > 0 ) d = 0;
764: else if ( t < 0 ) t += BASE, d = 1;
765: else {
766: while ( i < lu ) {
767: DMA( a, u[i], c, hi, t, not_used );
768: i++, c = hi;
769: if ( t != 0 ) break;
770: }
771: }
772: u += i, lu -= i;
773: TRAILINGZEROS( t, rsh );
774: l = i = 0;
775: if ( rsh == 0 ) {
776: *p++ = t;
777: goto no_more_v;
778: } else goto no_more_v_sh;
779: /**/
780: non0w:
781: u += i, lu -= i, v += i, lv -= i;
782: TRAILINGZEROS( t, rsh );
783: i = 0;
784: if ( rsh == 0 ) {
785: *p++ = t, l = 0;
786: for ( ; i < lv; *p++ = t ) {
787: DMA( a, u[i], c, hi, t, not_used );
788: c = hi;
789: DMA( b, v[i], d, hi, s, not_used );
790: i++, d = hi;
791: if ( (t -= s) == 0 ) continue;
792: else if ( t < 0 ) {
793: t += BASE;
794: if ( c != 0 ) c--;
795: else d++;
796: }
797: l = i;
798: }
799: no_more_v:
800: if ( i >= lu ) {
801: final_c_d:
802: if ( (c -= d) == 0 ) return( l + 1 );
803: else if ( c < 0 ) {
804: l = negate_nbd( ans, i + 1 ); /* <================ */
805: if ( (c = -c - 1) == 0 ) return l;
806: }
807: *p = c;
808: return( i + 2 );
809: }
810: /* a*u[i:lu-1] + c - d */
811: if ( c >= d ) c -= d;
812: else {
813: DMA( a, u[i], c, hi, t, not_used );
814: if ( i+1 >= lu ) {
815: if ( hi == 0 ) { c = t; goto final_c_d; }
816: if ( (t -= d) < 0 ) t += BASE, c = hi - 1;
817: else c = hi;
818: i++, *p++ = t;
819: goto final_c;
820: }
821: i++, c = hi;
822: if ( (t -= d) >= 0 ) *p++ = t;
823: else {
824: *p++ = t + BASE;
825: if ( c != 0 ) c--;
826: else {
827: for ( ; i < lu; *p++ = BMASK ) {
828: DMA( a, u[i], c, hi, t, not_used );
829: i++, c = hi;
830: if ( t == 0 ) continue;
831: *p++ = t - 1;
832: goto aU;
833: }
834: c--;
835: goto final_c;
836: }
837: }
838: }
839: /*** return( aV_plus_c( a, &u[i], lu - i, c, p ) + i + 1 ); ***/
840: aU:
841: for ( ; i < lu; i++, c = hi, *p++ = t ) { DMA( a, u[i], c, hi, t, not_used ); }
842: final_c:
843: if ( c == 0 ) return( i + 1 );
844: *p = c;
845: return( i + 2 );
846: } else { /* rsh != 0 */
847: for ( lsh = BSH - rsh; i < lv; t = w >> rsh ) {
848: DMA( a, u[i], c, hi, w, not_used );
849: c = hi;
850: DMA( b, v[i], d, hi, s, not_used );
851: i++, d = hi;
852: if ( (w -= s) < 0 ) {
853: w += BASE;
854: if ( c != 0 ) c--;
855: else d++;
856: }
857: if ( (*p++ = t | ((w << lsh) & BMASK)) != 0 ) l = i;
858: }
859: no_more_v_sh:
860: if ( i >= lu ) {
861: final_c_d_sh:
862: if ( (c -= d) == 0 )
863: if ( t == 0 ) return l;
864: else { *p = t; return( i + 1 ); }
865: else if ( c < 0 ) {
866: if ( i == 0 ) t = (BASE >> rsh) - t;
867: else {
868: l = negate_nbd( ans, i ); /* <================ */
869: /* t = (BASE >> rsh) - t;*/
870: t ^= ((BMASK >> rsh));
871: }
872: c = -c -1;
873: }
874: *p++ = t | ((c << lsh) & BMASK);
875: if ( (c >>= rsh) == 0 ) return( i + 1 );
876: *p = c;
877: return( i + 2 );
878: } else if ( c >= d ) c -= d;
879: else {
880: DMA( a, u[i], c, hi, w, not_used );
881: if ( i+1 >= lu ) {
882: if ( hi == 0 ) { c = w; goto final_c_d_sh; }
883: if ( (w -= d) < 0 ) w += BASE, c = hi - 1;
884: else c = hi;
885: i++, *p++ = t | ((w << lsh) & BMASK);
886: t = w >> rsh;
887: goto final_c_sh;
888: }
889: i++, c = hi;
890: if ( (w -= d) >= 0 ) { *p++ = t | ((w << lsh) & BMASK); t = w >> rsh; }
891: else {
892: w += BASE;
893: *p++ = t | ((w << lsh) & BMASK);
894: t = w >> rsh;
895: if ( c != 0 ) c--;
896: else {
897: while ( i < lu ) {
898: DMA( a, u[i], c, hi, w, not_used );
899: i++, c = hi;
900: if ( w == 0 ) {
901: *p++ = t | ((BMASK << lsh) & BMASK);
902: t = BMASK >> rsh;
903: continue;
904: } else {
905: w--;
906: *p++ = t | ((w << lsh) & BMASK);
907: t = w >> rsh;
908: goto aU_sh;
909: }
910: }
911: c--;
912: goto final_c_sh;
913: }
914: }
915: }
916: /*** return( aV_plus_c_sh( a, &u[i], lu - i, c, rsh, t ) + i ); ***/
917: aU_sh:
918: for ( ; i < lu; i++, c = hi, t = w >> rsh ) {
919: DMA( a, u[i], c, hi, w, not_used );
920: *p++ = t | ((w << lsh) & BMASK);
921: }
922: final_c_sh:
923: if ( c != 0 ) {
924: *p++ = t | ((c << lsh) & BMASK);
925: i++, t = c >> rsh;
926: }
927: if ( t == 0 ) return i;
928: *p = t;
929: return( i + 1 );
930: }
931: }
932: /****************/
933: /****************/
934: #ifdef FULLSET
935: static int UplusV_maxrshift(), aUplusV_maxrshift(), axUplusV_maxrshift(), aUplusbV_maxrshift();
936:
937: static int aU_plus_bV_maxrshift( a, u, lu, b, v, lv, ans )
938: int a, *u, lu, b, *v, lv, *ans;
939: {
940: if ( a == 1 )
941: return( b == 1 ? UplusV_maxrshift( u, lu, v, lv, ans ) :
942: aUplusV_maxrshift( b, v, lv, u, lu, ans ) );
943: else if ( b == 1 ) return aUplusV_maxrshift( a, u, lu, v, lv, ans );
944: else if ( a == b ) return axUplusV_maxrshift( a, u, lu, v, lv, ans );
945: return aUplusbV_maxrshift( a, u, lu, b, v, lv, ans );
946: }
947:
948: static int UplusV_maxrshift( u, lu, v, lv, ans )
949: int *u, lu, *v, lv, *ans;
950: /*
951: * Compute ( (U + V) >> (as much as possible) ) in ans[], where U=u[0:lu-1],
952: * V=v[0:lv-1], and 0 < a < BASE.
953: */
954: {
955: int i, t, c, rsh, *p = ans;
956: int hi;
957:
958: if ( lu < lv ) { SWAP( lu, lv, int ); SWAP( u, v, int * ); }
959: for ( c = i = 0; i < lv; ) {
960: t = (c += (u[i] + v[i])) & BMASK;
961: /* 0 <= c <= 2*(2^BSH -1) + 1 = 2^BSH + (2^BSH -1) */
962: i++, c >>= BSH;
963: if ( t != 0 ) goto non0w;
964: }
965: while ( i < lu ) {
966: t = (c += u[i]) & BMASK;
967: i++, c >>= BSH;
968: if ( t != 0 ) goto non0w;
969: }
970: *p = c;
971: return 1;
972: /**/
973: non0w:
974: u += i, lu -= i, v += i, lv -= i;
975: TRAILINGZEROS( t, rsh );
976: i = 0;
977: if ( rsh == 0 ) {
978: *p++ = t;
979: for ( ; i < lv; i++, c >>= BSH ) *p++ = (c += (u[i] + v[i])) & BMASK;
980: if ( c != 0 ) {
981: while ( i < lu )
982: if ( (c = u[i++] + 1) >= BASE ) *p++ = c & BMASK;
983: else {
984: *p++ = c;
985: goto cpu;
986: }
987: *p = 1; /* == c */
988: return( lu + 2 );
989: }
990: cpu:
991: for ( ; i < lu; i++ ) *p++ = u[i];
992: return( lu + 1 );
993: } else {
994: int lsh = BSH - rsh;
995:
996: for ( ; i < lv; i++, c >>= BSH ) {
997: c += (u[i] + v[i]);
998: *p++ = t | ((c << lsh) & BMASK);
999: t = (c & BMASK) >> rsh;
1000: }
1001: if ( c != 0 ) {
1002: while ( i < lu ) {
1003: c = u[i++] + 1;
1004: *p++ = t | ((c << lsh) & BMASK);
1005: t = (c & BMASK) >> rsh;
1006: if ( (c >>= BSH) == 0 ) goto cpu_sh;
1007: }
1008: *p = t | (1 << lsh);
1009: return( lu + 1 );
1010: }
1011: cpu_sh:
1012: for ( ; i < lu; t = c >> rsh ) *p++ = t | (((c = u[i++]) << lsh) & BMASK);
1013: if ( t == 0 ) return( lu );
1014: *p = t;
1015: return( lu + 1 );
1016: }
1017: }
1018:
1019: static int aUplusV_maxrshift( a, u, lu, v, lv, ans )
1020: int a, *u, lu, *v, lv, *ans;
1021: /*
1022: * Compute ( (a*U + V) >> (as much as possible) ) in ans[], where U=u[0:lu-1],
1023: * V=v[0:lv-1], and 1 < a < BASE.
1024: */
1025: {
1026: int i, l, t, c, rsh, *p = ans, w, lsh;
1027: int hi;
1028:
1029: for ( l = lu < lv ? lu : lv, c = i = 0; i < l; ) {
1030: /* Note that (carry in a*U[*]) <= 2^BSH-2 because {(2^j-1)}^2 + (2^j-2)
1031: * = 2^j*(2^j -2)+(2^j -1), and the following c is <= 2^j -1 because
1032: * {(2^j -1)}^2 + (2^j -1) + (2^j -1) = 2^j*(2^j -1) + (2^j -1).
1033: * ^^^^^^carry^^^^^^ pv[?].
1034: */
1035: c += v[i];
1036: DMA( a, u[i], c, hi, t, not_used );
1037: i++, c = hi;
1038: if ( t != 0 ) goto non0sum;
1039: }
1040: if ( i >= lu ) { /* compute (c + v[i:lv-1]) */
1041: while ( i < lv ) {
1042: t = (c += v[i++]) & BMASK;
1043: c >>= BSH;
1044: if ( t != 0 ) goto non0w_v;
1045: }
1046: *p = 1;
1047: return 1;
1048: non0w_v:
1049: v += i, lv -= i;
1050: i = 0;
1051: TRAILINGZEROS( t, rsh );
1052: if ( rsh == 0 ) {
1053: *p++ = t;
1054: L_addv:
1055: if ( c != 0 ) {
1056: while ( i < lv )
1057: if ( (c = u[i++] + 1) >= BASE ) *p++ = c & BMASK;
1058: else {
1059: *p++ = c;
1060: goto cpv;
1061: }
1062: *p = 1; /* == c */
1063: return( lv + 2 );
1064: }
1065: cpv:
1066: for ( ; i < lv; i++ ) *p++ = v[i];
1067: return( lv + 1 );
1068: } else {
1069: lsh = BSH - rsh;
1070: L_addv_sh:
1071: if ( c != 0 ) {
1072: while ( i < lv ) {
1073: c += v[i++];
1074: *p++ = t | ((c << lsh) & BMASK);
1075: t = (c & BMASK) >> rsh;
1076: if ( (c >>= BSH) == 0 ) goto cpv_sh;
1077: }
1078: *p = t | (1 << lsh);
1079: return( lv + 1 );
1080: }
1081: cpv_sh:
1082: for ( ; i < lv; t = c >> rsh ) *p++ = t | (((c = v[i++]) << lsh) & BMASK);
1083: if ( t == 0 ) return( lv );
1084: *p = t;
1085: return( lv + 1 );
1086: }
1087: } else { /* i >= lv, and compute (c + a*u[i:lu-1]) */
1088: while ( i < lu ) {
1089: DMA( a, u[i], c, hi, t, not_used );
1090: i++, c = hi;
1091: if ( t != 0 ) goto non0w_u;
1092: }
1093: TRAILINGZEROS( c, rsh );
1094: *p = c;
1095: return 1;
1096: /**/
1097: non0w_u:
1098: u += i, lu -= i;
1099: TRAILINGZEROS( t, rsh );
1100: i = 0;
1101: if ( rsh == 0 ) {
1102: *p++ = t;
1103: L_addu:
1104: #ifdef CALL
1105: return( aV_plus_c( a, u, lv, c, p ) + 1 );
1106: #else
1107: for ( ; i < lu; i++, *p++ = t, c = hi ) { DMA( a, u[i], c, hi, t, not_used ); }
1108: if ( c == 0 ) return( lu + 1 );
1109: *p = c;
1110: return( lu + 2 );
1111: #endif
1112: } else {
1113: lsh = BSH - rsh;
1114: L_addu_sh:
1115: #ifdef CALL
1116: return aV_plus_c_sh( a, u, lu, c, rsh, t, p );
1117: #else
1118: for ( ; i < lu; i++, c = hi, t = w >> rsh ) {
1119: DMA( a, u[i], c, hi, w, not_used );
1120: *p++ = t | ((w << lsh) & BMASK);
1121: }
1122: if ( c != 0 ) {
1123: *p++ = t | ((c << lsh) & BMASK);
1124: i++, t = c >> rsh;
1125: }
1126: if ( t == 0 ) return i;
1127: *p = t;
1128: return( i + 1 );
1129: #endif
1130: }
1131: }
1132: /**/
1133: non0sum:
1134: u += i, lu -= i, v += i, lv -= i, l -= i;
1135: TRAILINGZEROS( t, rsh );
1136: i = 0;
1137: if ( rsh == 0 ) {
1138: *p++ = t;
1139: for ( ; i < l; i++, *p++ = t, c = hi ) {
1140: c += v[i];
1141: DMA( a, u[i], c, hi, t, not_used );
1142: }
1143: if ( i >= lu ) /* compute (c + v[i:lv-1]). note i may be >= lv */ goto L_addv;
1144: else /* i >= lv, and compute (c + u[i:lu-1]) */ goto L_addu;
1145: } else {
1146: lsh = BSH - rsh;
1147: for ( ; i < l; i++, c = hi, t = w >> rsh ) {
1148: c += v[i];
1149: DMA( a, u[i], c, hi, w, not_used );
1150: *p++ = t | ((w << lsh) & BMASK);
1151: }
1152: if ( i >= lu ) /* compute (c + v[i:lv-1]) >> rsh. note i may be >= lv */ goto L_addv_sh;
1153: else /* i >= lv, and compute (c + u[i:lu-1]) >> rsh */ goto L_addu_sh;
1154: }
1155: }
1156:
1157:
1158: static int axUplusV_maxrshift( a, u, lu, v, lv, ans )
1159: int a, *u, lu, *v, lv, *ans;
1160: /*
1161: * Compute ( a*(U + V) >> (as much as possible) ) in ans[], where U=u[0:lu-1],
1162: * V=v[0:lv-1], and 1 < a < BASE.
1163: */
1164: {
1165: int i, t, c, d, rsh, w, *p = ans, lsh, x;
1166: int hi;
1167:
1168: if ( lu < lv ) { SWAP( lu, lv, int ); SWAP( u, v, int * ); }
1169: for ( i = c = d = 0; i < lv; ) {
1170: w = (d += (u[i] + v[i])) & BMASK;
1171: DMA( a, w, c, hi, t, not_used );
1172: i++, d >>= BSH, c = hi;
1173: if ( t != 0 ) goto non0w;
1174: }
1175: while ( i < lu ) {
1176: w = (d += u[i++]) & BMASK;
1177: DMA( a, w, c, hi, t, not_used );
1178: d >>= BSH, c = hi;
1179: if ( t != 0 ) goto non0w;
1180: }
1181: if ( d != 0 ) c += a;
1182: TRAILINGZEROS( c, rsh );
1183: if ( rsh != 0 || c <= BMASK ) { *p = c; return 1; }
1184: *p++ = c & BMASK;
1185: *p = c >> BSH;
1186: return 2;
1187: /**/
1188: non0w:
1189: u += i, lu -= i, v += i, lv -= i;
1190: TRAILINGZEROS( t, rsh );
1191: i = 0;
1192: if ( rsh == 0 ) {
1193: *p++ = t;
1194: for ( ; i < lv; i++, d >>= BSH, *p++ = t, c = hi ) {
1195: w = (d += (u[i] + v[i])) & BMASK;
1196: DMA( a, w, c, hi, t, not_used );
1197: }
1198: for ( ; i < lu; d >>= BSH, *p++ = t, c = hi ) {
1199: w = (d += u[i++]) & BMASK;
1200: DMA( a, w, c, hi, t, not_used );
1201: }
1202: if ( d != 0 ) c += a;
1203: if ( c == 0 ) return( lu + 1 );
1204: *p++ = c & BMASK;
1205: if ( (c >>= BSH) == 0 ) return( lu + 2 );
1206: *p = c;
1207: return( lu + 3 );
1208: } else {
1209: for ( lsh = BSH - rsh; i < lv; i++, d >>= BSH, t = x >> rsh, c = hi ) {
1210: w = (d += (u[i] + v[i])) & BMASK;
1211: DMA( a, w, c, hi, x, not_used );
1212: *p++ = t | ((x << lsh) & BMASK);
1213: }
1214: for ( ; i < lu; d >>= BSH, t = x >> rsh, c = hi ) {
1215: w = (d += u[i++]) & BMASK;
1216: DMA( a, w, c, hi, x, not_used );
1217: *p++ = t | ((x << lsh) & BMASK);
1218: }
1219: if ( d != 0 ) c += a;
1220: t |= ((c << lsh) & BMASK);
1221: c >>= rsh;
1222: if ( c != 0 ) {
1223: *p++ = t;
1224: *p = c;
1225: return( lu + 2 );
1226: } else if ( t == 0 ) return lu;
1227: *p = t;
1228: return( lu + 1 );
1229: }
1230: }
1231: #endif /* FULLSET */
1232:
1233: static int aUplusbV_maxrshift( a, u, lu, b, v, lv, ans )
1234: int a, *u, lu, b, *v, lv, *ans;
1235: /*
1236: * Compute ( (a*U + b*V) >> (as much as possible) ) in ans[], where U=u[0:lu-1],
1237: * V=v[0:lv-1], and 0 < a,b < BASE.
1238: */
1239: {
1240: int i, c, d, s, t, w, rsh, *p = ans, lsh;
1241: int hi;
1242:
1243: if ( lu < lv ) { SWAP( a, b, int ); SWAP( lu, lv, int ); SWAP( u, v, int * ); }
1244: for ( c = d = i = 0; i < lv; ) {
1245: DMA( a, u[i], c, hi, t, not_used ); /* 0 <= c <= 2^BSH -1 */
1246: c = hi;
1247: DMA( b, v[i], d, hi, s, not_used ); /* 0 <= d <= 2^BSH -2 */
1248: i++, d = hi, t += s;
1249: if ( t > BMASK ) c++, t &= BMASK;
1250: if ( t != 0 ) goto non0w;
1251: }
1252: c += d;
1253: for ( d = 0; i < lu; ) {
1254: DMA( a, u[i], c, hi, t, not_used );
1255: i++, c = hi;
1256: if ( t != 0 ) goto non0w;
1257: }
1258: TRAILINGZEROS( c, rsh );
1259: if ( rsh != 0 || c <= BMASK ) { *p = c; return 1; }
1260: *p++ = c & BMASK;
1261: *p = 1;
1262: return 2;
1263: /**/
1264: non0w:
1265: u += i, lu -= i, v += i, lv -= i;
1266: TRAILINGZEROS( t, rsh );
1267: i = 0;
1268: if ( rsh == 0 ) {
1269: *p++ = t;
1270: for ( ; i < lv; i++, *p++ = t ) {
1271: DMA( a, u[i], c, hi, t, not_used );
1272: c = hi;
1273: DMA( b, v[i], d, hi, s, not_used );
1274: d = hi, t += s;
1275: if ( t > BMASK ) c++, t &= BMASK;
1276: }
1277: c += d;
1278: for ( ; i < lu; i++, *p++ = t, c = hi ) {
1279: DMA( a, u[i], c, hi, t, not_used );
1280: }
1281: if ( c == 0 ) return( lu + 1 );
1282: else if ( c <= BMASK ) { *p = c; return( lu + 2 ); }
1283: *p++ = c & BMASK;
1284: *p = 1;
1285: return( lu + 3 );
1286: } else {
1287: for ( lsh = BSH - rsh; i < lv; i++, t = w >> rsh ) {
1288: DMA( a, u[i], c, hi, w, not_used );
1289: c = hi;
1290: DMA( b, v[i], d, hi, s, not_used );
1291: d = hi, w += s;
1292: if ( w > BMASK ) c++, w &= BMASK;
1293: *p++ = t | ((w << lsh) & BMASK);
1294: }
1295: c += d; /* <= 2^BSH + (2^BSH - 3) */
1296: for ( ; i < lu; i++, c = hi, t = w >> rsh ) {
1297: DMA( a, u[i], c, hi, w, not_used );
1298: *p++ = t | ((w << lsh) & BMASK);
1299: }
1300: if ( c == 0 ) {
1301: if ( t == 0 ) return lu;
1302: *p = t;
1303: return( lu + 1 );
1304: }
1305: *p++ = t | ((c << lsh) & BMASK);
1306: if ( (c >>= rsh) == 0 ) return( lu + 1 );
1307: *p = c;
1308: return( lu + 2 );
1309: }
1310: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to igcdhack.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine