do_fpar.c 7.4 KB

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