bammi@dsrgsun.ces.cwru.edu (Jwahar R. Bammi) (12/14/88)
#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# dflonum.c
# frexp.c
# ldexp.c
# modf.c
# This archive created: Wed Dec 14 04:30:47 1988
# By: Jwahar R. Bammi(Case Western Reserve University)
# Uucp: {decvax,sun,att}!cwjcc!dsrgsun!bammi
# Csnet: bammi@dsrgsun.ces.CWRU.edu
# Arpa: bammi@dsrgsun.ces.CWRU.edu
#
export PATH; PATH=/bin:$PATH
echo shar: extracting "'dflonum.c'" '(18187 characters)'
if test -f 'dflonum.c'
then
echo shar: over-writing existing file "'dflonum.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'dflonum.c'
X/* A working Gnulib68k, thanks to Scott McCauley for the origonal
X version, and John Dunning (jrd@STONY-BROOK.SCRC.Symbolics.COM)
X who got this working.
X
X Please not that only the double float. stuff is picked up from
X here, the other long-int and single float stuff come from
X fixnum.s and sflonum.s (see flonum.h for the appr. #defs).
X ++jrb
X */
X
X/* Subroutines needed by GCC output code for the 68000/20.
X
X Compile using -O flag with gcc.
X Use the -m68000 flag if you have a 68000
X
X This package includes 32 bit integer and 64-bit floating
X point routines.
X
X WARNING: alpha-test version (July 1988) -- probably full of bugs.
X If you find a bug, send bugs reports to jsm@phoenix.princeton.edu,
X or
X
X Scott McCauley
X PPPL P. O. Box 451
X Princeton NJ 08543
X
X Known bugs/features:
X
X 1) Depending on the version of GCC, this may produce a lot of
X warnings about conversion clashes, etc. It should still compile
X as is. Also, there appears to be a bug in the register-allocation
X parts of gcc that makes the compiler sometimes core dump
X or run into an infinite loop. This version works -- if you
X decide to change routines these compiler bugs can bite very hard....
X
X 2) all single-precision gets converted to double precision in the
X math routines (in accordance with K&R), so doing things in
X double precision is faster than single precision....
X
X 3) double precision floating point division may not be accurate to
X the last four or so bits. Most other routines round to the
X lsb.
X
X 4) no support of NaN and Inf
X
X 5) beware of undefined operations: i.e. a float to integer conversion
X when the float is larger than MAXIINT.
X
X */
X
X#include <stdio.h>
X/* #include <math.h> */
X
X#include "flonum.h"
X
X#define umultl _umulsi
X#define multl _mulsi3
X#define udivl _udivsi3
X#define divl _divsi3
X#define ddiv _divdf3
X#define qmult _qmult
X#define dmult _muldf3
X#define dneg _negdf2
X#define dadd _adddf3
X#define dcmp _cmpdf2
X#define dtoul _fixunsdfsi
X#define dtol _fixdfsi
X#define ltod _floatsidf
X
X#if 0 /* see eprintf.c in minix library */
X#ifdef L_eprintf
X#include <stdio.h>
X/* This is used by the `assert' macro. */
Xvoid
X _eprintf (string, line)
Xchar *string;
Xint line;
X{
X fprintf (stderr, string, line);
X}
X#endif
X#endif
X
X#ifdef L_umulsi3
X
X/*_umulsi3 (a, b) unsigned a, b; { return a * b; } */
X
Xunsigned long umultl(a,b)
Xunsigned long a, b;
X{
X register unsigned long d7, d6, d5, d4;
X
X d7 = a;
X d6 = b;
X d5 = d6;
X d4 = d6;
X /* without the next line, gcc may core dump. Gcc sometimes
X gets confused if you have too many registers */
X
X &a; &b;
X
X /*printf("a %u b %u\n", a, b);*/
X
X /* low word */
X MUL(d7, d6);
X SWAP(d5);
X MUL(d7, d5);
X SWAP(d7);
X MUL(d7, d4);
X d4 += d5;
X SWAP(d4);
X d4 &= 0xFFFF0000;
X d4 += d6;
X return(d4);
X}
X#endif
X
X#ifdef L_mulsi3
X/* _mulsi3 (a, b) int a, b; { return a * b; } */
X
Xlong multl(a, b)
Xlong a, b;
X{
X int sign = 0;
X long umultl();
X if ((a ^ b) < 0) sign = 1;
X if (a < 0) a = -a;
X if (b < 0) b = -b;
X /*printf("a is %d b is %d\n", a, b);*/
X if (sign) return(- umultl(a,b));
X else return(umultl(a,b));
X}
X
X#endif
X
X#ifdef L_udivsi3
X/*_udivsi3 (a, b) unsigned a, b; { return a / b; } */
X
X/*
X this routine based on one in the PD forth package for the sun by Mitch Bradley
X */
X
Xunsigned udivl(u, v)
Xregister unsigned long u, v;
X{
X register unsigned short um, ul, vm, vl;
X unsigned long ans;
X unsigned long u1, v1;
X long i;
X long rem;
X
X if (v == 0) {
X /* should cause an exception condition */
X DIV(u, v);
X fprintf(stderr, "division by zero\n");
X }
X if (v > u) return(0);
X
X ul = u; SWAP(u); um = u;
X vl = v; SWAP(v); vm = v;
X if (vm == 0) {
X u = vl; v = um;
X DIV(u, v);
X /* note -- if you delete the next line, gcc goes into
X an infinite loop */
X if (vm) printf("result is %d\n", v);
X vm = v & 0xFFFF; /* dividend is in low word */
X v &= 0xFFFF0000; /* remainder is in high word */
X v += ul;
X DIV(vl, v);
X v &= 0xFFFF; /* dividend is in low word */
X u = vm;
X SWAP(u);
X return(v + u);
X /*ans = ((um / vl) << 16) + ((((um % vl) << 16) + ul) / vl);
X return(ans);*/
X }
X
X if (vl == 0) return(um / vm);
X SWAP(u); SWAP(v);
X if ( (u >> 3) < v) {
X for(i = 0; u >= v; i++, u -= v);
X /*printf("lt 8\n");*/
X return(i);
X }
X u1 = u; v1 = v;
X
X /* scale divisor */
X v1--;
X for(i = 0; ((unsigned) v1) >= 0xFFFF; v1 >>= 1, i++);
X if (++v1 > 0xFFFF) {
X i++; v1 >>=1;
X }
X u1 >>= i;
X /*printf("i is %d, u1 is %x, v1 is %x\n", i, u1, v1);*/
X ans = u1 / v1;
X rem = u - (ans * v);
X if (rem > v) {ans++; rem -= v; }
X if (rem > v) {printf("oops\n");}
X return(ans);
X}
X#endif
X
X#ifdef L_divsi3
X
Xlong divl(a, b)
Xlong a, b;
X{
X int sign = 0;
X if ((a ^ b) < 0) sign = 1;
X if (a < 0) a = -a;
X if (b < 0) b = -b;
X if (sign) return(-udivl(a,b));
X else return(udivl(a,b));
X}
X#endif
X
X#ifdef L_umodsi3
X_umodsi3 (a, b)
Xunsigned a, b;
X{
X /*return a % b;*/
X return (a - ((a/b)*b));
X}
X#endif
X
X#ifdef L_modsi3
X_modsi3 (a, b)
Xint a, b;
X{
X /*return a % b;*/
X return( a - ((a/b) * b));
X}
X#endif
X
X#ifdef L_lshrsi3
X_lshrsi3 (a, b)
Xunsigned a, b;
X{
X return a >> b;
X}
X#endif
X
X#ifdef L_lshlsi3
X_lshlsi3 (a, b)
Xunsigned a, b;
X{
X return a << b;
X}
X#endif
X
X#ifdef L_ashrsi3
X_ashrsi3 (a, b)
Xint a, b;
X{
X return a >> b;
X}
X#endif
X
X#ifdef L_ashlsi3
X_ashlsi3 (a, b)
Xint a, b;
X{
X return a << b;
X}
X#endif
X
X
X/* this divide stuff is hopelessly broken, in addition to being much
X more complicated than in needs to be. The new version's at the
X end of this file */
X
X#if 0
X#ifdef L_divdf3
X
X/*double _divdf3 (a, b) double a, b; { return a / b; } */
X
Xdouble drecip1(f1)
Xdouble f1;
X{
X struct bitdouble *bdp = &f1;
X unsigned m1, m2;
X
X printf("drecip1(%X)", f1);
X if (bdp->exp == 0 ) return(0L);
X if (bdp->mant1 == 0L) {
X bdp->exp = 0x3ff + 0x3ff - bdp->exp;
X bdp->mant2 = 0L;
X return(f1);
X }
X bdp->exp = 0x3ff + 0x3ff - bdp->exp - 1;
X m1 = (0x00100000 + bdp->mant1) >> 5;
X m2 = (0x80000000 / m1);
X /*printf("m1 %x m2 %x\n", m1, m2);*/
X m2 <<= 5;
X m2 &= 0xFFFFF;
X /*printf("exp %x mant %x\n", bdp->exp, m2);*/
X bdp->mant1 = m2;
X bdp->mant2 = 0L;
X printf("drecip1->%X\n", f1);
X return(f1);
X}
X
Xdouble drecip(f)
Xdouble f;
X{
X struct bitdouble *bdp;
X double quotient, remainder;
X
X printf("drecip(%X)", f);
X quotient = drecip1(f);
X remainder = /* 1.0 */ ((double)one) - quotient * f;
X bdp = &remainder;
X for(; bdp->exp > 0x3ca; ) {
X printf("drc: exp=%X ", bdp->exp);
X remainder = /* 1.0 */ ((double)one) - (quotient*f);
X printf("rem=%X ", remainder);
X quotient = quotient + (drecip1(f)*remainder);
X }
X printf("drecip->%X\n", quotient);
X return(quotient);
X}
X
X
Xdouble ddiv(f1, f2)
Xdouble f1, f2;
X{
X return(f1 * drecip(f2));
X}
X#endif
X#endif /* commented out divide routines */
X
X#ifdef L_muldf3
X/*double _muldf3 (a, b) double a, b; { return a * b; } */
Xqmult(m11, m12, m21, m22, p1, p2)
Xunsigned long m11, m12, m21, m22, *p1, *p2;
X{
X/* register unsigned long d2 = m11; */
X register long d2 = m11;
X register unsigned long d3 = m12, d4 = m21, d5 = m22, d6 =0, d7 = 0;
X int i;
X /* guess what happens if you delete the next line.... */
X /* &i; */
X for (i = 0; i < 11; i++) ASL2(d2, d3);
X for (i = 0; i < 9; i++) ASL2(d4, d5);
X
X for (i = 0; i < 64; i++) {
X if (d2 < 0) { ADD2(d4, d5, d6, d7);}
X ASL2(d2, d3);
X ASR2(d4, d5);
X }
X d2 = d6;
X d3 = d7;
X for (i = 0; i < 9; i++) ASR2(d2, d3);
X *p1 = d2; *p2 = d3;
X}
X
Xdouble dmult(f1, f2)
Xdouble f1, f2;
X{
X register unsigned long d2;
X register unsigned d3, d4, d5, d6, d7;
X unsigned long p1, p2;
X
X struct bitdouble
X *bdp1 = (struct bitdouble *)&f1,
X *bdp2 = (struct bitdouble *)&f2;
X int exp1, exp2, i;
X
X exp1 = bdp1->exp; exp2 = bdp2->exp;
X /* check for zero */
X if (! exp1) return(0.0);
X if (! exp2) return(0.0);
X d2 = 0x00100000 + bdp1->mant1;
X d3 = bdp1->mant2;
X d4 = 0x00100000 + bdp2->mant1;
X d5 = bdp2->mant2;
X qmult(d2, d3, d4, d5, &p1, &p2);
X d6 = p1; d7 = p2;
X
X if (d6 & 0x00200000) {
X ASR2(d6, d7);
X exp1++;
X }
X
X if (bdp1->sign == bdp2->sign) bdp1->sign = 0;
X else bdp1->sign = 1;
X bdp1->exp = exp1 + exp2 - 0x3ff;
X bdp1->mant1 = d6;
X bdp1->mant2 = d7;
X return(f1);
X}
X#endif
X
X#ifdef L_negdf2
X/*double _negdf2 (a) double a; { return -a; } */
X
Xdouble dneg(num)
Xdouble num;
X{
X long *i = (long *)#
X *i ^= 0x80000000;
X return(num);
X}
X#endif
X
X#ifdef L_adddf3
X/*double _adddf3 (a, b) double a, b; { return a + b; } */
X
Xdouble dadd(f1, f2)
Xdouble f1, f2;
X{
X
X register long d4, d5, d6, d7;
X struct bitdouble
X *bdp1 = (struct bitdouble *)&f1,
X *bdp2 = (struct bitdouble *)&f2;
X short exp1, exp2, sign1, sign2, howmuch, temp;
X
X exp1 = bdp1->exp; exp2 = bdp2->exp;
X
X /* check for zero */
X if (! exp1) return(f2); if (! exp2) return(f1);
X
X /* align them */
X if (exp1 < exp2) {
X bdp1 = (struct bitdouble *)&f2; bdp2 = (struct bitdouble *)&f1;
X exp1 = bdp1->exp; exp2 = bdp2->exp;
X }
X howmuch = exp1 - exp2;
X if (howmuch > 53) return(f1);
X
X d7 = bdp2->mant1 + 0x00100000;
X d6 = bdp2->mant2;
X
X d5 = bdp1->mant1 + 0x00100000;
X d4 = bdp1->mant2;
X
X for (temp = 0; temp < howmuch; temp++) ASR2(d7, d6);
X
X /* take care of negative signs */
X if (bdp1->sign)
X {
X NEG(d5, d4);
X }
X if (bdp2->sign)
X {
X NEG(d7, d6);
X }
X
X ADD2(d7, d6, d5, d4);
X bdp1 = (struct bitdouble *)&f1;
X
X if (d5 < 0) {
X NEG(d5, d4);
X bdp1->sign = 1;
X } else bdp1->sign = 0;
X
X if (d5 == 0 && d4 == 0) return(0.0);
X
X for(howmuch = 0; d5 >= 0; howmuch++) ASL2(d5, d4);
X
X ASL2(d5, d4);
X for (temp = 0; temp < 12; temp++) ASR2(d5, d4);
X bdp1->mant1 = d5;
X bdp1->mant2 = d4;
X bdp1->exp = exp1 + 11 - howmuch;
X return(f1);
X}
X
X#endif
X
X#ifdef L_subdf3
Xdouble
X _subdf3 (a, b)
Xdouble a, b;
X{
X return a+(-b);
X}
X#endif
X
X#ifdef L_cmpdf2
X/*
X int _cmpdf2 (a, b) double a, b; { if (a > b) return 1;
X else if (a < b) return -1; return 0; }
X */
X
Xint dcmp(f1, f2)
Xdouble f1, f2;
X{
X struct bitdouble *bdp1, *bdp2;
X unsigned int s1, s2;
X bdp1 = (struct bitdouble *)&f1;
X bdp2 = (struct bitdouble *)&f2;
X s1 = bdp1->sign;
X s2 = bdp2->sign;
X if (s1 > s2) return(-1);
X if (s1 < s2) return(1);
X /* if sign of both negative, switch them */
X if (s1 != 0) {
X bdp1 = (struct bitdouble *)&f2;
X bdp2 = (struct bitdouble *)&f1;
X }
X s1 = bdp1->exp;
X s2 = bdp2->exp;
X if (s1 > s2) return(1);
X if (s1 < s2) return(-1);
X /* same exponent -- have to compare mantissa */
X s1 = bdp1->mant1;
X s2 = bdp2->mant1;
X if (s1 > s2) return(1);
X if (s1 < s2) return(-1);
X s1 = bdp1->mant2;
X s2 = bdp2->mant2;
X if (s1 > s2) return(1);
X if (s1 < s2) return(-1);
X return(0); /* the same! */
X}
X#endif
X
X#ifdef L_fixunsdfsi
X/*_fixunsdfsi (a) double a; { return (unsigned int) a; } */
X
X/* #ifdef L_fixdfsi _fixdfsi (a) double a; { return (int) a; } #endif */
X
Xunsigned long dtoul(f)
Xdouble f;
X{
X struct bitdouble *bdp;
X int si, ex, mant1, mant2;
X bdp = (struct bitdouble *)&f;
X si = bdp->sign;
X ex = bdp->exp;
X mant1 = bdp->mant1 + 0x00100000;
X mant2 = bdp->mant2;
X
X /* zero value */
X if (ex == 0) return(0L);
X /* less than 1 */
X if (ex < 0x3ff) return(0L);
X /* less than 0 */
X if (si ) return(0L);
X mant1 <<= 10;
X mant1 += (mant2 >> 22);
X mant1 >>= 30 + (0x3ff - ex);
X return(mant1);
X}
X
Xlong dtol(f)
Xdouble f;
X{
X struct bitdouble *bdp = (struct bitdouble *)&f;
X
X if (bdp->sign) {
X bdp->sign = 0;
X return( - dtoul(f));
X }
X return(dtoul(f));
X}
X#endif
X
X#ifdef L_fixunsdfdi
X
X/*double _fixunsdfdi (a) double a; { union double_di u;
X u.i[LOW] = (unsigned int) a; u.i[HIGH] = 0; return u.d; } */
X#endif
X
X
X#ifdef L_fixdfdi
Xdouble
X _fixdfdi (a)
Xdouble a;
X{
X union double_di u;
X u.i[LOW] = (int) a;
X u.i[HIGH] = (int) a < 0 ? -1 : 0;
X return u.d;
X}
X#endif
X
X#ifdef L_floatsidf
X/* double _floatsidf (a) int a; { return (double) a; } */
X
Xdouble ltod(i)
Xlong i;
X{
X int expo, shift;
X double retval;
X struct bitdouble *bdp = (struct bitdouble *)&retval;
X if (i == 0) {
X long *t = (long *)&retval;
X t[0] = 0L;
X t[1] = 0L;
X return(retval);
X }
X if (i < 0) {
X bdp->sign = 1;
X i = -i;
X } else bdp->sign = 0;
X shift = i;
X for (expo = 0x3ff + 31 ; shift > 0; expo--, shift <<= 1);
X shift <<= 1;
X bdp->exp = expo;
X bdp->mant1 = shift >> 12;
X bdp->mant2 = shift << 20;
X return(retval);
X}
X
X#endif
X
X#ifdef L_floatdidf
X/* ok as is -- will call other routines */
Xdouble
X _floatdidf (u)
Xunion double_di u;
X{
X register double hi
X = ((double) u.i[HIGH]) * (double) 0x10000 * (double) 0x10000;
X register double low = (unsigned int) u.i[LOW];
X return hi + low;
X}
X#endif
X
X#ifdef L_addsf3
XSFVALUE
X _addsf3 (a, b)
Xunion flt_or_int a, b;
X{
X union flt_or_int intify; return INTIFY ((double) a.f + (double) b.f);
X}
X#endif
X
X#ifdef L_negsf2
XSFVALUE
X _negsf2 (a)
Xunion flt_or_int a;
X{
X union flt_or_int intify;
X return INTIFY (-((double) a.f));
X}
X#endif
X
X#ifdef L_subsf3
XSFVALUE
X _subsf3 (a, b)
Xunion flt_or_int a, b;
X{
X union flt_or_int intify;
X return INTIFY (((double) a.f - (double) b.f));
X}
X#endif
X
X#ifdef L_cmpsf2
XSFVALUE
X _cmpsf2 (a, b)
Xunion flt_or_int a, b;
X{
X union flt_or_int intify;
X double a1, b1;
X a1 = a.f; b1 = b.f;
X if ( a1 > b1)
X return 1;
X else if (a1 < b1)
X return -1;
X return 0;
X}
X#endif
X
X#ifdef L_mulsf3
XSFVALUE
X _mulsf3 (a, b)
Xunion flt_or_int a, b;
X{
X union flt_or_int intify;
X return INTIFY (((double) a.f * (double) b.f));
X}
X#endif
X
X#ifdef L_divsf3
XSFVALUE
X _divsf3 (a, b)
Xunion flt_or_int a, b;
X{
X union flt_or_int intify;
X return INTIFY (((double) a.f / (double) b.f));
X}
X#endif
X
X#ifdef L_truncdfsf2
Xfloat dtof(d)
Xdouble d;
X{
X struct bitdouble *bdp = (struct bitdouble *)&d;
X float retval;
X struct bitfloat *bfp = (struct bitfloat *)&retval;
X int tempval;
X
X bfp->sign = bdp->sign;
X if (bdp->exp == 0) return ((float) 0.0);
X bfp->exp = bdp->exp - 0x400 + 0x80;
X tempval = (bdp->mant1 << 4 ) + ((0xF0000000 & bdp->mant2) >> 28);
X /* round */
X tempval++;
X if (tempval == 0x01000000) bfp->exp++;
X bfp->mant = tempval >> 1;
X return(retval);
X}
X
XSFVALUE
X _truncdfsf2 (a)
Xdouble a;
X{
X union flt_or_int intify;
X return INTIFY (dtof(a));
X}
X#endif
X
X#ifdef L_extendsfdf2
Xdouble ftod(f)
Xunion flt_or_int f;
X{
X double retval;
X struct bitfloat *bfp = (struct bitfloat *)&f.f;
X struct bitdouble *bdp = (struct bitdouble *)&retval;
X if (bfp->exp == 0) return(0.0);
X bdp->sign = bfp->sign;
X bdp->exp = 0x400 - 0x80 + bfp->exp;
X bdp->mant1 = bfp->mant >> 3;
X bdp->mant2 = (bfp->mant & 0x7) << 29;
X /*printf("returning %f from extendsfdf2\n", retval);*/
X return(retval);
X}
X
Xdouble
X _extendsfdf2 (a)
Xunion flt_or_int a;
X{
X union flt_or_int intify;
X double retval;
X return (ftod(a));
X}
X#endif
X
X#ifdef L_divdf3
X
X/* new double-float divide routine, by jrd */
X/* thanks jrd !! */
X
Xdouble _divdf3(num, denom)
Xdouble num, denom;
X{
X double local_num = num;
X double local_denom = denom;
X struct bitdouble * num_bdp = (struct bitdouble *)(&local_num);
X struct bitdouble * den_bdp = (struct bitdouble *)(&local_denom);
X short num_exp = num_bdp->exp,
X den_exp = den_bdp->exp;
X short result_sign = 0;
X /* register */ unsigned long num_m1, num_m2, num_m3, num_m4;
X register unsigned long den_m1, den_m2, den_m3, den_m4;
X unsigned long result_mant[2];
X unsigned long result_mask;
X short result_idx;
X
X if ((num_exp == 0) || (den_exp == 0)) /* zzz should really cause trap */
X return(0.0);
X
X /* deal with signs */
X result_sign = result_sign + num_bdp->sign - den_bdp->sign;
X
X /* unpack the numbers */
X num_m1 = num_bdp->mant1 | 0x00100000; /* hidden bit */
X num_m2 = num_bdp->mant2;
X num_m3 = /* ret_kludge(0); */ 0;
X num_m4 = num_m3;
X den_m1 = den_bdp->mant1 | 0x00100000; /* hidden bit */
X den_m2 = den_bdp->mant2;
X den_m3 = /* ret_kludge(0); */ 0;
X den_m4 = den_m3;
X
X#if 0
X /* buy a little extra accuracy by shifting num and denum up 10 bits */
X for (result_idx /* loop counter */ = 0 ; result_idx < 10 ; result_idx++)
X {
X ASL3(num_m1, num_m2, num_m3);
X ASL3(den_m1, den_m2, den_m3);
X }
X#endif
X
X /* hot wire result mask and index */
X result_mask = 0x00100000; /* start at top mantissa bit */
X result_idx = 0; /* of first word */
X result_mant[0] = num_m3;
X result_mant[1] = den_m3;
X
X /* if denom is greater than num here, shift denom right one and dec num expt */
X if (den_m1 < num_m1)
X goto kludge1; /* C is assembly language,
X remember? */
X if (den_m1 > num_m1)
X goto kludge0;
X if (den_m2 <= num_m2) /* first word eq, try 2nd */
X goto kludge1;
X
X kludge0:
X
X num_exp--;
X ASR4(den_m1, den_m2, den_m3, den_m4);
X
X kludge1:
X
X for ( ; result_mask ; )
X {
X /* if num >= den, subtract den from num and set bit in result */
X if (num_m1 > den_m1) goto kludge2;
X if (num_m1 < den_m1) goto kludge3;
X if (num_m2 > den_m2) goto kludge2;
X if (num_m2 < den_m2) goto kludge3;
X if (num_m3 > den_m3) goto kludge2;
X if (num_m3 < den_m3) goto kludge3;
X if (num_m4 < den_m4) goto kludge3;
X
X kludge2:
X result_mant[result_idx] |= result_mask;
X SUB4(den_m1, den_m2, den_m3, den_m4, num_m1, num_m2, num_m3, num_m4);
X kludge3:
X ASR4(den_m1, den_m2, den_m3, den_m4);
X result_mask >>= 1;
X if ((result_mask == 0) && (result_idx == 0))
X {
X result_mask = 0x80000000;
X result_idx++;
X }
X }
X
X /* compute the resultant exponent */
X num_exp = num_exp - den_exp + 0x3FF;
X
X /* reconstruct the result in local_num */
X num_bdp->sign = result_sign;
X num_bdp->exp = num_exp;
X num_bdp->mant1 = result_mant[0] & 0xFFFFF;
X num_bdp->mant2 = result_mant[1];
X
X /* done! */
X return(local_num);
X}
X#endif
SHAR_EOF
if test 18187 -ne "`wc -c 'dflonum.c'`"
then
echo shar: error transmitting "'dflonum.c'" '(should have been 18187 characters)'
fi
echo shar: extracting "'frexp.c'" '(488 characters)'
if test -f 'frexp.c'
then
echo shar: over-writing existing file "'frexp.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'frexp.c'
X/*
X * double frexp(value, eptr)
X * double value;
X * int *eptr;
X *
X * returns significand
X * in *eptr returns n such that value = significand * 2**n
X *
X * ++jrb bammi@dsrgsun.ces.cwru.edu
X *
X */
X#include "flonum.h"
X
X#define BIAS 0x3ff
X
Xdouble frexp(value, eptr)
Xdouble value;
X#ifdef SHORTLIB
Xshort *eptr;
X#else
Xint *eptr;
X#endif
X{
X double x = value;
X struct bitdouble *res = (struct bitdouble *) &x;
X
X *eptr = res->exp - BIAS + 1;
X res->exp = BIAS - 1;
X
X return x;
X}
SHAR_EOF
if test 488 -ne "`wc -c 'frexp.c'`"
then
echo shar: error transmitting "'frexp.c'" '(should have been 488 characters)'
fi
echo shar: extracting "'ldexp.c'" '(343 characters)'
if test -f 'ldexp.c'
then
echo shar: over-writing existing file "'ldexp.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'ldexp.c'
X/*
X * double ldexp(value, exp);
X * double value;
X * int exp;
X *
X * returns value * 2**exp
X *
X * ++jrb bammi@dsrgsun.ces.cwru.edu
X */
X#include "flonum.h"
X
Xdouble ldexp(value, ex)
Xdouble value;
X#ifdef SHORTLIB
Xshort ex;
X#else
Xint ex;
X#endif
X{
X double x = value;
X struct bitdouble *res = (struct bitdouble *) &x;
X
X res->exp += ex;
X return x;
X}
SHAR_EOF
if test 343 -ne "`wc -c 'ldexp.c'`"
then
echo shar: error transmitting "'ldexp.c'" '(should have been 343 characters)'
fi
echo shar: extracting "'modf.c'" '(2388 characters)'
if test -f 'modf.c'
then
echo shar: over-writing existing file "'modf.c'"
fi
sed 's/^X//' << \SHAR_EOF > 'modf.c'
X/*
X * double modf(value, iptr)
X * double value;
X * double *iptr;
X *
X * returns fractional part of value
X * in *iptr rets. the integral part such that (*iptr + modf) == value
X *
X * ++jrb bammi@dsrgsun.ces.cwru.edu
X */
X#include "flonum.h"
X
X/*
X * extract integral portion of an Ieee double -- return it as a double
X * ++jrb
X */
X
X#define BIAS 0x3ff
X#define MANT1 20
X#define MANT2 32
X#define MASK 0xFFFFFFFFL
X
Xstatic double extract_int(val)
Xdouble val;
X{
X struct bitdouble *bd = (struct bitdouble *)&val;
X register short ex = bd->exp - BIAS;
X
X /* trivial case */
X if(ex < 0) return 0.0;
X
X /* exponent <= # of bits in mant1 */
X if(ex <= MANT1)
X {
X bd->mant1 &= (MASK << (MANT1 - ex));
X bd->mant2 = 0L;
X }
X else /* if (ex <= (MANT1 + MANT2)) */
X {
X bd->mant2 &= (MASK << ((MANT1 + MANT2) - ex));
X }
X /* else the value is too big */
X
X return val;
X}
X
X
Xdouble modf(value, iptr)
Xdouble value, *iptr;
X{
X double integral = extract_int(value);
X
X *iptr = integral;
X return value - integral;
X}
X
X#ifdef TEST1
X#ifndef TEST
X#define TEST
X#endif
X#endif
X
X#ifdef TEST
X
X#define MAXRAND 0x7fffffffL
X
X#ifndef TEST1
X#define TESTSIZE 10240
Xdouble testin[TESTSIZE];
X
Xstatic void init()
X{
X register int i;
X extern long rand();
X
X for(i = 0; i < TESTSIZE; i++)
X testin[i] = ((double)rand()) + (((double)rand()) / ((double)MAXRAND));
X}
X
X#else
X
X#define TESTSIZE 23
Xdouble testin[TESTSIZE];
Xstatic void init()
X{
X register int i;
X double v;
X double inc = 0.1;
X
X for(v = -1.1, i = 0; i < TESTSIZE; v += inc, i++)
X testin[i] = v;
X}
X#endif /* TEST1 */
X
X#define ABS(x) (((x) < 0.0)? (-x) : (x))
X#define MARGIN 1.0e-6 /* 6 is arb */
X
X#include <stdio.h>
Xint main()
X{
X double frac, integ, e, r;
X register int i;
X register int errors = 0;
X extern double modf();
X
X init();
X for(i = 0; i < TESTSIZE; i++)
X {
X frac = modf(testin[i], &integ);
X if(frac >= 1.0)
X {
X printf("Error, frac >= 1, testin %f integ %f frac %f\n",
X testin[i], integ, frac);
X errors++;
X continue;
X }
X
X r = integ + frac;
X e = testin[i] - r;
X if(ABS(e) > MARGIN)
X {
X printf("In %f\tInt %f Frac %f Res %f\tError %f\n",
X testin[i], integ, frac, r, e);
X errors++;
X }
X#ifdef TEST1
X printf("%f\t%f %f\n", testin[i], integ, frac);
X#endif
X }
X
X printf("%d error(s)\n", errors);
X return errors;
X}
X#endif /* TEST */
SHAR_EOF
if test 2388 -ne "`wc -c 'modf.c'`"
then
echo shar: error transmitting "'modf.c'" '(should have been 2388 characters)'
fi
# End of shell archive
exit 0
usenet: {decvax,sun}!cwjcc!dsrgsun!bammi jwahar r. bammi
csnet: bammi@dsrgsun.ces.CWRU.edu
arpa: bammi@dsrgsun.ces.CWRU.edu
compuServe: 71515,155