[comp.sys.atari.st] GCC lib Update/Math lib Tos/MinixST shar 1/5

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 *)&num;
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