do_fpar.c 7.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459
  1. /*
  2. * Sources of the "FLOATING POINT ARITHMETIC" group instructions
  3. */
  4. /* $Id$ */
  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. #ifdef __STDC__
  17. #include <float.h>
  18. #define MAXDOUBLE DBL_MAX
  19. #else /* not __STDC__ */
  20. #if defined(vax) || defined(pdp) || defined(__vax) || defined(__pdp)
  21. #define MAXDOUBLE 1.701411834604692293e+38
  22. #else
  23. #define MAXDOUBLE 1.7976931348623157e+308
  24. #endif
  25. #endif /* not __STDC__ */
  26. #define SMALL (1.0/MAXDOUBLE)
  27. PRIVATE double adf(), sbf(), mlf(), dvf();
  28. PRIVATE double ttttp();
  29. PRIVATE double floor(), fabs();
  30. PRIVATE fef(), fif();
  31. #endif /* NOFLOAT */
  32. DoADF(l)
  33. register size l;
  34. {
  35. /* ADF w: Floating add (*) */
  36. #ifndef NOFLOAT
  37. double t = fpop(arg_wf(l));
  38. LOG(("@F6 DoADF(%ld)", l));
  39. spoilFRA();
  40. fpush(adf(fpop(l), t), l);
  41. #else /* NOFLOAT */
  42. nofloat();
  43. #endif /* NOFLOAT */
  44. }
  45. DoSBF(l)
  46. register size l;
  47. {
  48. /* SBF w: Floating subtract (*) */
  49. #ifndef NOFLOAT
  50. double t = fpop(arg_wf(l));
  51. LOG(("@F6 DoSBF(%ld)", l));
  52. spoilFRA();
  53. fpush(sbf(fpop(l), t), l);
  54. #else /* NOFLOAT */
  55. nofloat();
  56. #endif /* NOFLOAT */
  57. }
  58. DoMLF(l)
  59. register size l;
  60. {
  61. /* MLF w: Floating multiply (*) */
  62. #ifndef NOFLOAT
  63. double t = fpop(arg_wf(l));
  64. LOG(("@F6 DoMLF(%ld)", l));
  65. spoilFRA();
  66. fpush(mlf(fpop(l), t), l);
  67. #else /* NOFLOAT */
  68. nofloat();
  69. #endif /* NOFLOAT */
  70. }
  71. DoDVF(l)
  72. register size l;
  73. {
  74. /* DVF w: Floating divide (*) */
  75. #ifndef NOFLOAT
  76. double t = fpop(arg_wf(l));
  77. LOG(("@F6 DoDVF(%ld)", l));
  78. spoilFRA();
  79. fpush(dvf(fpop(l), t), l);
  80. #else /* NOFLOAT */
  81. nofloat();
  82. #endif /* NOFLOAT */
  83. }
  84. DoNGF(l)
  85. register size l;
  86. {
  87. /* NGF w: Floating negate (*) */
  88. #ifndef NOFLOAT
  89. double t = fpop(arg_wf(l));
  90. LOG(("@F6 DoNGF(%ld)", l));
  91. spoilFRA();
  92. fpush(-t, l);
  93. #else /* NOFLOAT */
  94. nofloat();
  95. #endif /* NOFLOAT */
  96. }
  97. DoFIF(l)
  98. register size l;
  99. {
  100. /* FIF w: Floating multiply and split integer and fraction part (*) */
  101. #ifndef NOFLOAT
  102. double t = fpop(arg_wf(l));
  103. LOG(("@F6 DoFIF(%ld)", l));
  104. spoilFRA();
  105. fif(fpop(l), t, l);
  106. #else /* NOFLOAT */
  107. nofloat();
  108. #endif /* NOFLOAT */
  109. }
  110. DoFEF(l)
  111. register size l;
  112. {
  113. /* FEF w: Split floating number in exponent and fraction part (*) */
  114. #ifndef NOFLOAT
  115. LOG(("@F6 DoFEF(%ld)", l));
  116. spoilFRA();
  117. fef(fpop(arg_wf(l)), l);
  118. #else /* NOFLOAT */
  119. nofloat();
  120. #endif /* NOFLOAT */
  121. }
  122. #ifndef NOFLOAT
  123. /* Service routines */
  124. PRIVATE double adf(f1, f2) /* returns f1 + f2 */
  125. double f1, f2;
  126. {
  127. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  128. if (f1 > 0.0 && f2 > 0.0) {
  129. if (MAXDOUBLE - f1 < f2) {
  130. trap(EFOVFL);
  131. return (0.0);
  132. }
  133. }
  134. else if (f1 < 0.0 && f2 < 0.0) {
  135. if (-(MAXDOUBLE + f1) > f2) {
  136. trap(EFOVFL);
  137. return (0.0);
  138. }
  139. }
  140. }
  141. return (f1 + f2);
  142. }
  143. PRIVATE double sbf(f1, f2) /* returns f1 - f2 */
  144. double f1, f2;
  145. {
  146. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  147. if (f2 < 0.0 && f1 > 0.0) {
  148. if (MAXDOUBLE - f1 < -f2) {
  149. trap(EFOVFL);
  150. return (0.0);
  151. }
  152. }
  153. else if (f2 > 0.0 && f1 < 0.0) {
  154. if (f2 - MAXDOUBLE > f1) {
  155. trap(EFOVFL);
  156. return (0.0);
  157. }
  158. }
  159. }
  160. return (f1 - f2);
  161. }
  162. PRIVATE double mlf(f1, f2) /* returns f1 * f2 */
  163. double f1, f2;
  164. {
  165. double ff1 = fabs(f1), ff2 = fabs(f2);
  166. if (f1 == 0.0 || f2 == 0.0)
  167. return (0.0);
  168. if ((ff1 >= 1.0 && ff2 <= 1.0) || (ff2 >= 1.0 && ff1 <= 1.0))
  169. return (f1 * f2);
  170. if (must_test && !(IgnMask&BIT(EFUNFL))) {
  171. if (ff1 < 1.0 && ff2 < 1.0) {
  172. if (SMALL / ff1 > ff2) {
  173. trap(EFUNFL);
  174. return (0.0);
  175. }
  176. return (f1 * f2);
  177. }
  178. }
  179. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  180. if (MAXDOUBLE / ff1 < ff2) {
  181. trap(EFOVFL);
  182. return (0.0);
  183. }
  184. }
  185. return (f1 * f2);
  186. }
  187. PRIVATE double dvf(f1, f2) /* returns f1 / f2 */
  188. double f1, f2;
  189. {
  190. double ff1 = fabs(f1), ff2 = fabs(f2);
  191. if (f2 == 0.0) {
  192. if (!(IgnMask&BIT(EFDIVZ))) {
  193. trap(EFDIVZ);
  194. }
  195. else return (0.0);
  196. }
  197. if (f1 == 0.0)
  198. return (0.0);
  199. if ((ff2 >= 1.0 && ff1 >= 1.0) || (ff1 <= 1.0 && ff2 <= 1.0))
  200. return (f1 / f2);
  201. if (must_test && !(IgnMask&BIT(EFUNFL))) {
  202. if (ff2 > 1.0 && ff1 < 1.0) {
  203. if (SMALL / ff2 > ff1) {
  204. trap(EFUNFL);
  205. return (0.0);
  206. }
  207. return (f1 / f2);
  208. }
  209. }
  210. if (must_test && !(IgnMask&BIT(EFOVFL))) {
  211. if (MAXDOUBLE * ff2 < ff1) {
  212. trap(EFOVFL);
  213. return (0.0);
  214. }
  215. }
  216. return (f1 / f2);
  217. }
  218. PRIVATE fif(f1, f2, n)
  219. double f1, f2;
  220. size n;
  221. {
  222. double f = mlf(f1, f2);
  223. double fl = floor(fabs(f));
  224. fpush(fabs(f) - fl, n); /* push fraction */
  225. fpush((f < 0.0) ? -fl : fl, n); /* push integer-part */
  226. }
  227. PRIVATE fef(f, n)
  228. double f;
  229. size n;
  230. {
  231. register long exponent, sign = (long) (f < 0.0);
  232. if (f == 0.0) {
  233. fpush(f, n);
  234. wpush(0L);
  235. return;
  236. }
  237. for (f = fabs(f), exponent = 0; f >= 1.0; exponent++)
  238. f /= 2.0;
  239. for (; f < 0.5; exponent--)
  240. f *= 2.0;
  241. fpush((sign) ? -f : f, n); /* push mantissa */
  242. wpush(exponent); /* push exponent */
  243. }
  244. /* floating point service routines, to avoid having to use -lm */
  245. PRIVATE double fabs(f)
  246. double f;
  247. {
  248. return (f < 0.0 ? -f : f);
  249. }
  250. PRIVATE double floor(f)
  251. double f;
  252. {
  253. double res, d;
  254. register int sign = 1;
  255. /* eliminate the sign */
  256. if (f < 0) {
  257. sign = -1, f = -f;
  258. }
  259. /* get the largest power of 2 <= f */
  260. d = 1.0;
  261. while (f - d >= d) {
  262. d *= 2.0;
  263. }
  264. /* reconstruct f by deminishing powers of 2 */
  265. res = 0.0;
  266. while (d >= 1.0) {
  267. if (res + d <= f)
  268. res += d;
  269. d /= 2.0;
  270. }
  271. /* undo the sign elimination */
  272. if (sign == -1) {
  273. res = -res, f = -f;
  274. if (res > f)
  275. res -= 1.0;
  276. }
  277. return res;
  278. }
  279. PRIVATE double ttttp(f, n) /* times ten to the power */
  280. double f;
  281. {
  282. while (n > 0) {
  283. f = mlf(f, 10.0);
  284. n--;
  285. }
  286. while (n < 0) {
  287. f = dvf(f, 10.0);
  288. n++;
  289. }
  290. return f;
  291. }
  292. /* Str2double is used to initialize the global data area with floats;
  293. we do not use, e.g., sscanf(), to be able to check the grammar of
  294. the string and to give warnings.
  295. */
  296. double str2double(str)
  297. char *str;
  298. {
  299. register char b;
  300. register int sign = 1; /* either +1 or -1 */
  301. register int frac = 0; /* how far in fraction part ? */
  302. register int ex; /* to store exponent */
  303. double mantissa = 0.0; /* to store mantissa */
  304. double d; /* double to be returned */
  305. b = *str++;
  306. if (b == '-') {
  307. sign = -1;
  308. b = *str++;
  309. }
  310. else if (b == '+') {
  311. sign = 1;
  312. b = *str++;
  313. }
  314. if ('0' <= b && b <= '9') {
  315. mantissa = (double) (b-'0');
  316. }
  317. else if (b == '.') {
  318. /* part before dot cannot be empty */
  319. warning(WBADFLOAT);
  320. frac = 1;
  321. }
  322. else {
  323. goto BadFloat;
  324. }
  325. LOG((" q9 str2double : (before while) mantissa = %20.20g", mantissa));
  326. while ((b = *str++) != 'e' && b != 'E' && b != '\0') {
  327. if (b == '.') {
  328. if (frac == 0) {
  329. frac++;
  330. }
  331. else { /* there already was a '.' in input */
  332. goto BadFloat;
  333. }
  334. }
  335. else if ('0' <= b && b <= '9') {
  336. double bval = b - '0';
  337. if (frac) {
  338. mantissa =
  339. adf(mantissa, ttttp(bval, -frac));
  340. frac++;
  341. }
  342. else {
  343. mantissa =
  344. adf(mlf(mantissa, 10.0), bval);
  345. }
  346. }
  347. else {
  348. goto BadFloat;
  349. }
  350. LOG((" q9 str2double : (inside while) mantissa = %20.20g",
  351. mantissa));
  352. }
  353. LOG((" q9 str2double : mantissa = %10.10g", mantissa));
  354. mantissa = sign * mantissa;
  355. if (b == '\0')
  356. return (mantissa);
  357. /* else we have b == 'e' or b== 'E' */
  358. /* Optional sign for exponent */
  359. b = *str++;
  360. if (b == '-') {
  361. sign = -1;
  362. b = *str++;
  363. }
  364. else if (b == '+') {
  365. sign = 1;
  366. b = *str++;
  367. }
  368. else {
  369. sign = 1;
  370. }
  371. ex = 0;
  372. do {
  373. if ('0' <= b && b <= '9') {
  374. ex = 10*ex + (b-'0');
  375. }
  376. else {
  377. goto BadFloat;
  378. }
  379. } while ((b = *str++) != '\0');
  380. LOG((" q9 str2double : exponent = %d", ex));
  381. /* Construct total value of float */
  382. ex = sign * ex;
  383. d = ttttp(mantissa, ex);
  384. return (d);
  385. BadFloat:
  386. fatal("Float garbled in loadfile");
  387. return (0.0);
  388. }
  389. #else /* NOFLOAT */
  390. nofloat() {
  391. fatal("attempt to execute a floating point instruction on an EM machine without FP");
  392. }
  393. #endif /* NOFLOAT */