Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 87 insertions, 0 deletions

diff --git a/libm/scalb.c b/libm/scalb.c
new file mode 100644
index 000000000..03d2de9dc
--- /dev/null
+++ b/libm/scalb.c
@@ -0,0 +1,87 @@
+#if defined(__ppc__)
+/***********************************************************************
+**      File:    scalb.c
+**      
+**      Contains: C source code for implementations of floating-point
+**                scalb functions defined in header <fp.h>.  In
+**                particular, this file contains implementations of
+**                functions scalb and scalbl for double and long double
+**                formats on PowerPC platforms.
+**                        
+**      Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
+**      
+**      Copyright: © 1992 by Apple Computer, Inc., all rights reserved
+**      
+**      Change History ( most recent first ):
+**
+**      28 May 97  ali   made an speed improvement for large n,
+**                       removed scalbl.
+**      12 Dec 92  JPO   First created.
+**                        
+***********************************************************************/
+
+typedef union           
+      { 
+      struct {
+#if defined(__BIG_ENDIAN__)
+        unsigned long int hi;
+        unsigned long int lo;
+#else
+        unsigned long int lo;
+        unsigned long int hi;
+#endif
+      } words;
+      double dbl;
+      } DblInHex;
+
+static const double twoTo1023  = 8.988465674311579539e307;   // 0x1p1023
+static const double twoToM1022 = 2.225073858507201383e-308;  // 0x1p-1022
+
+
+/***********************************************************************
+      double  scalb( double  x, long int n ) returns its argument x scaled
+      by the factor 2^m.  NaNs, signed zeros, and infinities are propagated
+      by this function regardless of the value of n.
+      
+      Exceptions:  OVERFLOW/INEXACT or UNDERFLOW inexact may occur;
+                         INVALID for signaling NaN inputs ( quiet NaN returned ).
+      
+      Calls:  none.
+***********************************************************************/
+
+double scalb ( double x, int n  )
+      {
+      DblInHex xInHex;
+      
+      xInHex.words.lo = 0UL;                     // init. low half of xInHex
+      
+      if ( n > 1023 ) 
+            {                                   // large positive scaling
+            if ( n > 2097 )                     // huge scaling
+            	return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023;
+            while ( n > 1023 ) 
+                  {                             // scale reduction loop
+                  x *= twoTo1023;               // scale x by 2^1023
+                  n -= 1023;                    // reduce n by 1023
+                  }
+            }
+      
+      else if ( n < -1022 ) 
+            {                                   // large negative scaling
+            if ( n < -2098 )                    // huge negative scaling
+            	return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022;
+            while ( n < -1022 ) 
+                  {                             // scale reduction loop
+                  x *= twoToM1022;              // scale x by 2^( -1022 )
+                  n += 1022;                    // incr n by 1022
+                  }
+            }
+
+/*******************************************************************************
+*      -1022 <= n <= 1023; convert n to double scale factor.                   *
+*******************************************************************************/
+
+      xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20;
+      return ( x * xInHex.dbl );
+      }
+#endif /* __ppc__ */