modules/src/flt_arith/flt_str2fl.c

   1 /*
   2   (c) copyright 1988 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: flt_str2fl.c,v 1.16 1994/06/24 11:16:07 ceriel Exp $ */
   7
   8 #include <ctype.h>
   9 #include "flt_misc.h"
  10
  11 /* The following tables can be computed with the following bc(1)
  12    program:
  13
  14 obase=16
  15 scale=0
  16 define t(x){
  17         auto a, b, c
  18         a=2;b=1;c=2^32;n=1
  19         while(a<x) {
  20                 b=a;n+=n;a*=a
  21         }
  22         n/=2
  23         a=b
  24         while(b<x) {
  25                 a=b;b*=c;n+=32
  26         }
  27         n-=32
  28         b=a
  29         while(a<x) {
  30                 b=a;a+=a;n+=1
  31         }
  32         n-=1
  33         x*=16^16
  34         b=x%a
  35         x/=a
  36         if(a<=(2*b)) x+=1
  37         obase=10
  38         n
  39         obase=16
  40         return(x)
  41 }
  42 for (i=1;i<28;i++) {
  43         t(10^i)
  44 }
  45 0
  46 for (i=1;i<20;i++) {
  47         t(10^(28*i))
  48 }
  49 0
  50 define r(x){
  51         auto a, b, c
  52         a=2;b=1;c=2^32;n=1
  53         while(a<x) {
  54                 b=a;n+=n;a*=a
  55         }
  56         n/=2
  57         a=b
  58         while(b<x) {
  59                 a=b;b*=c;n+=32
  60         }
  61         n-=32
  62         b=a
  63         while(a<x) {
  64                 b=a;a+=a;n+=1
  65         }
  66         a=b
  67         a*=16^16
  68         b=a%x
  69         a/=x
  70         if(x<=(2*b)) a+=1
  71         obase=10
  72         -n
  73         obase=16
  74         return(a)
  75 }
  76 for (i=1;i<28;i++) {
  77         r(10^i)
  78 }
  79 0
  80 for (i=1;i<20;i++) {
  81         r(10^(28*i))
  82 }
  83 0
  84
  85 */
  86 static flt_arith s10pow[] = {   /* representation of 10 ** i */
  87         { 0,    0,      0x80000000,     0 },
  88         { 0,    3,      0xA0000000,     0 },
  89         { 0,    6,      0xC8000000,     0 },
  90         { 0,    9,      0xFA000000,     0 },
  91         { 0,    13,     0x9C400000,     0 },
  92         { 0,    16,     0xC3500000,     0 },
  93         { 0,    19,     0xF4240000,     0 },
  94         { 0,    23,     0x98968000,     0 },
  95         { 0,    26,     0xBEBC2000,     0 },
  96         { 0,    29,     0xEE6B2800,     0 },
  97         { 0,    33,     0x9502F900,     0 },
  98         { 0,    36,     0xBA43B740,     0 },
  99         { 0,    39,     0xE8D4A510,     0 },
 100         { 0,    43,     0x9184E72A,     0 },
 101         { 0,    46,     0xB5E620F4,     0x80000000 },
 102         { 0,    49,     0xE35FA931,     0xA0000000 },
 103         { 0,    53,     0x8E1BC9BF,     0x04000000 },
 104         { 0,    56,     0xB1A2BC2E,     0xC5000000 },
 105         { 0,    59,     0xDE0B6B3A,     0x76400000 },
 106         { 0,    63,     0x8AC72304,     0x89E80000 },
 107         { 0,    66,     0xAD78EBC5,     0xAC620000 },
 108         { 0,    69,     0xD8D726B7,     0x177A8000 },
 109         { 0,    73,     0x87867832,     0x6EAC9000 },
 110         { 0,    76,     0xA968163F,     0x0A57B400 },
 111         { 0,    79,     0xD3C21BCE,     0xCCEDA100 },
 112         { 0,    83,     0x84595161,     0x401484A0 },
 113         { 0,    86,     0xA56FA5B9,     0x9019A5C8 },
 114         { 0,    89,     0xCECB8F27,     0xF4200F3A }
 115 };
 116 static flt_arith big_10pow[] = {  /* representation of 10 ** (28*i) */
 117         { 0,    0,      0x80000000,     0 },
 118         { 0,    93,     0x813F3978,     0xF8940984 },
 119         { 0,    186,    0x82818F12,     0x81ED44A0 },
 120         { 0,    279,    0x83C7088E,     0x1AAB65DB },
 121         { 0,    372,    0x850FADC0,     0x9923329E },
 122         { 0,    465,    0x865B8692,     0x5B9BC5C2 },
 123         { 0,    558,    0x87AA9AFF,     0x79042287 },
 124         { 0,    651,    0x88FCF317,     0xF22241E2 },
 125         { 0,    744,    0x8A5296FF,     0xE33CC930 },
 126         { 0,    837,    0x8BAB8EEF,     0xB6409C1A },
 127         { 0,    930,    0x8D07E334,     0x55637EB3 },
 128         { 0,    1023,   0x8E679C2F,     0x5E44FF8F },
 129         { 0,    1116,   0x8FCAC257,     0x558EE4E6 },
 130         { 0,    1209,   0x91315E37,     0xDB165AA9 },
 131         { 0,    1302,   0x929B7871,     0xDE7F22B9 },
 132         { 0,    1395,   0x940919BB,     0xD4620B6D },
 133         { 0,    1488,   0x957A4AE1,     0xEBF7F3D4 },
 134         { 0,    1581,   0x96EF14C6,     0x454AA840 },
 135         { 0,    1674,   0x98678061,     0x27ECE4F5 },
 136         { 0,    1767,   0x99E396C1,     0x3A3ACFF2 }
 137 };
 138
 139 static flt_arith r_10pow[] = { /* representation of 10 ** -i */
 140         { 0,    0,      0x80000000,     0 },
 141         { 0,    -4,     0xCCCCCCCC,     0xCCCCCCCD },
 142         { 0,    -7,     0xA3D70A3D,     0x70A3D70A },
 143         { 0,    -10,    0x83126E97,     0x8D4FDF3B },
 144         { 0,    -14,    0xD1B71758,     0xE219652C },
 145         { 0,    -17,    0xA7C5AC47,     0x1B478423 },
 146         { 0,    -20,    0x8637BD05,     0xAF6C69B6 },
 147         { 0,    -24,    0xD6BF94D5,     0xE57A42BC },
 148         { 0,    -27,    0xABCC7711,     0x8461CEFD },
 149         { 0,    -30,    0x89705F41,     0x36B4A597 },
 150         { 0,    -34,    0xDBE6FECE,     0xBDEDD5BF },
 151         { 0,    -37,    0xAFEBFF0B,     0xCB24AAFF },
 152         { 0,    -40,    0x8CBCCC09,     0x6F5088CC },
 153         { 0,    -44,    0xE12E1342,     0x4BB40E13 },
 154         { 0,    -47,    0xB424DC35,     0x095CD80F },
 155         { 0,    -50,    0x901D7CF7,     0x3AB0ACD9 },
 156         { 0,    -54,    0xE69594BE,     0xC44DE15B },
 157         { 0,    -57,    0xB877AA32,     0x36A4B449 },
 158         { 0,    -60,    0x9392EE8E,     0x921D5D07 },
 159         { 0,    -64,    0xEC1E4A7D,     0xB69561A5 },
 160         { 0,    -67,    0xBCE50864,     0x92111AEB },
 161         { 0,    -70,    0x971DA050,     0x74DA7BEF },
 162         { 0,    -74,    0xF1C90080,     0xBAF72CB1 },
 163         { 0,    -77,    0xC16D9A00,     0x95928A27 },
 164         { 0,    -80,    0x9ABE14CD,     0x44753B53 },
 165         { 0,    -84,    0xF79687AE,     0xD3EEC551 },
 166         { 0,    -87,    0xC6120625,     0x76589DDB },
 167         { 0,    -90,    0x9E74D1B7,     0x91E07E48 }
 168 };
 169
 170 static flt_arith r_big_10pow[] = { /* representation of 10 ** -(28*i) */
 171         { 0,    0,      0x80000000,     0 },
 172         { 0,    -94,    0xFD87B5F2,     0x8300CA0E },
 173         { 0,    -187,   0xFB158592,     0xBE068D2F },
 174         { 0,    -280,   0xF8A95FCF,     0x88747D94 },
 175         { 0,    -373,   0xF64335BC,     0xF065D37D },
 176         { 0,    -466,   0xF3E2F893,     0xDEC3F126 },
 177         { 0,    -559,   0xF18899B1,     0xBC3F8CA2 },
 178         { 0,    -652,   0xEF340A98,     0x172AACE5 },
 179         { 0,    -745,   0xECE53CEC,     0x4A314EBE },
 180         { 0,    -838,   0xEA9C2277,     0x23EE8BCB },
 181         { 0,    -931,   0xE858AD24,     0x8F5C22CA },
 182         { 0,    -1024,  0xE61ACF03,     0x3D1A45DF },
 183         { 0,    -1117,  0xE3E27A44,     0x4D8D98B8 },
 184         { 0,    -1210,  0xE1AFA13A,     0xFBD14D6E },
 185         { 0,    -1303,  0xDF82365C,     0x497B5454 },
 186         { 0,    -1396,  0xDD5A2C3E,     0xAB3097CC },
 187         { 0,    -1489,  0xDB377599,     0xB6074245 },
 188         { 0,    -1582,  0xD91A0545,     0xCDB51186 },
 189         { 0,    -1675,  0xD701CE3B,     0xD387BF48 },
 190         { 0,    -1768,  0xD4EEC394,     0xD6258BF8 }
 191 };
 192
 193 #define BIGSZ   (sizeof(big_10pow)/sizeof(big_10pow[0]))
 194 #define SMALLSZ (sizeof(s10pow)/sizeof(s10pow[0]))
 195
 196 static
 197 add_exponent(e, exp)
 198         register flt_arith *e;
 199 {
 200         int neg = exp < 0;
 201         int divsz, modsz;
 202         flt_arith x;
 203         register int status = 0;
 204
 205         if (neg) exp = -exp;
 206         divsz = exp / SMALLSZ;
 207         modsz = exp % SMALLSZ;
 208         if (!status) status = flt_status;
 209         flt_mul(e, (neg ? r_10pow : s10pow) + modsz, &x);
 210         if (!status) status = flt_status;
 211         while (divsz >= BIGSZ) {
 212                 flt_mul(&x, neg ? &r_big_10pow[BIGSZ-1] : &big_10pow[BIGSZ-1],&x);
 213                 if (!status) status = flt_status;
 214                 divsz -= BIGSZ-1;
 215         }
 216         flt_mul(&x, (neg ? r_big_10pow : big_10pow) + divsz, e);
 217         if (!status) status = flt_status;
 218         flt_status = status;
 219 }
 220
 221 void
 222 flt_str2flt(s, e)
 223         register char           *s;
 224         register flt_arith      *e;
 225 {
 226         register int    c;
 227         int             dotseen = 0;
 228         int             digitseen = 0;
 229         int             exp = 0;
 230
 231         while (isspace(*s)) s++;
 232
 233         flt_status = 0;
 234
 235         e->flt_sign = 0;
 236         e->flt_exp = 0;
 237         e->m1 = e->m2 = 0;
 238
 239         c = *s;
 240         switch(c) {
 241         case '-':
 242                 e->flt_sign = 1;
 243         case '+':
 244                 s++;
 245         }
 246         while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) {
 247                 if (c == '.') continue;
 248                 digitseen = 1;
 249                 if (e->m1 >= 0 && e->m1 <= 0x7FFFFFFF/5) {
 250                         struct flt_mantissa     a1;
 251
 252                         a1 = e->flt_mantissa;
 253                         flt_b64_sft(&(e->flt_mantissa), -3);
 254                         flt_b64_sft(&a1, -1);
 255                         flt_b64_add(&(e->flt_mantissa), &a1);
 256                         a1.flt_h_32 = 0;
 257                         a1.flt_l_32 = c - '0';
 258                         flt_b64_add(&(e->flt_mantissa), &a1);
 259                 }
 260                 else exp++;
 261                 if (dotseen) exp--;
 262         }
 263         if (! digitseen) {
 264                 flt_status = FLT_NOFLT;
 265                 return;
 266         }
 267
 268         if (c == 'E' || c == 'e') {
 269                 int     exp1 = 0;
 270                 int     sign = 1;
 271
 272                 switch(*s) {
 273                 case '-':
 274                         sign = -1;
 275                 case '+':
 276                         s++;
 277                 }
 278                 if (c = *s, isdigit(c)) {
 279                         do {
 280                                 exp1 = 10 * exp1 + (c - '0');
 281                         } while (c = *++s, isdigit(c));
 282                 }
 283                 exp += sign * exp1;
 284         }
 285         if (e->m1 == 0 && e->m2 == 0) return;
 286         e->flt_exp = 63;
 287         flt_nrm(e);
 288         add_exponent(e, exp);
 289         flt_chk(e);
 290 }
 291
 292 #define NDIG 18
 293
 294 static char *
 295 flt_ecvt(e, decpt, sign)
 296         register flt_arith *e;
 297         int *decpt, *sign;
 298 {
 299         /*      Like ecvt(), but for extended precision */
 300
 301         static char buf[NDIG+1];
 302         register char *p = buf;
 303         register char *pe;
 304         register int findex = 0;
 305
 306         pe = &buf[NDIG];
 307         buf[0] = '\0';
 308
 309         *sign = 0;
 310         if (e->flt_sign) {
 311                 *sign = 1;
 312                 e->flt_sign = 0;
 313         }
 314
 315         *decpt = 0;
 316         if (e->m1 != 0) {
 317                 register flt_arith *pp = &big_10pow[1];
 318
 319                 findex = 1;
 320                 while (flt_cmp(e, &big_10pow[BIGSZ-1]) >= 0) {
 321                         flt_mul(e,&r_big_10pow[BIGSZ-1],e);
 322                         *decpt += (BIGSZ-1)*SMALLSZ;
 323                 }
 324                 while (flt_cmp(e,pp) >= 0) {
 325                         pp++;
 326                         findex++;
 327                 }
 328                 findex--;
 329                 flt_mul(e,&r_big_10pow[findex],e);
 330                 *decpt += findex*SMALLSZ;
 331                 pp = &s10pow[1];
 332                 findex = 1;
 333                 while (pp < &s10pow[SMALLSZ] && flt_cmp(e, pp) >= 0) {
 334                         pp++;
 335                         findex++;
 336                 }
 337                 findex--;
 338                 *decpt += findex;
 339
 340                 if (flt_cmp(e, &s10pow[0]) < 0) {
 341                         while (flt_cmp(e, &r_big_10pow[BIGSZ-1]) < 0) {
 342                                 flt_mul(e,&big_10pow[BIGSZ-1],e);
 343                                 *decpt -= (BIGSZ-1)*SMALLSZ;
 344                         }
 345                         pp = &r_big_10pow[1];
 346                         findex = 1;
 347                         while(flt_cmp(e,pp) < 0) {
 348                                 pp++;
 349                                 findex++;
 350                         }
 351                         findex--;
 352                         flt_mul(e,&big_10pow[findex],e);
 353                         *decpt -= findex*SMALLSZ;
 354                         /* here, value >= 10 ** -28 */
 355                         flt_mul(e, &s10pow[1], e);
 356                         (*decpt)--;
 357                         pp = &r_10pow[0];
 358                         findex = 0;
 359                         while(flt_cmp(e, pp) < 0) {
 360                                 pp++;
 361                                 findex++;
 362                         }
 363                         flt_mul(e, &s10pow[findex], e);
 364                         *decpt -= findex;
 365                         findex = 0;
 366                 }
 367                 (*decpt)++;     /* because now value in [1.0, 10.0) */
 368         }
 369         while (p <= pe) {
 370                 if (findex) {
 371                         flt_arith tc, oldtc;
 372                         int count = 0;
 373
 374                         oldtc.flt_exp = 0;
 375                         oldtc.flt_sign = 0;
 376                         oldtc.m1 = 0;
 377                         oldtc.m2 = 0;
 378                         tc = s10pow[findex];
 379                         while (flt_cmp(e, &tc) >= 0) {
 380                                 oldtc = tc;
 381                                 flt_add(&tc, &s10pow[findex], &tc);
 382                                 count++;
 383                         }
 384                         *p++ = count + '0';
 385                         oldtc.flt_sign = 1;
 386                         flt_add(e, &oldtc, e);
 387                         findex--;
 388                         continue;
 389                 }
 390                 if (e->flt_exp >= 0 && e->m1 != 0) {
 391                         flt_arith x;
 392
 393                         x.m2 = 0; x.flt_exp = e->flt_exp;
 394                         x.flt_sign = 1;
 395                         x.m1 = (e->m1 >> 1) & 0x7FFFFFFF;
 396                         x.m1 = x.m1>>(30-e->flt_exp);
 397                         *p++ = (x.m1) + '0';
 398                         if (x.m1) {
 399                                 x.m1 = x.m1 << (31-e->flt_exp);
 400                                 flt_add(e, &x, e);
 401                         }
 402                 }
 403                 else *p++ = '0';
 404                 if (e->m1) flt_mul(e, &s10pow[1], e);
 405         }
 406         if (pe >= buf) {
 407                 p = pe;
 408                 *p += 5;        /* round of at the end */
 409                 while (*p > '9') {
 410                         *p = '0';
 411                         if (p > buf) ++*--p;
 412                         else {
 413                                 *p = '1';
 414                                 ++*decpt;
 415                         }
 416                 }
 417                 *pe = '\0';
 418                 while (--pe > buf && *pe == '0') *pe = '\0';
 419         }
 420         return buf;
 421 }
 422
 423 void
 424 flt_flt2str(e, buf, bufsize)
 425         flt_arith *e;
 426         char *buf;
 427 {
 428
 429         int sign, dp;
 430         register int i;
 431         register char *s1;
 432         char Xbuf[NDIG+12];
 433         register char *s = Xbuf;
 434         flt_arith e1;
 435
 436         e1 = *e;
 437         flt_status = 0;
 438         s1 = flt_ecvt(&e1,&dp,&sign);
 439         if (sign)
 440                 *s++ = '-';
 441         *s++ = *s1++;
 442         *s++ = '.';
 443         for (i = NDIG-1; i > 0; i--) {
 444                 if (*s1) *s++ = *s1++;
 445                 else {
 446                         if (i == NDIG-1) *s++ = '0';
 447                         break;
 448                 }
 449         }
 450         if (e->m1 | e->m2) {
 451                 --dp ;
 452         }
 453         if (dp != 0) {
 454                 *s++ = 'e';
 455                 if ( dp<0 ) {
 456                         *s++ = '-' ; dp= -dp ;
 457                 } else {
 458                         *s++ = '+' ;
 459                 }
 460                 s1 = &Xbuf[NDIG+12];
 461                 *--s1 = '\0';
 462                 do {
 463                         *--s1 = dp % 10 + '0';
 464                         dp /= 10;
 465                 } while (dp != 0);
 466                 while (*s1) *s++ = *s1++;
 467         }
 468         *s++ = '\0';
 469         if (s - Xbuf > bufsize) {
 470                 flt_status = FLT_BTSM;
 471                 return;
 472         }
 473         s = Xbuf;
 474         s1 = buf;
 475         do {
 476                 *s1++ = *s;
 477         } while (*s++);
 478 }