Return to BSD News archive
Path: sserve!manuel!munnari.oz.au!uunet!elroy.jpl.nasa.gov!usc!rpi!newsserver.pixel.kodak.com!laidbak!tellab5!vpnet!serveme!tallboy!news
From: kunst@prl.philips.nl (Pieter Kunst)
Newsgroups: comp.unix.bsd
Subject: pow bug fix
Message-ID: <9208170849.AA17863@hpas1.prl.philips.nl>
Date: Mon, 17 Aug 92 10:49:53 +0200
Sender: news@tallboy.chi.il.us
Reply-To: kunst@prl.philips.nl
Lines: 313
I fixed the pow(-2,3) bug in DJGPP version 1.07 (and below).
As the diff -c3 file is almost as long as the fixed file,
I send both; use whichever is of most convenience...
Pieter Kunst (kunst@prl.philips.nl)
------------ demo program (test.c) below this line --------------------
#include <stdio.h>
extern double pow(double,double);
main()
{
printf ("pow(-2,2) = %g\n", pow(-2,2)); /* 4 */
printf ("pow(-2,3) = %g\n", pow(-2,3)); /* -8 */
printf ("pow(2,3.1) = %g\n", pow(2,3.1)); /* 8.574188 */
printf ("pow(-2,3.1) = %g\n", pow(-2,3.1)); /* DOMAIN error */
}
---------------- demo makefile below this line -----------------------
#
PROG = test
OBJ = test.o exp.o
..c.o:
gcc -c $<
..s.o:
gcc -c $<
$(PROG): $(OBJ)
gcc $(OBJ) -o $(PROG)
---------------- fixed exp.s below this line -------------------------
/* This is file EXP.S */
/*
** Copyright (C) 1991 DJ Delorie, 24 Kirsten Ave, Rochester NH 03867-2954
**
** Modified O.ROBERT 24, Avenue de Verdun 92170 VANVES, FRANCE
**
** E-mail: roberto@germinal.ibp.fr
**
** This file is distributed under the terms listed in the document
** "copying.dj", available from DJ Delorie at the address above.
** A copy of "copying.dj" should accompany this file; if not, a copy
** should be available from where this file was obtained. This file
** may not be distributed without a verbatim copy of "copying.dj".
**
** This file is distributed WITHOUT ANY WARRANTY; without even the implied
** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
**
** or 03-Apr-1991 corrected bug about argument zero to pow
** fyl2x didn't like it
**
** 17-08-92: modified to fix pow(-2,3) bug by kunst@prl.philips.nl
** pow(x,y) x^y. A domain error occurs if x=0 and y<=0,
** or if x<0 and y is not an integer.
**
*/
/* History:15,24 */
.data
yint:
.word 0,0
LCW1:
.word 0
LCW2:
.word 0
.text
LC0:
.double 0d1.0e+00
frac:
fstcww LCW1
fstcww LCW2
fwait
andw $0xf3ff,LCW2
orw $0x0400,LCW2
fldcww LCW2
fldl %st(0)
frndint
fldcww LCW1
fxch %st(1)
fsub %st(1),%st
ret
.globl _pow2
_pow2:
fldl 4(%esp)
Lpow2:
call frac
f2xm1
faddl LC0
fscale
fstp %st(1)
ret
.globl _exp
_exp:
fldl 4(%esp)
fldl2e
fmulp
jmp Lpow2
.globl _pow10
_pow10:
fldl 4(%esp)
fldl2t
fmulp
jmp Lpow2
.globl _pow
_pow:
fldl 12(%esp)
fldl 4(%esp)
ftst
fnstsww %ax
sahf
jbe xltez
fyl2x
jmp Lpow2
xltez:
jb xltz
fstp %st(0)
ftst
fnstsww %ax
sahf
ja ygtz
jb error
fstp %st(0)
fld1l
fchs
error:
fsqrt
ret
ygtz:
fstp %st(0)
fldzl
ret
xltz:
fabs
fxch %st(1)
call frac
ftst
fnstsww %ax
fstp %st(0)
sahf
je yisint
fstp %st(0)
fchs
jmp error
yisint:
fistl yint
fxch %st(1)
fyl2x
call Lpow2
andl $1,yint
jz yeven
fchs
yeven:
ret
---------------- diff -c3 starts below this line ---------------------
*** exp.org Mon Aug 17 10:26:53 1992
--- exp.s Mon Aug 17 10:33:45 1992
***************
*** 17,26 ****
--- 17,33 ----
**
** or 03-Apr-1991 corrected bug about argument zero to pow
** fyl2x didn't like it
+ **
+ ** 17-08-92: modified to fix pow(-2,3) bug by kunst@prl.philips.nl
+ ** pow(x,y) x^y. A domain error occurs if x=0 and y<=0,
+ ** or if x<0 and y is not an integer.
+ **
*/
/* History:15,24 */
.data
+ yint:
+ .word 0,0
LCW1:
.word 0
LCW2:
***************
*** 27,40 ****
.word 0
.text
-
LC0:
.double 0d1.0e+00
! .globl _pow2
! _pow2:
! fldl 4(%esp)
! Lpow2:
fstcww LCW1
fstcww LCW2
fwait
--- 34,43 ----
.word 0
.text
LC0:
.double 0d1.0e+00
! frac:
fstcww LCW1
fstcww LCW2
fwait
***************
*** 46,51 ****
--- 49,61 ----
fldcww LCW1
fxch %st(1)
fsub %st(1),%st
+ ret
+
+ .globl _pow2
+ _pow2:
+ fldl 4(%esp)
+ Lpow2:
+ call frac
f2xm1
faddl LC0
fscale
***************
*** 68,92 ****
.globl _pow
_pow:
fldl 4(%esp)
ftst
! fnstsww %ax
sahf
! je zero1
! fldl 12(%esp)
ftst
! fnstsww %ax
sahf
! je zero2
! fxch %st(1)
fyl2x
! jmp Lpow2
! zero1:
! fstpl 4(%esp)
! fldzl
ret
! zero2:
! fstpl 4(%esp)
! fstpl 4(%esp)
! fld1l
! ret
--- 78,129 ----
.globl _pow
_pow:
+ fldl 12(%esp)
fldl 4(%esp)
+ ftst
+ fnstsww %ax
+ sahf
+ jbe xltez
+ fyl2x
+ jmp Lpow2
+ xltez:
+ jb xltz
+ fstp %st(0)
ftst
! fnstsww %ax
sahf
! ja ygtz
! jb error
! fstp %st(0)
! fld1l
! fchs
! error:
! fsqrt
! ret
! ygtz:
! fstp %st(0)
! fldzl
! ret
! xltz:
! fabs
! fxch %st(1)
! call frac
ftst
! fnstsww %ax
! fstp %st(0)
sahf
! je yisint
! fstp %st(0)
! fchs
! jmp error
! yisint:
! fistl yint
! fxch %st(1)
fyl2x
! call Lpow2
! andl $1,yint
! jz yeven
! fchs
! yeven:
ret
!