ack/util/int/do_fpar.c
George Koehler 777d0abb00 Add util/int/switch.h to declare do_instr() and Do???()
Use switch.h to reduce warnings from clang about implicit declarations
of functions.  I used `grep ... do_*.c | sed ... | sort`, and some
manual editing, to make the big list of Do???() functions.
2019-10-30 18:36:55 -04:00

448 lines
7.9 KiB
C

/** @file
* Sources of the "FLOATING POINT ARITHMETIC" group instructions
*/
/* $Id$ */
#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 "io.h"
#include "warn.h"
#include "switch.h"
#ifndef NOFLOAT
extern double fpop();
#ifdef __STDC__
#include <float.h>
#define MAXDOUBLE DBL_MAX
#else /* not __STDC__ */
#if defined(vax) || defined(pdp) || defined(__vax) || defined(__pdp)
#define MAXDOUBLE 1.701411834604692293e+38
#else
#define MAXDOUBLE 1.7976931348623157e+308
#endif
#endif /* not __STDC__ */
#define SMALL (1.0/MAXDOUBLE)
PRIVATE double adf(double, double), sbf(double, double), mlf(double, double), dvf(double, double);
PRIVATE double ttttp(double, int);
PRIVATE double floor(double), fabs(double);
PRIVATE void fef(double, size), fif(double, double, size);
#endif /* NOFLOAT */
void DoADF(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 */
}
void DoSBF(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 */
}
void DoMLF(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 */
}
void DoDVF(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 */
}
void DoNGF(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 */
}
void DoFIF(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 */
}
void DoFEF(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 */
/** Returns "f1" + "f2" */
PRIVATE double adf(double f1, double 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);
}
/** Returns "f1" - "f2" */
PRIVATE double sbf(double f1, double 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);
}
/** Returns "f1" * "f2" */
PRIVATE double mlf(double f1, double 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);
}
/** Returns "f1" / "f2" */
PRIVATE double dvf(double f1, double 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 void fif(double f1, double 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 void fef(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(double f)
{
return (f < 0.0 ? -f : f);
}
PRIVATE double floor(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;
}
/** Times ten to the power. */
PRIVATE double ttttp(double f, int n)
{
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(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 */
void nofloat(void)
{
fatal("attempt to execute a floating point instruction on an EM machine without FP");
}
#endif /* NOFLOAT */