plan9port

fork of plan9port with libvec, libstr and libsdb
Log | Files | Refs | README | LICENSE

ieee.c (3265B)


      1 #include <u.h>
      2 #include <libc.h>
      3 #include <bio.h>
      4 #include <mach.h>
      5 
      6 /*
      7  * These routines assume that if the number is representable
      8  * in IEEE floating point, it will be representable in the native
      9  * double format.  Naive but workable, probably.
     10  */
     11 int
     12 ieeeftoa64(char *buf, uint n, u32int h, u32int l)
     13 {
     14 	double fr;
     15 	int exp;
     16 
     17 	if (n <= 0)
     18 		return 0;
     19 
     20 
     21 	if(h & (1UL<<31)){
     22 		*buf++ = '-';
     23 		h &= ~(1UL<<31);
     24 	}else
     25 		*buf++ = ' ';
     26 	n--;
     27 	if(l == 0 && h == 0)
     28 		return snprint(buf, n, "0.");
     29 	exp = (h>>20) & ((1L<<11)-1L);
     30 	if(exp == 0)
     31 		return snprint(buf, n, "DeN(%.8lux%.8lux)", h, l);
     32 	if(exp == ((1L<<11)-1L)){
     33 		if(l==0 && (h&((1L<<20)-1L)) == 0)
     34 			return snprint(buf, n, "Inf");
     35 		else
     36 			return snprint(buf, n, "NaN(%.8lux%.8lux)", h&((1L<<20)-1L), l);
     37 	}
     38 	exp -= (1L<<10) - 2L;
     39 	fr = l & ((1L<<16)-1L);
     40 	fr /= 1L<<16;
     41 	fr += (l>>16) & ((1L<<16)-1L);
     42 	fr /= 1L<<16;
     43 	fr += (h & (1L<<20)-1L) | (1L<<20);
     44 	fr /= 1L<<21;
     45 	fr = ldexp(fr, exp);
     46 	return snprint(buf, n, "%.18g", fr);
     47 }
     48 
     49 int
     50 ieeeftoa32(char *buf, uint n, u32int h)
     51 {
     52 	double fr;
     53 	int exp;
     54 
     55 	if (n <= 0)
     56 		return 0;
     57 
     58 	if(h & (1UL<<31)){
     59 		*buf++ = '-';
     60 		h &= ~(1UL<<31);
     61 	}else
     62 		*buf++ = ' ';
     63 	n--;
     64 	if(h == 0)
     65 		return snprint(buf, n, "0.");
     66 	exp = (h>>23) & ((1L<<8)-1L);
     67 	if(exp == 0)
     68 		return snprint(buf, n, "DeN(%.8lux)", h);
     69 	if(exp == ((1L<<8)-1L)){
     70 		if((h&((1L<<23)-1L)) == 0)
     71 			return snprint(buf, n, "Inf");
     72 		else
     73 			return snprint(buf, n, "NaN(%.8lux)", h&((1L<<23)-1L));
     74 	}
     75 	exp -= (1L<<7) - 2L;
     76 	fr = (h & ((1L<<23)-1L)) | (1L<<23);
     77 	fr /= 1L<<24;
     78 	fr = ldexp(fr, exp);
     79 	return snprint(buf, n, "%.9g", fr);
     80 }
     81 
     82 int
     83 beieeeftoa32(char *buf, uint n, void *s)
     84 {
     85 	return ieeeftoa32(buf, n, beswap4(*(u32int*)s));
     86 }
     87 
     88 int
     89 beieeeftoa64(char *buf, uint n, void *s)
     90 {
     91 	return ieeeftoa64(buf, n, beswap4(*(u32int*)s), beswap4(((u32int*)(s))[1]));
     92 }
     93 
     94 int
     95 leieeeftoa32(char *buf, uint n, void *s)
     96 {
     97 	return ieeeftoa32(buf, n, leswap4(*(u32int*)s));
     98 }
     99 
    100 int
    101 leieeeftoa64(char *buf, uint n, void *s)
    102 {
    103 	return ieeeftoa64(buf, n, leswap4(((u32int*)(s))[1]), leswap4(*(u32int*)s));
    104 }
    105 
    106 /* packed in 12 bytes, with s[2]==s[3]==0; mantissa starts at s[4]*/
    107 int
    108 beieeeftoa80(char *buf, uint n, void *s)
    109 {
    110 	uchar *reg = (uchar*)s;
    111 	int i;
    112 	ulong x;
    113 	uchar ieee[8+8];	/* room for slop */
    114 	uchar *p, *q;
    115 
    116 	memset(ieee, 0, sizeof(ieee));
    117 	/* sign */
    118 	if(reg[0] & 0x80)
    119 		ieee[0] |= 0x80;
    120 
    121 	/* exponent */
    122 	x = ((reg[0]&0x7F)<<8) | reg[1];
    123 	if(x == 0)		/* number is ±0 */
    124 		goto done;
    125 	if(x == 0x7FFF){
    126 		if(memcmp(reg+4, ieee+1, 8) == 0){ /* infinity */
    127 			x = 2047;
    128 		}else{				/* NaN */
    129 			x = 2047;
    130 			ieee[7] = 0x1;		/* make sure */
    131 		}
    132 		ieee[0] |= x>>4;
    133 		ieee[1] |= (x&0xF)<<4;
    134 		goto done;
    135 	}
    136 	x -= 0x3FFF;		/* exponent bias */
    137 	x += 1023;
    138 	if(x >= (1<<11) || ((reg[4]&0x80)==0 && x!=0))
    139 		return snprint(buf, n, "not in range");
    140 	ieee[0] |= x>>4;
    141 	ieee[1] |= (x&0xF)<<4;
    142 
    143 	/* mantissa */
    144 	p = reg+4;
    145 	q = ieee+1;
    146 	for(i=0; i<56; i+=8, p++, q++){	/* move one byte */
    147 		x = (p[0]&0x7F) << 1;
    148 		if(p[1] & 0x80)
    149 			x |= 1;
    150 		q[0] |= x>>4;
    151 		q[1] |= (x&0xF)<<4;
    152 	}
    153     done:
    154 	return beieeeftoa64(buf, n, (void*)ieee);
    155 }
    156 
    157 
    158 int
    159 leieeeftoa80(char *buf, uint n, void *s)
    160 {
    161 	int i;
    162 	char *cp;
    163 	char b[12];
    164 
    165 	cp = (char*) s;
    166 	for(i=0; i<12; i++)
    167 		b[11-i] = *cp++;
    168 	return beieeeftoa80(buf, n, b);
    169 }