[9front] [PATCH] APE: add ceilf, fabsf and lroundf

[9front] [PATCH] APE: add ceilf, fabsf and lroundf

[telephil9]

Hi,

Here is a patch to add ceilf, fabsf and lroundf to APE math library.
These functions are needed for the netsurf port.

Source code was copied from sourceware newlib. Looking at implementations from various OSes, I see that OpenBSD uses these exact same versions copied from the same sources.

Thanks
--phil

diff -r ce98610ce572 sys/include/ape/math.h
--- a/sys/include/ape/math.h	Wed Feb 10 15:42:18 2021 -0800
+++ b/sys/include/ape/math.h	Sun Feb 14 10:02:40 2021 +0100
@@ -30,10 +30,13 @@
 extern double modf(double, double *);
 extern double pow(double, double);
 extern double sqrt(double);
+extern float ceilf(float);
 extern double ceil(double);
 extern double fabs(double);
+extern float fabsf(float);
 extern double floor(double);
 extern double fmod(double, double);
+extern long int lroundf(float);
 extern double NaN(void);
 extern int isNaN(double);
 extern double Inf(int);
diff -r ce98610ce572 sys/src/ape/lib/ap/math/ceilf.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sys/src/ape/lib/ap/math/ceilf.c	Sun Feb 14 10:02:40 2021 +0100
@@ -0,0 +1,46 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <math.h>
+#include <stdint.h>
+#include "math_private.h"
+
+static const float huge = 1.0e30;
+
+float
+ceilf(float x)
+{
+	int32_t i0,jj0;
+	uint32_t i;
+
+	GET_FLOAT_WORD(i0,x);
+	jj0 = ((i0>>23)&0xff)-0x7f;
+	if(jj0<23) {
+	    if(jj0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0<0) {i0=0x80000000;} 
+		    else if(i0!=0) { i0=0x3f800000;}
+		}
+	    } else {
+		i = (0x007fffff)>>jj0;
+		if((i0&i)==0) return x; /* x is integral */
+		if(huge+x>(float)0.0) {	/* raise inexact flag */
+		    if(i0>0) i0 += (0x00800000)>>jj0;
+		    i0 &= (~i);
+		}
+	    }
+	} else {
+	    if(jj0==0x80) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	}
+	SET_FLOAT_WORD(x,i0);
+	return x;
+}
diff -r ce98610ce572 sys/src/ape/lib/ap/math/fabsf.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sys/src/ape/lib/ap/math/fabsf.c	Sun Feb 14 10:02:40 2021 +0100
@@ -0,0 +1,27 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * fabsf(x) returns the absolute value of x.
+ */
+
+#include <math.h>
+#include <stdint.h>
+#include "math_private.h"
+
+float
+fabsf(float x)
+{
+	uint32_t ix;
+	GET_FLOAT_WORD(ix,x);
+	SET_FLOAT_WORD(x,ix&0x7fffffff);
+        return x;
+}
diff -r ce98610ce572 sys/src/ape/lib/ap/math/lroundf.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sys/src/ape/lib/ap/math/lroundf.c	Sun Feb 14 10:02:40 2021 +0100
@@ -0,0 +1,42 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+#include <math.h>
+#include <stdint.h>
+#include "math_private.h"
+
+long int
+lroundf(float x)
+{
+  int32_t exponent_less_127;
+  uint32_t w;
+  long int result;
+  int32_t sign;
+  GET_FLOAT_WORD (w, x);
+  exponent_less_127 = ((w & 0x7f800000) >> 23) - 127;
+  sign = (w & 0x80000000) != 0 ? -1 : 1;
+  w &= 0x7fffff;
+  w |= 0x800000;
+  if (exponent_less_127 < (int)((8 * sizeof (long int)) - 1))
+    {
+      if (exponent_less_127 < 0)
+        return exponent_less_127 < -1 ? 0 : sign;
+      else if (exponent_less_127 >= 23)
+        result = (long int) w << (exponent_less_127 - 23);
+      else
+        {
+          w += 0x400000 >> exponent_less_127;
+          result = w >> (23 - exponent_less_127);
+        }
+    }
+  else
+      return (long int) x;
+  return sign * result;
+}
diff -r ce98610ce572 sys/src/ape/lib/ap/math/math_private.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sys/src/ape/lib/ap/math/math_private.h	Sun Feb 14 10:02:40 2021 +0100
@@ -0,0 +1,42 @@
+#ifndef __MATH_PRIVATE
+#define __MATH_PRIVATE
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* A union which permits us to convert between a float and a 32 bit
+   int.  */
+
+typedef union
+{
+  float value;
+  uint32_t word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)					\
+do {								\
+  ieee_float_shape_type gf_u;					\
+  gf_u.value = (d);						\
+  (i) = gf_u.word;						\
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)					\
+do {								\
+  ieee_float_shape_type sf_u;					\
+  sf_u.word = (i);						\
+  (d) = sf_u.value;						\
+} while (0)
+
+
+#endif /* __MATH_PRIVATE */
diff -r ce98610ce572 sys/src/ape/lib/ap/math/mkfile
--- a/sys/src/ape/lib/ap/math/mkfile	Wed Feb 10 15:42:18 2021 -0800
+++ b/sys/src/ape/lib/ap/math/mkfile	Sun Feb 14 10:02:40 2021 +0100
@@ -5,9 +5,11 @@
 	asin.$O\
 	atan.$O\
 	atan2.$O\
+	ceilf.$O\
 	erf.$O\
 	exp.$O\
 	fabs.$O\
+	fabsf.$O\
 	floor.$O\
 	fmin.$O\
 	fmod.$O\
@@ -17,6 +19,7 @@
 	j1.$O\
 	jn.$O\
 	log.$O\
+	lroundf.$O\
 	pow.$O\
 	sin.$O\
 	sinh.$O\

Re: [9front] [PATCH] APE: add ceilf, fabsf and lroundf

[ori]

Quoth telephil9@gmail.com:
> i,
> 
> Here is a patch to add ceilf, fabsf and lroundf to APE math library.
> These functions are needed for the netsurf port.
> 
> Source code was copied from sourceware newlib. Looking at implementations from various OSes, I see that OpenBSD uses these exact same versions copied from the same sources.
> 
> 
> +
> +float
> +ceilf(float x)
> 
> +float
> +fabsf(float x)

These two have a lot more bit twiddling than I'd like;
we've got versios that work on the floating point values
for doubles; it should be straightforward to convert
them over to floats: mostly a matter of suffixing the
functions in those files with an 'f'.

> diff -r ce98610ce572 sys/src/ape/lib/ap/math/lroundf.c
> --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
> +++ b/sys/src/ape/lib/ap/math/lroundf.c	Sun Feb 14 10:02:40 2021 +0100

This one is.. still more twiddling than I'd like, but it looks
ok; while we're here, can we add lround(double)?