Return to divbig.c CVS log

Up to [CSRG BSD Unix] / 42BSD / ucb / lisp / franz

Annotation of 42BSD/ucb/lisp/franz/divbig.c, revision 1.1

1.1     ! root        1: #ifndef lint
        !             2: static char *rcsid =
        !             3:    "$Header: divbig.c,v 1.3 83/09/12 14:17:07 sklower Exp $";
        !             4: #endif
        !             5: 
        !             6: /*                                     -[Sat Jan 29 12:22:36 1983 by jkf]-
        !             7:  *     divbig.c                                $Locker:  $
        !             8:  * bignum division
        !             9:  *
        !            10:  * (c) copyright 1982, Regents of the University of California
        !            11:  */
        !            12: 
        !            13: 
        !            14: #include "global.h"
        !            15: 
        !            16: #define b 0x40000000
        !            17: #define toint(p) ((int) (p))
        !            18: 
        !            19: divbig(dividend, divisor, quotient, remainder)
        !            20: lispval dividend, divisor, *quotient, *remainder;
        !            21: {
        !            22:        register *ujp, *vip;
        !            23:        int *alloca(), d, negflag = 0, m, n, carry, rem, qhat, j;
        !            24:        int borrow, negrem = 0;
        !            25:        long *utop = sp(), *ubot, *vbot, *qbot;
        !            26:        register lispval work; lispval export();
        !            27:        Keepxs();
        !            28: 
        !            29:        /* copy dividend */
        !            30:        for(work = dividend; work; work = work ->s.CDR)
        !            31:                stack(work->s.I);
        !            32:        ubot = sp();
        !            33:        if(*ubot < 0) {         /* knuth's division alg works only for pos
        !            34:                                        bignums                         */
        !            35:                negflag ^= 1;
        !            36:                negrem = 1;
        !            37:                dsmult(utop-1,ubot,-1);
        !            38:        }
        !            39:        stack(0);
        !            40:        ubot = sp();
        !            41: 
        !            42:        
        !            43:        /*copy divisor */
        !            44:        for(work = divisor; work; work = work->s.CDR)
        !            45:                stack(work->s.I);
        !            46: 
        !            47:        vbot = sp();
        !            48:        stack(0);
        !            49:        if(*vbot < 0) {
        !            50:                negflag ^= 1;
        !            51:                dsmult(ubot-1,vbot,-1);
        !            52:        }
        !            53: 
        !            54:        /* check validity of data */
        !            55:        n = ubot - vbot;
        !            56:        m = utop - ubot - n - 1;
        !            57:        if (n == 1) {
        !            58:                /* do destructive division by  a single. */
        !            59:                rem = dsdiv(utop-1,ubot,*vbot);
        !            60:                if(negrem)
        !            61:                        rem = -rem;
        !            62:                if(negflag)
        !            63:                        dsmult(utop-1,ubot,-1);
        !            64:                if(remainder)
        !            65:                        *remainder = inewint(rem);
        !            66:                if(quotient)
        !            67:                        *quotient = export(utop,ubot);
        !            68:                Freexs();
        !            69:                return;
        !            70:        }
        !            71:        if (m < 0) {
        !            72:                if (remainder)
        !            73:                        *remainder = dividend;
        !            74:                if(quotient)
        !            75:                        *quotient = inewint(0);
        !            76:                Freexs();
        !            77:                return;
        !            78:        }
        !            79:        qbot = alloca(toint(utop) + toint(vbot) - 2 * toint(ubot));
        !            80: d1:
        !            81:        d = b /(*vbot +1);
        !            82:        dsmult(utop-1,ubot,d);
        !            83:        dsmult(ubot-1,vbot,d);
        !            84: 
        !            85: d2:    for(j=0,ujp=ubot; j <= m; j++,ujp++) {
        !            86: 
        !            87:        d3:     
        !            88:                qhat = calqhat(ujp,vbot);
        !            89:        d4:
        !            90:                if((borrow = mlsb(ujp + n, ujp, ubot, -qhat)) < 0) {
        !            91:                        adback(ujp + n, ujp, ubot);
        !            92:                        qhat--;
        !            93:                }
        !            94:                qbot[j] = qhat;
        !            95:        }
        !            96: d8:    if(remainder) {
        !            97:                dsdiv(utop-1, utop - n, d);
        !            98:                if(negrem) dsmult(utop-1,utop-n,-1);
        !            99:                *remainder = export(utop,utop-n);
        !           100:        }
        !           101:        if(quotient) {
        !           102:                if(negflag)
        !           103:                        dsmult(qbot+m,qbot,-1);
        !           104:                *quotient = export(qbot + m + 1, qbot);
        !           105:        }
        !           106:        Freexs();
        !           107: }
        !           108: /*
        !           109:  * asm code commented out due to optimizer bug
        !           110:  * also, this file is now shared with the 68k version!
        !           111: calqhat(ujp,v1p)
        !           112: register int *ujp, *v1p;
        !           113: {
        !           114: asm("  cmpl    (r10),(r11)             # v[1] == u[j] ??");
        !           115: asm("  beql    2f                      ");
        !           116: asm("  # calculate qhat and rhat simultaneously,");
        !           117: asm("  #  qhat in r0");
        !           118: asm("  #  rhat in r1");
        !           119: asm("  emul    (r11),$0x40000000,4(r11),r4 # u[j]b+u[j+1] into r4,r5");
        !           120: asm("  ediv    (r10),r4,r0,r1          # qhat = ((u[j]b+u[j+1])/v[1]) into r0");
        !           121: asm("                                  # (u[j]b+u[j+1] -qhat*v[1]) into r1");
        !           122: asm("                                  # called rhat");
        !           123: asm("1:");
        !           124: asm("  # check if v[2]*qhat > rhat*b+u[j+2]");
        !           125: asm("  emul    r0,4(r10),$0,r2         # qhat*v[2] into r3,r2");
        !           126: asm("  emul    r1,$0x40000000,8(r11),r8 #rhat*b + u[j+2] into r9,r8");
        !           127: asm("  # give up if r3,r2 <= r9,r8, otherwise iterate");
        !           128: asm("  subl2   r8,r2                   # perform r3,r2 - r9,r8");
        !           129: asm("  sbwc    r9,r3");
        !           130: asm("  bleq    3f                      # give up if negative or equal");
        !           131: asm("  decl    r0                      # otherwise, qhat = qhat - 1");
        !           132: asm("  addl2   (r10),r1                # since dec'ed qhat, inc rhat by v[1]");
        !           133: asm("  jbr     1b");
        !           134: asm("2:        ");
        !           135: asm("  # get here if v[1]==u[j]");
        !           136: asm("  # set qhat to b-1");
        !           137: asm("  # rhat is easily calculated since if we substitute b-1 for qhat in");
        !           138: asm("  # the formula, then it simplifies to (u[j+1] + v[1])");
        !           139: asm("  # ");
        !           140: asm("  addl3   4(r11),(r10),r1         # rhat = u[j+1] + v[1]");
        !           141: asm("  movl    $0x3fffffff,r0          # qhat = b-1");
        !           142: asm("  jbr     1b");
        !           143: asm("3:");
        !           144: }
        !           145: mlsb(utop,ubot,vtop,nqhat)
        !           146: register int *utop, *ubot, *vtop;
        !           147: register int nqhat;
        !           148: {
        !           149: asm("  clrl    r0");
        !           150: asm("loop2:    addl2   (r11),r0");
        !           151: asm("  emul    r8,-(r9),r0,r2");
        !           152: asm("  extzv   $0,$30,r2,(r11)");
        !           153: asm("  extv    $30,$32,r2,r0");
        !           154: asm("  acbl    r10,$-4,r11,loop2");
        !           155: }
        !           156: adback(utop,ubot,vtop)
        !           157: register int *utop, *ubot, *vtop;
        !           158: {
        !           159: asm("  clrl    r0");
        !           160: asm("loop3:    addl2   -(r9),r0");
        !           161: asm("  addl2   (r11),r0");
        !           162: asm("  extzv   $0,$30,r0,(r11)");
        !           163: asm("  extv    $30,$2,r0,r0");
        !           164: asm("  acbl    r10,$-4,r11,loop3");
        !           165: }
        !           166: dsdiv(top,bot,div)
        !           167: register int* bot;
        !           168: {
        !           169: asm("  clrl    r0");
        !           170: asm("loop4:    emul    r0,$0x40000000,(r11),r1");
        !           171: asm("  ediv    12(ap),r1,(r11),r0");
        !           172: asm("  acbl    4(ap),$4,r11,loop4");
        !           173: }
        !           174: dsmult(top,bot,mult)
        !           175: register int* top;
        !           176: {
        !           177: asm("  clrl    r0");
        !           178: asm("loop5:    emul    12(ap),(r11),r0,r1");
        !           179: asm("  extzv   $0,$30,r1,(r11)");
        !           180: asm("  extv    $30,$32,r1,r0");
        !           181: asm("  acbl    8(ap),$-4,r11,loop5");
        !           182: asm("  movl    r1,4(r11)");
        !           183: }
        !           184: */
        !           185: lispval
        !           186: export(top,bot)
        !           187: register long *top, *bot;
        !           188: {
        !           189:        register lispval p;
        !           190:        lispval result;
        !           191: 
        !           192:        top--; /* screwey convention matches original
        !           193:                  vax assembler convenience */
        !           194:        while(bot < top)
        !           195:        {
        !           196:                if(*bot==0)
        !           197:                        bot++;
        !           198:                else if(*bot==-1)
        !           199:                        *++bot |= 0xc0000000;
        !           200:                else break;
        !           201:        }
        !           202:        if(bot==top) return(inewint(*bot));
        !           203:        result = p = newsdot();
        !           204:        protect(p);
        !           205:        p->s.I = *top--;
        !           206:        while(top >= bot) {
        !           207:                p = p->s.CDR = newdot();
        !           208:                p->s.I = *top--;
        !           209:        }
        !           210:        p->s.CDR = 0;
        !           211:        np--;
        !           212:        return(result);
        !           213: }
        !           214: 
        !           215: #define MAXINT 0x80000000L
        !           216: 
        !           217: Ihau(fix)
        !           218: register int fix;
        !           219: {
        !           220:        register count;
        !           221:        if(fix==MAXINT)
        !           222:                return(32);
        !           223:        if(fix < 0)
        !           224:                fix = -fix;
        !           225:        for(count = 0; fix; count++)
        !           226:                fix /= 2;
        !           227:        return(count);
        !           228: }
        !           229: lispval
        !           230: Lhau()
        !           231: {
        !           232:        register count;
        !           233:        register lispval handy;
        !           234:        register dum1,dum2;
        !           235:        lispval Labsval();
        !           236: 
        !           237:        handy = lbot->val;
        !           238: top:
        !           239:        switch(TYPE(handy)) {
        !           240:        case INT:
        !           241:                count = Ihau(handy->i);
        !           242:                break;
        !           243:        case SDOT:
        !           244:                handy = Labsval();
        !           245:                for(count = 0; handy->s.CDR!=((lispval) 0); handy = handy->s.CDR)
        !           246:                        count += 30;
        !           247:                count += Ihau(handy->s.I);
        !           248:                break;
        !           249:        default:
        !           250:                handy = errorh1(Vermisc,"Haulong: bad argument",nil,
        !           251:                               TRUE,997,handy);
        !           252:                goto top;
        !           253:        }
        !           254:        return(inewint(count));
        !           255: }
        !           256: lispval
        !           257: Lhaipar()
        !           258: {
        !           259:        int *sp();
        !           260:        register lispval work;
        !           261:        register n;
        !           262:        register int *top = sp() - 1;
        !           263:        register int *bot;
        !           264:        int mylen;
        !           265: 
        !           266:        /*chkarg(2);*/
        !           267:        work = lbot->val;
        !           268:                                        /* copy data onto stack */
        !           269: on1:
        !           270:        switch(TYPE(work)) {
        !           271:        case INT:
        !           272:                stack(work->i);
        !           273:                break;
        !           274:        case SDOT:
        !           275:                for(; work!=((lispval) 0); work = work->s.CDR)
        !           276:                        stack(work->s.I);
        !           277:                break;
        !           278:        default:
        !           279:                work = errorh1(Vermisc,"Haipart: bad first argument",nil,
        !           280:                                TRUE,996,work);
        !           281:                goto on1;
        !           282:        }
        !           283:        bot = sp();
        !           284:        if(*bot < 0) {
        !           285:                stack(0);
        !           286:                dsmult(top,bot,-1);
        !           287:                bot--;
        !           288:        }
        !           289:        for(; *bot==0 && bot < top; bot++);
        !           290:                                /* recalculate haulong internally */
        !           291:        mylen = (top - bot) * 30 + Ihau(*bot);
        !           292:                                /* get second argument            */
        !           293:        work = lbot[1].val;
        !           294:        while(TYPE(work)!=INT)
        !           295:                work = errorh1(Vermisc,"Haipart: 2nd arg not int",nil,
        !           296:                                TRUE,995,work);
        !           297:        n = work->i;
        !           298:        if(n >= mylen || -n >= mylen)
        !           299:                goto done;
        !           300:        if(n==0) return(inewint(0));
        !           301:        if(n > 0) {
        !           302:                                /* Here we want n most significant bits
        !           303:                                   so chop off mylen - n bits */
        !           304:                stack(0);
        !           305:                n = mylen - n;
        !           306:                for(n; n >= 30; n -= 30)
        !           307:                        top--;
        !           308:                if(top < bot)
        !           309:                        error("Internal error in haipart #1",FALSE);
        !           310:                dsdiv(top,bot,1<<n);
        !           311: 
        !           312:        } else {
        !           313:                                /* here we want abs(n) low order bits */
        !           314:                stack(0);
        !           315:                bot = top + 1;
        !           316:                for(; n <= 0; n += 30)
        !           317:                        bot--;
        !           318:                n = 30 - n;
        !           319:                *bot &= ~ (-1<<n);
        !           320:        }
        !           321: done:
        !           322:        return(export(top + 1,bot));
        !           323: }
        !           324: #define STICKY 1
        !           325: #define TOEVEN 2
        !           326: lispval
        !           327: Ibiglsh(bignum,count,mode)
        !           328: lispval bignum, count;
        !           329: {
        !           330:        int *sp();
        !           331:        register lispval work;
        !           332:        register n;
        !           333:        register int *top = sp() - 1;
        !           334:        register int *bot;
        !           335:        int mylen, guard = 0, sticky = 0, round = 0;
        !           336:        lispval export();
        !           337: 
        !           338:                                /* get second argument            */
        !           339:        work = count;
        !           340:        while(TYPE(work)!=INT)
        !           341:                work = errorh1(Vermisc,"Bignum-shift: 2nd arg not int",nil,
        !           342:                                TRUE,995,work);
        !           343:        n = work->i;
        !           344:        if(n==0) return(bignum);
        !           345:        for(; n >= 30; n -= 30) {/* Here we want to multiply by 2^n
        !           346:                                   so start by copying n/30 zeroes
        !           347:                                   onto stack */
        !           348:                stack(0);
        !           349:        }
        !           350: 
        !           351:        work = bignum;          /* copy data onto stack */
        !           352: on1:
        !           353:        switch(TYPE(work)) {
        !           354:        case INT:
        !           355:                stack(work->i);
        !           356:                break;
        !           357:        case SDOT:
        !           358:                for(; work!=((lispval) 0); work = work->s.CDR)
        !           359:                        stack(work->s.I);
        !           360:                break;
        !           361:        default:
        !           362:                work = errorh1(Vermisc,"Bignum-shift: bad bignum argument",nil,
        !           363:                                TRUE,996,work);
        !           364:                goto on1;
        !           365:        }
        !           366:        bot = sp();
        !           367:        if(n >= 0) {
        !           368:                stack(0);
        !           369:                bot--;
        !           370:                dsmult(top,bot,1<<n);
        !           371:        } else {
        !           372:                        /* Trimming will only work without leading
        !           373:                           zeroes without my having to think
        !           374:                           a lot harder about it, if the inputs
        !           375:                           are canonical */
        !           376:                for(n = -n; n > 30; n -= 30) {
        !           377:                        if(guard) sticky |= 1;
        !           378:                        guard = round;
        !           379:                        if(top > bot) {
        !           380:                                round = *top;
        !           381:                                top --;
        !           382:                        } else  {
        !           383:                                round = *top;
        !           384:                                *top >>= 30;
        !           385:                        }
        !           386:                }
        !           387:                if(n > 0) {
        !           388:                        if(guard) sticky |= 1;
        !           389:                        guard = round;
        !           390:                        round = dsrsh(top,bot,-n,-1<<n);
        !           391:                }
        !           392:                stack(0); /*so that dsadd1 will work;*/
        !           393:                if (mode==STICKY) {
        !           394:                        if(((*top&1)==0) && (round | guard | sticky))
        !           395:                                dsadd1(top,bot);
        !           396:                } else if (mode==TOEVEN) {
        !           397:                        int mask;
        !           398: 
        !           399:                        if(n==0) n = 30;
        !           400:                        mask = (1<<(n-1));
        !           401:                        if(! (round & mask) ) goto chop;
        !           402:                        mask -= 1;
        !           403:                        if(  ((round&mask)==0)
        !           404:                          && guard==0
        !           405:                          && sticky==0
        !           406:                          && (*top&1)==0 ) goto chop;
        !           407:                        dsadd1(top,bot);
        !           408:                }
        !           409:                chop:;
        !           410:        }
        !           411:        work = export(top + 1,bot);
        !           412:        return(work);
        !           413: }
        !           414: 
        !           415: /*From drb  Mon Jul 27 01:25:56 1981
        !           416: To: sklower
        !           417: 
        !           418: The idea is that the answer/2
        !           419: is equal to the exact answer/2 rounded towards - infinity.  The final bit
        !           420: of the answer is the "or" of the true final bit, together with all true
        !           421: bits after the binary point.  In other words, the 1's bit of the answer
        !           422: is almost always 1.  THE FINAL BIT OF THE ANSWER IS 0 IFF n*2^i = THE
        !           423: ANSWER RETURNED EXACTLY, WITH A 0 FINAL BIT.
        !           424: 
        !           425: 
        !           426: To try again, more succintly:  the answer is correct to within 1, and
        !           427: the 1's bit of the answer will be 0 only if the answer is exactly
        !           428: correct. */
        !           429: 
        !           430: lispval
        !           431: Lsbiglsh()
        !           432: {
        !           433:        register struct argent *mylbot = lbot;
        !           434:        chkarg(2,"sticky-bignum-leftshift");
        !           435:        return(Ibiglsh(lbot->val,lbot[1].val,STICKY));
        !           436: }
        !           437: lispval
        !           438: Lbiglsh()
        !           439: {
        !           440:        register struct argent *mylbot = lbot;
        !           441:        chkarg(2,"bignum-leftshift");
        !           442:        return(Ibiglsh(lbot->val,lbot[1].val,TOEVEN));
        !           443: }
        !           444: lispval
        !           445: HackHex() /* this is a one minute function so drb and kls can debug biglsh */
        !           446: /* (HackHex i) returns a string which is the result of printing i in hex */
        !           447: {
        !           448:        register struct argent *mylbot = lbot;
        !           449:        char buf[32];
        !           450:        sprintf(buf,"%lx",lbot->val->i);
        !           451:        return((lispval)inewstr(buf));
        !           452: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.