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

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

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