Alternative 80386 math lib v2.0 part04/06
Glenn Geers
glenn at qed.physics.su.OZ.AU
Mon Dec 17 08:05:41 AEST 1990
Submitted-by: root at trantor
Archive-name: mathlib2.0/part04
---- Cut Here and feed the following to sh ----
#!/bin/sh
# this is mathlib.04 (part 4 of mathlib2.0)
# do not concatenate these parts, unpack them in order with /bin/sh
# file asinf.s continued
#
if test ! -r _shar_seq_.tmp; then
echo 'Please unpack part 1 first!'
exit 1
fi
(read Scheck
if test "$Scheck" != 4; then
echo Please unpack part "$Scheck" next!
exit 1
else
exit 0
fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
echo 'x - still skipping asinf.s'
else
echo 'x - continuing file asinf.s'
sed 's/^X//' << 'SHAR_EOF' >> 'asinf.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 asinf
asinf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fst %st(1)
X fst %st(2)
X fmulp
X fld1
X fsubp
X fsqrt
X fld1
X fdivp
X fmulp
X fld1
X fpatan
X
X leave
X ret
SHAR_EOF
echo 'File asinf.s is complete' &&
chmod 0644 asinf.s ||
echo 'restore of asinf.s failed'
Wc_c="`wc -c < 'asinf.s'`"
test 403 -eq "$Wc_c" ||
echo 'asinf.s: original size 403, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= coshf.s ==============
if test -f 'coshf.s' -a X"$1" != X"-c"; then
echo 'x - skipping coshf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting coshf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'coshf.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
LC2:
X .double 0.500
X
X .align 4
X .globl coshf
coshf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X ftst
X fstsw %ax
X sahf
X ja .Lpos
X fchs
X
Lpos:
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X fst %st(1)
X
X fld1
X fdivp
X faddp
X
X fldl .LC2
X fmulp
X
X leave
X ret
SHAR_EOF
chmod 0644 coshf.s ||
echo 'restore of coshf.s failed'
Wc_c="`wc -c < 'coshf.s'`"
test 535 -eq "$Wc_c" ||
echo 'coshf.s: original size 535, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= floorf.s ==============
if test -f 'floorf.s' -a X"$1" != X"-c"; then
echo 'x - skipping floorf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting floorf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'floorf.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 floorf
floorf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X flds 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 floorf.s ||
echo 'restore of floorf.s failed'
Wc_c="`wc -c < 'floorf.s'`"
test 630 -eq "$Wc_c" ||
echo 'floorf.s: original size 630, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= log2f.s ==============
if test -f 'log2f.s' -a X"$1" != X"-c"; then
echo 'x - skipping log2f.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting log2f.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'log2f.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 log2f
log2f:
X pushl %ebp
X movl %esp,%ebp
X
X fld1
X flds 8(%ebp)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 log2f.s ||
echo 'restore of log2f.s failed'
Wc_c="`wc -c < 'log2f.s'`"
test 331 -eq "$Wc_c" ||
echo 'log2f.s: original size 331, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= sqrtf.s ==============
if test -f 'sqrtf.s' -a X"$1" != X"-c"; then
echo 'x - skipping sqrtf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting sqrtf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'sqrtf.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
**
** sqrt in the prevailing precision
*/
X
X .align 4
X .globl sqrtf
sqrtf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fsqrt
X
X leave
X ret
SHAR_EOF
chmod 0644 sqrtf.s ||
echo 'restore of sqrtf.s failed'
Wc_c="`wc -c < 'sqrtf.s'`"
test 361 -eq "$Wc_c" ||
echo 'sqrtf.s: original size 361, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= asinhf.s ==============
if test -f 'asinhf.s' -a X"$1" != X"-c"; then
echo 'x - skipping asinhf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting asinhf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'asinhf.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 asinhf
asinhf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X
X fmuls 8(%ebp)
X fld1
X faddp
X fsqrt
X fadds 8(%ebp)
X fldln2
X fxch %st(1)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 asinhf.s ||
echo 'restore of asinhf.s failed'
Wc_c="`wc -c < 'asinhf.s'`"
test 399 -eq "$Wc_c" ||
echo 'asinhf.s: original size 399, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= dremf.s ==============
if test -f 'dremf.s' -a X"$1" != X"-c"; then
echo 'x - skipping dremf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting dremf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'dremf.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 drem
drem:
X pushl %ebp
X movl %esp,%ebp
X
X fldl 16(%ebp)
X fldl 8(%ebp)
Lnotred:
X fprem1
X
X fstsw %ax
X sahf
X jp .Lnotred
X
X leave
X ret
SHAR_EOF
chmod 0644 dremf.s ||
echo 'restore of dremf.s failed'
Wc_c="`wc -c < 'dremf.s'`"
test 381 -eq "$Wc_c" ||
echo 'dremf.s: original size 381, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= fmodf.s ==============
if test -f 'fmodf.s' -a X"$1" != X"-c"; then
echo 'x - skipping fmodf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting fmodf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'fmodf.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 fmodf
fmodf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 12(%ebp)
X flds 8(%ebp)
Lnotred:
X fprem
X
X fstsw %ax
X sahf
X jp .Lnotred
X
X leave
X ret
SHAR_EOF
chmod 0644 fmodf.s ||
echo 'restore of fmodf.s failed'
Wc_c="`wc -c < 'fmodf.s'`"
test 382 -eq "$Wc_c" ||
echo 'fmodf.s: original size 382, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= logbf.s ==============
if test -f 'logbf.s' -a X"$1" != X"-c"; then
echo 'x - skipping logbf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting logbf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'logbf.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 logbf
logbf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fxtract
X fldl %st(1)
X
X leave
X ret
SHAR_EOF
chmod 0644 logbf.s ||
echo 'restore of logbf.s failed'
Wc_c="`wc -c < 'logbf.s'`"
test 341 -eq "$Wc_c" ||
echo 'logbf.s: original size 341, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= sqrtpf.s ==============
if test -f 'sqrtpf.s' -a X"$1" != X"-c"; then
echo 'x - skipping sqrtpf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting sqrtpf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'sqrtpf.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
**
** sqrt in 64 bit precision always
*/
X
X .align 4
X .globl sqrtpf
sqrtpf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X fnstcw -12(%ebp) /* store current control word */
X fnstcw -16(%ebp) /* store again */
X andw $0xf0ff, -12(%ebp) /* enable 32 bit precision */
X fldcw -12(%ebp)
X
X flds 8(%ebp)
X fsqrt
X
X fldcw -16(%ebp) /* restore original control word */
X
X leave
X ret
SHAR_EOF
chmod 0644 sqrtpf.s ||
echo 'restore of sqrtpf.s failed'
Wc_c="`wc -c < 'sqrtpf.s'`"
test 607 -eq "$Wc_c" ||
echo 'sqrtpf.s: original size 607, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= atan2f.s ==============
if test -f 'atan2f.s' -a X"$1" != X"-c"; then
echo 'x - skipping atan2f.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting atan2f.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atan2f.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
Lpi:
X .double 3.14159265358979323846
X
X .align 4
Lmpi:
X .double -3.14159265358979323846
X
X .align 4
Lhalfpi:
X .double 1.57079632679489661923
X
X .align 4
Lmhalfpi:
X .double -1.57079632679489661923
X
X .align 4
X .globl atan2f
atan2f:
X pushl %ebp
X movl %esp,%ebp
X
X flds 12(%ebp)
X ftst
X fnstsw %ax
X sahf
X flds 8(%ebp)
X jz .Lgotzero
X jc .Lgotneg
X
X fdivp
X fld1
X fpatan
X
X leave
X ret
X
Lgotneg:
X ftst
X fnstsw %ax
X sahf
X jc .Lneg1
X
X fdivp
X fld1
X fpatan
X flds .Lmpi
X fsubrp
X
X leave
X ret
X
Lneg1:
X fdivp
X fld1
X fpatan
X flds .Lpi
X fsubrp
X
X leave
X ret
X
Lgotzero:
X ftst
X fnstsw %ax
X sahf
X jz .Lzero
X jc .Lneg
X
X flds .Lhalfpi
X
X leave
X ret
X
Lzero:
X fldz
X
X leave
X ret
X
Lneg:
X flds .Lmhalfpi
X
X leave
X ret
SHAR_EOF
chmod 0644 atan2f.s ||
echo 'restore of atan2f.s failed'
Wc_c="`wc -c < 'atan2f.s'`"
test 933 -eq "$Wc_c" ||
echo 'atan2f.s: original size 933, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= exp10f.s ==============
if test -f 'exp10f.s' -a X"$1" != X"-c"; then
echo 'x - skipping exp10f.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting exp10f.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'exp10f.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 exp10f
exp10f:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fldl2t
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 exp10f.s ||
echo 'restore of exp10f.s failed'
Wc_c="`wc -c < 'exp10f.s'`"
test 406 -eq "$Wc_c" ||
echo 'exp10f.s: original size 406, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= hypotf.s ==============
if test -f 'hypotf.s' -a X"$1" != X"-c"; then
echo 'x - skipping hypotf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting hypotf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'hypotf.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 hypotf
hypotf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fmuls 8(%ebp)
X flds 12(%ebp)
X fmuls 12(%ebp)
X faddp
X fsqrt
X
X leave
X ret
SHAR_EOF
chmod 0644 hypotf.s ||
echo 'restore of hypotf.s failed'
Wc_c="`wc -c < 'hypotf.s'`"
test 379 -eq "$Wc_c" ||
echo 'hypotf.s: original size 379, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= logf.s ==============
if test -f 'logf.s' -a X"$1" != X"-c"; then
echo 'x - skipping logf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting logf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'logf.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 logf
logf:
X pushl %ebp
X movl %esp,%ebp
X
X fldln2
X flds 8(%ebp)
X fyl2x
X
X leave
X ret
SHAR_EOF
chmod 0644 logf.s ||
echo 'restore of logf.s failed'
Wc_c="`wc -c < 'logf.s'`"
test 331 -eq "$Wc_c" ||
echo 'logf.s: original size 331, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= tanf.s ==============
if test -f 'tanf.s' -a X"$1" != X"-c"; then
echo 'x - skipping tanf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting tanf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'tanf.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 tanf
tanf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fptan
X fstp %st(0)
X
X leave
X ret
SHAR_EOF
chmod 0644 tanf.s ||
echo 'restore of tanf.s failed'
Wc_c="`wc -c < 'tanf.s'`"
test 336 -eq "$Wc_c" ||
echo 'tanf.s: original size 336, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= atanf.s ==============
if test -f 'atanf.s' -a X"$1" != X"-c"; then
echo 'x - skipping atanf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting atanf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atanf.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 atanf
atanf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fld1
X fpatan
X
X leave
X ret
SHAR_EOF
chmod 0644 atanf.s ||
echo 'restore of atanf.s failed'
Wc_c="`wc -c < 'atanf.s'`"
test 333 -eq "$Wc_c" ||
echo 'atanf.s: original size 333, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= exp2f.s ==============
if test -f 'exp2f.s' -a X"$1" != X"-c"; then
echo 'x - skipping exp2f.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting exp2f.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'exp2f.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 exp2f
exp2f:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 exp2f.s ||
echo 'restore of exp2f.s failed'
Wc_c="`wc -c < 'exp2f.s'`"
test 389 -eq "$Wc_c" ||
echo 'exp2f.s: original size 389, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= ieee_extf.s ==============
if test -f 'ieee_extf.s' -a X"$1" != X"-c"; then
echo 'x - skipping ieee_extf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting ieee_extf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ieee_extf.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 isnanf
isnanf:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%ebp), %eax
X andl $0x7f800000, %eax
X cmpl $0x7f800000, %eax
X jne .Lnotnan
X movl 8(%ebp), %eax
X andl $0x7ffff, %eax
X cmpl $0x0, %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 isinff
isinff:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%ebp), %eax
X andl $0x7fffffff, %eax
X cmpl $0x7f800000, %eax
X je .Lgotinf
X
X movl $0, %eax
X leave
X ret
X
Lgotinf:
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl iszerof
iszerof:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%ebp), %eax
X cmpl $0x0, %eax
X je .Lgotzero
X
X movl $0, %eax
X leave
X ret
X
Lgotzero:
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl signbitf
signbitf:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%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 issubnormalf
issubnormalf:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%ebp), %eax
X andl $0x7f800000, %eax
X cmpl $0x0, %eax
X je .Lcouldbesub
X
Lnotsubnorm:
X movl $0, %eax
X leave
X ret
X
Lcouldbesub:
X movl 8(%ebp), %eax
X andl $0x7ffff, %eax
X cmpl $0x0, %eax
X je .Lnotsubnorm
X
X movl $1, %eax
X leave
X ret
X
X .align 4
X .globl isnormalf
isnormalf:
X pushl %ebp
X movl %esp,%ebp
X
X movl 8(%ebp), %eax
X andl $0x7f800000, %eax /* mask sign bit */
X xorl $0x7f800000, %eax
X cmpl $0x0, %eax
X je .Lnotnorm
X cmpl $0x7f800000, %eax
X je .Lnotnorm
X
Lnorm:
X movl $1, %eax
X leave
X ret
X
Lnotnorm:
X movl $0, %eax
X leave
X ret
SHAR_EOF
chmod 0644 ieee_extf.s ||
echo 'restore of ieee_extf.s failed'
Wc_c="`wc -c < 'ieee_extf.s'`"
test 1738 -eq "$Wc_c" ||
echo 'ieee_extf.s: original size 1738, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= powf.s ==============
if test -f 'powf.s' -a X"$1" != X"-c"; then
echo 'x - skipping powf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting powf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'powf.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.
**
** This file started life as a C implemenation of pow, it then evolved
** into an embeded asm version (with some C) and finally, all assembler.
** It's actually hacked assembler produced by gcc.
**
** Copyright 1990 G. Geers
**
*/
X
text
LC0:
X .ascii "pow: DOMAIN ERROR\12\0"
X .align 4
LC1:
X .long 0x0
X .align 4
LC2:
X .long 0x3f800000
X .align 4
globl powf
powf:
X pushl %ebp
X movl %esp,%ebp
X pushl %esi
X pushl %ebx
X xorl %esi,%esi
X movl $1,%ebx
X flds 8(%ebp)
X ftst
X fstp %st(0)
X fnstsw %ax
X sahf
X jae .L2
X movl %ebx,%esi
X pushl 12(%ebp)
X
/**/
X flds 12(%ebp)
X
X fstcw -8(%ebp)
X fstcw -12(%ebp)
X orw $0x0400, -12(%ebp)
X fldcw -12(%ebp)
X
X frndint
X fldcw -8(%ebp)
/**/
X fsubrs 12(%ebp)
X addl $4,%esp
X ftst
X fstp %st(0)
X fnstsw %ax
X sahf
X je .L3
X movl $33,errno
X pushl $.LC0
X pushl $_iob+32
X call fprintf
X flds .LC1
X jmp .L1
X .align 4
L3:
X flds 12(%ebp)
/**/
X fstpl -8(%ebp)
X
X fldl -8(%ebp)
X fistpl -8(%ebp)
X movl -8(%ebp),%eax
/**/
X movl %eax,%edx
X testl %eax,%eax
X jge .L5
X incl %eax
L5:
X andl $-2,%eax
X subl %eax,%edx
X movl %edx,%eax
X testl %eax,%eax
X je .L2
X xorl %ebx,%ebx
L2:
X flds 8(%ebp)
X ftst
X fstp %st(0)
X fnstsw %ax
X sahf
X jne .L6
X flds 12(%ebp)
X ftst
X fstp %st(0)
X fnstsw %ax
X sahf
X je .L6
X flds .LC1
X jmp .L1
X .align 4
L6:
X flds 12(%ebp)
X ftst
X fstp %st(0)
X fnstsw %ax
X sahf
X jne .L8
X flds .LC2
X jmp .L1
X .align 4
L8:
/APP
X flds 12(%ebp)
X flds 8(%ebp)
/NO_APP
X testl %esi,%esi
X je .L10
/APP
X fchs
/NO_APP
L10:
/APP
X fyl2x
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubr %st(1)
X f2xm1
X fld1
X faddp
X fld %st(2)
X fstp %st(2)
X fscale
/NO_APP
X testl %ebx,%ebx
X jne .L1
/APP
X fchs
/NO_APP
L1:
X leal -8(%ebp),%esp
X popl %ebx
X popl %esi
X leave
X ret
SHAR_EOF
chmod 0644 powf.s ||
echo 'restore of powf.s failed'
Wc_c="`wc -c < 'powf.s'`"
test 1850 -eq "$Wc_c" ||
echo 'powf.s: original size 1850, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= tanhf.s ==============
if test -f 'tanhf.s' -a X"$1" != X"-c"; then
echo 'x - skipping tanhf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting tanhf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'tanhf.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 tanhf
tanhf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X flds 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X fstl %st(1)
X fsts -16(%ebp)
X
X fld1
X fdivp
X
X fsubr
X fadds -16(%ebp)
X fdivrp
X
X leave
X ret
SHAR_EOF
chmod 0644 tanhf.s ||
echo 'restore of tanhf.s failed'
Wc_c="`wc -c < 'tanhf.s'`"
test 495 -eq "$Wc_c" ||
echo 'tanhf.s: original size 495, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= atanhf.s ==============
if test -f 'atanhf.s' -a X"$1" != X"-c"; then
echo 'x - skipping atanhf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting atanhf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'atanhf.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
Lhalf:
X .double 0.500
X
X .align 4
X .globl atanhf
atanhf:
X pushl %ebp
X movl %esp,%ebp
X
X fld1
X fadds 8(%ebp)
X fld1
X fsubs 8(%ebp)
X fdivrp
X
X fldln2
X fxch %st(1)
X fyl2x
X
X flds .Lhalf
X fmulp
X
X leave
X ret
SHAR_EOF
chmod 0644 atanhf.s ||
echo 'restore of atanhf.s failed'
Wc_c="`wc -c < 'atanhf.s'`"
test 440 -eq "$Wc_c" ||
echo 'atanhf.s: original size 440, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= expf.s ==============
if test -f 'expf.s' -a X"$1" != X"-c"; then
echo 'x - skipping expf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting expf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'expf.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 expf
expf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 expf.s ||
echo 'restore of expf.s failed'
Wc_c="`wc -c < 'expf.s'`"
test 402 -eq "$Wc_c" ||
echo 'expf.s: original size 402, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= ieee_valuesf.s ==============
if test -f 'ieee_valuesf.s' -a X"$1" != X"-c"; then
echo 'x - skipping ieee_valuesf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting ieee_valuesf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ieee_valuesf.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_normalf
X
max_normalf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7f7fffff, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
X
X .align 4
X .globl min_normalf
X
min_normalf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x00800001, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
X
X .align 4
X .globl min_subnormalf
X
min_subnormalf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x00000001, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
X
X .align 4
X .globl max_subnormalf
X
max_subnormalf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x007fffff, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
X
X .align 4
X .globl quiet_nanf
X
quiet_nanf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7fffffff, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
X
X .align 4
X .globl signaling_nanf
X
signaling_nanf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7f800001, -8(%ebp)
X
X fldl -8(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 ieee_valuesf.s ||
echo 'restore of ieee_valuesf.s failed'
Wc_c="`wc -c < 'ieee_valuesf.s'`"
test 1128 -eq "$Wc_c" ||
echo 'ieee_valuesf.s: original size 1128, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= rintf.s ==============
if test -f 'rintf.s' -a X"$1" != X"-c"; then
echo 'x - skipping rintf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting rintf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'rintf.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 rintf
rintf:
X pushl %ebp
X movl %esp,%ebp
X
X flds 8(%ebp)
X frndint
X
X leave
X ret
SHAR_EOF
chmod 0644 rintf.s ||
echo 'restore of rintf.s failed'
Wc_c="`wc -c < 'rintf.s'`"
test 328 -eq "$Wc_c" ||
echo 'rintf.s: original size 328, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= ceilf.s ==============
if test -f 'ceilf.s' -a X"$1" != X"-c"; then
echo 'x - skipping ceilf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting ceilf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'ceilf.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.
**
** Could use ceil(x) = -floor(-x) but this is quicker.
**
** Copyright 1990 G. Geers
**
*/
X
X .align 4
X .globl ceilf
ceilf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X flds 8(%ebp) /* load data */
X
X fstcw -12(%ebp) /* store control word */
X fstcw -16(%ebp) /* store it again */
X orw $0x0800, -16(%ebp) /* round toward +inf */
X fldcw -16(%ebp) /* store new control word */
X frndint /* rounding gives ceil(x) */
X fldcw -12(%ebp) /* restore original control word */
X
X leave
X ret
SHAR_EOF
chmod 0644 ceilf.s ||
echo 'restore of ceilf.s failed'
Wc_c="`wc -c < 'ceilf.s'`"
test 684 -eq "$Wc_c" ||
echo 'ceilf.s: original size 684, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= expm1f.s ==============
if test -f 'expm1f.s' -a X"$1" != X"-c"; then
echo 'x - skipping expm1f.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting expm1f.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'expm1f.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 expm1
expm1:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X flds 8(%ebp)
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X fstpl -8(%ebp)
X
X flds 8(%ebp)
X fchs
X fldl2e
X fmulp
X fstl %st(1)
X frndint
X fstl %st(2)
X fsubrp
X f2xm1
X fld1
X faddp
X fscale
X
X fld1
X fsubp
X fmull -8(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 expm1f.s ||
echo 'restore of expm1f.s failed'
Wc_c="`wc -c < 'expm1f.s'`"
test 573 -eq "$Wc_c" ||
echo 'expm1f.s: original size 573, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= infinityf.s ==============
if test -f 'infinityf.s' -a X"$1" != X"-c"; then
echo 'x - skipping infinityf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting infinityf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'infinityf.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 infinityf
infinityf:
X pushl %ebp
X movl %esp,%ebp
X subl $8, %esp
X
X movl $0x7f800000, -8(%ebp)
X
X flds -8(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 infinityf.s ||
echo 'restore of infinityf.s failed'
Wc_c="`wc -c < 'infinityf.s'`"
test 372 -eq "$Wc_c" ||
echo 'infinityf.s: original size 372, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= scalbf.s ==============
if test -f 'scalbf.s' -a X"$1" != X"-c"; then
echo 'x - skipping scalbf.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting scalbf.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'scalbf.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 scalbf
scalbf:
X pushl %ebp
X movl %esp,%ebp
X
X fildl 12(%ebp)
X flds 8(%ebp)
X fscale
X
X leave
X ret
SHAR_EOF
chmod 0644 scalbf.s ||
echo 'restore of scalbf.s failed'
Wc_c="`wc -c < 'scalbf.s'`"
test 345 -eq "$Wc_c" ||
echo 'scalbf.s: original size 345, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
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)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
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
extern short _getsw();
X
ieee_retrospective(f)
FILE *f;
{
X short 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 879 -eq "$Wc_c" ||
echo 'ieee_retro.c: original size 879, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= _getsw.s ==============
if test -f '_getsw.s' -a X"$1" != X"-c"; then
echo 'x - skipping _getsw.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting _getsw.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > '_getsw.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 _getsw
X
_getsw:
X pushl %ebp
X movl %esp,%ebp
X
X fstsw %ax
X
X leave
X ret
SHAR_EOF
chmod 0644 _getsw.s ||
echo 'restore of _getsw.s failed'
Wc_c="`wc -c < '_getsw.s'`"
test 320 -eq "$Wc_c" ||
echo '_getsw.s: original size 320, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= setcont.s ==============
if test -f 'setcont.s' -a X"$1" != X"-c"; then
echo 'x - skipping setcont.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting setcont.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'setcont.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 setcont
X
setcont:
X pushl %ebp
X movl %esp, %ebp
X subl $4, %esp
X
X fnstcw -8(%ebp)
X movw -8(%ebp), %ax
X fldcw 8(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 setcont.s ||
echo 'restore of setcont.s failed'
Wc_c="`wc -c < 'setcont.s'`"
test 377 -eq "$Wc_c" ||
echo 'setcont.s: original size 377, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= setinternal.s ==============
if test -f 'setinternal.s' -a X"$1" != X"-c"; then
echo 'x - skipping setinternal.s (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting setinternal.s (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'setinternal.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 setinternal
X
setinternal:
X pushl %ebp
X movl %esp, %ebp
X subl $4, %esp
X
X fnstcw -8(%ebp)
X movw -8(%ebp), %ax
X orw $0x007e, -8(%ebp) /* All exceptions except IOP masked */
X fldcw -8(%ebp)
X
X leave
X ret
SHAR_EOF
chmod 0644 setinternal.s ||
echo 'restore of setinternal.s failed'
Wc_c="`wc -c < 'setinternal.s'`"
test 449 -eq "$Wc_c" ||
echo 'setinternal.s: original size 449, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= paranoia.c ==============
if test -f 'paranoia.c' -a X"$1" != X"-c"; then
echo 'x - skipping paranoia.c (File already exists)'
rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting paranoia.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'paranoia.c' &&
/* A C version of Kahan's Floating Point Test "Paranoia"
X
X Thos Sumner, UCSF, Feb. 1985
X David Gay, BTL, Jan. 1986
X
X This is a rewrite from the Pascal version by
X
X B. A. Wichmann, 18 Jan. 1985
X
X (and does NOT exhibit good C programming style).
X
(C) Apr 19 1983 in BASIC version by:
X Professor W. M. Kahan,
X 567 Evans Hall
X Electrical Engineering & Computer Science Dept.
X University of California
X Berkeley, California 94720
X USA
X
converted to Pascal by:
X B. A. Wichmann
X National Physical Laboratory
X Teddington Middx
X TW11 OLW
X UK
X
converted to C by:
X
X David M. Gay and Thos Sumner
X AT&T Bell Labs Computer Center, Rm. U-76
X 600 Mountain Avenue University of California
X Murray Hill, NJ 07974 San Francisco, CA 94143
X USA USA
X
with simultaneous corrections to the Pascal source (reflected
in the Pascal source available over netlib).
[A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
X
Reports of results on various systems from all the versions
of Paranoia are being collected by Richard Karpinski at the
same address as Thos Sumner. This includes sample outputs,
bug reports, and criticisms.
X
You may copy this program freely if you acknowledge its source.
Comments on the Pascal version to NPL, please.
X
X
The C version catches signals from floating-point exceptions.
If signal(SIGFPE,...) is unavailable in your environment, you may
#define NOSIGNAL to comment out the invocations of signal.
X
This source file is too big for some C compilers, but may be split
into pieces. Comments containing "SPLIT" suggest convenient places
for this splitting. At the end of these comments is an "ed script"
(for the UNIX(tm) editor ed) that will do this splitting.
X
By #defining Single when you compile this source, you may obtain
a single-precision C version of Paranoia.
X
X
The following is from the introductory commentary from Wichmann's work:
X
The BASIC program of Kahan is written in Microsoft BASIC using many
facilities which have no exact analogy in Pascal. The Pascal
version below cannot therefore be exactly the same. Rather than be
a minimal transcription of the BASIC program, the Pascal coding
follows the conventional style of block-structured languages. Hence
the Pascal version could be useful in producing versions in other
structured languages.
X
Rather than use identifiers of minimal length (which therefore have
little mnemonic significance), the Pascal version uses meaningful
identifiers as follows [Note: A few changes have been made for C]:
X
X
BASIC C BASIC C BASIC C
X
X A J S StickyBit
X A1 AInverse J0 NoErrors T
X B Radix [Failure] T0 Underflow
X B1 BInverse J1 NoErrors T2 ThirtyTwo
X B2 RadixD2 [SeriousDefect] T5 OneAndHalf
X B9 BMinusU2 J2 NoErrors T7 TwentySeven
X C [Defect] T8 TwoForty
X C1 CInverse J3 NoErrors U OneUlp
X D [Flaw] U0 UnderflowThreshold
X D4 FourD K PageNo U1
X E0 L Milestone U2
X E1 M V
X E2 Exp2 N V0
X E3 N1 V8
X E5 MinSqEr O Zero V9
X E6 SqEr O1 One W
X E7 MaxSqEr O2 Two X
X E8 O3 Three X1
X E9 O4 Four X8
X F1 MinusOne O5 Five X9 Random1
X F2 Half O8 Eight Y
X F3 Third O9 Nine Y1
X F6 P Precision Y2
X F9 Q Y9 Random2
X G1 GMult Q8 Z
X G2 GDiv Q9 Z0 PseudoZero
X G3 GAddSub R Z1
X H R1 RMult Z2
X H1 HInverse R2 RDiv Z9
X I R3 RAddSub
X IO NoTrials R4 RSqrt
X I3 IEEE R9 Random9
X
X SqRWrng
X
All the variables in BASIC are true variables and in consequence,
the program is more difficult to follow since the "constants" must
be determined (the glossary is very helpful). The Pascal version
uses Real constants, but checks are added to ensure that the values
are correctly converted by the compiler.
X
The major textual change to the Pascal version apart from the
identifiersis that named procedures are used, inserting parameters
wherehelpful. New procedures are also introduced. The
correspondence is as follows:
X
X
BASIC Pascal
lines
X
X 90- 140 Pause
X 170- 250 Instructions
X 380- 460 Heading
X 480- 670 Characteristics
X 690- 870 History
2940-2950 Random
3710-3740 NewD
4040-4080 DoesYequalX
4090-4110 PrintIfNPositive
4640-4850 TestPartialUnderflow
X
=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
X
Below is an "ed script" that splits para.c into 10 files
of the form part[1-8].c, subs.c, and msgs.c, plus a header
file, paranoia.h, that these files require.
X
r paranoia.c
$
?SPLIT
+,$w msgs.c
X .,$d
?SPLIT
X .d
+d
-,$w subs.c
-,$d
?part8
+d
?include
X .,$w part8.c
X .,$d
-d
?part7
+d
?include
X .,$w part7.c
X .,$d
-d
?part6
+d
?include
X .,$w part6.c
X .,$d
-d
?part5
+d
?include
X .,$w part5.c
X .,$d
-d
?part4
+d
?include
X .,$w part4.c
X .,$d
-d
?part3
+d
?include
X .,$w part3.c
X .,$d
-d
?part2
+d
?include
X .,$w part2.c
X .,$d
?SPLIT
X .d
1,/^#include/-1d
1,$w part1.c
/Computed constants/,$d
1,$s/^int/extern &/
1,$s/^FLOAT/extern &/
1,$s/^char/extern &/
1,$s! = .*!;!
/^Guard/,/^Round/s/^/extern /
/^jmp_buf/s/^/extern /
/^Sig_type/s/^/extern /
s/$/\
extern void sigfpe();/
w paranoia.h
q
X
*/
X
#include <stdio.h>
#ifndef NOSIGNAL
#include <signal.h>
#endif
#include <setjmp.h>
X
#include "fpumath.h"
X
#ifdef Single
#define FLOAT float
#define FABS(x) fabsf(x)
#define FLOOR(x) floorf(x)
#define LOG(x) logf(x)
#define POW(x,y) powf(x,y)
#ifdef TEST
#define SQRT(x) sqrtpf(x)
#else
#define SQRT(x) sqrtf(x)
#endif
#else
#define FLOAT double
#define FABS(x) fabs(x)
#define FLOOR(x) floor(x)
#define LOG(x) log(x)
#define POW(x,y) pow(x,y)
#ifdef TEST
#define SQRT(x) sqrtp(x)
#else
#define SQRT(x) sqrt(x)
#endif
#endif
X
jmp_buf ovfl_buf;
typedef void (*Sig_type)();
Sig_type sigsave;
X
#define KEYBOARD 0
X
FLOAT Radix, BInvrse, RadixD2, BMinusU2;
FLOAT Sign(), Random();
X
/*Small floating point constants.*/
FLOAT Zero = 0.0;
FLOAT Half = 0.5;
FLOAT One = 1.0;
FLOAT Two = 2.0;
FLOAT Three = 3.0;
FLOAT Four = 4.0;
FLOAT Five = 5.0;
FLOAT Eight = 8.0;
FLOAT Nine = 9.0;
FLOAT TwentySeven = 27.0;
FLOAT ThirtyTwo = 32.0;
FLOAT TwoForty = 240.0;
FLOAT MinusOne = -1.0;
FLOAT OneAndHalf = 1.5;
/*Integer constants*/
int NoTrials = 20; /*Number of tests for commutativity. */
#define False 0
#define True 1
X
/* Definitions for declared types
X Guard == (Yes, No);
X Rounding == (Chopped, Rounded, Other);
X Message == packed array [1..40] of char;
X Class == (Flaw, Defect, Serious, Failure);
X */
#define Yes 1
#define No 0
#define Chopped 2
#define Rounded 1
#define Other 0
#define Flaw 3
#define Defect 2
#define Serious 1
#define Failure 0
typedef int Guard, Rounding, Class;
typedef char Message;
X
/* Declarations of Variables */
int Indx;
char ch[8];
FLOAT AInvrse, A1;
FLOAT C, CInvrse;
FLOAT D, FourD;
FLOAT E0, E1, Exp2, E3, MinSqEr;
FLOAT SqEr, MaxSqEr, E9;
FLOAT Third;
FLOAT F6, F9;
FLOAT H, HInvrse;
int I;
FLOAT StickyBit, J;
FLOAT MyZero;
FLOAT Precision;
FLOAT Q, Q9;
FLOAT R, Random9;
FLOAT T, Underflow, S;
FLOAT OneUlp, UfThold, U1, U2;
FLOAT V, V0, V9;
FLOAT W;
FLOAT X, X1, X2, X8, Random1;
FLOAT Y, Y1, Y2, Random2;
FLOAT Z, PseudoZero, Z1, Z2, Z9;
int ErrCnt[4];
int fpecount;
int Milestone;
int PageNo;
int M, N, N1;
Guard GMult, GDiv, GAddSub;
Rounding RMult, RDiv, RAddSub, RSqrt;
int Break, Done, NotMonot, Monot, Anomaly, IEEE,
X SqRWrng, UfNGrad;
/* Computed constants. */
/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
X
/* floating point exception receiver */
X void
sigfpe()
{
X fpecount++;
X printf("\n* * * FLOATING-POINT ERROR * * *\n");
X fflush(stdout);
X if (sigsave) {
#ifndef NOSIGNAL
X signal(SIGFPE, sigsave);
#endif
X sigsave = 0;
X longjmp(ovfl_buf, 1);
X }
X abort();
}
X
main()
{
/*
** Modified by G. Geers - glenn at qed.physics.su.oz.au
**
** Define TEST if you want to include my code
**
*/
X
#ifdef TEST
#ifdef Single
X setcont(MASK_ALL0);
#else
X setcont(MASK_ALL);
#endif
#endif
X
X /* First two assignments use integer right-hand sides. */
X Zero = 0;
X One = 1;
X Two = One + One;
X Three = Two + One;
X Four = Three + One;
X Five = Four + One;
X Eight = Four + Four;
X Nine = Three * Three;
X TwentySeven = Nine * Three;
X ThirtyTwo = Four * Eight;
X TwoForty = Four * Five * Three * Four;
X MinusOne = -One;
X Half = One / Two;
X OneAndHalf = One + Half;
X ErrCnt[Failure] = 0;
X ErrCnt[Serious] = 0;
X ErrCnt[Defect] = 0;
X ErrCnt[Flaw] = 0;
X PageNo = 1;
X /*=============================================*/
X Milestone = 0;
X /*=============================================*/
#ifndef NOSIGNAL
X signal(SIGFPE, sigfpe);
#endif
X Instructions();
X Pause();
X Heading();
X Pause();
X Characteristics();
X Pause();
X History();
X Pause();
X /*=============================================*/
X Milestone = 7;
X /*=============================================*/
X printf("Program is now RUNNING tests on small integers:\n");
X
X TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
X && (One > Zero) && (One + One == Two),
X "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
X Z = - Zero;
X if (Z != 0.0) {
X ErrCnt[Failure] = ErrCnt[Failure] + 1;
X printf("Comparison alleges that -0.0 is Non-zero!\n");
X U1 = 0.001;
X Radix = 1;
X TstPtUf();
X }
X TstCond (Failure, (Three == Two + One) && (Four == Three + One)
X && (Four + Two * (- Two) == Zero)
X && (Four - Three - One == Zero),
X "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
X TstCond (Failure, (MinusOne == (0 - One))
X && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
X && (MinusOne + FABS(One) == Zero)
X && (MinusOne + MinusOne * MinusOne == Zero),
X "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
X TstCond (Failure, Half + MinusOne + Half == Zero,
X "1/2 + (-1) + 1/2 != 0");
X /*=============================================*/
X /*SPLIT
X part2();
X part3();
X part4();
X part5();
X part6();
X part7();
X part8();
X
X }
#include "paranoia.h"
part2(){
*/
X Milestone = 10;
X /*=============================================*/
X TstCond (Failure, (Nine == Three * Three)
X && (TwentySeven == Nine * Three) && (Eight == Four + Four)
X && (ThirtyTwo == Eight * Four)
X && (ThirtyTwo - TwentySeven - Four - One == Zero),
X "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
X TstCond (Failure, (Five == Four + One) &&
X (TwoForty == Four * Five * Three * Four)
X && (TwoForty / Three - Four * Four * Five == Zero)
X && ( TwoForty / Four - Five * Three * Four == Zero)
X && ( TwoForty / Five - Four * Three * Four == Zero),
X "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
X if (ErrCnt[Failure] == 0) {
X printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
X printf("\n");
X }
X printf("Searching for Radix and Precision.\n");
X W = One;
X do {
X W = W + W;
X Y = W + One;
X Z = Y - W;
X Y = Z - One;
X } while (MinusOne + FABS(Y) < Zero);
X /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
X Precision = Zero;
X Y = One;
X do {
X Radix = W + Y;
X Y = Y + Y;
X Radix = Radix - W;
X } while ( Radix == Zero);
X if (Radix < Two) Radix = One;
X printf("Radix = %f .\n", Radix);
X if (Radix != 1) {
X W = One;
X do {
X Precision = Precision + One;
X W = W * Radix;
X Y = W + One;
X } while ((Y - W) == One);
X }
X /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
X ...*/
X U1 = One / W;
X U2 = Radix * U1;
X printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
X printf("Recalculating radix and precision\n ");
X
X /*save old values*/
X E0 = Radix;
X E1 = U1;
X E9 = U2;
X E3 = Precision;
X
X X = Four / Three;
X Third = X - One;
X F6 = Half - Third;
X X = F6 + F6;
X X = FABS(X - Third);
X if (X < U2) X = U2;
X
X /*... now X = (unknown no.) ulps of 1+...*/
X do {
X U2 = X;
X Y = Half * U2 + ThirtyTwo * U2 * U2;
X Y = One + Y;
X X = Y - One;
X } while ( ! ((U2 <= X) || (X <= Zero)));
X
X /*... now U2 == 1 ulp of 1 + ... */
X X = Two / Three;
X F6 = X - Half;
X Third = F6 + F6;
X X = Third - Half;
X X = FABS(X + F6);
X if (X < U1) X = U1;
X
X /*... now X == (unknown no.) ulps of 1 -... */
X do {
X U1 = X;
X Y = Half * U1 + ThirtyTwo * U1 * U1;
X Y = Half - Y;
X X = Half + Y;
X Y = Half - X;
X X = Half + Y;
X } while ( ! ((U1 <= X) || (X <= Zero)));
X /*... now U1 == 1 ulp of 1 - ... */
X if (U1 == E1) printf("confirms closest relative separation U1 .\n");
X else printf("gets better closest relative separation U1 = %.7e .\n", U1);
X W = One / U1;
X F9 = (Half - U1) + Half;
X Radix = FLOOR(0.01 + U2 / U1);
X if (Radix == E0) printf("Radix confirmed.\n");
X else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
X TstCond (Defect, Radix <= Eight + Eight,
X "Radix is too big: roundoff problems");
X TstCond (Flaw, (Radix == Two) || (Radix == 10)
X || (Radix == One), "Radix is not as good as 2 or 10");
X /*=============================================*/
X Milestone = 20;
X /*=============================================*/
X TstCond (Failure, F9 - Half < Half,
X "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
X X = F9;
X I = 1;
X Y = X - Half;
X Z = Y - Half;
X TstCond (Failure, (X != One)
X || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
X X = One + U2;
X I = 0;
X /*=============================================*/
X Milestone = 25;
X /*=============================================*/
X /*... BMinusU2 = nextafter(Radix, 0) */
X BMinusU2 = Radix - One;
X BMinusU2 = (BMinusU2 - U2) + One;
X /* Purify Integers */
X if (Radix != One) {
X X = - TwoForty * LOG(U1) / LOG(Radix);
X Y = FLOOR(Half + X);
X if (FABS(X - Y) * Four < One) X = Y;
X Precision = X / TwoForty;
X Y = FLOOR(Half + Precision);
X if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
X }
X if ((Precision != FLOOR(Precision)) || (Radix == One)) {
X printf("Precision cannot be characterized by an Integer number\n");
X printf("of significant digits but, by itself, this is a minor flaw.\n");
X }
X if (Radix == One)
X printf("logarithmic encoding has precision characterized solely by U1.\n");
X else printf("The number of significant digits of the Radix is %f .\n",
X Precision);
X TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
X "Precision worse than 5 decimal figures ");
X /*=============================================*/
X Milestone = 30;
X /*=============================================*/
X /* Test for extra-precise subepressions */
X X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
X do {
X Z2 = X;
X X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
X } while ( ! ((Z2 <= X) || (X <= Zero)));
X X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
X do {
X Z1 = Z;
X Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
X + One / Two)) + One / Two;
X } while ( ! ((Z1 <= Z) || (Z <= Zero)));
X do {
X do {
X Y1 = Y;
X Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
X )) + Half;
X } while ( ! ((Y1 <= Y) || (Y <= Zero)));
X X1 = X;
X X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
X } while ( ! ((X1 <= X) || (X <= Zero)));
X if ((X1 != Y1) || (X1 != Z1)) {
X BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
X printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
X printf("are symptoms of inconsistencies introduced\n");
X printf("by extra-precise evaluation of arithmetic subexpressions.\n");
X notify("Possibly some part of this");
X if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
X "That feature is not tested further by this program.\n") ;
X }
X else {
X if ((Z1 != U1) || (Z2 != U2)) {
X if ((Z1 >= U1) || (Z2 >= U2)) {
X BadCond(Failure, "");
X notify("Precision");
X printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
X printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
X }
X else {
X if ((Z1 <= Zero) || (Z2 <= Zero)) {
X printf("Because of unusual Radix = %f", Radix);
X printf(", or exact rational arithmetic a result\n");
X printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
X notify("of an\nextra-precision");
X }
X if (Z1 != Z2 || Z1 > Zero) {
X X = Z1 / U1;
X Y = Z2 / U2;
X if (Y > X) X = Y;
X Q = - LOG(X);
X printf("Some subexpressions appear to be calculated extra\n");
X printf("precisely with about %g extra B-digits, i.e.\n",
X (Q / LOG(Radix)));
X printf("roughly %g extra significant decimals.\n",
X Q / LOG(10.));
X }
X printf("That feature is not tested further by this program.\n");
X }
X }
X }
X Pause();
X /*=============================================*/
X /*SPLIT
X }
#include "paranoia.h"
part3(){
*/
X Milestone = 35;
X /*=============================================*/
X if (Radix >= Two) {
X X = W / (Radix * Radix);
X Y = X + One;
X Z = Y - X;
X T = Z + U2;
X X = T - Z;
X TstCond (Failure, X == U2,
X "Subtraction is not normalized X=Y,X+Z != Y+Z!");
X if (X == U2) printf(
X "Subtraction appears to be normalized, as it should be.");
X }
X printf("\nChecking for guard digit in *, /, and -.\n");
X Y = F9 * One;
X Z = One * F9;
X X = F9 - Half;
X Y = (Y - Half) - X;
X Z = (Z - Half) - X;
X X = One + U2;
X T = X * Radix;
X R = Radix * X;
X X = T - Radix;
X X = X - Radix * U2;
X T = R - Radix;
X T = T - Radix * U2;
X X = X * (Radix - One);
X T = T * (Radix - One);
X if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
X else {
X GMult = No;
X TstCond (Serious, False,
X "* lacks a Guard Digit, so 1*X != X");
X }
X Z = Radix * U2;
X X = One + Z;
X Y = FABS((X + Z) - X * X) - U2;
X X = One - U2;
X Z = FABS((X - U2) - X * X) - U1;
X TstCond (Failure, (Y <= Zero)
X && (Z <= Zero), "* gets too many final digits wrong.\n");
X Y = One - U2;
X X = One + U2;
X Z = One / Y;
X Y = Z - X;
X X = One / Three;
X Z = Three / Nine;
X X = X - Z;
X T = Nine / TwentySeven;
X Z = Z - T;
X TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
X "Division lacks a Guard Digit, so error can exceed 1 ulp\n\
or 1/3 and 3/9 and 9/27 may disagree");
X Y = F9 / One;
X X = F9 - Half;
X Y = (Y - Half) - X;
X X = One + U2;
X T = X / One;
X X = T - X;
X if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
X else {
X GDiv = No;
X TstCond (Serious, False,
X "Division lacks a Guard Digit, so X/1 != X");
X }
X X = One / (One + U2);
X Y = X - Half - Half;
X TstCond (Serious, Y < Zero,
X "Computed value of 1/1.000..1 >= 1");
X X = One - U2;
X Y = One + Radix * U2;
X Z = X * Radix;
X T = Y * Radix;
X R = Z / Radix;
X StickyBit = T / Radix;
X X = R - X;
X Y = StickyBit - Y;
X TstCond (Failure, X == Zero && Y == Zero,
X "* and/or / gets too many last digits wrong");
X Y = One - U1;
X X = One - F9;
X Y = One - Y;
X T = Radix - U2;
X Z = Radix - BMinusU2;
X T = Radix - T;
X if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
X else {
X GAddSub = No;
X TstCond (Serious, False,
X "- lacks Guard Digit, so cancellation is obscured");
X }
X if (F9 != One && F9 - One >= Zero) {
X BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
X printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
X printf(" such precautions against division by zero as\n");
X printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
X }
X if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
X " *, /, and - appear to have guard digits, as they should.\n");
X /*=============================================*/
X Milestone = 40;
X /*=============================================*/
X Pause();
X printf("Checking rounding on multiply, divide and add/subtract.\n");
X RMult = Other;
X RDiv = Other;
X RAddSub = Other;
X RadixD2 = Radix / Two;
X A1 = Two;
X Done = False;
X do {
X AInvrse = Radix;
X do {
X X = AInvrse;
X AInvrse = AInvrse / A1;
X } while ( ! (FLOOR(AInvrse) != AInvrse));
X Done = (X == One) || (A1 > Three);
X if (! Done) A1 = Nine + One;
X } while ( ! (Done));
X if (X == One) A1 = Radix;
X AInvrse = One / A1;
X X = A1;
X Y = AInvrse;
X Done = False;
X do {
X Z = X * Y - Half;
X TstCond (Failure, Z == Half,
X "X * (1/X) differs from 1");
X Done = X == Radix;
X X = Radix;
X Y = One / X;
X } while ( ! (Done));
X Y2 = One + U2;
X Y1 = One - U2;
SHAR_EOF
true || echo 'restore of paranoia.c failed'
fi
echo 'End of mathlib2.0 part 4'
echo 'File paranoia.c is continued in part 5'
echo 5 > _shar_seq_.tmp
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