flt_str2fl.c 9.7 KB

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