80386 alternative math library part02/04
Glenn Geers
glenn at extro.ucc.su.OZ.AU
Sat Dec 8 08:56:54 AEST 1990
Submitted-by: glenn at trantor
Archive-name: libfpu/part02
#!/bin/sh
# This is part 02 of libfpu
# ============= fabs.s ==============
if test -f 'fabs.s' -a X"$1" != X"-c"; then
echo 'x - skipping fabs.s (File already exists)'
else
echo 'x - extracting fabs.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fabs.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl fabs
fabs:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fabs
X
X leave
X ret
SHAR_EOF
chmod 0644 fabs.s ||
echo 'restore of fabs.s failed'
Wc_c="`wc -c < 'fabs.s'`"
test 322 -eq "$Wc_c" ||
echo 'fabs.s: original size 322, current size' "$Wc_c"
fi
# ============= finite.s ==============
if test -f 'finite.s' -a X"$1" != X"-c"; then
echo 'x - skipping finite.s (File already exists)'
else
echo 'x - extracting finite.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'finite.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl finite
finite:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x7ff00000, %eax
X cmpl $0x7ff00000, %eax
X je .Lnotfinite
X
X movl $1, %eax
X leave
X ret
X
.Lnotfinite:
X movl $0, %eax
X leave
X ret
SHAR_EOF
chmod 0644 finite.s ||
echo 'restore of finite.s failed'
Wc_c="`wc -c < 'finite.s'`"
test 447 -eq "$Wc_c" ||
echo 'finite.s: original size 447, current size' "$Wc_c"
fi
# ============= floor.s ==============
if test -f 'floor.s' -a X"$1" != X"-c"; then
echo 'x - skipping floor.s (File already exists)'
else
echo 'x - extracting floor.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'floor.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl floor
floor:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X fldl 8(%ebp) /* load data */
X
X fstcw -12(%ebp) /* store control word */
X fstcw -16(%ebp) /* store it again */
X orw $0x0400, -16(%ebp) /* round toward -inf */
X fldcw -16(%ebp) /* store new control word */
X
X frndint /* rounding gives floor(x) */
X fldcw -12(%ebp) /* restore original control word */
X
X leave
X ret
SHAR_EOF
chmod 0644 floor.s ||
echo 'restore of floor.s failed'
Wc_c="`wc -c < 'floor.s'`"
test 628 -eq "$Wc_c" ||
echo 'floor.s: original size 628, current size' "$Wc_c"
fi
# ============= fmod.s ==============
if test -f 'fmod.s' -a X"$1" != X"-c"; then
echo 'x - skipping fmod.s (File already exists)'
else
echo 'x - extracting fmod.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fmod.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl fmod
fmod:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 16(%ebp)
X fldl 8(%ebp)
.Lnotred:
X fprem
X
X fstsw %ax
X sahf
X jp .Lnotred
X
X leave
X ret
SHAR_EOF
chmod 0644 fmod.s ||
echo 'restore of fmod.s failed'
Wc_c="`wc -c < 'fmod.s'`"
test 380 -eq "$Wc_c" ||
echo 'fmod.s: original size 380, current size' "$Wc_c"
fi
# ============= fpumath.h ==============
if test -f 'fpumath.h' -a X"$1" != X"-c"; then
echo 'x - skipping fpumath.h (File already exists)'
else
echo 'x - extracting fpumath.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fpumath.h' &&
/*
** This file is covered by the GNU General Public Licence
**
** This file should be installed somewhere in your standard include
** search path and given the name fpumath.h.
*/
X
/*
** This is a good way of flagging whether
** you are using the alternate stuff
*/
#define __FPU__
X
extern double j0(), j1(), jn(), y0(), y1(), yn();
extern double erf(), erfc();
extern double exp(), log(), log10(), pow(), sqrt();
extern double expm1(), log1p(), exp10();
extern double exp2(), log2();
extern double floor(), ceil(), fabs();
extern double lgamma(), gamma();
extern double sinh(), cosh(), tanh();
extern double asinh(), acosh(), atanh();
extern double sin(), cos(), tan(), asin();
extern double acos(), atan(), atan2(), hypot();
extern double sqrtp();
X
/*
** IEEE functions
*/
X
extern double rint(), drem(), scalb(), logb(), copysign();
extern double infinity(), nextafter();
extern int isnan(), iszero(), isinf(), isnormal(), issubnormal();
extern int signbit(), finite();
extern void ieee_retrospective();
extern double max_normal(), min_normal(), min_subnormal(), max_subnormal();
extern double quiet_nan(),signaling_nan();
extern double nextafter();
X
/*
** Extensions
*/
X
extern int setinternal(), setcont();
X
X
/*
** Useful defines for setcont
*/
X
#define MASK_ALL 0x127f
#define MASK_ALL1 0x137f /* extra precision */
X
#define __infinity infinity
X
#define M_E 2.7182818284590452354
#define M_LOG2E 1.4426950408889634074
#define M_LOG10E 0.43429448190325182765
#define M_LN2 0.69314718055994530942
#define M_LN10 2.30258509299404568402
#define M_PI 3.14159265358979323846
#define M_PI_2 1.57079632679489661923
#define M_PI_4 0.78539816339744830962
#define M_1_PI 0.31830988618379067154
#define M_2_PI 0.63661977236758134308
#define M_2_SQRTPI 1.12837916709551257390
#define M_SQRT2 1.41421356237309504880
#define M_SQRT1_2 0.70710678118654752440
X
#define MAXFLOAT ((float)3.40282346638528860e+38)
#define MINFLOAT ((float)1.40129846432481707e-45)
X
#ifndef IEEE
#define MAXDOUBLE 1.79769313486231570e+308 /* wrong in math.h */
#define MINDOUBLE 4.94065645841246544e-324
#else
#define MAXDOUBLE (max_normal())
#define MINDOUBLE (min_subnormal())
#endif
X
#define ULP 1.1102230246251565e-16
X
#define HUGE_VAL 3.40282346638528860e+38
X
#ifdef IEEE
#define HUGE (__infinity())
#else
#define HUGE MAXDOUBLE
#endif
SHAR_EOF
chmod 0644 fpumath.h ||
echo 'restore of fpumath.h failed'
Wc_c="`wc -c < 'fpumath.h'`"
test 2314 -eq "$Wc_c" ||
echo 'fpumath.h: original size 2314, current size' "$Wc_c"
fi
# ============= gamma.c ==============
if test -f 'gamma.c' -a X"$1" != X"-c"; then
echo 'x - skipping gamma.c (File already exists)'
else
echo 'x - extracting gamma.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'gamma.c' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** This implementation takes no notice of the fact that gamma(z) is undefined
** if z is 0 or a negative integer - beware.
**
** Copyright 1990 G. Geers
**
*/
X
/*
** The limits are approximate
*/
#define U_LIM 100.0
#define L_LIM -100.0
X
extern double lgamma();
extern double exp();
extern int signgam;
X
double gamma(x)
double x;
{
X if (x >= L_LIM && x <= U_LIM)
X return((double)signgam*exp(lgamma(x)));
}
SHAR_EOF
chmod 0644 gamma.c ||
echo 'restore of gamma.c failed'
Wc_c="`wc -c < 'gamma.c'`"
test 604 -eq "$Wc_c" ||
echo 'gamma.c: original size 604, current size' "$Wc_c"
fi
# ============= hypot.s ==============
if test -f 'hypot.s' -a X"$1" != X"-c"; then
echo 'x - skipping hypot.s (File already exists)'
else
echo 'x - extracting hypot.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'hypot.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
.globl hypot
hypot:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fmull 8(%ebp)
X fldl 16(%ebp)
X fmull 16(%ebp)
X faddp
X fsqrt
X
X leave
X ret
SHAR_EOF
chmod 0644 hypot.s ||
echo 'restore of hypot.s failed'
Wc_c="`wc -c < 'hypot.s'`"
test 377 -eq "$Wc_c" ||
echo 'hypot.s: original size 377, current size' "$Wc_c"
fi
# ============= ieee_ext.s ==============
if test -f 'ieee_ext.s' -a X"$1" != X"-c"; then
echo 'x - skipping ieee_ext.s (File already exists)'
else
echo 'x - extracting ieee_ext.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ieee_ext.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl isnan
isnan:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x7ff00000, %eax
X cmpl $0x7ff00000, %eax
X jne .Lnotnan
X movl 12(%ebp), %eax
X andl $0xfffff, %eax
X orl 8(%ebp), %eax
X je .Lnotnan
X
X movl $1, %eax
X leave
X ret
X
.Lnotnan:
X movl $0, %eax
X
.Ldone:
X leave
X ret
X
X .align 4
X .globl isinf
isinf:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x7ff00000, %eax
X cmpl $0x7ff00000, %eax
X je .Lcouldbeinf
X
.Lnotinf:
X movl $0, %eax
X leave
X ret
X
.Lcouldbeinf:
X andl $0xfffff, %eax
X orl 8(%ebp), %eax
X jne .Lnotinf
X
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl iszero
iszero:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X cmpl $0x0, %eax
X je .Lcouldbezero
.Lnotzero:
X movl $0, %eax
X leave
X ret
X
.Lcouldbezero:
X andl $0xfffff, %eax
X orl 8(%ebp), %eax
X jne .Lnotzero
X
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl signbit
signbit:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x80000000, %eax
X cmpl $0x80000000, %eax
X jne .Lpos
X movl $1, %eax
X leave
X ret
X
.Lpos:
X movl $0, %eax
X leave
X ret
X
X .align 4
X .globl issubnormal
issubnormal:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x7ff00000, %eax
X cmpl $0x0, %eax
X je .Lcouldbesub
X
.Lnotsubnorm:
X movl $0, %eax
X leave
X ret
X
.Lcouldbesub:
X movl 12(%ebp), %eax
X andl $0xfffff, %eax
X orl 8(%ebp), %eax
X je .Lnotsubnorm
X
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl isnormal
isnormal:
X pushl %ebp
X movl %esp,%ebp
X
X movl 12(%ebp), %eax
X andl $0x7ff00000, %eax /* mask sign bit */
X xorl $0x7ff00000, %eax
X cmpl $0x0, %eax
X je .Lnotnorm
X cmpl $0x7ff00000, %eax
X je .Lcouldbenorm
X
.Lnorm:
X movl $1, %eax
X leave
X ret
X
.Lcouldbenorm:
X movl 12(%ebp), %eax
X andl $0xfffff, %eax
X orl 8(%ebp), %eax
X je .Lnorm
X
.Lnotnorm:
X movl $0, %eax
X leave
X ret
SHAR_EOF
chmod 0644 ieee_ext.s ||
echo 'restore of ieee_ext.s failed'
Wc_c="`wc -c < 'ieee_ext.s'`"
test 1977 -eq "$Wc_c" ||
echo 'ieee_ext.s: original size 1977, current size' "$Wc_c"
fi
# ============= ieee_retro.c ==============
if test -f 'ieee_retro.c' -a X"$1" != X"-c"; then
echo 'x - skipping ieee_retro.c (File already exists)'
else
echo 'x - extracting ieee_retro.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ieee_retro.c' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
#include <stdio.h>
X
ieee_retrospective(f)
FILE *f;
{
X int status;
X
X if (f == (FILE *)NULL)
X f = stderr;
X
X status = _getsw();
X
X if (status & 0xdf) {
X fprintf(f,"\n");
X fprintf(f,"Warning: the following IEEE floating-point");
X fprintf(f," arithmetic exceptions\n");
X fprintf(f,"occurred in this program and were never cleared:\n");
X if (status & 0x0001) fprintf(f,"Invalid Operation;\n");
X if (status & 0x0002) fprintf(f,"Denormal;\n");
X if (status & 0x0004) fprintf(f,"Division by Zero;\n");
X if (status & 0x0008) fprintf(f,"Overflow;\n");
X if (status & 0x0010) fprintf(f,"Underflow;\n");
X fprintf(f,"\n");
X }
}
SHAR_EOF
chmod 0644 ieee_retro.c ||
echo 'restore of ieee_retro.c failed'
Wc_c="`wc -c < 'ieee_retro.c'`"
test 853 -eq "$Wc_c" ||
echo 'ieee_retro.c: original size 853, current size' "$Wc_c"
fi
# ============= ieee_values.s ==============
if test -f 'ieee_values.s' -a X"$1" != X"-c"; then
echo 'x - skipping ieee_values.s (File already exists)'
else
echo 'x - extracting ieee_values.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ieee_values.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl max_normal
X
max_normal:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7fefffff, -12(%ebp)
X movl $0xffffffff, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
X
X .align 4
X .globl min_normal
X
min_normal:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x00100000, -12(%ebp)
X movl $0x00000001, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
X
X .align 4
X .globl min_subnormal
X
min_subnormal:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x0, -12(%ebp)
X movl $0x00000001, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
X
X .align 4
X .globl max_subnormal
X
max_subnormal:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x000fffff, -12(%ebp)
X movl $0xffffffff, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
X
X .align 4
X .globl quiet_nan
X
quiet_nan:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7fffffff, -12(%ebp)
X movl $0xffffffff, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
X
X .align 4
X .globl signaling_nan
X
signaling_nan:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7ff00000, -12(%ebp)
X movl $0x00000001, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 ieee_values.s ||
echo 'restore of ieee_values.s failed'
Wc_c="`wc -c < 'ieee_values.s'`"
test 1295 -eq "$Wc_c" ||
echo 'ieee_values.s: original size 1295, current size' "$Wc_c"
fi
# ============= infinity.s ==============
if test -f 'infinity.s' -a X"$1" != X"-c"; then
echo 'x - skipping infinity.s (File already exists)'
else
echo 'x - extracting infinity.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'infinity.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl infinity
infinity:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7ff00000, -12(%ebp)
X movl $0x0, -16(%ebp)
X
X fldl -16(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 infinity.s ||
echo 'restore of infinity.s failed'
Wc_c="`wc -c < 'infinity.s'`"
test 394 -eq "$Wc_c" ||
echo 'infinity.s: original size 394, current size' "$Wc_c"
fi
# ============= j0.c ==============
if test -f 'j0.c' -a X"$1" != X"-c"; then
echo 'x - skipping j0.c (File already exists)'
else
echo 'x - extracting j0.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'j0.c' &&
/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved. The Berkeley software License Agreement
X * specifies the terms and conditions for redistribution.
X */
X
#ifndef lint
static char sccsid[] = "@(#)j0.c 5.3 (Berkeley) 9/22/88";
#endif /* not lint */
X
/*
X floating point Bessel's function
X of the first and second kinds
X of order zero
X
X j0(x) returns the value of J0(x)
X for all real values of x.
X
X There are no error returns.
X Calls sin, cos, sqrt.
X
X There is a niggling bug in J0 which
X causes errors up to 2e-16 for x in the
X interval [-8,8].
X The bug is caused by an inappropriate order
X of summation of the series. rhm will fix it
X someday.
X
X Coefficients are from Hart & Cheney.
X #5849 (19.22D)
X #6549 (19.25D)
X #6949 (19.41D)
X
X y0(x) returns the value of Y0(x)
X for positive real values of x.
X For x<=0, if on the VAX, error number EDOM is set and
X the reserved operand fault is generated;
X otherwise (an IEEE machine) an invalid operation is performed.
X
X Calls sin, cos, sqrt, log, j0.
X
X The values of Y0 have not been checked
X to more than ten places.
X
X Coefficients are from Hart & Cheney.
X #6245 (18.78D)
X #6549 (19.25D)
X #6949 (19.41D)
*/
X
#include "mathimpl.h"
X
extern double sin();
extern double cos();
extern double sqrt();
extern double log();
X
double j0();
X
#if defined(vax)||defined(tahoe)
#include <errno.h>
#else /* defined(vax)||defined(tahoe) */
static const double zero = 0.e0;
#endif /* defined(vax)||defined(tahoe) */
X
static double pzero, qzero;
X
static const double tpi = .6366197723675813430755350535e0;
static const double pio4 = .7853981633974483096156608458e0;
static const double p1[] = {
X 0.4933787251794133561816813446e21,
X -.1179157629107610536038440800e21,
X 0.6382059341072356562289432465e19,
X -.1367620353088171386865416609e18,
X 0.1434354939140344111664316553e16,
X -.8085222034853793871199468171e13,
X 0.2507158285536881945555156435e11,
X -.4050412371833132706360663322e8,
X 0.2685786856980014981415848441e5,
};
static const double q1[] = {
X 0.4933787251794133562113278438e21,
X 0.5428918384092285160200195092e19,
X 0.3024635616709462698627330784e17,
X 0.1127756739679798507056031594e15,
X 0.3123043114941213172572469442e12,
X 0.6699987672982239671814028660e9,
X 0.1114636098462985378182402543e7,
X 0.1363063652328970604442810507e4,
X 1.0
};
static const double p2[] = {
X 0.5393485083869438325262122897e7,
X 0.1233238476817638145232406055e8,
X 0.8413041456550439208464315611e7,
X 0.2016135283049983642487182349e7,
X 0.1539826532623911470917825993e6,
X 0.2485271928957404011288128951e4,
X 0.0,
};
static const double q2[] = {
X 0.5393485083869438325560444960e7,
X 0.1233831022786324960844856182e8,
X 0.8426449050629797331554404810e7,
X 0.2025066801570134013891035236e7,
X 0.1560017276940030940592769933e6,
X 0.2615700736920839685159081813e4,
X 1.0,
};
static const double p3[] = {
X -.3984617357595222463506790588e4,
X -.1038141698748464093880530341e5,
X -.8239066313485606568803548860e4,
X -.2365956170779108192723612816e4,
X -.2262630641933704113967255053e3,
X -.4887199395841261531199129300e1,
X 0.0,
};
static const double q3[] = {
X 0.2550155108860942382983170882e6,
X 0.6667454239319826986004038103e6,
X 0.5332913634216897168722255057e6,
X 0.1560213206679291652539287109e6,
X 0.1570489191515395519392882766e5,
X 0.4087714673983499223402830260e3,
X 1.0,
};
static const double p4[] = {
X -.2750286678629109583701933175e20,
X 0.6587473275719554925999402049e20,
X -.5247065581112764941297350814e19,
X 0.1375624316399344078571335453e18,
X -.1648605817185729473122082537e16,
X 0.1025520859686394284509167421e14,
X -.3436371222979040378171030138e11,
X 0.5915213465686889654273830069e8,
X -.4137035497933148554125235152e5,
};
static const double q4[] = {
X 0.3726458838986165881989980e21,
X 0.4192417043410839973904769661e19,
X 0.2392883043499781857439356652e17,
X 0.9162038034075185262489147968e14,
X 0.2613065755041081249568482092e12,
X 0.5795122640700729537480087915e9,
X 0.1001702641288906265666651753e7,
X 0.1282452772478993804176329391e4,
X 1.0,
};
X
static void asympt();
X
double
j0(arg) double arg;{
X double argsq, n, d;
X int i;
X
X if(arg < 0.) arg = -arg;
X if(arg > 8.){
X asympt(arg);
X n = arg - pio4;
X return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
X }
X argsq = arg*arg;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*argsq + p1[i];
X d = d*argsq + q1[i];
X }
X return(n/d);
}
X
double
y0(arg) double arg;{
X double argsq, n, d;
X int i;
X
X if(arg <= 0.){
#if defined(vax)||defined(tahoe)
X return(infnan(EDOM)); /* NaN */
#else /* defined(vax)||defined(tahoe) */
X return(zero/zero); /* IEEE machines: invalid operation */
#endif /* defined(vax)||defined(tahoe) */
X }
X if(arg > 8.){
X asympt(arg);
X n = arg - pio4;
X return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
X }
X argsq = arg*arg;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*argsq + p4[i];
X d = d*argsq + q4[i];
X }
X return(n/d + tpi*j0(arg)*log(arg));
}
X
static void
asympt(arg) double arg;{
X double zsq, n, d;
X int i;
X zsq = 64./(arg*arg);
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p2[i];
X d = d*zsq + q2[i];
X }
X pzero = n/d;
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p3[i];
X d = d*zsq + q3[i];
X }
X qzero = (8./arg)*(n/d);
}
SHAR_EOF
chmod 0644 j0.c ||
echo 'restore of j0.c failed'
Wc_c="`wc -c < 'j0.c'`"
test 5099 -eq "$Wc_c" ||
echo 'j0.c: original size 5099, current size' "$Wc_c"
fi
# ============= j1.c ==============
if test -f 'j1.c' -a X"$1" != X"-c"; then
echo 'x - skipping j1.c (File already exists)'
else
echo 'x - extracting j1.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'j1.c' &&
/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved. The Berkeley software License Agreement
X * specifies the terms and conditions for redistribution.
X */
X
#ifndef lint
static char sccsid[] = "@(#)j1.c 5.3 (Berkeley) 9/22/88";
#endif /* not lint */
X
/*
X floating point Bessel's function
X of the first and second kinds
X of order one
X
X j1(x) returns the value of J1(x)
X for all real values of x.
X
X There are no error returns.
X Calls sin, cos, sqrt.
X
X There is a niggling bug in J1 which
X causes errors up to 2e-16 for x in the
X interval [-8,8].
X The bug is caused by an inappropriate order
X of summation of the series. rhm will fix it
X someday.
X
X Coefficients are from Hart & Cheney.
X #6050 (20.98D)
X #6750 (19.19D)
X #7150 (19.35D)
X
X y1(x) returns the value of Y1(x)
X for positive real values of x.
X For x<=0, if on the VAX, error number EDOM is set and
X the reserved operand fault is generated;
X otherwise (an IEEE machine) an invalid operation is performed.
X
X Calls sin, cos, sqrt, log, j1.
X
X The values of Y1 have not been checked
X to more than ten places.
X
X Coefficients are from Hart & Cheney.
X #6447 (22.18D)
X #6750 (19.19D)
X #7150 (19.35D)
*/
X
#include "mathimpl.h"
X
extern double sin();
extern double cos();
extern double sqrt();
extern double log();
X
double j1();
X
#if defined(vax)||defined(tahoe)
#include <errno.h>
#else /* defined(vax)||defined(tahoe) */
static const double zero = 0.e0;
#endif /* defined(vax)||defined(tahoe) */
X
static double pzero, qzero;
X
static const double tpi = .6366197723675813430755350535e0;
static const double pio4 = .7853981633974483096156608458e0;
static const double p1[] = {
X 0.581199354001606143928050809e21,
X -.6672106568924916298020941484e20,
X 0.2316433580634002297931815435e19,
X -.3588817569910106050743641413e17,
X 0.2908795263834775409737601689e15,
X -.1322983480332126453125473247e13,
X 0.3413234182301700539091292655e10,
X -.4695753530642995859767162166e7,
X 0.2701122710892323414856790990e4,
};
static const double q1[] = {
X 0.1162398708003212287858529400e22,
X 0.1185770712190320999837113348e20,
X 0.6092061398917521746105196863e17,
X 0.2081661221307607351240184229e15,
X 0.5243710262167649715406728642e12,
X 0.1013863514358673989967045588e10,
X 0.1501793594998585505921097578e7,
X 0.1606931573481487801970916749e4,
X 1.0,
};
static const double p2[] = {
X -.4435757816794127857114720794e7,
X -.9942246505077641195658377899e7,
X -.6603373248364939109255245434e7,
X -.1523529351181137383255105722e7,
X -.1098240554345934672737413139e6,
X -.1611616644324610116477412898e4,
X 0.0,
};
static const double q2[] = {
X -.4435757816794127856828016962e7,
X -.9934124389934585658967556309e7,
X -.6585339479723087072826915069e7,
X -.1511809506634160881644546358e7,
X -.1072638599110382011903063867e6,
X -.1455009440190496182453565068e4,
X 1.0,
};
static const double p3[] = {
X 0.3322091340985722351859704442e5,
X 0.8514516067533570196555001171e5,
X 0.6617883658127083517939992166e5,
X 0.1849426287322386679652009819e5,
X 0.1706375429020768002061283546e4,
X 0.3526513384663603218592175580e2,
X 0.0,
};
static const double q3[] = {
X 0.7087128194102874357377502472e6,
X 0.1819458042243997298924553839e7,
X 0.1419460669603720892855755253e7,
X 0.4002944358226697511708610813e6,
X 0.3789022974577220264142952256e5,
X 0.8638367769604990967475517183e3,
X 1.0,
};
static const double p4[] = {
X -.9963753424306922225996744354e23,
X 0.2655473831434854326894248968e23,
X -.1212297555414509577913561535e22,
X 0.2193107339917797592111427556e20,
X -.1965887462722140658820322248e18,
X 0.9569930239921683481121552788e15,
X -.2580681702194450950541426399e13,
X 0.3639488548124002058278999428e10,
X -.2108847540133123652824139923e7,
X 0.0,
};
static const double q4[] = {
X 0.5082067366941243245314424152e24,
X 0.5435310377188854170800653097e22,
X 0.2954987935897148674290758119e20,
X 0.1082258259408819552553850180e18,
X 0.2976632125647276729292742282e15,
X 0.6465340881265275571961681500e12,
X 0.1128686837169442121732366891e10,
X 0.1563282754899580604737366452e7,
X 0.1612361029677000859332072312e4,
X 1.0,
};
X
static void asympt();
X
double
j1(arg) double arg;{
X double xsq, n, d, x;
X int i;
X
X x = arg;
X if(x < 0.) x = -x;
X if(x > 8.){
X asympt(x);
X n = x - 3.*pio4;
X n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
X if(arg <0.) n = -n;
X return(n);
X }
X xsq = x*x;
X for(n=0,d=0,i=8;i>=0;i--){
X n = n*xsq + p1[i];
X d = d*xsq + q1[i];
X }
X return(arg*n/d);
}
X
double
y1(arg) double arg;{
X double xsq, n, d, x;
X int i;
X
X x = arg;
X if(x <= 0.){
#if defined(vax)||defined(tahoe)
X return(infnan(EDOM)); /* NaN */
#else /* defined(vax)||defined(tahoe) */
X return(zero/zero); /* IEEE machines: invalid operation */
#endif /* defined(vax)||defined(tahoe) */
X }
X if(x > 8.){
X asympt(x);
X n = x - 3*pio4;
X return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
X }
X xsq = x*x;
X for(n=0,d=0,i=9;i>=0;i--){
X n = n*xsq + p4[i];
X d = d*xsq + q4[i];
X }
X return(x*n/d + tpi*(j1(x)*log(x)-1./x));
}
X
static void
asympt(arg) double arg;{
X double zsq, n, d;
X int i;
X zsq = 64./(arg*arg);
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p2[i];
X d = d*zsq + q2[i];
X }
X pzero = n/d;
X for(n=0,d=0,i=6;i>=0;i--){
X n = n*zsq + p3[i];
X d = d*zsq + q3[i];
X }
X qzero = (8./arg)*(n/d);
}
SHAR_EOF
chmod 0644 j1.c ||
echo 'restore of j1.c failed'
Wc_c="`wc -c < 'j1.c'`"
test 5169 -eq "$Wc_c" ||
echo 'j1.c: original size 5169, current size' "$Wc_c"
fi
# ============= jn.c ==============
if test -f 'jn.c' -a X"$1" != X"-c"; then
echo 'x - skipping jn.c (File already exists)'
else
echo 'x - extracting jn.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'jn.c' &&
/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved. The Berkeley software License Agreement
X * specifies the terms and conditions for redistribution.
X */
X
#ifndef lint
static char sccsid[] = "@(#)jn.c 5.3 (Berkeley) 9/22/88";
#endif /* not lint */
X
/*
X floating point Bessel's function of
X the first and second kinds and of
X integer order.
X
X int n;
X double x;
X jn(n,x);
X
X returns the value of Jn(x) for all
X integer values of n and all real values
X of x.
X
X There are no error returns.
X Calls j0, j1.
X
X For n=0, j0(x) is called,
X for n=1, j1(x) is called,
X for n<x, forward recursion us used starting
X from values of j0(x) and j1(x).
X for n>x, a continued fraction approximation to
X j(n,x)/j(n-1,x) is evaluated and then backward
X recursion is used starting from a supposed value
X for j(n,x). The resulting value of j(0,x) is
X compared with the actual value to correct the
X supposed value of j(n,x).
X
X yn(n,x) is similar in all respects, except
X that forward recursion is used for all
X values of n>1.
*/
X
#include "mathimpl.h"
X
extern double j0();
extern double j1();
X
#if defined(vax)||defined(tahoe)
#include <errno.h>
#else /* defined(vax)||defined(tahoe) */
static double zero = 0.e0;
#endif /* defined(vax)||defined(tahoe) */
X
double
jn(n,x) int n; double x;{
X int i;
X double a, b, temp;
X double xsq, t;
X
X if(n<0){
X n = -n;
X x = -x;
X }
X if(n==0) return(j0(x));
X if(n==1) return(j1(x));
X if(x == 0.) return(0.);
X if(n>x) goto recurs;
X
X a = j0(x);
X b = j1(x);
X for(i=1;i<n;i++){
X temp = b;
X b = (2.*i/x)*b - a;
X a = temp;
X }
X return(b);
X
recurs:
X xsq = x*x;
X for(t=0,i=n+16;i>n;i--){
X t = xsq/(2.*i - t);
X }
X t = x/(2.*n-t);
X
X a = t;
X b = 1;
X for(i=n-1;i>0;i--){
X temp = b;
X b = (2.*i/x)*b - a;
X a = temp;
X }
X return(t*j0(x)/b);
}
X
double
yn(n,x) int n; double x;{
X int i;
X int sign;
X double a, b, temp;
X
X if (x <= 0) {
#if defined(vax)||defined(tahoe)
X return(infnan(EDOM)); /* NaN */
#else /* defined(vax)||defined(tahoe) */
X return(zero/zero); /* IEEE machines: invalid operation */
#endif /* defined(vax)||defined(tahoe) */
X }
X sign = 1;
X if(n<0){
X n = -n;
X if(n%2 == 1) sign = -1;
X }
X if(n==0) return(y0(x));
X if(n==1) return(sign*y1(x));
X
X a = y0(x);
X b = y1(x);
X for(i=1;i<n;i++){
X temp = b;
X b = (2.*i/x)*b - a;
X a = temp;
X }
X return(sign*b);
}
SHAR_EOF
chmod 0644 jn.c ||
echo 'restore of jn.c failed'
Wc_c="`wc -c < 'jn.c'`"
test 2310 -eq "$Wc_c" ||
echo 'jn.c: original size 2310, current size' "$Wc_c"
fi
# ============= lgamma.c ==============
if test -f 'lgamma.c' -a X"$1" != X"-c"; then
echo 'x - skipping lgamma.c (File already exists)'
else
echo 'x - extracting lgamma.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'lgamma.c' &&
/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved. The Berkeley software License Agreement
X * specifies the terms and conditions for redistribution.
X */
X
#ifndef lint
static char sccsid[] = "@(#)lgamma.c 5.3 (Berkeley) 9/22/88";
#endif /* not lint */
X
/*
X C program for floating point log Gamma function
X
X lgamma(x) computes the log of the absolute
X value of the Gamma function.
X The sign of the Gamma function is returned in the
X external quantity signgam.
X
X The coefficients for expansion around zero
X are #5243 from Hart & Cheney; for expansion
X around infinity they are #5404.
X
X Calls log, floor and sin.
*/
X
#include "mathimpl.h"
X
extern double log();
extern double floor();
extern double sin();
X
#if defined(vax)||defined(tahoe)
#include <errno.h>
#endif /* defined(vax)||defined(tahoe) */
X
int signgam = 0;
static const double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */
static const double pi = 3.1415926535897932384626434;
X
#define M 6
#define N 8
static const double p1[] = {
X 0.83333333333333101837e-1,
X -.277777777735865004e-2,
X 0.793650576493454e-3,
X -.5951896861197e-3,
X 0.83645878922e-3,
X -.1633436431e-2,
};
static const double p2[] = {
X -.42353689509744089647e5,
X -.20886861789269887364e5,
X -.87627102978521489560e4,
X -.20085274013072791214e4,
X -.43933044406002567613e3,
X -.50108693752970953015e2,
X -.67449507245925289918e1,
X 0.0,
};
static const double q2[] = {
X -.42353689509744090010e5,
X -.29803853309256649932e4,
X 0.99403074150827709015e4,
X -.15286072737795220248e4,
X -.49902852662143904834e3,
X 0.18949823415702801641e3,
X -.23081551524580124562e2,
X 0.10000000000000000000e1,
};
X
static double pos(), neg(), asym();
X
double
lgamma(arg)
double arg;
{
X
X signgam = 1.;
X if(arg <= 0.) return(neg(arg));
X if(arg > 8.) return(asym(arg));
X return(log(pos(arg)));
}
X
static double
asym(arg)
double arg;
{
X double n, argsq;
X int i;
X
X argsq = 1./(arg*arg);
X for(n=0,i=M-1; i>=0; i--){
X n = n*argsq + p1[i];
X }
X return((arg-.5)*log(arg) - arg + goobie + n/arg);
}
X
static double
neg(arg)
double arg;
{
X double t;
X
X arg = -arg;
X /*
X * to see if arg were a true integer, the old code used the
X * mathematically correct observation:
X * sin(n*pi) = 0 <=> n is an integer.
X * but in finite precision arithmetic, sin(n*PI) will NEVER
X * be zero simply because n*PI is a rational number. hence
X * it failed to work with our newer, more accurate sin()
X * which uses true pi to do the argument reduction...
X * temp = sin(pi*arg);
X */
X t = floor(arg);
X if (arg - t > 0.5e0)
X t += 1.e0; /* t := integer nearest arg */
#if defined(vax)||defined(tahoe)
X if (arg == t) {
X return(infnan(ERANGE)); /* +INF */
X }
#endif /* defined(vax)||defined(tahoe) */
X signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */
X /* 0 if t was even */
X signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */
X /* -1 if t was even */
X t = arg - t; /* -0.5 <= t <= 0.5 */
X if (t < 0.e0) {
X t = -t;
X signgam = -signgam;
X }
X return(-log(arg*pos(arg)*sin(pi*t)/pi));
}
X
static double
pos(arg)
double arg;
{
X double n, d, s;
X register i;
X
X if(arg < 2.) return(pos(arg+1.)/arg);
X if(arg > 3.) return((arg-1.)*pos(arg-1.));
X
X s = arg - 2.;
X for(n=0,d=0,i=N-1; i>=0; i--){
X n = n*s + p2[i];
X d = d*s + q2[i];
X }
X return(n/d);
}
SHAR_EOF
chmod 0644 lgamma.c ||
echo 'restore of lgamma.c failed'
Wc_c="`wc -c < 'lgamma.c'`"
test 3382 -eq "$Wc_c" ||
echo 'lgamma.c: original size 3382, current size' "$Wc_c"
fi
# ============= log.s ==============
if test -f 'log.s' -a X"$1" != X"-c"; then
echo 'x - skipping log.s (File already exists)'
else
echo 'x - extracting log.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'log.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl log
log:
X pushl %ebp
X movl %esp,%ebp
X
X fldln2
X fldl 8(%ebp)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 log.s ||
echo 'restore of log.s failed'
Wc_c="`wc -c < 'log.s'`"
test 329 -eq "$Wc_c" ||
echo 'log.s: original size 329, current size' "$Wc_c"
fi
# ============= log10.s ==============
if test -f 'log10.s' -a X"$1" != X"-c"; then
echo 'x - skipping log10.s (File already exists)'
else
echo 'x - extracting log10.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'log10.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl log10
log10:
X pushl %ebp
X movl %esp,%ebp
X
X fldlg2
X fldl 8(%ebp)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 log10.s ||
echo 'restore of log10.s failed'
Wc_c="`wc -c < 'log10.s'`"
test 333 -eq "$Wc_c" ||
echo 'log10.s: original size 333, current size' "$Wc_c"
fi
# ============= log1p.s ==============
if test -f 'log1p.s' -a X"$1" != X"-c"; then
echo 'x - skipping log1p.s (File already exists)'
else
echo 'x - extracting log1p.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'log1p.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl log1p
log1p:
X pushl %ebp
X movl %esp,%ebp
X
X fldln2
X fldl 8(%ebp)
X fld1
X faddp
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 log1p.s ||
echo 'restore of log1p.s failed'
Wc_c="`wc -c < 'log1p.s'`"
test 346 -eq "$Wc_c" ||
echo 'log1p.s: original size 346, current size' "$Wc_c"
fi
# ============= log2.s ==============
if test -f 'log2.s' -a X"$1" != X"-c"; then
echo 'x - skipping log2.s (File already exists)'
else
echo 'x - extracting log2.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'log2.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl log2
log2:
X pushl %ebp
X movl %esp,%ebp
X
X fld1
X fldl 8(%ebp)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 log2.s ||
echo 'restore of log2.s failed'
Wc_c="`wc -c < 'log2.s'`"
test 329 -eq "$Wc_c" ||
echo 'log2.s: original size 329, current size' "$Wc_c"
fi
# ============= logb.s ==============
if test -f 'logb.s' -a X"$1" != X"-c"; then
echo 'x - skipping logb.s (File already exists)'
else
echo 'x - extracting logb.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'logb.s' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl logb
logb:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 8(%ebp)
X fxtract
X fldl %st(1)
X
X leave
X ret
SHAR_EOF
chmod 0644 logb.s ||
echo 'restore of logb.s failed'
Wc_c="`wc -c < 'logb.s'`"
test 339 -eq "$Wc_c" ||
echo 'logb.s: original size 339, current size' "$Wc_c"
fi
# ============= mathimpl.h ==============
if test -f 'mathimpl.h' -a X"$1" != X"-c"; then
echo 'x - skipping mathimpl.h (File already exists)'
else
echo 'x - extracting mathimpl.h (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'mathimpl.h' &&
/*
X * Copyright (c) 1988 The Regents of the University of California.
X * All rights reserved.
X *
X * Redistribution and use in source and binary forms are permitted
X * provided that: (1) source distributions retain this entire copyright
X * notice and comment, and (2) distributions including binaries display
X * the following acknowledgement: ``This product includes software
X * developed by the University of California, Berkeley and its contributors''
X * in the documentation or other materials provided with the distribution
X * and in all advertising materials mentioning features or use of this
X * software. Neither the name of the University nor the names of its
X * contributors may be used to endorse or promote products derived
X * from this software without specific prior written permission.
X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
X * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
X *
X * All recipients should regard themselves as participants in an ongoing
X * research project and hence should feel obligated to report their
X * experiences (good or bad) with these elementary function codes, using
X * the sendbug(8) program, to the authors.
X *
X * @(#)mathimpl.h 5.2 (Berkeley) 6/1/90
X */
X
#include <math.h>
X
#if defined(__STDC__) || defined(__GNUC__)
#define const const
#else
#define const /**/
#endif
X
#if defined(vax)||defined(tahoe)
X
/* Deal with different ways to concatenate in cpp */
# if defined(__STDC__) || defined(__GNUC__)
# define cat3(a,b,c) a ## b ## c
# else
# define cat3(a,b,c) a/**/b/**/c
# endif
X
/* Deal with vax/tahoe byte order issues */
# ifdef vax
# define cat3t(a,b,c) cat3(a,b,c)
# else
# define cat3t(a,b,c) cat3(a,c,b)
# endif
X
# define vccast(name) (*(const double *)(cat3(name,,x)))
X
X /*
X * Define a constant to high precision on a Vax or Tahoe.
X *
X * Args are the name to define, the decimal floating point value,
X * four 16-bit chunks of the float value in hex
X * (because the vax and tahoe differ in float format!), the power
X * of 2 of the hex-float exponent, and the hex-float mantissa.
X * Most of these arguments are not used at compile time; they are
X * used in a post-check to make sure the constants were compiled
X * correctly.
X *
X * People who want to use the constant will have to do their own
X * #define foo vccast(foo)
X * since CPP cannot do this for them from inside another macro (sigh).
X * We define "vccast" if this needs doing.
X */
# define vc(name, value, x1,x2,x3,x4, bexp, xval) \
X const static long cat3(name,,x)[] = {cat3t(0x,x1,x2), cat3t(0x,x3,x4)};
X
# define ic(name, value, bexp, xval) ;
X
#else /* vax or tahoe */
X
X /* Hooray, we have an IEEE machine */
# undef vccast
# define vc(name, value, x1,x2,x3,x4, bexp, xval) ;
X
# define ic(name, value, bexp, xval) \
X const static double name = value;
X
#endif /* defined(vax)||defined(tahoe) */
SHAR_EOF
chmod 0644 mathimpl.h ||
echo 'restore of mathimpl.h failed'
Wc_c="`wc -c < 'mathimpl.h'`"
test 2996 -eq "$Wc_c" ||
echo 'mathimpl.h: original size 2996, current size' "$Wc_c"
fi
# ============= nextafter.c ==============
if test -f 'nextafter.c' -a X"$1" != X"-c"; then
echo 'x - skipping nextafter.c (File already exists)'
else
echo 'x - extracting nextafter.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'nextafter.c' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
** A mix of C and assembler - well I've got the functions so I might
** as well use them!
**
*/
X
#include "fpumath.h"
X
asm(".align 4");
asm(".Lulp:");
asm(".double 1.1102230246251565e-16");
X
asm(".align 4");
asm(".Lulpup:");
asm(".double 2.2204460492503131e-16");
X
double
nextafter(x, y)
double x, y;
{
X asm("sub $8, %esp");
X
X if (isnan(x) || isnan(y))
X return(quiet_nan());
X
X if (isinf(x))
X if (y > x)
X return(-max_normal());
X else
X if (y < x)
X return(max_normal());
X
X if (x == 0.0)
X if (y > 0.0)
X return(min_subnormal());
X else
X return(-min_subnormal());
X
X if (isnormal(x)) {
X asm("movl 12(%ebp), %eax");
X asm("andl $0x7ff00000, %eax");
X asm("movl %eax, -12(%ebp)");
X asm("movl $0x0, -16(%ebp)");
X asm("fldl -16(%ebp)");
X
X if (fabs(x) <= 1.0)
X asm("fldl .Lulp");
X else
X asm("fldl .Lulpup");
X
X asm("fmulp");
X
X if (y > x) {
X asm("fldl 8(%ebp)");
X asm("faddp");
X asm("leave");
X asm("ret");
X }
X if (y < x) {
X asm("fldl 8(%ebp)");
X asm("fsubp");
X asm("leave");
X asm("ret");
X }
X }
X else
X if (issubnormal(x)) {
X if (y > x)
X return(x + min_subnormal());
X
X if (y < x)
X return(x - min_subnormal());
X }
X
X return(x);
}
SHAR_EOF
chmod 0644 nextafter.c ||
echo 'restore of nextafter.c failed'
Wc_c="`wc -c < 'nextafter.c'`"
test 1392 -eq "$Wc_c" ||
echo 'nextafter.c: original size 1392, current size' "$Wc_c"
fi
true || echo 'restore of paranoia.c failed'
echo End of part 2, continue with part 3
exit 0
--
Glenn Geers | "So when it's over, we're back to people.
Department of Theoretical Physics | Just to prove that human touch can have
The University of Sydney | no equal."
Sydney NSW 2006 Australia | - Basia Trzetrzelewska, 'Prime Time TV'
More information about the Alt.sources
mailing list