1 #include "cc.h" 2 3 enum 4 { 5 Mpscale = 29, /* safely smaller than bits in a long */ 6 Mpprec = 36, /* Mpscale*Mpprec sb > largest fp exp */ 7 Mpbase = 1L<<Mpscale, 8 }; 9 10 typedef 11 struct 12 { 13 long a[Mpprec]; 14 char ovf; 15 } Mp; 16 17 int mpatof(char*, double*); 18 int mpatov(char *s, vlong *v); 19 void mpint(Mp*, int); 20 void mppow(Mp*, int, int); 21 void mpmul(Mp*, int); 22 void mpadd(Mp*, Mp*); 23 int mptof(Mp*, double*); 24 25 /* 26 * convert a string, s, to floating in *d 27 * return conversion overflow. 28 * required syntax is [+-]d*[.]d*[e[+-]d*] 29 */ 30 int 31 mpatof(char *s, double *d) 32 { 33 Mp a, b; 34 int dp, c, f, ef, ex, zer; 35 double d1, d2; 36 37 dp = 0; /* digits after decimal point */ 38 f = 0; /* sign */ 39 ex = 0; /* exponent */ 40 zer = 1; /* zero */ 41 memset(&a, 0, sizeof(a)); 42 for(;;) { 43 switch(c = *s++) { 44 default: 45 goto bad; 46 case '-': 47 f = 1; 48 case ' ': 49 case '\t': 50 case '+': 51 continue; 52 case '.': 53 dp = 1; 54 continue; 55 case '1': 56 case '2': 57 case '3': 58 case '4': 59 case '5': 60 case '6': 61 case '7': 62 case '8': 63 case '9': 64 zer = 0; 65 case '0': 66 mpint(&b, c-'0'); 67 mpmul(&a, 10); 68 mpadd(&a, &b); 69 if(dp) 70 dp++; 71 continue; 72 case 'E': 73 case 'e': 74 ex = 0; 75 ef = 0; 76 for(;;) { 77 c = *s++; 78 if(c == '+' || c == ' ' || c == '\t') 79 continue; 80 if(c == '-') { 81 ef = 1; 82 continue; 83 } 84 if(c >= '0' && c <= '9') { 85 ex = ex*10 + (c-'0'); 86 continue; 87 } 88 break; 89 } 90 if(ef) 91 ex = -ex; 92 case 0: 93 break; 94 } 95 break; 96 } 97 if(a.ovf) 98 goto bad; 99 if(zer) { 100 *d = 0; 101 return 0; 102 } 103 if(dp) 104 dp--; 105 dp -= ex; 106 if(dp > 0) { 107 /* 108 * must divide by 10**dp 109 */ 110 if(mptof(&a, &d1)) 111 goto bad; 112 113 /* 114 * trial exponent of 8**dp 115 * 8 (being between 5 and 10) 116 * should pick up all underflows 117 * in the division of 5**dp. 118 */ 119 d2 = frexp(d1, &ex); 120 d2 = ldexp(d2, ex-3*dp); 121 if(d2 == 0) 122 goto bad; 123 124 /* 125 * decompose each 10 into 5*2. 126 * create 5**dp in fixed point 127 * and then play with the exponent 128 * for the remaining 2**dp. 129 * note that 5**dp will overflow 130 * with as few as 134 input digits. 131 */ 132 mpint(&a, 1); 133 mppow(&a, 5, dp); 134 if(mptof(&a, &d2)) 135 goto bad; 136 d1 = frexp(d1/d2, &ex); 137 d1 = ldexp(d1, ex-dp); 138 if(d1 == 0) 139 goto bad; 140 } else { 141 /* 142 * must multiply by 10**|dp| -- 143 * just do it in fixed point. 144 */ 145 mppow(&a, 10, -dp); 146 if(mptof(&a, &d1)) 147 goto bad; 148 } 149 if(f) 150 d1 = -d1; 151 *d = d1; 152 return 0; 153 154 bad: 155 return 1; 156 } 157 158 /* 159 * convert a to floating in *d 160 * return conversion overflow 161 */ 162 int 163 mptof(Mp *a, double *d) 164 { 165 double f, g; 166 long x, *a1; 167 int i; 168 169 if(a->ovf) 170 return 1; 171 a1 = a->a; 172 f = ldexp(*a1++, 0); 173 for(i=Mpscale; i<Mpprec*Mpscale; i+=Mpscale) 174 if(x = *a1++) { 175 g = ldexp(x, i); 176 /* 177 * NOTE: the test (g==0) is plan9 178 * specific. ansi compliant overflow 179 * is signaled by HUGE and errno==ERANGE. 180 * change this for your particular ldexp. 181 */ 182 if(g == 0) 183 return 1; 184 f += g; /* this could bomb! */ 185 } 186 *d = f; 187 return 0; 188 } 189 190 /* 191 * return a += b 192 */ 193 void 194 mpadd(Mp *a, Mp *b) 195 { 196 int i, c; 197 long x, *a1, *b1; 198 199 if(b->ovf) 200 a->ovf = 1; 201 if(a->ovf) 202 return; 203 c = 0; 204 a1 = a->a; 205 b1 = b->a; 206 for(i=0; i<Mpprec; i++) { 207 x = *a1 + *b1++ + c; 208 c = 0; 209 if(x >= Mpbase) { 210 x -= Mpbase; 211 c = 1; 212 } 213 *a1++ = x; 214 } 215 a->ovf = c; 216 } 217 218 /* 219 * return a = c 220 */ 221 void 222 mpint(Mp *a, int c) 223 { 224 225 memset(a, 0, sizeof(*a)); 226 a->a[0] = c; 227 } 228 229 /* 230 * return a *= c 231 */ 232 void 233 mpmul(Mp *a, int c) 234 { 235 Mp p; 236 int b; 237 238 memmove(&p, a, sizeof(p)); 239 if(!(c & 1)) 240 memset(a, 0, sizeof(*a)); 241 c &= ~1; 242 for(b=2; c; b<<=1) { 243 mpadd(&p, &p); 244 if(c & b) { 245 mpadd(a, &p); 246 c &= ~b; 247 } 248 } 249 } 250 251 /* 252 * return a *= b**e 253 */ 254 void 255 mppow(Mp *a, int b, int e) 256 { 257 int b1; 258 259 b1 = b*b; 260 b1 = b1*b1; 261 while(e >= 4) { 262 mpmul(a, b1); 263 e -= 4; 264 if(a->ovf) 265 return; 266 } 267 while(e > 0) { 268 mpmul(a, b); 269 e--; 270 } 271 } 272