source: trunk/minix/lib/ansi/ext_comp.c@ 16

Last change on this file since 16 was 9, checked in by Mattia Monga, 14 years ago

Minix 3.1.2a
File size: 15.5 KB
RevLine 
[9]1/*
2 (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
3 See the copyright notice in the ACK home directory, in the file "Copyright".
4*/
5
6/* $Id: ext_comp.c,v 1.1.1.1 2005/04/21 14:56:05 beng Exp $ */
7
8/* extended precision arithmetic for the strtod() and cvt() routines */
9
10/* This may require some more work when long doubles get bigger than 8
11 bytes. In this case, these routines may become obsolete. ???
12*/
13
14#include "ext_fmt.h"
15#include <float.h>
16#include <errno.h>
17#include <ctype.h>
18
19static int b64_add(struct mantissa *e1, struct mantissa *e2);
20static b64_sft(struct mantissa *e1, int n);
21
22static
23mul_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
24{
25 /* Multiply the extended numbers e1 and e2, and put the
26 result in e3.
27 */
28 register int i,j; /* loop control */
29 unsigned short mp[4];
30 unsigned short mc[4];
31 unsigned short result[8]; /* result */
32
33 register unsigned short *pres;
34
35 /* first save the sign (XOR) */
36 e3->sign = e1->sign ^ e2->sign;
37
38 /* compute new exponent */
39 e3->exp = e1->exp + e2->exp + 1;
40
41 /* check for overflow/underflow ??? */
42
43 /* 128 bit multiply of mantissas */
44
45 /* assign unknown long formats */
46 /* to known unsigned word formats */
47 mp[0] = e1->m1 >> 16;
48 mp[1] = (unsigned short) e1->m1;
49 mp[2] = e1->m2 >> 16;
50 mp[3] = (unsigned short) e1->m2;
51 mc[0] = e2->m1 >> 16;
52 mc[1] = (unsigned short) e2->m1;
53 mc[2] = e2->m2 >> 16;
54 mc[3] = (unsigned short) e2->m2;
55 for (i = 8; i--;) {
56 result[i] = 0;
57 }
58 /*
59 * fill registers with their components
60 */
61 for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
62 unsigned short k = 0;
63 unsigned long mpi = mp[i];
64 for(j=4;j--;) {
65 unsigned long tmp = (unsigned long)pres[j] + k;
66 if (mc[j]) tmp += mpi * mc[j];
67 pres[j] = tmp;
68 k = tmp >> 16;
69 }
70 pres[-1] = k;
71 }
72
73 if (! (result[0] & 0x8000)) {
74 e3->exp--;
75 for (i = 0; i <= 3; i++) {
76 result[i] <<= 1;
77 if (result[i+1]&0x8000) result[i] |= 1;
78 }
79 result[4] <<= 1;
80 }
81 /*
82 * combine the registers to a total
83 */
84 e3->m1 = ((unsigned long)(result[0]) << 16) + result[1];
85 e3->m2 = ((unsigned long)(result[2]) << 16) + result[3];
86 if (result[4] & 0x8000) {
87 if (++e3->m2 == 0) {
88 if (++e3->m1 == 0) {
89 e3->m1 = 0x80000000;
90 e3->exp++;
91 }
92 }
93 }
94}
95
96static
97add_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
98{
99 /* Add two extended numbers e1 and e2, and put the result
100 in e3
101 */
102 struct EXTEND ce2;
103 int diff;
104
105 if ((e2->m1 | e2->m2) == 0L) {
106 *e3 = *e1;
107 return;
108 }
109 if ((e1->m1 | e1->m2) == 0L) {
110 *e3 = *e2;
111 return;
112 }
113 ce2 = *e2;
114 *e3 = *e1;
115 e1 = &ce2;
116
117 /* adjust mantissas to equal power */
118 diff = e3->exp - e1->exp;
119 if (diff < 0) {
120 diff = -diff;
121 e3->exp += diff;
122 b64_sft(&(e3->mantissa), diff);
123 }
124 else if (diff > 0) {
125 e1->exp += diff;
126 b64_sft(&(e1->mantissa), diff);
127 }
128 if (e1->sign != e3->sign) {
129 /* e3 + e1 = e3 - (-e1) */
130 if (e1->m1 > e3->m1 ||
131 (e1->m1 == e3->m1 && e1->m2 > e3->m2)) {
132 /* abs(e1) > abs(e3) */
133 if (e3->m2 > e1->m2) {
134 e1->m1 -= 1; /* carry in */
135 }
136 e1->m1 -= e3->m1;
137 e1->m2 -= e3->m2;
138 *e3 = *e1;
139 }
140 else {
141 if (e1->m2 > e3->m2)
142 e3->m1 -= 1; /* carry in */
143 e3->m1 -= e1->m1;
144 e3->m2 -= e1->m2;
145 }
146 }
147 else {
148 if (b64_add(&e3->mantissa,&e1->mantissa)) {/* addition carry */
149 b64_sft(&e3->mantissa,1);/* shift mantissa one bit RIGHT */
150 e3->m1 |= 0x80000000L; /* set max bit */
151 e3->exp++; /* increase the exponent */
152 }
153 }
154 if ((e3->m2 | e3->m1) != 0L) {
155 /* normalize */
156 if (e3->m1 == 0L) {
157 e3->m1 = e3->m2; e3->m2 = 0L; e3->exp -= 32;
158 }
159 if (!(e3->m1 & 0x80000000)) {
160 unsigned long l = 0x40000000;
161 int cnt = -1;
162
163 while (! (l & e3->m1)) {
164 l >>= 1; cnt--;
165 }
166 e3->exp += cnt;
167 b64_sft(&(e3->mantissa), cnt);
168 }
169 }
170}
171
172static int
173cmp_ext(struct EXTEND *e1, struct EXTEND *e2)
174{
175 struct EXTEND tmp;
176
177 e2->sign = ! e2->sign;
178 add_ext(e1, e2, &tmp);
179 e2->sign = ! e2->sign;
180 if (tmp.m1 == 0 && tmp.m2 == 0) return 0;
181 if (tmp.sign) return -1;
182 return 1;
183}
184
185static
186b64_sft(struct mantissa *e1, int n)
187{
188 if (n > 0) {
189 if (n > 63) {
190 e1->l_32 = 0;
191 e1->h_32 = 0;
192 return;
193 }
194 if (n >= 32) {
195 e1->l_32 = e1->h_32;
196 e1->h_32 = 0;
197 n -= 32;
198 }
199 if (n > 0) {
200 e1->l_32 >>= n;
201 if (e1->h_32 != 0) {
202 e1->l_32 |= (e1->h_32 << (32 - n));
203 e1->h_32 >>= n;
204 }
205 }
206 return;
207 }
208 n = -n;
209 if (n > 0) {
210 if (n > 63) {
211 e1->l_32 = 0;
212 e1->h_32 = 0;
213 return;
214 }
215 if (n >= 32) {
216 e1->h_32 = e1->l_32;
217 e1->l_32 = 0;
218 n -= 32;
219 }
220 if (n > 0) {
221 e1->h_32 <<= n;
222 if (e1->l_32 != 0) {
223 e1->h_32 |= (e1->l_32 >> (32 - n));
224 e1->l_32 <<= n;
225 }
226 }
227 }
228}
229
230static int
231b64_add(struct mantissa *e1, struct mantissa *e2)
232 /*
233 * pointers to 64 bit 'registers'
234 */
235{
236 register int overflow;
237 int carry;
238
239 /* add higher pair of 32 bits */
240 overflow = ((unsigned long) 0xFFFFFFFF - e1->h_32 < e2->h_32);
241 e1->h_32 += e2->h_32;
242
243 /* add lower pair of 32 bits */
244 carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
245 e1->l_32 += e2->l_32;
246 if ((carry) && (++e1->h_32 == 0))
247 return(1); /* had a 64 bit overflow */
248 else
249 return(overflow); /* return status from higher add */
250}
251
252/* The following tables can be computed with the following bc(1)
253 program:
254
255obase=16
256scale=0
257define t(x){
258 auto a, b, c
259 a=2;b=1;c=2^32;n=1
260 while(a<x) {
261 b=a;n+=n;a*=a
262 }
263 n/=2
264 a=b
265 while(b<x) {
266 a=b;b*=c;n+=32
267 }
268 n-=32
269 b=a
270 while(a<x) {
271 b=a;a+=a;n+=1
272 }
273 n-=1
274 x*=16^16
275 b=x%a
276 x/=a
277 if(a<=(2*b)) x+=1
278 obase=10
279 n
280 obase=16
281 return(x)
282}
283for (i=1;i<28;i++) {
284 t(10^i)
285}
2860
287for (i=1;i<20;i++) {
288 t(10^(28*i))
289}
2900
291define r(x){
292 auto a, b, c
293 a=2;b=1;c=2^32;n=1
294 while(a<x) {
295 b=a;n+=n;a*=a
296 }
297 n/=2
298 a=b
299 while(b<x) {
300 a=b;b*=c;n+=32
301 }
302 n-=32
303 b=a
304 while(a<x) {
305 b=a;a+=a;n+=1
306 }
307 a=b
308 a*=16^16
309 b=a%x
310 a/=x
311 if(x<=(2*b)) a+=1
312 obase=10
313 -n
314 obase=16
315 return(a)
316}
317for (i=1;i<28;i++) {
318 r(10^i)
319}
3200
321for (i=1;i<20;i++) {
322 r(10^(28*i))
323}
3240
325
326*/
327static struct EXTEND ten_powers[] = { /* representation of 10 ** i */
328 { 0, 0, 0x80000000, 0 },
329 { 0, 3, 0xA0000000, 0 },
330 { 0, 6, 0xC8000000, 0 },
331 { 0, 9, 0xFA000000, 0 },
332 { 0, 13, 0x9C400000, 0 },
333 { 0, 16, 0xC3500000, 0 },
334 { 0, 19, 0xF4240000, 0 },
335 { 0, 23, 0x98968000, 0 },
336 { 0, 26, 0xBEBC2000, 0 },
337 { 0, 29, 0xEE6B2800, 0 },
338 { 0, 33, 0x9502F900, 0 },
339 { 0, 36, 0xBA43B740, 0 },
340 { 0, 39, 0xE8D4A510, 0 },
341 { 0, 43, 0x9184E72A, 0 },
342 { 0, 46, 0xB5E620F4, 0x80000000 },
343 { 0, 49, 0xE35FA931, 0xA0000000 },
344 { 0, 53, 0x8E1BC9BF, 0x04000000 },
345 { 0, 56, 0xB1A2BC2E, 0xC5000000 },
346 { 0, 59, 0xDE0B6B3A, 0x76400000 },
347 { 0, 63, 0x8AC72304, 0x89E80000 },
348 { 0, 66, 0xAD78EBC5, 0xAC620000 },
349 { 0, 69, 0xD8D726B7, 0x177A8000 },
350 { 0, 73, 0x87867832, 0x6EAC9000 },
351 { 0, 76, 0xA968163F, 0x0A57B400 },
352 { 0, 79, 0xD3C21BCE, 0xCCEDA100 },
353 { 0, 83, 0x84595161, 0x401484A0 },
354 { 0, 86, 0xA56FA5B9, 0x9019A5C8 },
355 { 0, 89, 0xCECB8F27, 0xF4200F3A }
356};
357static struct EXTEND big_ten_powers[] = { /* representation of 10 ** (28*i) */
358 { 0, 0, 0x80000000, 0 },
359 { 0, 93, 0x813F3978, 0xF8940984 },
360 { 0, 186, 0x82818F12, 0x81ED44A0 },
361 { 0, 279, 0x83C7088E, 0x1AAB65DB },
362 { 0, 372, 0x850FADC0, 0x9923329E },
363 { 0, 465, 0x865B8692, 0x5B9BC5C2 },
364 { 0, 558, 0x87AA9AFF, 0x79042287 },
365 { 0, 651, 0x88FCF317, 0xF22241E2 },
366 { 0, 744, 0x8A5296FF, 0xE33CC930 },
367 { 0, 837, 0x8BAB8EEF, 0xB6409C1A },
368 { 0, 930, 0x8D07E334, 0x55637EB3 },
369 { 0, 1023, 0x8E679C2F, 0x5E44FF8F },
370 { 0, 1116, 0x8FCAC257, 0x558EE4E6 },
371 { 0, 1209, 0x91315E37, 0xDB165AA9 },
372 { 0, 1302, 0x929B7871, 0xDE7F22B9 },
373 { 0, 1395, 0x940919BB, 0xD4620B6D },
374 { 0, 1488, 0x957A4AE1, 0xEBF7F3D4 },
375 { 0, 1581, 0x96EF14C6, 0x454AA840 },
376 { 0, 1674, 0x98678061, 0x27ECE4F5 },
377 { 0, 1767, 0x99E396C1, 0x3A3ACFF2 }
378};
379
380static struct EXTEND r_ten_powers[] = { /* representation of 10 ** -i */
381 { 0, 0, 0x80000000, 0 },
382 { 0, -4, 0xCCCCCCCC, 0xCCCCCCCD },
383 { 0, -7, 0xA3D70A3D, 0x70A3D70A },
384 { 0, -10, 0x83126E97, 0x8D4FDF3B },
385 { 0, -14, 0xD1B71758, 0xE219652C },
386 { 0, -17, 0xA7C5AC47, 0x1B478423 },
387 { 0, -20, 0x8637BD05, 0xAF6C69B6 },
388 { 0, -24, 0xD6BF94D5, 0xE57A42BC },
389 { 0, -27, 0xABCC7711, 0x8461CEFD },
390 { 0, -30, 0x89705F41, 0x36B4A597 },
391 { 0, -34, 0xDBE6FECE, 0xBDEDD5BF },
392 { 0, -37, 0xAFEBFF0B, 0xCB24AAFF },
393 { 0, -40, 0x8CBCCC09, 0x6F5088CC },
394 { 0, -44, 0xE12E1342, 0x4BB40E13 },
395 { 0, -47, 0xB424DC35, 0x095CD80F },
396 { 0, -50, 0x901D7CF7, 0x3AB0ACD9 },
397 { 0, -54, 0xE69594BE, 0xC44DE15B },
398 { 0, -57, 0xB877AA32, 0x36A4B449 },
399 { 0, -60, 0x9392EE8E, 0x921D5D07 },
400 { 0, -64, 0xEC1E4A7D, 0xB69561A5 },
401 { 0, -67, 0xBCE50864, 0x92111AEB },
402 { 0, -70, 0x971DA050, 0x74DA7BEF },
403 { 0, -74, 0xF1C90080, 0xBAF72CB1 },
404 { 0, -77, 0xC16D9A00, 0x95928A27 },
405 { 0, -80, 0x9ABE14CD, 0x44753B53 },
406 { 0, -84, 0xF79687AE, 0xD3EEC551 },
407 { 0, -87, 0xC6120625, 0x76589DDB },
408 { 0, -90, 0x9E74D1B7, 0x91E07E48 }
409};
410
411static struct EXTEND r_big_ten_powers[] = { /* representation of 10 ** -(28*i) */
412 { 0, 0, 0x80000000, 0 },
413 { 0, -94, 0xFD87B5F2, 0x8300CA0E },
414 { 0, -187, 0xFB158592, 0xBE068D2F },
415 { 0, -280, 0xF8A95FCF, 0x88747D94 },
416 { 0, -373, 0xF64335BC, 0xF065D37D },
417 { 0, -466, 0xF3E2F893, 0xDEC3F126 },
418 { 0, -559, 0xF18899B1, 0xBC3F8CA2 },
419 { 0, -652, 0xEF340A98, 0x172AACE5 },
420 { 0, -745, 0xECE53CEC, 0x4A314EBE },
421 { 0, -838, 0xEA9C2277, 0x23EE8BCB },
422 { 0, -931, 0xE858AD24, 0x8F5C22CA },
423 { 0, -1024, 0xE61ACF03, 0x3D1A45DF },
424 { 0, -1117, 0xE3E27A44, 0x4D8D98B8 },
425 { 0, -1210, 0xE1AFA13A, 0xFBD14D6E },
426 { 0, -1303, 0xDF82365C, 0x497B5454 },
427 { 0, -1396, 0xDD5A2C3E, 0xAB3097CC },
428 { 0, -1489, 0xDB377599, 0xB6074245 },
429 { 0, -1582, 0xD91A0545, 0xCDB51186 },
430 { 0, -1675, 0xD701CE3B, 0xD387BF48 },
431 { 0, -1768, 0xD4EEC394, 0xD6258BF8 }
432};
433
434#define TP (int)(sizeof(ten_powers)/sizeof(ten_powers[0]))
435#define BTP (int)(sizeof(big_ten_powers)/sizeof(big_ten_powers[0]))
436#define MAX_EXP (TP * BTP - 1)
437
438static
439add_exponent(struct EXTEND *e, int exp)
440{
441 int neg = exp < 0;
442 int divsz, modsz;
443 struct EXTEND x;
444
445 if (neg) exp = -exp;
446 divsz = exp / TP;
447 modsz = exp % TP;
448 if (neg) {
449 mul_ext(e, &r_ten_powers[modsz], &x);
450 mul_ext(&x, &r_big_ten_powers[divsz], e);
451 }
452 else {
453 mul_ext(e, &ten_powers[modsz], &x);
454 mul_ext(&x, &big_ten_powers[divsz], e);
455 }
456}
457
458_str_ext_cvt(const char *s, char **ss, struct EXTEND *e)
459{
460 /* Like strtod, but for extended precision */
461 register int c;
462 int dotseen = 0;
463 int digitseen = 0;
464 int exp = 0;
465
466 if (ss) *ss = (char *)s;
467 while (isspace(*s)) s++;
468
469 e->sign = 0;
470 e->exp = 0;
471 e->m1 = e->m2 = 0;
472
473 c = *s;
474 switch(c) {
475 case '-':
476 e->sign = 1;
477 case '+':
478 s++;
479 }
480 while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) {
481 if (c == '.') continue;
482 digitseen = 1;
483 if (e->m1 <= (unsigned long)(0xFFFFFFFF)/10) {
484 struct mantissa a1;
485
486 a1 = e->mantissa;
487 b64_sft(&(e->mantissa), -3);
488 b64_sft(&a1, -1);
489 b64_add(&(e->mantissa), &a1);
490 a1.h_32 = 0;
491 a1.l_32 = c - '0';
492 b64_add(&(e->mantissa), &a1);
493 }
494 else exp++;
495 if (dotseen) exp--;
496 }
497 if (! digitseen) return;
498
499 if (ss) *ss = (char *)s - 1;
500
501 if (c == 'E' || c == 'e') {
502 int exp1 = 0;
503 int sign = 1;
504 int exp_overflow = 0;
505
506 switch(*s) {
507 case '-':
508 sign = -1;
509 case '+':
510 s++;
511 }
512 if (c = *s, isdigit(c)) {
513 do {
514 int tmp;
515
516 exp1 = 10 * exp1 + (c - '0');
517 if ((tmp = sign * exp1 + exp) > MAX_EXP ||
518 tmp < -MAX_EXP) {
519 exp_overflow = 1;
520 }
521 } while (c = *++s, isdigit(c));
522 if (ss) *ss = (char *)s;
523 }
524 exp += sign * exp1;
525 if (exp_overflow) {
526 exp = sign * MAX_EXP;
527 if (e->m1 != 0 || e->m2 != 0) errno = ERANGE;
528 }
529 }
530 if (e->m1 == 0 && e->m2 == 0) return;
531 e->exp = 63;
532 while (! (e->m1 & 0x80000000)) {
533 b64_sft(&(e->mantissa),-1);
534 e->exp--;
535 }
536 add_exponent(e, exp);
537}
538
539#include <math.h>
540
541static
542ten_mult(struct EXTEND *e)
543{
544 struct EXTEND e1 = *e;
545
546 e1.exp++;
547 e->exp += 3;
548 add_ext(e, &e1, e);
549}
550
551#define NDIGITS 128
552#define NSIGNIFICANT 19
553
554char *
555_ext_str_cvt(struct EXTEND *e, int ndigit, int *decpt, int *sign, int ecvtflag)
556{
557 /* Like cvt(), but for extended precision */
558
559 static char buf[NDIGITS+1];
560 struct EXTEND m;
561 register char *p = buf;
562 register char *pe;
563 int findex = 0;
564
565 if (ndigit < 0) ndigit = 0;
566 if (ndigit > NDIGITS) ndigit = NDIGITS;
567 pe = &buf[ndigit];
568 buf[0] = '\0';
569
570 *sign = 0;
571 if (e->sign) {
572 *sign = 1;
573 e->sign = 0;
574 }
575
576 *decpt = 0;
577 if (e->m1 != 0) {
578 register struct EXTEND *pp = &big_ten_powers[1];
579
580 while(cmp_ext(e,pp) >= 0) {
581 pp++;
582 findex = pp - big_ten_powers;
583 if (findex >= BTP) break;
584 }
585 pp--;
586 findex = pp - big_ten_powers;
587 mul_ext(e,&r_big_ten_powers[findex],e);
588 *decpt += findex * TP;
589 pp = &ten_powers[1];
590 while(pp < &ten_powers[TP] && cmp_ext(e, pp) >= 0) pp++;
591 pp--;
592 findex = pp - ten_powers;
593 *decpt += findex;
594
595 if (cmp_ext(e, &ten_powers[0]) < 0) {
596 pp = &r_big_ten_powers[1];
597 while(cmp_ext(e,pp) < 0) pp++;
598 pp--;
599 findex = pp - r_big_ten_powers;
600 mul_ext(e, &big_ten_powers[findex], e);
601 *decpt -= findex * TP;
602 /* here, value >= 10 ** -28 */
603 ten_mult(e);
604 (*decpt)--;
605 pp = &r_ten_powers[0];
606 while(cmp_ext(e, pp) < 0) pp++;
607 findex = pp - r_ten_powers;
608 mul_ext(e, &ten_powers[findex], e);
609 *decpt -= findex;
610 findex = 0;
611 }
612 (*decpt)++; /* because now value in [1.0, 10.0) */
613 }
614 if (! ecvtflag) {
615 /* for fcvt() we need ndigit digits behind the dot */
616 pe += *decpt;
617 if (pe > &buf[NDIGITS]) pe = &buf[NDIGITS];
618 }
619 m.exp = -62;
620 m.sign = 0;
621 m.m1 = 0xA0000000;
622 m.m2 = 0;
623 while (p <= pe) {
624 struct EXTEND oneminm;
625
626 if (p - pe > NSIGNIFICANT) {
627 findex = 0;
628 e->m1 = 0;
629 }
630 if (findex) {
631 struct EXTEND tc, oldtc;
632 int count = 0;
633
634 oldtc.exp = 0;
635 oldtc.sign = 0;
636 oldtc.m1 = 0;
637 oldtc.m2 = 0;
638 tc = ten_powers[findex];
639 while (cmp_ext(e, &tc) >= 0) {
640 oldtc = tc;
641 add_ext(&tc, &ten_powers[findex], &tc);
642 count++;
643 }
644 *p++ = count + '0';
645 oldtc.sign = 1;
646 add_ext(e, &oldtc, e);
647 findex--;
648 continue;
649 }
650 if (e->m1) {
651 m.sign = 1;
652 add_ext(&ten_powers[0], &m, &oneminm);
653 m.sign = 0;
654 if (e->exp >= 0) {
655 struct EXTEND x;
656
657 x.m2 = 0; x.exp = e->exp;
658 x.sign = 1;
659 x.m1 = e->m1>>(31-e->exp);
660 *p++ = (x.m1) + '0';
661 x.m1 = x.m1 << (31-e->exp);
662 add_ext(e, &x, e);
663 }
664 else *p++ = '0';
665 /* Check that remainder is still significant */
666 if (cmp_ext(&m, e) > 0 || cmp_ext(e, &oneminm) > 0) {
667 if (e->m1 && e->exp >= -1) *(p-1) += 1;
668 e->m1 = 0;
669 continue;
670 }
671 ten_mult(&m);
672 ten_mult(e);
673 }
674 else *p++ = '0';
675 }
676 if (pe >= buf) {
677 p = pe;
678 *p += 5; /* round of at the end */
679 while (*p > '9') {
680 *p = '0';
681 if (p > buf) ++*--p;
682 else {
683 *p = '1';
684 ++*decpt;
685 if (! ecvtflag) {
686 /* maybe add another digit at the end,
687 because the point was shifted right
688 */
689 if (pe > buf) *pe = '0';
690 pe++;
691 }
692 }
693 }
694 *pe = '\0';
695 }
696 return buf;
697}
698
699_dbl_ext_cvt(double value, struct EXTEND *e)
700{
701 /* Convert double to extended
702 */
703 int exponent;
704
705 value = frexp(value, &exponent);
706 e->sign = value < 0.0;
707 if (e->sign) value = -value;
708 e->exp = exponent - 1;
709 value *= 4294967296.0;
710 e->m1 = value;
711 value -= e->m1;
712 value *= 4294967296.0;
713 e->m2 = value;
714}
715
716static struct EXTEND max_d;
717
718double
719_ext_dbl_cvt(struct EXTEND *e)
720{
721 /* Convert extended to double
722 */
723 double f;
724 int sign = e->sign;
725
726 e->sign = 0;
727 if (e->m1 == 0 && e->m2 == 0) {
728 return 0.0;
729 }
730 if (max_d.exp == 0) {
731 _dbl_ext_cvt(DBL_MAX, &max_d);
732 }
733 if (cmp_ext(&max_d, e) < 0) {
734 f = HUGE_VAL;
735 errno = ERANGE;
736 }
737 else f = ldexp((double)e->m1*4294967296.0 + (double)e->m2, e->exp-63);
738 if (sign) f = -f;
739 if (f == 0.0 && (e->m1 != 0 || e->m2 != 0)) {
740 errno = ERANGE;
741 }
742 return f;
743}
Note: See TracBrowser for help on using the repository browser.