do_fpar.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639
  1. /*
  2. * Sources of the "FLOATING POINT ARITHMETIC" group instructions
  3. */
  4. /* $Header$ */
  5. #include <em_abs.h>
  6. #include "nofloat.h"
  7. #include "global.h"
  8. #include "log.h"
  9. #include "mem.h"
  10. #include "trap.h"
  11. #include "text.h"
  12. #include "fra.h"
  13. #include "warn.h"
  14. #ifndef NOFLOAT
  15. extern double fpop();
  16. #define MAXDOUBLE 99.e999 /* IEEE infinity */ /*???*/
  17. #define SMALL (1.0/MAXDOUBLE)
  18. PRIVATE double adf(), sbf(), mlf(), dvf();
  19. PRIVATE double ttttp();
  20. PRIVATE double floor(), fabs();
  21. PRIVATE fef(), fif();
  22. #endif NOFLOAT
  23. DoADFl2(arg)
  24. size arg;
  25. {
  26. /* ADF w: Floating add (*) */
  27. #ifndef NOFLOAT
  28. register size l = (L_arg_2() * arg);
  29. double t = fpop(arg_wf(l));
  30. LOG(("@F6 DoADFl2(%ld)", l));
  31. spoilFRA();
  32. fpush(adf(fpop(l), t), l);
  33. #else NOFLOAT
  34. arg = arg;
  35. nofloat();
  36. #endif NOFLOAT
  37. }
  38. DoADFs(hob, wfac)
  39. long hob;
  40. size wfac;
  41. {
  42. /* ADF w: Floating add (*) */
  43. #ifndef NOFLOAT
  44. register size l = (S_arg(hob) * wfac);
  45. double t = fpop(arg_wf(l));
  46. LOG(("@F6 DoADFs(%ld)", l));
  47. spoilFRA();
  48. fpush(adf(fpop(l), t), l);
  49. #else NOFLOAT
  50. hob = hob;
  51. wfac = wfac;
  52. nofloat();
  53. #endif NOFLOAT
  54. }
  55. DoADFz()
  56. {
  57. /* ADF w: Floating add (*) */
  58. #ifndef NOFLOAT
  59. register size l = upop(wsize);
  60. double t = fpop(arg_wf(l));
  61. LOG(("@F6 DoADFz(%ld)", l));
  62. spoilFRA();
  63. fpush(adf(fpop(l), t), l);
  64. #else NOFLOAT
  65. nofloat();
  66. #endif NOFLOAT
  67. }
  68. DoSBFl2(arg)
  69. size arg;
  70. {
  71. /* SBF w: Floating subtract (*) */
  72. #ifndef NOFLOAT
  73. register size l = (L_arg_2() * arg);
  74. double t = fpop(arg_wf(l));
  75. LOG(("@F6 DoSBFl2(%ld)", l));
  76. spoilFRA();
  77. fpush(sbf(fpop(l), t), l);
  78. #else NOFLOAT
  79. arg = arg;
  80. nofloat();
  81. #endif NOFLOAT
  82. }
  83. DoSBFs(hob, wfac)
  84. long hob;
  85. size wfac;
  86. {
  87. /* SBF w: Floating subtract (*) */
  88. #ifndef NOFLOAT
  89. register size l = (S_arg(hob) * wfac);
  90. double t = fpop(arg_wf(l));
  91. LOG(("@F6 DoSBFs(%ld)", l));
  92. spoilFRA();
  93. fpush(sbf(fpop(l), t), l);
  94. #else NOFLOAT
  95. hob = hob;
  96. wfac = wfac;
  97. nofloat();
  98. #endif NOFLOAT
  99. }
  100. DoSBFz()
  101. {
  102. /* SBF w: Floating subtract (*) */
  103. #ifndef NOFLOAT
  104. register size l = upop(wsize);
  105. double t = fpop(arg_wf(l));
  106. LOG(("@F6 DoSBFz(%ld)", l));
  107. spoilFRA();
  108. fpush(sbf(fpop(l), t), l);
  109. #else NOFLOAT
  110. nofloat();
  111. #endif NOFLOAT
  112. }
  113. DoMLFl2(arg)
  114. size arg;
  115. {
  116. /* MLF w: Floating multiply (*) */
  117. #ifndef NOFLOAT
  118. register size l = (L_arg_2() * arg);
  119. double t = fpop(arg_wf(l));
  120. LOG(("@F6 DoMLFl2(%ld)", l));
  121. spoilFRA();
  122. fpush(mlf(fpop(l), t), l);
  123. #else NOFLOAT
  124. arg = arg;
  125. nofloat();
  126. #endif NOFLOAT
  127. }
  128. DoMLFs(hob, wfac)
  129. long hob;
  130. size wfac;
  131. {
  132. /* MLF w: Floating multiply (*) */
  133. #ifndef NOFLOAT
  134. register size l = (S_arg(hob) * wfac);
  135. double t = fpop(arg_wf(l));
  136. LOG(("@F6 DoMLFs(%ld)", l));
  137. spoilFRA();
  138. fpush(mlf(fpop(l), t), l);
  139. #else NOFLOAT
  140. hob = hob;
  141. wfac = wfac;
  142. nofloat();
  143. #endif NOFLOAT
  144. }
  145. DoMLFz()
  146. {
  147. /* MLF w: Floating multiply (*) */
  148. #ifndef NOFLOAT
  149. register size l = upop(wsize);
  150. double t = fpop(arg_wf(l));
  151. LOG(("@F6 DoMLFz(%ld)", l));
  152. spoilFRA();
  153. fpush(mlf(fpop(l), t), l);
  154. #else NOFLOAT
  155. nofloat();
  156. #endif NOFLOAT
  157. }
  158. DoDVFl2(arg)
  159. size arg;
  160. {
  161. /* DVF w: Floating divide (*) */
  162. #ifndef NOFLOAT
  163. register size l = (L_arg_2() * arg);
  164. double t = fpop(arg_wf(l));
  165. LOG(("@F6 DoDVFl2(%ld)", l));
  166. spoilFRA();
  167. fpush(dvf(fpop(l), t), l);
  168. #else NOFLOAT
  169. arg = arg;
  170. nofloat();
  171. #endif NOFLOAT
  172. }
  173. DoDVFs(hob, wfac)
  174. long hob;
  175. size wfac;
  176. {
  177. /* DVF w: Floating divide (*) */
  178. #ifndef NOFLOAT
  179. register size l = (S_arg(hob) * wfac);
  180. double t = fpop(arg_wf(l));
  181. LOG(("@F6 DoDVFs(%ld)", l));
  182. spoilFRA();
  183. fpush(dvf(fpop(l), t), l);
  184. #else NOFLOAT
  185. hob = hob;
  186. wfac = wfac;
  187. nofloat();
  188. #endif NOFLOAT
  189. }
  190. DoDVFz()
  191. {
  192. /* DVF w: Floating divide (*) */
  193. #ifndef NOFLOAT
  194. register size l = upop(wsize);
  195. double t = fpop(arg_wf(l));
  196. LOG(("@F6 DoDVFz(%ld)", l));
  197. spoilFRA();
  198. fpush(dvf(fpop(l), t), l);
  199. #else NOFLOAT
  200. nofloat();
  201. #endif NOFLOAT
  202. }
  203. DoNGFl2(arg)
  204. size arg;
  205. {
  206. /* NGF w: Floating negate (*) */
  207. #ifndef NOFLOAT
  208. register size l = (L_arg_2() * arg);
  209. double t = fpop(arg_wf(l));
  210. LOG(("@F6 DoNGFl2(%ld)", l));
  211. spoilFRA();
  212. fpush(-t, l);
  213. #else NOFLOAT
  214. arg = arg;
  215. nofloat();
  216. #endif NOFLOAT
  217. }
  218. DoNGFz()
  219. {
  220. /* NGF w: Floating negate (*) */
  221. #ifndef NOFLOAT
  222. register size l = upop(wsize);
  223. double t = fpop(arg_wf(l));
  224. LOG(("@F6 DoNGFz(%ld)", l));
  225. spoilFRA();
  226. fpush(-t, l);
  227. #else NOFLOAT
  228. nofloat();
  229. #endif NOFLOAT
  230. }
  231. DoFIFl2(arg)
  232. size arg;
  233. {
  234. /* FIF w: Floating multiply and split integer and fraction part (*) */
  235. #ifndef NOFLOAT
  236. register size l = (L_arg_2() * arg);
  237. double t = fpop(arg_wf(l));
  238. LOG(("@F6 DoFIFl2(%ld)", l));
  239. spoilFRA();
  240. fif(fpop(l), t, l);
  241. #else NOFLOAT
  242. arg = arg;
  243. nofloat();
  244. #endif NOFLOAT
  245. }
  246. DoFIFz()
  247. {
  248. /* FIF w: Floating multiply and split integer and fraction part (*) */
  249. #ifndef NOFLOAT
  250. register size l = upop(wsize);
  251. double t = fpop(arg_wf(l));
  252. LOG(("@F6 DoFIFz(%ld)", l));
  253. spoilFRA();
  254. fif(fpop(l), t, l);
  255. #else NOFLOAT
  256. nofloat();
  257. #endif NOFLOAT
  258. }
  259. DoFEFl2(arg)
  260. size arg;
  261. {
  262. /* FEF w: Split floating number in exponent and fraction part (*) */
  263. #ifndef NOFLOAT
  264. register size l = (L_arg_2() * arg);
  265. LOG(("@F6 DoFEFl2(%ld)", l));
  266. spoilFRA();
  267. fef(fpop(arg_wf(l)), l);
  268. #else NOFLOAT
  269. arg = arg;
  270. nofloat();
  271. #endif NOFLOAT
  272. }
  273. DoFEFz()
  274. {
  275. /* FEF w: Split floating number in exponent and fraction part (*) */
  276. #ifndef NOFLOAT
  277. register size l = upop(wsize);
  278. LOG(("@F6 DoFEFz(%ld)", l));
  279. spoilFRA();
  280. fef(fpop(arg_wf(l)), l);
  281. #else NOFLOAT
  282. nofloat();
  283. #endif NOFLOAT
  284. }
  285. #ifndef NOFLOAT
  286. /* Service routines */
  287. PRIVATE double adf(f1, f2) /* returns f1 + f2 */
  288. double f1, f2;
  289. {
  290. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  291. if (f1 > 0.0 && f2 > 0.0) {
  292. if (MAXDOUBLE - f1 < f2) {
  293. trap(EFOVFL);
  294. return (0.0);
  295. }
  296. }
  297. else if (f1 < 0.0 && f2 < 0.0) {
  298. if (-(MAXDOUBLE + f1) > f2) {
  299. trap(EFOVFL);
  300. return (0.0);
  301. }
  302. }
  303. }
  304. return (f1 + f2);
  305. }
  306. PRIVATE double sbf(f1, f2) /* returns f1 - f2 */
  307. double f1, f2;
  308. {
  309. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  310. if (f2 < 0.0 && f1 > 0.0) {
  311. if (MAXDOUBLE - f1 < -f2) {
  312. trap(EFOVFL);
  313. return (0.0);
  314. }
  315. }
  316. else if (f2 > 0.0 && f1 < 0.0) {
  317. if (f2 - MAXDOUBLE > f1) {
  318. trap(EFOVFL);
  319. return (0.0);
  320. }
  321. }
  322. }
  323. return (f1 - f2);
  324. }
  325. PRIVATE double mlf(f1, f2) /* returns f1 * f2 */
  326. double f1, f2;
  327. {
  328. double ff1 = fabs(f1), ff2 = fabs(f2);
  329. if (f1 == 0.0 || f2 == 0.0)
  330. return (0.0);
  331. if ((ff1 >= 1.0 && ff2 <= 1.0) || (ff2 >= 1.0 && ff1 <= 1.0))
  332. return (f1 * f2);
  333. if (must_test && !(IgnMask&BIT(EFUNFL))) {
  334. if (ff1 < 1.0 && ff2 < 1.0) {
  335. if (SMALL / ff1 > ff2) {
  336. trap(EFUNFL);
  337. return (0.0);
  338. }
  339. return (f1 * f2);
  340. }
  341. }
  342. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  343. if (MAXDOUBLE / ff1 < ff2) {
  344. trap(EFOVFL);
  345. return (0.0);
  346. }
  347. }
  348. return (f1 * f2);
  349. }
  350. PRIVATE double dvf(f1, f2) /* returns f1 / f2 */
  351. double f1, f2;
  352. {
  353. double ff1 = fabs(f1), ff2 = fabs(f2);
  354. if (f2 == 0.0) {
  355. if (!(IgnMask&BIT(EFDIVZ))) {
  356. trap(EFDIVZ);
  357. }
  358. else return (0.0);
  359. }
  360. if (f1 == 0.0)
  361. return (0.0);
  362. if ((ff2 >= 1.0 && ff1 >= 1.0) || (ff1 <= 1.0 && ff2 <= 1.0))
  363. return (f1 / f2);
  364. if (must_test && !(IgnMask&BIT(EFUNFL))) {
  365. if (ff2 > 1.0 && ff1 < 1.0) {
  366. if (SMALL / ff2 > ff1) {
  367. trap(EFUNFL);
  368. return (0.0);
  369. }
  370. return (f1 / f2);
  371. }
  372. }
  373. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  374. if (MAXDOUBLE * ff2 < ff1) {
  375. trap(EFOVFL);
  376. return (0.0);
  377. }
  378. }
  379. return (f1 / f2);
  380. }
  381. PRIVATE fif(f1, f2, n)
  382. double f1, f2;
  383. size n;
  384. {
  385. double f = mlf(f1, f2);
  386. double fl = floor(fabs(f));
  387. fpush(fabs(f) - fl, n); /* push fraction */
  388. fpush((f < 0.0) ? -fl : fl, n); /* push integer-part */
  389. }
  390. PRIVATE fef(f, n)
  391. double f;
  392. size n;
  393. {
  394. register long exponent, sign = (long) (f < 0.0);
  395. for (f = fabs(f), exponent = 0; f >= 1.0; exponent++)
  396. f /= 2.0;
  397. for (; f < 0.5; exponent--)
  398. f *= 2.0;
  399. fpush((sign) ? -f : f, n); /* push mantissa */
  400. npush(exponent, wsize); /* push exponent */
  401. }
  402. /* floating point service routines, to avoid having to use -lm */
  403. PRIVATE double fabs(f)
  404. double f;
  405. {
  406. return (f < 0.0 ? -f : f);
  407. }
  408. PRIVATE double floor(f)
  409. double f;
  410. {
  411. double res, d;
  412. register int sign = 1;
  413. /* eliminate the sign */
  414. if (f < 0) {
  415. sign = -1, f = -f;
  416. }
  417. /* get the largest power of 2 <= f */
  418. d = 1.0;
  419. while (f - d >= d) {
  420. d *= 2.0;
  421. }
  422. /* reconstruct f by deminishing powers of 2 */
  423. res = 0.0;
  424. while (d >= 1.0) {
  425. if (res + d <= f)
  426. res += d;
  427. d /= 2.0;
  428. }
  429. /* undo the sign elimination */
  430. if (sign == -1) {
  431. res = -res, f = -f;
  432. if (res > f)
  433. res -= 1.0;
  434. }
  435. return res;
  436. }
  437. PRIVATE double ttttp(f, n) /* times ten to the power */
  438. double f;
  439. {
  440. while (n > 0) {
  441. f = mlf(f, 10.0);
  442. n--;
  443. }
  444. while (n < 0) {
  445. f = dvf(f, 10.0);
  446. n++;
  447. }
  448. return f;
  449. }
  450. /* Str2double is used to initialize the global data area with floats;
  451. we do not use, e.g., sscanf(), to be able to check the grammar of
  452. the string and to give warnings.
  453. */
  454. double str2double(str)
  455. char *str;
  456. {
  457. register char b;
  458. register int sign = 1; /* either +1 or -1 */
  459. register int frac = 0; /* how far in fraction part ? */
  460. register int ex; /* to store exponent */
  461. double mantissa = 0.0; /* to store mantissa */
  462. double d; /* double to be returned */
  463. b = *str++;
  464. if (b == '-') {
  465. sign = -1;
  466. b = *str++;
  467. }
  468. else if (b == '+') {
  469. sign = 1;
  470. b = *str++;
  471. }
  472. if ('0' <= b && b <= '9') {
  473. mantissa = (double) (b-'0');
  474. }
  475. else if (b == '.') {
  476. /* part before dot cannot be empty */
  477. warning(WBADFLOAT);
  478. frac = 1;
  479. }
  480. else {
  481. goto BadFloat;
  482. }
  483. LOG((" q9 str2double : (before while) mantissa = %20.20g", mantissa));
  484. while ((b = *str++) != 'e' && b != 'E' && b != '\0') {
  485. if (b == '.') {
  486. if (frac == 0) {
  487. frac++;
  488. }
  489. else { /* there already was a '.' in input */
  490. goto BadFloat;
  491. }
  492. }
  493. else if ('0' <= b && b <= '9') {
  494. double bval = b - '0';
  495. if (frac) {
  496. mantissa =
  497. adf(mantissa, ttttp(bval, -frac));
  498. frac++;
  499. }
  500. else {
  501. mantissa =
  502. adf(mlf(mantissa, 10.0), bval);
  503. }
  504. }
  505. else {
  506. goto BadFloat;
  507. }
  508. LOG((" q9 str2double : (inside while) mantissa = %20.20g",
  509. mantissa));
  510. }
  511. LOG((" q9 str2double : mantissa = %10.10g", mantissa));
  512. mantissa = sign * mantissa;
  513. if (b == '\0')
  514. return (mantissa);
  515. /* else we have b == 'e' or b== 'E' */
  516. /* Optional sign for exponent */
  517. b = *str++;
  518. if (b == '-') {
  519. sign = -1;
  520. b = *str++;
  521. }
  522. else if (b == '+') {
  523. sign = 1;
  524. b = *str++;
  525. }
  526. else {
  527. sign = 1;
  528. }
  529. ex = 0;
  530. do {
  531. if ('0' <= b && b <= '9') {
  532. ex = 10*ex + (b-'0');
  533. }
  534. else {
  535. goto BadFloat;
  536. }
  537. } while ((b = *str++) != '\0');
  538. LOG((" q9 str2double : exponent = %d", ex));
  539. /* Construct total value of float */
  540. ex = sign * ex;
  541. d = ttttp(mantissa, ex);
  542. return (d);
  543. BadFloat:
  544. fatal("Float garbled in loadfile");
  545. return (0.0);
  546. }
  547. #else NOFLOAT
  548. nofloat() {
  549. fatal("attempt to execute a floating point instruction on an EM machine without FP");
  550. }
  551. #endif NOFLOAT