Return to fpu_bcd.c CVS log

Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe

Annotation of nono/fpe/fpu_bcd.c, revision 1.1

1.1     ! root        1: //
        !             2: // XM6i
        !             3: // Copyright (c) 2013 Y.Sugahara
        !             4: //
        !             5: // [ fpn → BCD 変換 ]
        !             6: //
        !             7: 
        !             8: #include "fpu_emulate.h"
        !             9: #include "fpu_bcd.h"
        !            10: #include <string.h>
        !            11: #include <stdint.h>
        !            12: #include <stdio.h>
        !            13: #include <math.h>
        !            14: #include <stdlib.h>
        !            15: 
        !            16: //#define FPU_BCD_DEBUG
        !            17: 
        !            18: #if defined(FPU_BCD_DEBUG)
        !            19: #define DPRINTF(msg, ...)      printf(msg, ## __VA_ARGS__)
        !            20: #define dump(s, a)     debug_dump(s, a)
        !            21: #else
        !            22: #define DPRINTF(msg, ...)      /**/
        !            23: #define dump(s, a)     /**/
        !            24: #endif
        !            25: 
        !            26: static int fpu_bcd_iszero(const BCD *a);
        !            27: static uint8_t fpu_bcd_shr(BCD *a);
        !            28: static uint8_t fpu_bcd_shr_n(BCD *a, int n);
        !            29:        uint8_t fpu_bcd_shl(BCD *a);
        !            30: static uint8_t fpu_bcd_shl_n(BCD *a, int n);
        !            31: static uint8_t fpu_bcd_add_d(BCD *a, const BCD *b);
        !            32: static void fpu_bcd_add_last_n(BCD *a, uint8_t n);
        !            33: static void fpu_bcd_round_cy(BCD *a, uint8_t cy);
        !            34: static void fpu_bcd_round(BCD *a);
        !            35: static void fpu_bcd_add(BCD *a, const BCD *b);
        !            36: static void fpu_bcd_mul(BCD *a, const BCD *b);
        !            37: static void fpu_bcd_exp2(BCD *a, int n);
        !            38: static int getbit(struct fpn *fp, int n);
        !            39: static void fpu_itobcd(BCD *a, int x);
        !            40: static void fpu_fpntobcd(BCD *res, struct fpn *fp);
        !            41: 
        !            42: #if defined(FPU_BCD_DEBUG)
        !            43: static void __attribute__((__unused__))
        !            44: debug_dump(const char *s, const BCD *a)
        !            45: {
        !            46:        int i;
        !            47: 
        !            48:        printf("%10s ", s);
        !            49:        printf("e=%+5d ", a->e);
        !            50:        printf("%d.", a->d[0]);
        !            51:        for (i = 1; i < FPU_BCD_DIGIT; i++) {
        !            52:                printf("%d", a->d[i]);
        !            53:        }
        !            54:        printf("\n");
        !            55: }
        !            56: #endif
        !            57: 
        !            58: /*
        !            59:  * (a == 0)
        !            60:  */
        !            61: static int
        !            62: fpu_bcd_iszero(const BCD *a)
        !            63: {
        !            64:        int i;
        !            65:        for (i = 0; i < FPU_BCD_DIGIT; i++) {
        !            66:                if (a->d[i]) return 0;
        !            67:        }
        !            68:        return 1;
        !            69: }
        !            70: 
        !            71: /*
        !            72:  * a >>= 1
        !            73:  * return: carry digit
        !            74:  */
        !            75: static uint8_t
        !            76: fpu_bcd_shr(BCD *a)
        !            77: {
        !            78:        int i;
        !            79:        uint8_t cy;
        !            80: 
        !            81:        cy = a->d[FPU_BCD_DIGIT - 1];
        !            82:        for (i = FPU_BCD_DIGIT - 1; i >= 1; i--) {
        !            83:                a->d[i] = a->d[i - 1];
        !            84:        }
        !            85:        a->d[0] = 0;
        !            86: dump("shr", a);
        !            87:        return cy;
        !            88: }
        !            89: 
        !            90: /*
        !            91:  * a >>= n
        !            92:  */
        !            93: static uint8_t
        !            94: fpu_bcd_shr_n(BCD *a, int n)
        !            95: {
        !            96:        uint8_t cy;
        !            97:        int count;
        !            98: 
        !            99:        if (n <= 0) return 0;
        !           100:        count = FPU_BCD_DIGIT - n;
        !           101:        if (count < 0) {
        !           102:                memset(a->d, 0, FPU_BCD_DIGIT * sizeof(a->d[0]));
        !           103:                return 0;
        !           104:        }
        !           105:        cy = a->d[count];
        !           106: 
        !           107:        memmove(&a->d[n], &a->d[0], count * sizeof(a->d[0]));
        !           108:        memset(&a->d[0], 0, n * sizeof(a->d[0]));
        !           109:        return cy;
        !           110: }
        !           111: 
        !           112: /*
        !           113:  * a <<= 1
        !           114:  * return: carry digit
        !           115:  */
        !           116: /* 使われてないけど対称性のため残してある。どうすべ */
        !           117: uint8_t
        !           118: fpu_bcd_shl(BCD *a)
        !           119: {
        !           120:        int i;
        !           121:        uint8_t cy;
        !           122: 
        !           123:        cy = a->d[0];
        !           124:        for (i = 0; i < FPU_BCD_DIGIT - 1; i++) {
        !           125:                a->d[i] = a->d[i + 1];
        !           126:        }
        !           127:        a->d[FPU_BCD_DIGIT - 1] = 0;
        !           128: dump("shl", a);
        !           129:        return cy;
        !           130: }
        !           131: 
        !           132: /*
        !           133:  * a <<= n
        !           134:  * return: carry digit
        !           135:  */
        !           136: static uint8_t
        !           137: fpu_bcd_shl_n(BCD *a, int n)
        !           138: {
        !           139:        uint8_t cy;
        !           140: 
        !           141:        if (n <= 0)
        !           142:                return 0;
        !           143:        if (n >= FPU_BCD_DIGIT) {
        !           144:                memset(a->d, 0, FPU_BCD_DIGIT * sizeof(a->d[0]));
        !           145:                return 0;
        !           146:        }
        !           147:        cy = a->d[n - 1];
        !           148: 
        !           149:        memmove(&a->d[0], &a->d[n], (FPU_BCD_DIGIT - n) * sizeof(a->d[0]));
        !           150:        memset(&a->d[FPU_BCD_DIGIT - n], 0, n * sizeof(a->d[0]));
        !           151:        return cy;
        !           152: }
        !           153: 
        !           154: /*
        !           155:  * a = a + b
        !           156:  * add digit only
        !           157:  */
        !           158: static uint8_t
        !           159: fpu_bcd_add_d(BCD *a, const BCD *b)
        !           160: {
        !           161:        uint8_t cy = 0;
        !           162:        int i;
        !           163: 
        !           164:        for (i = FPU_BCD_DIGIT - 1; i >= 0; i--) {
        !           165:                a->d[i] = a->d[i] + b->d[i] + cy;
        !           166:                if (a->d[i] >= 10) {
        !           167:                        a->d[i] -= 10;
        !           168:                        cy = 1;
        !           169:                } else {
        !           170:                        cy = 0;
        !           171:                }
        !           172:        }
        !           173:        return cy;
        !           174: }
        !           175: 
        !           176: /*
        !           177:  * a = a[last digit] + n
        !           178:  * for rounding
        !           179:  */
        !           180: static void
        !           181: fpu_bcd_add_last_n(BCD *a, uint8_t n)
        !           182: {
        !           183:        uint8_t cy = n;
        !           184:        int i;
        !           185: 
        !           186:        for (i = FPU_BCD_DIGIT - 1; i >= 0; i--) {
        !           187:                a->d[i] += cy;
        !           188:                if (a->d[i] >= 10) {
        !           189:                        a->d[i] -= 10;
        !           190:                        cy = 1;
        !           191:                } else {
        !           192:                        cy = 0;
        !           193:                        /* breakable because no more carry */
        !           194:                        break;
        !           195:                }
        !           196:        }
        !           197: 
        !           198:        if (cy) {
        !           199:                cy = fpu_bcd_shr(a);
        !           200:                a->d[0] = 1;
        !           201:                a->e++;
        !           202:                fpu_bcd_round_cy(a, cy);
        !           203:        }
        !           204: dump("add_last_n", a);
        !           205: }
        !           206: 
        !           207: /*
        !           208:  * banker's rounding by carry
        !           209:  */
        !           210: static void
        !           211: fpu_bcd_round_cy(BCD *a, uint8_t cy)
        !           212: {
        !           213:        if (cy > 5
        !           214:         || (cy == 5 && a->d[FPU_BCD_DIGIT - 1] % 2 != 0)) {
        !           215:                fpu_bcd_add_last_n(a, 1);
        !           216:        }
        !           217: }
        !           218: 
        !           219: /*
        !           220:  * rounding at last digit
        !           221:  */
        !           222: static void
        !           223: fpu_bcd_round(BCD *a)
        !           224: {
        !           225:        if (a->d[FPU_BCD_DIGIT - 1] >= 5) {
        !           226:                fpu_bcd_add_last_n(a, 10);
        !           227:                a->d[FPU_BCD_DIGIT - 1] = 0;
        !           228:        }
        !           229: }
        !           230: 
        !           231: /*
        !           232:  * a = a + b
        !           233:  */
        !           234: static void
        !           235: fpu_bcd_add(BCD *a, const BCD *b)
        !           236: {
        !           237:        uint8_t cy = 0;
        !           238:        BCD t0;
        !           239:        BCD *t = &t0;
        !           240: 
        !           241:        if (fpu_bcd_iszero(a)) {
        !           242:                *a = *b;
        !           243:                return;
        !           244:        }
        !           245:        if (fpu_bcd_iszero(b)) {
        !           246:                return;
        !           247:        }
        !           248: 
        !           249:        *t = *b;
        !           250: 
        !           251:        while (a->e != t->e) {
        !           252:                if (a->e > t->e) {
        !           253:                        t->e++;
        !           254:                        cy = fpu_bcd_shr(t);
        !           255:                } else if (a->e < t->e) {
        !           256:                        a->e++;
        !           257:                        cy = fpu_bcd_shr(a);
        !           258:                }
        !           259:        }
        !           260:        fpu_bcd_round_cy(a, cy);
        !           261: 
        !           262: dump("add:a", a);
        !           263: dump("add:t", t);
        !           264:        cy = fpu_bcd_add_d(a, t);
        !           265: 
        !           266:        if (cy) {
        !           267:                // last digit rounding
        !           268:                cy = fpu_bcd_shr(a);
        !           269:                a->d[0] = 1;
        !           270:                a->e++;
        !           271:                fpu_bcd_round_cy(a, cy);
        !           272:        }
        !           273: dump("add", a);
        !           274: }
        !           275: 
        !           276: /*
        !           277:  * a = a * b
        !           278:  */
        !           279: static void
        !           280: fpu_bcd_mul(BCD *a, const BCD *b)
        !           281: {
        !           282:        /* need DIGIT * 2 + 1 for mul + cy */
        !           283:        uint8_t t0[FPU_BCD_DIGIT * 2 + 1 + 1];
        !           284:        uint8_t *t = &t0[1];
        !           285:        int i, j, k;
        !           286: 
        !           287: dump("mul:a", a);
        !           288: dump("mul:b", b);
        !           289:        memset(t0, 0, sizeof(t0));
        !           290: 
        !           291:        for (i = FPU_BCD_DIGIT - 1; i >= 0; i--) {
        !           292:                for (j = FPU_BCD_DIGIT - 1; j >= 0; j--) {
        !           293:                        t[i + j] += a->d[i] * b->d[j];
        !           294:                        for (k = i + j; k >= 0; k--) {
        !           295:                                if (t[k] >= 10) {
        !           296:                                        /* t[-1] available */
        !           297:                                        t[k - 1] += t[k] / 10;
        !           298:                                        t[k] %= 10;
        !           299:                                } else {
        !           300:                                        break;
        !           301:                                }
        !           302:                        }
        !           303:                }
        !           304:        }
        !           305: 
        !           306:        a->e = a->e + b->e;
        !           307:        if (t[-1]) {
        !           308:                a->e++;
        !           309:                /* copy 1 digit shifted */
        !           310:                memcpy(a->d, &t[-1], sizeof(a->d));
        !           311:                fpu_bcd_round_cy(a, t[FPU_BCD_DIGIT - 1]);
        !           312:        } else {
        !           313:                memcpy(a->d, t, sizeof(a->d));
        !           314:                fpu_bcd_round_cy(a, t[FPU_BCD_DIGIT]);
        !           315:        }
        !           316: dump("mul", a);
        !           317: }
        !           318: 
        !           319: /*
        !           320:  * a = pow(2, n)
        !           321:  */
        !           322: static void
        !           323: fpu_bcd_exp2(BCD *a, int n)
        !           324: {
        !           325:        BCD r0;
        !           326:        BCD *r = &r0;
        !           327:        int i;
        !           328: #if defined(FPU_BCD_DEBUG)
        !           329: int save_n = n;
        !           330: #endif
        !           331: DPRINTF("exp2 n=%d\n", n);
        !           332:        memset(a, 0, sizeof(*a));
        !           333:        a->d[0] = 1;
        !           334: 
        !           335:        if (n == 0) {
        !           336:        } else {
        !           337:                memset(r, 0, sizeof(*r));
        !           338:                if (n > 0) {
        !           339:                        /* r = 2 */
        !           340:                        r->d[0] = 2;
        !           341:                } else {
        !           342:                        n = -n;
        !           343:                        /* r = 0.5 */
        !           344:                        r->e = -1;
        !           345:                        r->d[0] = 5;
        !           346:                }
        !           347: 
        !           348:                for (i = 0; i < 32; i++) {
        !           349:                        if (n & (1 << i)) {
        !           350:                                n ^= (1 << i);
        !           351:                                fpu_bcd_mul(a, r);
        !           352:                        }
        !           353:                        if (n == 0) break;
        !           354:                        fpu_bcd_mul(r, r);
        !           355:                }
        !           356:        }
        !           357: dump("exp2", a);
        !           358: #if defined(__NetBSD__)        /* とりあえず */
        !           359: DPRINTF("%18s %.*e\n", "pow", FPU_BCD_DIGIT-1, pow(2, save_n));
        !           360: #else
        !           361: DPRINTF("%18s %.*Le\n", "powl", FPU_BCD_DIGIT-1, powl(2, save_n));
        !           362: #endif
        !           363: }
        !           364: 
        !           365: /* relative mantissa bit */
        !           366: int
        !           367: getbit(struct fpn *fp, int n)
        !           368: {
        !           369:        int t;
        !           370:        uint32_t x;
        !           371:        int b;
        !           372: 
        !           373:        t = n + (32 - FP_LG);
        !           374:        x = fp->fp_mant[(t - 1) / 32];
        !           375:        b = x & (1 << (31 - ((t - 1) % 32)));
        !           376: DPRINTF("n=%d %d\n", n, b);
        !           377:        return b;
        !           378: }
        !           379: 
        !           380: static void
        !           381: fpu_itobcd(BCD *a, int x)
        !           382: {
        !           383:        int i;
        !           384:        BCD r0;
        !           385:        BCD *r = &r0;
        !           386: 
        !           387:        memset(a, 0, sizeof(*a));
        !           388: 
        !           389:        x = abs(x);
        !           390: 
        !           391:        for (i = 0; i < 32; i++) {
        !           392:                if (x == 0) break;
        !           393:                if (x & (1 << i)) {
        !           394:                        x ^= (1 << i);
        !           395:                        fpu_bcd_exp2(r, i);
        !           396:                        fpu_bcd_add(a, r);
        !           397:                }
        !           398:        }
        !           399: }
        !           400: 
        !           401: /* Convert fpn to internal BCD structure */
        !           402: static void
        !           403: fpu_fpntobcd(BCD *res, struct fpn *fp)
        !           404: {
        !           405:        BCD bcd_r0;
        !           406:        BCD *bcd_r = &bcd_r0;
        !           407:        BCD bcd_a0;
        !           408:        BCD *bcd_a = &bcd_a0;
        !           409:        BCD bcd_b0;
        !           410:        BCD *bcd_b = &bcd_b0;
        !           411:        int i;
        !           412: 
        !           413:        memset(bcd_r, 0, sizeof(*bcd_r));
        !           414:        memset(bcd_a, 0, sizeof(*bcd_a));
        !           415:        memset(bcd_b, 0, sizeof(*bcd_b));
        !           416: 
        !           417:        for (i = FP_NMANT - 1; i >= 0; i--) {
        !           418:                // 最下位ビットから取り出していく
        !           419:                if (getbit(fp, i)) {
        !           420: 
        !           421:                        // 調べるビットの相当する 2^(exp-i) を計算する
        !           422:                        fpu_bcd_exp2(bcd_r, fp->fp_exp - i);
        !           423: dump("r", bcd_r);
        !           424: 
        !           425:                        fpu_bcd_add(bcd_b, bcd_r);
        !           426: dump("b", bcd_b);
        !           427:                }
        !           428:        }
        !           429: dump("b", bcd_b);
        !           430: 
        !           431:        fpu_bcd_round(bcd_b);
        !           432: dump("b", bcd_b);
        !           433: 
        !           434:        /* とりあえず */
        !           435:        *res = *bcd_b;
        !           436:        return;
        !           437: }
        !           438: 
        !           439: /*
        !           440:  * fpn -> 96-bit packed BCD(space)
        !           441:  */
        !           442: void
        !           443: fpu_ftop(struct fpemu *fe, struct fpn *fp, uint32_t *space, int kfactor)
        !           444: {
        !           445:        BCD b0, *b;
        !           446:        BCD e0, *e;
        !           447:        int i;
        !           448:        int shift;
        !           449: 
        !           450:        b = &b0;
        !           451:        e = &e0;
        !           452: 
        !           453:        space[0] = 0;
        !           454:        space[1] = 0;
        !           455:        space[2] = 0;
        !           456: 
        !           457:        if (fp->fp_sign) {
        !           458:                space[0] |= 0x80000000;
        !           459:        }
        !           460: 
        !           461:        if (ISINF(fp)) {
        !           462:                space[0] |= 0x7fff0000;
        !           463:                return;
        !           464:        }
        !           465: 
        !           466:        /* fp をアンパックド BCD に変換 */
        !           467:        fpu_fpntobcd(b, fp);
        !           468: 
        !           469:        /*
        !           470:         * k-factor
        !           471:         */
        !           472: DPRINTF("kfactor=%d\n", kfactor);
        !           473:        if (kfactor <= 0) {
        !           474:                /*
        !           475:                 * FORTRAN "F" format.
        !           476:                 * 小数以下 -k 桁まで出力する。
        !           477:                 */
        !           478:                if (b->e >= 17) {
        !           479:                        /* 10^17以上の値なら17桁で小数に到達しないので k=+17 と等価 */
        !           480:                        kfactor = 17;
        !           481:                } else {
        !           482:                        shift = FPU_BCD_DIGIT - b->e + kfactor - 2;
        !           483:                }
        !           484:        }
        !           485:        if (kfactor > 0) {
        !           486:                if (kfactor > 17) {
        !           487:                        fe->fe_fpsr |= FPSR_OPERR;
        !           488:                        kfactor = 17;
        !           489:                }
        !           490:                /*
        !           491:                 * FORTRAN "E" format.
        !           492:                 * k は小数点位置に関わらず出力する桁数(文字数) を表す。
        !           493:                 * 1 なら 1桁だけ(つまり整数部のみ)、17 で全桁となる。
        !           494:                 */
        !           495:                shift = FPU_BCD_DIGIT - kfactor - 1;
        !           496:        }
        !           497: dump("k:0", b);
        !           498:        fpu_bcd_shr_n(b, shift);
        !           499: dump("k:1", b);
        !           500:        fpu_bcd_round(b);
        !           501:        /* 四捨五入に使った最終桁はもう不要 */
        !           502:        b->d[FPU_BCD_DIGIT - 1] = 0;
        !           503: dump("k:2", b);
        !           504:        fpu_bcd_shl_n(b, shift);
        !           505: dump("k:3", b);
        !           506: 
        !           507:        /*
        !           508:         * ここから 68881 の パックド BCD に変換していく
        !           509:         */
        !           510:        /* exp の変換 */
        !           511:        /* P 形式の EXP フィールドは 10 基数 */
        !           512:        if (b->e < 0) {
        !           513:                space[0] |= 0x40000000;
        !           514:        }
        !           515: 
        !           516:        fpu_itobcd(e, abs(b->e));
        !           517:        fpu_bcd_shr_n(e, 3 - e->e);
        !           518:        space[0] |= e->d[0] << 12;      /* EXP3 */
        !           519:        space[0] |= e->d[1] << 24;      /* EXP2 */
        !           520:        space[0] |= e->d[2] << 20;      /* EXP1 */
        !           521:        space[0] |= e->d[3] << 16;      /* EXP0 */
        !           522: 
        !           523:        space[0] |= b->d[0];
        !           524: DPRINTF("BCD[0]=%08x\n", space[0]);
        !           525: 
        !           526:        /* 小数部 */
        !           527:        for (i = 1; i < 9; i++) {
        !           528:                space[1] <<= 4;
        !           529:                space[1] |= b->d[i];
        !           530:        }
        !           531: 
        !           532:        for (i = 9; i < 17; i++) {
        !           533:                space[2] <<= 4;
        !           534:                space[2] |= b->d[i];
        !           535:        }
        !           536: }

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