Source Code 1992 March

home *** CD-ROM | disk | FTP | other *** search

/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2333 < prev next >

Wrap

Internet Message Format | 1990-12-28 | 54.9 KB

From: glenn@qed.physics.su.OZ.AU (Glenn Geers) Newsgroups: alt.sources Subject: Alternative 80386 math lib v2.0 part04/06 Message-ID: <1990Dec16.210541.28363@metro.ucc.su.OZ.AU> Date: 16 Dec 90 21:05:41 GMT Submitted-by: root@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@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'