lib/ansi/ext_comp.c

   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$ */
   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
  19 static int b64_add(struct mantissa *e1, struct mantissa *e2);
  20 static b64_sft(struct mantissa *e1, int n);
  21
  22 static
  23 mul_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
  96 static
  97 add_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
 172 static int
 173 cmp_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
 185 static
 186 b64_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
 230 static int
 231 b64_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
 255 obase=16
 256 scale=0
 257 define 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 }
 283 for (i=1;i<28;i++) {
 284         t(10^i)
 285 }
 286 0
 287 for (i=1;i<20;i++) {
 288         t(10^(28*i))
 289 }
 290 0
 291 define 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 }
 317 for (i=1;i<28;i++) {
 318         r(10^i)
 319 }
 320 0
 321 for (i=1;i<20;i++) {
 322         r(10^(28*i))
 323 }
 324 0
 325
 326 */
 327 static 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 };
 357 static 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
 380 static 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
 411 static 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
 438 static
 439 add_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
 541 static
 542 ten_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
 554 char *
 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
 716 static struct EXTEND max_d;
 717
 718 double
 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 }