451 lines
7.4 KiB
C
451 lines
7.4 KiB
C
/*
|
|
* Sources of the "FLOATING POINT ARITHMETIC" group instructions
|
|
*/
|
|
|
|
/* $Header$ */
|
|
|
|
#include <em_abs.h>
|
|
#include "nofloat.h"
|
|
#include "global.h"
|
|
#include "log.h"
|
|
#include "mem.h"
|
|
#include "trap.h"
|
|
#include "text.h"
|
|
#include "fra.h"
|
|
#include "warn.h"
|
|
|
|
#ifndef NOFLOAT
|
|
|
|
extern double fpop();
|
|
|
|
#define MAXDOUBLE 99.e999 /* IEEE infinity */ /*???*/
|
|
#define SMALL (1.0/MAXDOUBLE)
|
|
|
|
PRIVATE double adf(), sbf(), mlf(), dvf();
|
|
PRIVATE double ttttp();
|
|
PRIVATE double floor(), fabs();
|
|
PRIVATE fef(), fif();
|
|
|
|
#endif NOFLOAT
|
|
|
|
DoADF(l)
|
|
register size l;
|
|
{
|
|
/* ADF w: Floating add (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoADF(%ld)", l));
|
|
spoilFRA();
|
|
fpush(adf(fpop(l), t), l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoSBF(l)
|
|
register size l;
|
|
{
|
|
/* SBF w: Floating subtract (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoSBF(%ld)", l));
|
|
spoilFRA();
|
|
fpush(sbf(fpop(l), t), l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoMLF(l)
|
|
register size l;
|
|
{
|
|
/* MLF w: Floating multiply (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoMLF(%ld)", l));
|
|
spoilFRA();
|
|
fpush(mlf(fpop(l), t), l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoDVF(l)
|
|
register size l;
|
|
{
|
|
/* DVF w: Floating divide (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoDVF(%ld)", l));
|
|
spoilFRA();
|
|
fpush(dvf(fpop(l), t), l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoNGF(l)
|
|
register size l;
|
|
{
|
|
/* NGF w: Floating negate (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoNGF(%ld)", l));
|
|
spoilFRA();
|
|
fpush(-t, l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoFIF(l)
|
|
register size l;
|
|
{
|
|
/* FIF w: Floating multiply and split integer and fraction part (*) */
|
|
#ifndef NOFLOAT
|
|
double t = fpop(arg_wf(l));
|
|
|
|
LOG(("@F6 DoFIF(%ld)", l));
|
|
spoilFRA();
|
|
fif(fpop(l), t, l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
DoFEF(l)
|
|
register size l;
|
|
{
|
|
/* FEF w: Split floating number in exponent and fraction part (*) */
|
|
#ifndef NOFLOAT
|
|
LOG(("@F6 DoFEF(%ld)", l));
|
|
spoilFRA();
|
|
fef(fpop(arg_wf(l)), l);
|
|
#else NOFLOAT
|
|
nofloat();
|
|
#endif NOFLOAT
|
|
}
|
|
|
|
#ifndef NOFLOAT
|
|
|
|
/* Service routines */
|
|
|
|
PRIVATE double adf(f1, f2) /* returns f1 + f2 */
|
|
double f1, f2;
|
|
{
|
|
if (must_test && !(IgnMask&BIT(EFOVFL))) {
|
|
if (f1 > 0.0 && f2 > 0.0) {
|
|
if (MAXDOUBLE - f1 < f2) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
else if (f1 < 0.0 && f2 < 0.0) {
|
|
if (-(MAXDOUBLE + f1) > f2) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
}
|
|
return (f1 + f2);
|
|
}
|
|
|
|
PRIVATE double sbf(f1, f2) /* returns f1 - f2 */
|
|
double f1, f2;
|
|
{
|
|
if (must_test && !(IgnMask&BIT(EFOVFL))) {
|
|
if (f2 < 0.0 && f1 > 0.0) {
|
|
if (MAXDOUBLE - f1 < -f2) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
else if (f2 > 0.0 && f1 < 0.0) {
|
|
if (f2 - MAXDOUBLE > f1) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
}
|
|
return (f1 - f2);
|
|
}
|
|
|
|
PRIVATE double mlf(f1, f2) /* returns f1 * f2 */
|
|
double f1, f2;
|
|
{
|
|
double ff1 = fabs(f1), ff2 = fabs(f2);
|
|
|
|
if (f1 == 0.0 || f2 == 0.0)
|
|
return (0.0);
|
|
|
|
if ((ff1 >= 1.0 && ff2 <= 1.0) || (ff2 >= 1.0 && ff1 <= 1.0))
|
|
return (f1 * f2);
|
|
|
|
if (must_test && !(IgnMask&BIT(EFUNFL))) {
|
|
if (ff1 < 1.0 && ff2 < 1.0) {
|
|
if (SMALL / ff1 > ff2) {
|
|
trap(EFUNFL);
|
|
return (0.0);
|
|
}
|
|
return (f1 * f2);
|
|
}
|
|
}
|
|
|
|
if (must_test && !(IgnMask&BIT(EFOVFL))) {
|
|
if (MAXDOUBLE / ff1 < ff2) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
return (f1 * f2);
|
|
}
|
|
|
|
PRIVATE double dvf(f1, f2) /* returns f1 / f2 */
|
|
double f1, f2;
|
|
{
|
|
double ff1 = fabs(f1), ff2 = fabs(f2);
|
|
|
|
if (f2 == 0.0) {
|
|
if (!(IgnMask&BIT(EFDIVZ))) {
|
|
trap(EFDIVZ);
|
|
}
|
|
else return (0.0);
|
|
}
|
|
|
|
if (f1 == 0.0)
|
|
return (0.0);
|
|
|
|
if ((ff2 >= 1.0 && ff1 >= 1.0) || (ff1 <= 1.0 && ff2 <= 1.0))
|
|
return (f1 / f2);
|
|
|
|
if (must_test && !(IgnMask&BIT(EFUNFL))) {
|
|
if (ff2 > 1.0 && ff1 < 1.0) {
|
|
if (SMALL / ff2 > ff1) {
|
|
trap(EFUNFL);
|
|
return (0.0);
|
|
}
|
|
return (f1 / f2);
|
|
}
|
|
}
|
|
|
|
if (must_test && !(IgnMask&BIT(EFOVFL))) {
|
|
if (MAXDOUBLE * ff2 < ff1) {
|
|
trap(EFOVFL);
|
|
return (0.0);
|
|
}
|
|
}
|
|
return (f1 / f2);
|
|
}
|
|
|
|
PRIVATE fif(f1, f2, n)
|
|
double f1, f2;
|
|
size n;
|
|
{
|
|
double f = mlf(f1, f2);
|
|
double fl = floor(fabs(f));
|
|
|
|
fpush(fabs(f) - fl, n); /* push fraction */
|
|
fpush((f < 0.0) ? -fl : fl, n); /* push integer-part */
|
|
}
|
|
|
|
PRIVATE fef(f, n)
|
|
double f;
|
|
size n;
|
|
{
|
|
register long exponent, sign = (long) (f < 0.0);
|
|
|
|
if (f == 0.0) {
|
|
fpush(f, n);
|
|
wpush(0L);
|
|
return;
|
|
}
|
|
|
|
for (f = fabs(f), exponent = 0; f >= 1.0; exponent++)
|
|
f /= 2.0;
|
|
|
|
for (; f < 0.5; exponent--)
|
|
f *= 2.0;
|
|
|
|
fpush((sign) ? -f : f, n); /* push mantissa */
|
|
wpush(exponent); /* push exponent */
|
|
}
|
|
|
|
/* floating point service routines, to avoid having to use -lm */
|
|
|
|
PRIVATE double fabs(f)
|
|
double f;
|
|
{
|
|
return (f < 0.0 ? -f : f);
|
|
}
|
|
|
|
PRIVATE double floor(f)
|
|
double f;
|
|
{
|
|
double res, d;
|
|
register int sign = 1;
|
|
|
|
/* eliminate the sign */
|
|
if (f < 0) {
|
|
sign = -1, f = -f;
|
|
}
|
|
|
|
/* get the largest power of 2 <= f */
|
|
d = 1.0;
|
|
while (f - d >= d) {
|
|
d *= 2.0;
|
|
}
|
|
|
|
/* reconstruct f by deminishing powers of 2 */
|
|
res = 0.0;
|
|
while (d >= 1.0) {
|
|
if (res + d <= f)
|
|
res += d;
|
|
d /= 2.0;
|
|
}
|
|
|
|
/* undo the sign elimination */
|
|
if (sign == -1) {
|
|
res = -res, f = -f;
|
|
if (res > f)
|
|
res -= 1.0;
|
|
}
|
|
|
|
return res;
|
|
}
|
|
|
|
PRIVATE double ttttp(f, n) /* times ten to the power */
|
|
double f;
|
|
{
|
|
while (n > 0) {
|
|
f = mlf(f, 10.0);
|
|
n--;
|
|
}
|
|
while (n < 0) {
|
|
f = dvf(f, 10.0);
|
|
n++;
|
|
}
|
|
return f;
|
|
}
|
|
|
|
/* Str2double is used to initialize the global data area with floats;
|
|
we do not use, e.g., sscanf(), to be able to check the grammar of
|
|
the string and to give warnings.
|
|
*/
|
|
|
|
double str2double(str)
|
|
char *str;
|
|
{
|
|
register char b;
|
|
register int sign = 1; /* either +1 or -1 */
|
|
register int frac = 0; /* how far in fraction part ? */
|
|
register int ex; /* to store exponent */
|
|
double mantissa = 0.0; /* to store mantissa */
|
|
double d; /* double to be returned */
|
|
|
|
b = *str++;
|
|
if (b == '-') {
|
|
sign = -1;
|
|
b = *str++;
|
|
}
|
|
else if (b == '+') {
|
|
sign = 1;
|
|
b = *str++;
|
|
}
|
|
|
|
if ('0' <= b && b <= '9') {
|
|
mantissa = (double) (b-'0');
|
|
}
|
|
else if (b == '.') {
|
|
/* part before dot cannot be empty */
|
|
warning(WBADFLOAT);
|
|
frac = 1;
|
|
}
|
|
else {
|
|
goto BadFloat;
|
|
}
|
|
|
|
LOG((" q9 str2double : (before while) mantissa = %20.20g", mantissa));
|
|
|
|
while ((b = *str++) != 'e' && b != 'E' && b != '\0') {
|
|
if (b == '.') {
|
|
if (frac == 0) {
|
|
frac++;
|
|
}
|
|
else { /* there already was a '.' in input */
|
|
goto BadFloat;
|
|
}
|
|
}
|
|
else if ('0' <= b && b <= '9') {
|
|
double bval = b - '0';
|
|
|
|
if (frac) {
|
|
mantissa =
|
|
adf(mantissa, ttttp(bval, -frac));
|
|
frac++;
|
|
}
|
|
else {
|
|
mantissa =
|
|
adf(mlf(mantissa, 10.0), bval);
|
|
}
|
|
}
|
|
else {
|
|
goto BadFloat;
|
|
}
|
|
LOG((" q9 str2double : (inside while) mantissa = %20.20g",
|
|
mantissa));
|
|
}
|
|
LOG((" q9 str2double : mantissa = %10.10g", mantissa));
|
|
mantissa = sign * mantissa;
|
|
if (b == '\0')
|
|
return (mantissa);
|
|
/* else we have b == 'e' or b== 'E' */
|
|
|
|
/* Optional sign for exponent */
|
|
b = *str++;
|
|
if (b == '-') {
|
|
sign = -1;
|
|
b = *str++;
|
|
}
|
|
else if (b == '+') {
|
|
sign = 1;
|
|
b = *str++;
|
|
}
|
|
else {
|
|
sign = 1;
|
|
}
|
|
|
|
ex = 0;
|
|
do {
|
|
if ('0' <= b && b <= '9') {
|
|
ex = 10*ex + (b-'0');
|
|
}
|
|
else {
|
|
goto BadFloat;
|
|
}
|
|
} while ((b = *str++) != '\0');
|
|
LOG((" q9 str2double : exponent = %d", ex));
|
|
|
|
/* Construct total value of float */
|
|
ex = sign * ex;
|
|
d = ttttp(mantissa, ex);
|
|
return (d);
|
|
|
|
BadFloat:
|
|
fatal("Float garbled in loadfile");
|
|
return (0.0);
|
|
}
|
|
|
|
#else NOFLOAT
|
|
|
|
nofloat() {
|
|
fatal("attempt to execute a floating point instruction on an EM machine without FP");
|
|
}
|
|
|
|
#endif NOFLOAT
|
|
|