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

Annotation of OpenXM_contrib2/asir2000/asm/ddN.c, Revision 1.2

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
        !            26:  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification
        !            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:  *
        !            48:  * $OpenXM: OpenXM_contrib2/asir2000/asm/ddN.c,v 1.1.1.1 1999/12/03 07:39:06 noro Exp $
        !            49: */
1.1       noro       50: #ifndef FBASE
                     51: #define FBASE
                     52: #endif
                     53:
                     54: #include "ca.h"
                     55: #include "base.h"
                     56: #include "inline.h"
                     57:
                     58: void bxprintn(N),bprintn(N);
                     59: void divnmain_special(int,int,unsigned int *,unsigned int *,unsigned int *);
                     60: void dupn(N,N);
                     61:
                     62: void divn(n1,n2,nq,nr)
                     63: N n1,n2,*nq,*nr;
                     64: {
                     65:        int tmp,b;
                     66:        int i,j,d1,d2,dd;
                     67:        unsigned int *m1,*m2;
                     68:        unsigned int uq,ur,msw;
                     69:        N q,r,t1,t2;
                     70:
                     71:        if ( !n2 )
                     72:                error("divn: divsion by 0");
                     73:        else if ( !n1 ) {
                     74:                *nq = 0; *nr = 0;
                     75:        } else if ( UNIN(n2) ) {
                     76:                COPY(n1,*nq); *nr = 0;
                     77:        } else if ( (d1 = PL(n1)) < (d2 = PL(n2)) ) {
                     78:                COPY(n1,*nr); *nq = 0;
                     79:        } else if ( d2 == 1 ) {
                     80:                if ( d1 == 1 ) {
                     81:                        DQR(BD(n1)[0],BD(n2)[0],uq,ur)
                     82:                        STON(uq,*nq); STON(ur,*nr);
                     83:                } else {
                     84:                        ur = divin(n1,BD(n2)[0],nq); STON(ur,*nr);
                     85:                }
                     86:                return;
                     87:        } else if ( (d1 == d2) && ((tmp = cmpn(n1,n2)) <= 0) )
                     88:                if ( !tmp ) {
                     89:                        *nr = 0; COPY(ONEN,*nq);
                     90:                } else {
                     91:                        COPY(n1,*nr); *nq = 0;
                     92:                }
                     93:        else {
                     94:                msw = BD(n2)[d2-1];
                     95:                for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                     96:                b = BSH-i;
                     97:                if ( b ) {
                     98:                        bshiftn(n1,-b,&t1); bshiftn(n2,-b,&t2);
                     99:                } else {
                    100:                        COPY(n1,t1); COPY(n2,t2);
                    101:                }
                    102:                m1 = (unsigned int *)ALLOCA((PL(t1)+1)*sizeof(unsigned int));
                    103:                bcopy(BD(t1),m1,PL(t1)*sizeof(unsigned int)); m1[PL(t1)] = 0;
                    104:                m2 = BD(t2); dd = d1-d2;
                    105:                *nq = q = NALLOC(dd+1); INITRC(q);
                    106:                divnmain(d1,d2,m1,m2,BD(q)); FREEN(t1); FREEN(t2);
                    107:                PL(q) = (BD(q)[dd]?dd+1:dd);
                    108:                if ( !PL(q) ) {
                    109:                        FREEN(q);
                    110:                }
                    111:                for ( j = d2-1; (j >= 0) && (m1[j] == 0); j--);
                    112:                if ( j == -1 )
                    113:                        *nr = 0;
                    114:                else {
                    115:                        r = NALLOC(j+1); INITRC(r); PL(r) = j+1;
                    116:                        bcopy(m1,BD(r),(j+1)*sizeof(unsigned int));
                    117:                        if ( b ) {
                    118:                                bshiftn(r,b,nr); FREEN(r);
                    119:                        } else
                    120:                                *nr = r;
                    121:                }
                    122:        }
                    123: }
                    124:
                    125: void divsn(n1,n2,nq)
                    126: N n1,n2,*nq;
                    127: {
                    128:        int d1,d2,dd,i,b;
                    129:        unsigned int *m1,*m2;
                    130:        unsigned int uq,msw;
                    131:        N q,t1,t2;
                    132:
                    133:        if ( !n1 )
                    134:                *nq = 0;
                    135:        else if ( UNIN(n2) ) {
                    136:                COPY(n1,*nq);
                    137:        } else if ( (d1 = PL(n1)) < (d2 = PL(n2)) )
                    138:                error("divsn: cannot happen");
                    139:        else if ( d2 == 1 ) {
                    140:                if ( d1 == 1 ) {
                    141:                        uq = BD(n1)[0] / BD(n2)[0]; STON(uq,*nq);
                    142:                } else
                    143:                        divin(n1,BD(n2)[0],nq);
                    144:        } else {
                    145:                msw = BD(n2)[d2-1];
                    146:                for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                    147:                b = BSH-i;
                    148:                if ( b ) {
                    149:                        bshiftn(n1,-b,&t1); bshiftn(n2,-b,&t2);
                    150:                } else {
                    151:                        COPY(n1,t1); COPY(n2,t2);
                    152:                }
                    153:                m1 = (unsigned int *)ALLOCA((PL(t1)+1)*sizeof(unsigned int));
                    154:                bcopy(BD(t1),m1,PL(t1)*sizeof(unsigned int)); m1[PL(t1)] = 0;
                    155:                m2 = BD(t2); dd = d1-d2;
                    156:                *nq = q = NALLOC(dd+1); INITRC(q);
                    157:                divnmain(d1,d2,m1,m2,BD(q));
                    158:                FREEN(t1); FREEN(t2);
                    159:                PL(q) = (BD(q)[dd]?dd+1:dd);
                    160:                if ( !PL(q) )
                    161:                        FREEN(q);
                    162:        }
                    163: }
                    164:
                    165: void remn(n1,n2,nr)
                    166: N n1,n2,*nr;
                    167: {
                    168:        int d1,d2,tmp;
                    169:        unsigned int uq,ur;
                    170:        N q;
                    171:
                    172:        if ( !n2 )
                    173:                error("remn: divsion by 0");
                    174:        else if ( !n1 || UNIN(n2) )
                    175:                *nr = 0;
                    176:        else if ( (d1 = PL(n1)) < (d2 = PL(n2)) ) {
                    177:                COPY(n1,*nr);
                    178:        } else if ( d2 == 1 ) {
                    179:                if ( d1 == 1 ) {
                    180:                        DQR(BD(n1)[0],BD(n2)[0],uq,ur)
                    181:                        STON(ur,*nr);
                    182:                } else {
                    183:                        ur = rem(n1,BD(n2)[0]); UTON(ur,*nr);
                    184:                }
                    185:                return;
                    186:        } else if ( (d1 == d2) && ((tmp = cmpn(n1,n2)) <= 0) )
                    187:                if ( !tmp ) {
                    188:                        *nr = 0;
                    189:                } else {
                    190:                        COPY(n1,*nr);
                    191:                }
                    192:        else
                    193:                divn(n1,n2,&q,nr);
                    194: }
                    195:
                    196: /* d = 2^(32*words)-lower */
                    197:
                    198: void remn_special(a,d,bits,lower,b)
                    199: N a,d;
                    200: int bits;
                    201: unsigned int lower;
                    202: N *b;
                    203: {
                    204:        int words;
                    205:        N r;
                    206:        int rl,i;
                    207:        unsigned int c,tmp;
                    208:        unsigned int *rb,*src,*dst;
                    209:
                    210:        if ( cmpn(a,d) < 0 ) {
                    211:                *b = a;
                    212:                return;
                    213:        }
                    214:        words = bits>>5;
                    215:        r = NALLOC(PL(a)); dupn(a,r);
                    216:        rl = PL(r); rb = BD(r);
                    217:        while ( rl > words ) {
                    218:                src = rb+words;
                    219:                dst = rb;
                    220:                i = rl-words;
                    221:                for ( c = 0; i > 0; i--, src++, dst++ ) {
                    222:                        DMA2(*src,lower,c,*dst,c,*dst)
                    223:                        *src = 0;
                    224:                }
                    225:                for ( ; c; dst++ ) {
                    226:                        tmp = *dst + c;
                    227:                        c = tmp < c ? 1 : 0;
                    228:                        *dst = tmp;
                    229:                }
                    230:                rl = dst-rb;
                    231:        }
                    232:        for ( i = words-1; i >= 0 && !rb[i]; i-- );
                    233:        PL(r) = i + 1;
                    234:        if ( cmpn(r,d) >= 0 ) {
                    235:                tmp = rb[0]-BD(d)[0];
                    236:                UTON(tmp,*b);
                    237:        } else
                    238:                *b = r;
                    239: }
                    240: void mulin(n,d,p)
                    241: N n;
                    242: unsigned int d;
                    243: unsigned int *p;
                    244: {
                    245:        unsigned int carry;
                    246:        unsigned int *m,*me;
                    247:
                    248:        bzero((char *)p,(int)((PL(n)+1)*sizeof(int)));
                    249:        for ( carry = 0, m = BD(n), me = m+PL(n); m < me; p++, m++ ) {
                    250:                DMA2(*m,d,*p,carry,carry,*p)
                    251:        }
                    252:        *p = carry;
                    253: }
                    254:
                    255: unsigned int divin(n,dvr,q)
                    256: N n;
                    257: unsigned int dvr;
                    258: N *q;
                    259: {
                    260:        int d;
                    261:        unsigned int up;
                    262:        unsigned int *m,*mq,*md;
                    263:        N nq;
                    264:
                    265:        d = PL(n); m = BD(n);
                    266:        if ( ( d == 1 ) && ( dvr > *m ) ) {
                    267:                *q = 0;
                    268:                return *m;
                    269:        }
                    270:        *q = nq = NALLOC(d); INITRC(nq);
                    271:        for ( md = m+d-1, mq = BD(nq)+d-1, up = 0; md >= m; md--, mq-- ) {
                    272:                DSAB(dvr,up,*md,*mq,up)
                    273:        }
                    274:        PL(nq) = (BD(nq)[d-1]?d:d-1);
                    275:        if ( !PL(nq) )
                    276:                FREEN(nq);
                    277:        return ( up );
                    278: }
                    279:
                    280: void bprintn(n)
                    281: N n;
                    282: {
                    283:        int l,i;
                    284:        unsigned int *b;
                    285:
                    286:        if ( !n )
                    287:                printf("0");
                    288:        else {
                    289:                l = PL(n); b = BD(n);
                    290:                for ( i = l-1; i >= 0; i-- )
                    291:                        printf("%010u|",b[i]);
                    292:        }
                    293: }
                    294:
                    295: void bxprintn(n)
                    296: N n;
                    297: {
                    298:        int l,i;
                    299:        unsigned int *b;
                    300:
                    301:        if ( !n )
                    302:                printf("0");
                    303:        else {
                    304:                l = PL(n); b = BD(n);
                    305:                for ( i = l-1; i >= 0; i-- )
                    306:                        printf("%08x|",b[i]);
                    307:        }
                    308: }
                    309:
                    310: #if defined(VISUAL) || defined(i386)
                    311: void muln(n1,n2,nr)
                    312: N n1,n2,*nr;
                    313: {
                    314:        unsigned int tmp,carry,mul;
                    315:        unsigned int *p1,*m1,*m2;
                    316:        int i,d1,d2,d;
                    317:        N r;
                    318:
                    319:        if ( !n1 || !n2 )
                    320:                *nr = 0;
                    321:        else if ( UNIN(n1) )
                    322:                COPY(n2,*nr);
                    323:        else if ( UNIN(n2) )
                    324:                COPY(n1,*nr);
                    325:        else if ( (PL(n1) == 1) && (PL(n2) == 1) ) {
                    326:                DM(*BD(n1),*BD(n2),carry,tmp)
                    327:                if ( carry ) {
                    328:                        *nr = r = NALLOC(2); INITRC(r);
                    329:                        PL(r) = 2; p1 = BD(r); *p1++ = tmp; *p1 = carry;
                    330:                } else
                    331:                        STON(tmp,*nr);
                    332:        } else {
                    333:                d1 = PL(n1); d2 = PL(n2);
                    334:                d = d1+d2;
                    335:                *nr = r = NALLOC(d); INITRC(r);
                    336:                bzero((char *)BD(r),(int)((d1+d2)*sizeof(int)));
                    337:                for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
                    338:                        if ( mul = *m2 )
                    339:                                muln_1(m1,d1,mul,BD(r)+i);
                    340:                PL(r) = (BD(r)[d-1]?d:d-1);
                    341:        }
                    342: }
                    343:
                    344: void _muln(n1,n2,nr)
                    345: N n1,n2,nr;
                    346: {
                    347:        unsigned int mul;
                    348:        unsigned int *m1,*m2;
                    349:        int i,d1,d2,d;
                    350:
                    351:        if ( !n1 || !PL(n1) || !n2 || !PL(n2) )
                    352:                PL(nr) = 0;
                    353:        else if ( UNIN(n1) )
                    354:                dupn(n2,nr);
                    355:        else if ( UNIN(n2) )
                    356:                dupn(n1,nr);
                    357:        else {
                    358:                d1 = PL(n1); d2 = PL(n2);
                    359:                d = d1+d2;
                    360:                bzero((char *)BD(nr),(int)((d1+d2)*sizeof(int)));
                    361:                for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
                    362:                        if ( mul = *m2 )
                    363:                                muln_1(m1,d1,mul,BD(nr)+i);
                    364:                PL(nr) = (BD(nr)[d-1]?d:d-1);
                    365:        }
                    366: }
                    367:
                    368: void muln_1(p,s,d,r)
                    369: unsigned int *p;
                    370: int s;
                    371: unsigned int d;
                    372: unsigned int *r;
                    373: {
                    374:        /* esi : p, edi : r, carry : ebx, s : ecx */
                    375: #if defined(VISUAL)
                    376:        __asm {
                    377:        push esi
                    378:        push edi
                    379:        mov esi,p
                    380:        mov edi,r
                    381:        mov ecx,s
                    382:        xor ebx,ebx
                    383:        Lstart_muln:
                    384:        mov eax,DWORD PTR [esi]
                    385:        mul d
                    386:        add eax,DWORD PTR [edi]
                    387:        adc edx,0
                    388:        add eax,ebx
                    389:        adc edx,0
                    390:        mov DWORD PTR[edi],eax
                    391:        mov ebx,edx
                    392:        lea     esi,DWORD PTR [esi+4]
                    393:        lea     edi,DWORD PTR [edi+4]
                    394:        dec ecx
                    395:        jnz Lstart_muln
                    396:        mov DWORD PTR [edi],ebx
                    397:        pop edi
                    398:        pop esi
                    399:        }
                    400: #else
                    401:        asm volatile("\
                    402:        movl    %0,%%esi;\
                    403:        movl    %1,%%edi;\
                    404:        movl    $0,%%ebx;\
                    405:        Lstart_muln:;\
                    406:        movl    (%%esi),%%eax;\
                    407:        mull    %2;\
                    408:        addl    (%%edi),%%eax;\
                    409:        adcl    $0,%%edx;\
                    410:        addl    %%ebx,%%eax;\
                    411:        adcl    $0,%%edx;\
                    412:        movl    %%eax,(%%edi);\
                    413:        movl    %%edx,%%ebx;\
                    414:        leal    4(%%esi),%%esi;\
                    415:        leal    4(%%edi),%%edi;\
                    416:        decl    %3;\
                    417:        jnz             Lstart_muln;\
                    418:        movl    %%ebx,(%%edi)"\
                    419:        :\
                    420:        :"m"(p),"m"(r),"m"(d),"m"(s)\
                    421:        :"eax","ebx","edx","esi","edi");
                    422: #endif
                    423: }
                    424:
                    425: void divnmain(d1,d2,m1,m2,q)
                    426: int d1,d2;
                    427: unsigned int *m1,*m2,*q;
                    428: {
                    429:        int i,j;
                    430:        UL r,ltmp;
                    431:        unsigned int l,ur;
                    432:        unsigned int *n1,*n2;
                    433:        unsigned int u,qhat;
                    434:        unsigned int v1,v2;
                    435:
                    436:        v1 = m2[d2-1]; v2 = m2[d2-2];
                    437: #if 0
                    438:        if ( v1 == (1<<31) && !v2 ) {
                    439:                divnmain_special(d1,d2,m1,m2,q);
                    440:                return;
                    441:        }
                    442: #endif
                    443:        for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
                    444:                n1 = m1+d2; n2 = m1+d2-1;
                    445:                if ( *n1 == v1 ) {
                    446:                        qhat = ULBASE - 1;
                    447:                        r = (UL)*n1 + (UL)*n2;
                    448:                } else {
                    449:                        DSAB(v1,*n1,*n2,qhat,ur)
                    450:                        r = (UL)ur;
                    451:                }
                    452:                DM(v2,qhat,u,l)
                    453:                while ( 1 ) {
                    454:                        if ((r > (UL)u) || ((r == (UL)u) && (*(n1-2) >= l)))
                    455:                                break;
                    456:                        if ( l >= v2 )
                    457:                                l -= v2;
                    458:                        else {
                    459:                                l = (unsigned int)((UL)l+(ULBASE-(UL)v2)); u--;
                    460:                        }
                    461:                r += v1; qhat--;
                    462:                }
                    463:                if ( qhat ) {
                    464:                        u = divn_1(m2,d2,qhat,m1);
                    465:                        if ( *(m1+d2) < u ) {
                    466:                                for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
                    467:                                        i > 0; i--, n1++, n2++ ) {
                    468:                                        ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
                    469:                                        *n1 = (unsigned int)(ltmp & BMASK);
                    470:                                        ur = (unsigned int)(ltmp >> BSH);
                    471:                                }
                    472:                        }
                    473:                        *n1 = 0;
                    474:                }
                    475:                *q = qhat;
                    476:        }
                    477: }
                    478:
                    479: void divnmain_special(d1,d2,m1,m2,q)
                    480: int d1,d2;
                    481: unsigned int *m1,*m2,*q;
                    482: {
                    483:        int i,j;
                    484:        UL ltmp;
                    485:        unsigned int ur;
                    486:        unsigned int *n1,*n2;
                    487:        unsigned int u,qhat;
                    488:        unsigned int v1,v2;
                    489:
                    490:        v1 = m2[d2-1]; v2 = 0;
                    491:        for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
                    492:                n1 = m1+d2; n2 = m1+d2-1;
                    493:                qhat = ((*n1)<<1)|((*n2)>>31);
                    494:                if ( qhat ) {
                    495:                        u = divn_1(m2,d2,qhat,m1);
                    496:                        if ( *(m1+d2) < u ) {
                    497:                                for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
                    498:                                        i > 0; i--, n1++, n2++ ) {
                    499:                                        ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
                    500:                                        *n1 = (unsigned int)(ltmp & BMASK);
                    501:                                        ur = (unsigned int)(ltmp >> BSH);
                    502:                                }
                    503:                        }
                    504:                        *n1 = 0;
                    505:                }
                    506:                *q = qhat;
                    507:        }
                    508: }
                    509:
                    510: unsigned int divn_1(p,s,d,r)
                    511: unsigned int *p;
                    512: int s;
                    513: unsigned int d;
                    514: unsigned int *r;
                    515: {
                    516: /*
                    517:        unsigned int borrow,l;
                    518:
                    519:        for ( borrow = 0; s--; p++, r++ ) {
                    520:                DMA(*p,d,borrow,borrow,l)
                    521:                if ( *r >= l )
                    522:                        *r -= l;
                    523:                else {
                    524:                        *r = (unsigned int)((UL)*r+(ULBASE-(UL)l)); borrow++;
                    525:                }
                    526:        }
                    527:        return borrow;
                    528: */
                    529:        /* esi : p, edi : r, borrow : ebx, s : ecx */
                    530: #if defined(VISUAL)
                    531:        __asm {
                    532:        push esi
                    533:        push edi
                    534:        mov esi,p
                    535:        mov edi,r
                    536:        mov ecx,s
                    537:        xor ebx,ebx
                    538:        Lstart_divn:
                    539:        mov eax,DWORD PTR [esi]
                    540:        mul d
                    541:        add eax,ebx
                    542:        adc edx,0
                    543:        sub DWORD PTR [edi],eax
                    544:        adc edx,0
                    545:        mov ebx,edx
                    546:        lea     esi,DWORD PTR [esi+4]
                    547:        lea     edi,DWORD PTR [edi+4]
                    548:        dec ecx
                    549:        jnz Lstart_divn
                    550:        mov eax,ebx
                    551:        pop edi
                    552:        pop esi
                    553:        }
                    554:        /* return value is in eax. */
                    555: #else
                    556:        unsigned int borrow;
                    557:
                    558:        asm volatile("\
                    559:        movl    %1,%%esi;\
                    560:        movl    %2,%%edi;\
                    561:        movl    $0,%%ebx;\
                    562:        Lstart_divn:;\
                    563:        movl    (%%esi),%%eax;\
                    564:        mull    %3;\
                    565:        addl    %%ebx,%%eax;\
                    566:        adcl    $0,%%edx;\
                    567:        subl    %%eax,(%%edi);\
                    568:        adcl    $0,%%edx;\
                    569:        movl    %%edx,%%ebx;\
                    570:        leal    4(%%esi),%%esi;\
                    571:        leal    4(%%edi),%%edi;\
                    572:        decl    %4;\
                    573:        jnz             Lstart_divn;\
                    574:        movl    %%ebx,%0"\
                    575:        :"=m"(borrow)\
                    576:        :"m"(p),"m"(r),"m"(d),"m"(s)\
                    577:        :"eax","ebx","edx","esi","edi");
                    578:
                    579:        return borrow;
                    580: #endif
                    581: }
                    582:
                    583: #else
                    584:
                    585: void muln(n1,n2,nr)
                    586: N n1,n2,*nr;
                    587: {
                    588:        unsigned int tmp,carry,mul;
                    589:        unsigned int *p1,*pp,*m1,*m2;
                    590:        int i,j,d1,d2;
                    591:        N r;
                    592:
                    593:        if ( !n1 || !n2 )
                    594:                *nr = 0;
                    595:        else if ( UNIN(n1) )
                    596:                COPY(n2,*nr);
                    597:        else if ( UNIN(n2) )
                    598:                COPY(n1,*nr);
                    599:        else if ( (PL(n1) == 1) && (PL(n2) == 1) ) {
                    600:                DM(*BD(n1),*BD(n2),carry,tmp)
                    601:                if ( carry ) {
                    602:                        *nr = r = NALLOC(2); INITRC(r);
                    603:                        PL(r) = 2; p1 = BD(r); *p1++ = tmp; *p1 = carry;
                    604:                } else
                    605:                        STON(tmp,*nr);
                    606:        } else {
                    607:                d1 = PL(n1); d2 = PL(n2);
                    608:                *nr = r = NALLOC(d1+d2); INITRC(r);
                    609:                bzero((char *)BD(r),(int)((d1+d2)*sizeof(int)));
                    610:                for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
                    611:                        if ( mul = *m2 ) {
                    612:                                for ( j = d1, carry = 0, p1 = m1, pp = BD(r)+i;
                    613:                                        j--; pp++ ) {
                    614:                                        DMA2(*p1++,mul,*pp,carry,carry,*pp)
                    615:                                }
                    616:                                *pp = carry;
                    617:                        }
                    618:                PL(r) = (carry?d1+d2:d1+d2-1);
                    619:        }
                    620: }
                    621:
                    622: void _muln(n1,n2,nr)
                    623: N n1,n2,nr;
                    624: {
                    625:        unsigned int tmp,carry,mul;
                    626:        unsigned int *p1,*pp,*m1,*m2;
                    627:        int i,j,d1,d2;
                    628:
                    629:        if ( !n1 || !PL(n1) || !n2 || !PL(n2) )
                    630:                PL(nr) = 0;
                    631:        else if ( UNIN(n1) )
                    632:                dupn(n2,nr);
                    633:        else if ( UNIN(n2) )
                    634:                dupn(n1,nr);
                    635:        else {
                    636:                d1 = PL(n1); d2 = PL(n2);
                    637:                bzero((char *)BD(nr),(int)((d1+d2)*sizeof(int)));
                    638:                for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
                    639:                        if ( mul = *m2 ) {
                    640:                                for ( j = d1, carry = 0, p1 = m1, pp = BD(nr)+i;
                    641:                                        j--; pp++ ) {
                    642:                                        DMA2(*p1++,mul,*pp,carry,carry,*pp)
                    643:                                }
                    644:                                *pp = carry;
                    645:                        }
                    646:                PL(nr) = (carry?d1+d2:d1+d2-1);
                    647:        }
                    648: }
                    649:
                    650: /* r[0...s] = p[0...s-1]*d */
                    651:
                    652: void muln_1(p,s,d,r)
                    653: unsigned int *p;
                    654: int s;
                    655: unsigned int d;
                    656: unsigned int *r;
                    657: {
                    658:        unsigned int carry;
                    659:
                    660:        for ( carry = 0; s--; p++, r++ ) {
                    661:                DMA2(*p,d,carry,*r,carry,*r)
                    662:        }
                    663:        *r = carry;
                    664: }
                    665:
                    666: void divnmain(d1,d2,m1,m2,q)
                    667: int d1,d2;
                    668: unsigned int *m1,*m2,*q;
                    669: {
                    670:        int i,j;
                    671:        UL r,ltmp;
                    672:        unsigned int l,ur,tmp;
                    673:        unsigned int *n1,*n2;
                    674:        unsigned int u,qhat;
                    675:        unsigned int v1,v2;
                    676:
                    677:        v1 = m2[d2-1]; v2 = m2[d2-2];
                    678:        for ( j = d1-d2, m1 += j, q += j; j >= 0; j--, q--, m1-- ) {
                    679:                n1 = m1+d2; n2 = m1+d2-1;
                    680:                if ( *n1 == v1 ) {
                    681:                        qhat = ULBASE - 1;
                    682:                        r = (UL)*n1 + (UL)*n2;
                    683:                } else {
                    684:                        DSAB(v1,*n1,*n2,qhat,ur)
                    685:                        r = (UL)ur;
                    686:                }
                    687:                DM(v2,qhat,u,l)
                    688:                while ( 1 ) {
                    689:                        if ((r > (UL)u) || ((r == (UL)u) && (*(n1-2) >= l)))
                    690:                                break;
                    691:                        if ( l >= v2 )
                    692:                                l -= v2;
                    693:                        else {
                    694:                                l = (unsigned int)((UL)l+(ULBASE-(UL)v2)); u--;
                    695:                        }
                    696:                r += v1; qhat--;
                    697:                }
                    698:                if ( qhat ) {
                    699:                        u = divn_1(m2,d2,qhat,m1);
                    700:                        if ( *(m1+d2) < u ) {
                    701:                                for ( i = d2, qhat--, ur = 0, n1 = m1, n2 = m2;
                    702:                                        i > 0; i--, n1++, n2++ ) {
                    703:                                        ltmp = (UL)*n1 + (UL)*n2 + (UL)ur;
                    704:                                        *n1 = (unsigned int)(ltmp & BMASK);
                    705:                                        ur = (unsigned int)(ltmp >> BSH);
                    706:                                }
                    707:                        }
                    708:                        *n1 = 0;
                    709:                }
                    710:                *q = qhat;
                    711:        }
                    712: }
                    713:
                    714: unsigned int divn_1(p,s,d,r)
                    715: unsigned int *p;
                    716: int s;
                    717: unsigned int d;
                    718: unsigned int *r;
                    719: {
                    720:        unsigned int borrow,l;
                    721:
                    722:        for ( borrow = 0; s--; p++, r++ ) {
                    723:                DMA(*p,d,borrow,borrow,l)
                    724:                if ( *r >= l )
                    725:                        *r -= l;
                    726:                else {
                    727:                        *r = (unsigned int)((UL)*r+(ULBASE-(UL)l)); borrow++;
                    728:                }
                    729:        }
                    730:        return borrow;
                    731: }
                    732: #endif

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