From d5a94b37224cfb9f0f3de79330ba5cf9d10bae9d Mon Sep 17 00:00:00 2001
From: Lev Kujawski <int21h@mailbox.org>
Date: Thu, 11 Feb 2021 15:25:51 -0700
Subject: [PATCH] Produce IEEE compliant output from pow() despite platform
 deviations

src/cmd/ksh93/features/math.sh:
- Specify ast_float.h within iffehdrs instead of math.h, so that iffe
  will pick up on macro substitutions within libast. This should make
  any future efforts to remedy floating point behavior easier as well.
- Always include ast_float.h within the generated math header file,
  not just on IA64 platforms.

src/cmd/ksh93/tests/arith.sh:
- Test pow(1.0,-Inf) and pow(1.0,NaN) for IEEE compliance as well.
- Test the exponentiation operator (**) in addition, as streval.c,
  which processes the same, calls pow() separately.

src/lib/libast/features/float:
- Test the IEEE compliance of the underlying math library's pow()
  function and substitute macros producing compliant behavior if
  necessary.
---
 src/cmd/ksh93/features/math.sh |  4 ++--
 src/cmd/ksh93/tests/arith.sh   | 10 ++++++++++
 src/lib/libast/features/float  | 13 +++++++++++++
 3 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/src/cmd/ksh93/features/math.sh b/src/cmd/ksh93/features/math.sh
index 3d9b851ef..e2913d4f7 100644
--- a/src/cmd/ksh93/features/math.sh
+++ b/src/cmd/ksh93/features/math.sh
@@ -29,7 +29,7 @@ esac
 
 command=$0
 iffeflags="-n -v"
-iffehdrs="math.h"
+iffehdrs="ast_float.h"
 iffelibs="-lm"
 table=/dev/null
 
@@ -136,11 +136,11 @@ case $_hdr_ieeefp in
 1)	echo "#include <ieeefp.h>" ;;
 esac
 cat <<!
+#include <ast_float.h>
 #if defined(__ia64__) && defined(signbit)
 # if defined __GNUC__ && __GNUC__ >= 4
 #  define __signbitl(f)		__builtin_signbitl(f)
 # else
-#  include <ast_float.h>
 #  if _lib_copysignl
 #   define __signbitl(f)	(int)(copysignl(1.0,(f))<0.0)
 #  endif
diff --git a/src/cmd/ksh93/tests/arith.sh b/src/cmd/ksh93/tests/arith.sh
index 146baa441..9a0eae0ba 100755
--- a/src/cmd/ksh93/tests/arith.sh
+++ b/src/cmd/ksh93/tests/arith.sh
@@ -297,6 +297,12 @@ fi
 if	(( (4**3)**2 != pow(pow(4,3),2) ))
 then	err_exit '(4**3)**2 not working'
 fi
+if	(( 1**Inf != pow(1,Inf) ))
+then	err_exit '1**Inf not working'
+fi
+if	(( 1**NaN != pow(1,NaN) ))
+then	err_exit '1**NaN not working'
+fi
 typeset -Z3 x=11
 typeset -i x
 if	(( x != 11 ))
@@ -465,8 +471,12 @@ then	set \
 	(( NaN != NaN )) || err_exit 'NaN == NaN'
 	(( -5*Inf == -Inf )) || err_exit '-5*Inf != -Inf'
 	[[ $(print -- $((sqrt(-1.0)))) == ?(-)nan ]]|| err_exit 'sqrt(-1.0) != NaN'
+	(( pow(1.0,-Inf) == 1.0 )) || err_exit 'pow(1.0,-Inf) != 1.0'
+	(( pow(-Inf,0.0) == 1.0 )) || err_exit 'pow(-Inf,0.0) != 1.0'
 	(( pow(1.0,Inf) == 1.0 )) || err_exit 'pow(1.0,Inf) != 1.0'
 	(( pow(Inf,0.0) == 1.0 )) || err_exit 'pow(Inf,0.0) != 1.0'
+	(( pow(1.0,NaN) == 1.0 )) || err_exit 'pow(1.0,NaN) != 1.0'
+	(( pow(Nan,0.0) == 1.0 )) || err_exit 'pow(Nan,0.0) != 1.0'
 	[[ $(print -- $((NaN/Inf))) == ?(-)nan ]] || err_exit 'NaN/Inf != NaN'
 	(( 4.0/Inf == 0.0 )) || err_exit '4.0/Inf != 0.0'
 else	err_exit 'Inf and NaN not working'
diff --git a/src/lib/libast/features/float b/src/lib/libast/features/float
index 781f76e0f..0a65200e7 100644
--- a/src/lib/libast/features/float
+++ b/src/lib/libast/features/float
@@ -1230,3 +1230,16 @@ tst	- -DSCAN=1 - -lm -DSTRTO=1 - -DMAC=1 - -DDIV=1 - -DEXP=1 - -DADD=1 - -DMPY=1
 		return 0;
 	}
 }end
+
+tst	ast_pow_ieee -lm note{ pow(1,inf) is IEEE compliant }end execute{
+	#include <sys/types.h>
+	#include <sys/stat.h>
+	#include <stdlib.h>
+	#include <unistd.h>
+	#include <math.h>
+	int main() { return pow(1.0, 1.0 / 0.0) != 1.0; }
+}end fail{
+	echo '#define powf(x,y)	(((x)==1.0)?1.0:powf((x),(y)))'
+	echo '#define pow(x,y)	(((x)==1.0)?1.0:pow((x),(y)))'
+	echo '#define powl(x,y)	(((x)==1.0)?1.0:powl((x),(y)))'
+}end