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