ext_comp.c 15 KB

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