koreth@ssyx.ucsc.edu (Steven Grimm) (03/20/89)
Submitted-by: brooks@csss-a.prime.com (David Brooks)
Posting-number: Volume 2, Issue 33
Archive-name: fplib10
[These are the sources to some floating-point library routines for Sozobon C.
Apparently this is not the complete floating point library; the whole thing,
in binary-only form, has been posted to the binaries group. -sg]
#!/bin/sh
# shar: Shell Archiver (v1.22)
#
# Run the following text with /bin/sh to create:
# ATOF.C
# EXP.C
# FABS.C
# FMOD.C
# FPCMP.S
# ITRIG.C
# LOG.C
# MATH.H
# POW.C
# README
# SQRT.C
# TRIG.C
#
sed 's/^X//' << 'SHAR_EOF' > ATOF.C &&
X/* ATOF function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X/* Hack warning. Hack warning. Hack warning... */
X/* Somehow I feel scanf should call atof, not the other way around. */
X/* Guys??? */
X
Xfloat atof(s)
Xchar *s;
X{ float res;
X
X sscanf(s, "%g", &res);
X return res;
X}
X
X_scandum()
X{
X _scanf(); /* Force floating version to load */
X}
SHAR_EOF
chmod 0600 ATOF.C || echo "restore of ATOF.C fails"
sed 's/^X//' << 'SHAR_EOF' > EXP.C &&
X/* EXP and hyperbolic trig functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <errno.h>
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* EXP: */
X/* */
X/* Uses the synthesis (2**n) * (2**[m/8]) ** b**g */
X/* where b = 2**(1/8) and 8*n+m+g = arg/ln(b). b**g has a cubic */
X/* approximation: source and accuracy unknown. */
X/* */
X/* Beware a Sozobon C bug: can't compare two negative floating numbers. */
X
Xfloat exp(a) register float a;
X{ fstruct res;
X register int aint;
X
X static fstruct huge = 0xFFFFFF7FL;
X static float powtab[] = {1.0, 1.09050773, 1.1892071, 1.2968396,
X 1.41421356, 1.5422108, 1.68179283, 1.8340081};
X
X if (a > 43.5)
X { errno = EDOM;
X return huge.f;
X }
X
X if ((a + 43.5) < 0.0)
X return 0.0; /* Underflow */
X
X if ((aint = (int)(a *= 11.5415603)) < 0) /* 8/ln(2) */
X --aint; /* Correct mathematically */
X a = (float)aint - a; /* -(frac part) */
X res.f = 1.0 + a * (-0.0866439378 + a * (0.003750577 +
X a * -0.11321783e-3));
X res.sc[3] += aint >> 3; /* Mult by 2**n */
X return res.f * powtab[aint&7]; /* Mult by 2**[m/8] */
X}
X
X/* HYPERBOLIC FUNCTIONS: virtually free */
X
Xfloat sinh(a) float a;
X{ float exp();
X register float ea = exp(a);
X
X return 0.5 * (ea - (1.0 / ea));
X}
X
Xfloat cosh(a) float a;
X{ float exp();
X register float ea = exp(a);
X
X return 0.5 * (ea + (1.0 / ea));
X}
X
Xfloat tanh(a) float a;
X{ float exp();
X register float e2a;
X
X e2a = exp(a);
X e2a = e2a * e2a; /* exp-squared-a */
X return (e2a - 1.0) / (e2a + 1.0);
X}
X
SHAR_EOF
chmod 0600 EXP.C || echo "restore of EXP.C fails"
sed 's/^X//' << 'SHAR_EOF' > FABS.C &&
X/* FABS, LDEXP, and FREXP functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <errno.h>
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* FABS: directly manipulate the sign bit, punning the type. Actually */
X/* the argument and return are both float. */
X
Xunsigned long fabs(a)
Xunsigned long a;
X{ return a & 0xFFFFFF7FL; }
X
X/* LDEXP: again direct manipulation. */
X
Xfloat ldexp(a, n)
Xfstruct a;
Xint n;
X{ int exp = (a.sc[3] & 0x7F) + n;
X
X if (exp < 0)
X return 0L; /* Underflow */
X
X if (exp > 0x7F)
X { errno = EDOM; /* Overflow */
X return a.ul | 0xFFFFFF7FL; /* Preserve sign */
X }
X
X a.sc[3] = (a.sc[3] & 0x80) | exp;
X return a.f;
X}
X
X/* FREXP: split into convenient mantissa and exponent */
X
Xfloat frexp(a, ip)
Xfstruct a;
Xint *ip;
X{ if (a.ul == 0L)
X *ip = 0;
X else
X { *ip = (a.sc[3] & 0x7f) - 0x40;
X a.sc[3] = (a.sc[3] & 0x80) | 0x40; /* preserve sign,
X force exponent */
X }
X
X return a.f;
X}
SHAR_EOF
chmod 0600 FABS.C || echo "restore of FABS.C fails"
sed 's/^X//' << 'SHAR_EOF' > FMOD.C &&
X/* CEIL, FLOOR, MODF, and FMOD functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* This macro returns true if the number is so large that is has no */
X/* representable fraction. */
X
X#define ISLGINT(a) (((a).sc[3] & 0x7F) > (0x40 + 24))
X
X/* CEIL: smallest integral value not less than a. */
X
Xfloat ceil(a)
Xfstruct a;
X{ register float aint;
X
X if (ISLGINT(a))
X return a.f;
X
X if ((aint = (float)((long)a.f)) != a.f && a.sc[3] >= 0)
X ++aint;
X
X return aint;
X}
X
X/* FLOOR: largest integral value not greater than a. */
X
Xfloat floor(a)
Xfstruct a;
X{ register float aint;
X
X if (ISLGINT(a))
X return a.f;
X
X if ((aint = (float)((long)a.f)) != a.f && a.sc[3] < 0)
X --aint;
X
X return aint;
X}
X
X/* Just for once, we luck out: negative numbers behave as needed... */
X
Xfloat modf(a, ip)
Xfstruct a;
Xfloat *ip;
X{ long ipart;
X
X *ip = (ISLGINT(a))?a.f:(float)(long)a.f;
X
X return a.f - *ip;
X}
X
Xfloat fmod(n, d)
Xfloat n, d;
X{ fstruct quot;
X long iquot;
X
X quot.f = n / d;
X if (ISLGINT(quot))
X return 0.0;
X
X iquot = (long)quot.f;
X return n - (float)iquot * d;
X}
SHAR_EOF
chmod 0600 FMOD.C || echo "restore of FMOD.C fails"
sed 's/^X//' << 'SHAR_EOF' > FPCMP.S &&
X* Copyright (c) 1988 by Sozobon, Limited. Author: Johann Ruegg
X*
X* Permission is granted to anyone to use this software for any purpose
X* on any computer system, and to redistribute it freely, with the
X* following restrictions:
X* 1) No charge may be made other than reasonable charges for reproduction.
X* 2) Modified versions must be clearly marked as such.
X* 3) The authors are not responsible for any harmful consequences
X* of using this software, even if they result from defects in it.
X*
X* MODIFICATIONS:
X*
X* D.W.Brooks 3/5/89 Fixed fpcmp (when both -ve) and fpneg (when zero)
X*
X
X .globl fpcmp
X .globl _fpcmp
Xfpcmp:
X_fpcmp:
X move.l 4(sp),d0
X move.l 8(sp),d1
X* Start DWB
X move.b d0,d2 * See if both signs are negative
X and.b d1,d2
X bpl L0
X move.l d0,d2 * They are: switch arguments so the rest works
X move.l d1,d0
X move.l d2,d1
XL0:
X* End DWB
X cmp.b d1,d0
X bne L1
X cmp.l d1,d0
XL1:
X rts
X
X .globl fpneg
X .globl _fpneg
Xfpneg:
X_fpneg:
X move.l 4(sp),d0
X
X* DWB removed: tst.b d0
X* This would have not negated very small numbers.
X
X beq L2
X eor.b #$80,d0
XL2:
X rts
SHAR_EOF
chmod 0600 FPCMP.S || echo "restore of FPCMP.S fails"
sed 's/^X//' << 'SHAR_EOF' > ITRIG.C &&
X/* ASIN, ACOS, ATAN, and ATAN2 functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <errno.h>
X
X#define PI 3.14159265
X#define PIBY2 1.57079633
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* ARCTAN: */
X/* */
X/* Single precision is 24 bits, or better than 7 decimal digits. The */
X/* method is from Hart et al: "Computer Approximations". It uses a */
X/* divide-and-conquer algorithm. The argument range is divided into */
X/* five using the partition points tan({1,3,5,7}*pi/16), a polynomial */
X/* valid over +/- pi/16 is calculated, and other magic is used to */
X/* reposition the result. Precision is >8 places. */
X
Xfloat atan(a) fstruct a;
X{ fstruct absval;
X register float tx0, tsq, temp;
X register int part;
X register char sign;
X
X static float adj[] = {0.0, 0.414213562, 1.0, 2.414213562};
X static float atof[] = {0.0, 0.39269908, 0.78539816,
X 1.178097245, 1.57079633};
X
X sign = a.sc[3] & 0x80;
X a.sc[3] &= 0x7f; /* get fabs(a) */
X tx0 = a.f;
X
X/* Figure out the partition */
X
X part = tx0<0.66817864?(tx0<0.19891237?0:1):(tx0<1.49660576?2:
X tx0<5.0273395?3:4);
X if (part == 4)
X tx0 = -1.0 / tx0;
X else if (part != 0)
X { temp = adj[part];
X tx0 = (tx0 - temp) / (tx0 * temp + 1.0);
X }
X
X/* Here's the calculation */
X
X tsq = tx0 * tx0 + 1.67784279;
X a.f = (0.93833093 / tsq + 0.44075154) * tx0 + atof[part];
X
X a.sc[3] ^= sign; /* Negate if negative */
X return a.f;
X}
X
X/* ARCSIN and ARCCOS use the standard identities. There's some less */
X/* optimal code here because Sozobon C can't properly compare two */
X/* negative floating numbers. */
X
Xfloat asin(a) float a;
X{ float atan(), sqrt();
X
X if (a > 1.0 || (a+1.0) < 0.0)
X { errno = EDOM;
X return 0.0;
X }
X
X return atan(a / sqrt(1.0 - a*a));
X}
X
Xfloat acos(a) float a;
X{ float atan(), sqrt();
X float temp;
X
X if (a > 1.0 || (a+1.0) < 0.0)
X { errno = EDOM;
X return 0.0;
X }
X
X temp = atan(sqrt(1.0 - a*a) / a);
X if (a >= 0.0)
X return temp;
X else
X return PI + temp;
X}
X
X/* ATAN2: */
X/* */
X/* Computes atan(quotient), and returns that for positive cos, else */
X/* extends the range depending on the sin. */
X
Xfloat atan2(s, c) fstruct s, c;
X{ register float r;
X float atan();
X
X if (c.sc[3] == 0L) /* Infinite argument */
X return (s.sc[3]<0)?-PIBY2:PIBY2;
X
X r = atan(s.f / c.f);
X
X if (c.sc[3] >= 0)
X return r; /* Range -pi/2..pi/2 */
X if (s.sc[3] >= 0)
X return PI + r;
X return -PI + r;
X}
SHAR_EOF
chmod 0600 ITRIG.C || echo "restore of ITRIG.C fails"
sed 's/^X//' << 'SHAR_EOF' > LOG.C &&
X/* LOG functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <errno.h>
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* LOG: */
X/* */
X/* This yanks out the exponent, performs a polynomial valid over +/- */
X/* sqrt(2), and adds the scaled exponent back in. Domain error on */
X/* arguments <= 0. */
X
Xfloat log(a) fstruct a;
X{ register float t1, t2, t3;
X static fstruct nhuge = 0xFFFFFFFFL; /* Largest negative */
X
X if (a.f <= 0.0)
X { errno = EDOM; /* Impossible really */
X return nhuge.f;
X }
X
X if (a.f < 1.41421356 && a.f > 0.70710678)
X { t3 = 0.0; /* arg between sqrt(2) and sqrt(2)/2 */
X t2 = (t1 = a.f - 1.0) + 2.0;
X }
X else
X { t3 = (float)(((int)a.sc[3] << 1) - 0x81) * 0.34657359;
X /* ...(log2(a) * 2 + 1) * ln(2)/2 */
X a.sc[3] = 0x40; /* Scale back to range */
X t2 = a.f + 0.70710678;
X t1 = a.f - 0.70710678;
X }
X
X t1 /= t2;
X t2 = t1 * t1;
X
X return t3 + t1 * (2.0000008 + t2 * (0.66644078 + t2 * 0.41517739));
X}
X
X/* LOG10 is easy. */
X
Xfloat log10(a) float a;
X{ float log();
X return 0.43429448 * log(a); /* ln(a) * log10(e) */
X}
SHAR_EOF
chmod 0600 LOG.C || echo "restore of LOG.C fails"
sed 's/^X//' << 'SHAR_EOF' > MATH.H &&
X/*
X * MATH.H Math function declarations
X */
X
X#ifndef MATH_H
X#define MATH_H
X
Xextern float sin();
Xextern float cos();
Xextern float tan();
Xextern float asin();
Xextern float acos();
Xextern float atan();
Xextern float atan2();
Xextern float sinh();
Xextern float cosh();
Xextern float tanh();
Xextern float exp();
Xextern float log();
Xextern float log10();
Xextern float pow();
Xextern float sqrt();
Xextern float ceil();
Xextern float floor();
Xextern float fabs();
Xextern float ldexp();
Xextern float frexp();
Xextern float modf();
Xextern float fmod();
Xextern float atof();
X
X#endif MATH_H
SHAR_EOF
chmod 0600 MATH.H || echo "restore of MATH.H fails"
sed 's/^X//' << 'SHAR_EOF' > POW.C &&
X/* POW function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X/* */
X/* This uses exp and log in the obvious way. A domain error occurs if */
X/* a=0 and p<=0, or a<0 and p is not integral. */
X
X#include <errno.h>
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
Xfloat pow(a, p)
Xfstruct a;
Xregister float p;
X{ fstruct r;
X float exp(), log();
X register long pint;
X register int sign = 0;
X
X if (a.ul == 0L && p <= 0.0)
X { errno = EDOM;
X return 0.0;
X }
X
X if (a.sc[3] < 0) /* ...if a.f < 0 */
X { pint = (long)p;
X if ((float)pint != p) /* ... if nonintegral */
X errno = EDOM;
X sign = pint & 1; /* ... carry on anyway */
X a.sc[3] &= 0x7f;
X }
X
X r.f = exp(log(a.f) * p);
X if (sign)
X r.sc[3] |= 0x80;
X return r.f;
X}
SHAR_EOF
chmod 0600 POW.C || echo "restore of POW.C fails"
sed 's/^X//' << 'SHAR_EOF' > README &&
X FLOATING POINT LIBRARY FOR SOZOBON C
X
X Rev 1.0 -- prepared March 6, 1989
X
XThe files contained in this archive are an initial release of some
Xfloating point functions for Sozobon C. We provide:
X
X sqrt exp fabs
X sin log ldexp
X cos log10 frexp
X tan pow atof
X atan ceil
X atan2 floor math.h
X asin modf
X acos fmod
X
XThese are available in source. There is also a new LIBM.A that has the
Xnew routines included (and the other routines re-arranged to remove the
Xneed to scan the library twice).
X
XThese routines may be used without restriction, but I retain the copyright
Xon them, and the Only True Upgrades will be those I make. If you have any
Ximprovements, please send them to me.
X
XDavid W. Brooks
X36 Elm Street
XMedfield, MA 02052
X
XBROOKS@CSSS-A.PRIME.COM
X{uunet,mit-eddie}!csss-a.prime.com!brooks
SHAR_EOF
chmod 0600 README || echo "restore of README fails"
sed 's/^X//' << 'SHAR_EOF' > SQRT.C &&
X/* SQRT function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X/* */
X/* Single precision is 24 bits, or better than 7 decimal digits. The */
X/* method is from Hart et al: "Computer Approximations". We reduce */
X/* the range to [0.25..1.0], then apply the polynomial: */
X/* 0.25927688 + 1.0520212n - 0.31632214n^2, giving 2.3 digits */
X/* approximation. Then we do two Newton's iterations (giving ~9 digits */
X/* precision. */
X/* */
X/* A negative argument gives a return value of 0.0 and sets EDOM. */
X
X#include <errno.h>
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* Sadly we can't use register for objects of type fstruct. */
X
Xfloat sqrt(n)
Xfstruct n;
X{ register int exp; /* Separate unbiased exponent. */
X fstruct s;
X
X if ((n.ul & 0xFFFFFF00L) == 0L) /* Check for unnormalized 0. */
X return 0.0;
X
X if ((exp = n.sc[3]) < 0) /* Test argument's sign. */
X { errno = EDOM;
X return 0.0;
X }
X
X exp -= 0x40; /* Get unbiased exponent. */
X if (exp & 1)
X { n.sc[3] = 0x3F; /* Convert to 0.25..0.5. */
X ++exp;
X }
X else
X n.sc[3] = 0x40; /* Or convert to 0.5..1. */
X
X s.f = 0.25927688 + n.f * (1.0520212 + n.f * -0.31632214);
X
X s.f = s.f + n.f / s.f; /* One Newton. */
X s.sc[3] -= 1; /* (here's the divide by 2). */
X s.f = s.f + n.f / s.f; /* The other Newton. */
X s.sc[3] += ((unsigned int)exp >> 1) - 1;
X /* Divide by 2 and insert
X half the original exponent. */
X return s.f;
X}
SHAR_EOF
chmod 0600 SQRT.C || echo "restore of SQRT.C fails"
sed 's/^X//' << 'SHAR_EOF' > TRIG.C &&
X/* SIN, COS, and TAN functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <errno.h>
X
X#define PIBY2 1.57079633
X#define TWOBYPI 0.636619772
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
X
X/* SIN:
X/* */
X/* Single precision is 24 bits, or better than 7 decimal digits. The */
X/* method is from Hart et al: "Computer Approximations". It reduces */
X/* the range to 0..pi/2, scaled to 0..1, and applies: */
X/* 1.57079629t - 0.64596336t^3 + 0.0796884805t^5 - 0.0046722279t^7 */
X/* + 0.00015082056t^9 (precision of 8.5 digits) and */
X/* finally fixes the sign. */
X
Xfloat sin(a) fstruct a;
X{ register float t, ipart, tsq;
X fstruct x;
X register unsigned long tint;
X register char sign;
X
X sign = a.sc[3] & 0x80;
X a.sc[3] ^= sign; /* sin(negative x) = -sin(-x) */
X
X/* If the argument is huge, the result is rather arbitrary -- and, */
X/* besides, the following code will break. */
X
X if (a.sc[3] >= (64+24))
X { errno = EDOM;
X return 0.0;
X }
X
X t = a.f * TWOBYPI; /* Express in quarter-turns */
X
X/* Get integer part and fraction. Optimise the common case, in the
X first quadrant */
X
X if (t <= 1.0)
X tint = 0L;
X else
X { tint = t; /* tint = whole quarterturns */
X ipart = (float)tint;
X t = (tint & 1)?ipart - t + 1.0:t - ipart; /* Get fraction */
X }
X tsq = t * t;
X x.f = t * (1.57079629 + tsq * (-0.64596336 + tsq * (0.0796884805
X + tsq * (-0.46722279e-2 + tsq * 0.15082056e-3))));
X
X if ((int)tint & 2) /* Quadrants 3 and 4 */
X x.sc[3] |= 0x80; /* Set negative */
X if (sign != 0)
X x.sc[3] ^= 0x80; /* Negate */
X
X return x.f;
X}
X
X/* COS: */
X/* */
X/* Although approximations are available, this will do fine. */
X
Xfloat cos(a) float a;
X{ float sin();
X return sin(PIBY2 - a);
X}
X
X/* TAN: */
X/* (v-e-r-y s-l-o-w-l-y) */
X
Xfloat tan(a) float a;
X{ float sin();
X return sin(a) / sin(PIBY2 - a);
X}
SHAR_EOF
chmod 0600 TRIG.C || echo "restore of TRIG.C fails"
exit 0