Annotation of OpenXM_contrib2/asir2000/engine/igcdhack.c, Revision 1.4
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.4 ! murao 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/igcdhack.c,v 1.3 2000/08/22 05:04:05 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: 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.1 noro 872: static int bw_int32();
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