Date: Mon, 17 Aug 92 10:49:53 +0200
From: kunst AT prl DOT philips DOT nl (Pieter Kunst)
Subject: pow bug fix
Reply-To: kunst AT prl DOT philips DOT nl
To: djgpp AT sun DOT soe DOT clarkson DOT edu, dj AT ctron DOT com
Cc: ronis AT ronis DOT chem DOT mcgill DOT ca, roberto AT germinal DOT ibp DOT fr

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 AT prl DOT philips DOT 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 AT germinal DOT ibp DOT 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 AT prl DOT philips DOT 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 AT prl DOT philips DOT 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
!