![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fpu_bcd.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe|
Return to fpu_bcd.c CVS log | Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe |
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.