[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.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>