[BACK]Return to igcdhack.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Annotation of OpenXM_contrib2/asir2000/engine/igcdhack.c, Revision 1.1.1.1

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>