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

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

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