jepler/RTAI
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Compare commits

merge into: jepler:master
Branches Tags
jepler:master
jepler:remove-kernel-standard
jepler:libm-no-userspace-headers
jepler:remove-nonfree-code-harder
jepler:remove-nonfree-code
...
pull from: jepler:remove-nonfree-code
Branches Tags
jepler:remove-kernel-standard
jepler:libm-no-userspace-headers
jepler:remove-nonfree-code-harder
jepler:remove-nonfree-code
jepler:master

1 commit
master ... remove-non

Author SHA1 Message Date
Jeff Epler
04c0b06102 Remove files that do not have a distributable license 
These files do not carry a licences that makes them distributable.
However, all of them appear to be either unused or not required by
linuxcnc, based on my analysis at
    http://mid.gmane.org/20131101024943.GA31516%40unpythonic.net
 2014-08-20 07:04:33 -05:00
9 changed files with 5 additions and 1489 deletions
Show stats Expand all files Collapse all files

10
base/math/GNUmakefile.am
View file

@ -7,15 +7,15 @@ CROSS_COMPILE = @CROSS_COMPILE@
if CONFIG_RTAI_MATH_C99
libmath_a_SOURCES = \
ceilfloor.c e_acos.c e_acosh.c e_asin.c e_atan2.c e_atanh.c e_cosh.c \
e_acos.c e_acosh.c e_asin.c e_atan2.c e_atanh.c e_cosh.c \
e_exp.c e_fmod.c e_gamma.c e_gamma_r.c e_hypot.c e_j0.c e_j1.c e_jn.c \
e_lgamma.c e_lgamma_r.c e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
e_scalb.c e_sinh.c e_sqrt.c fpmacros.c frexpldexp.c k_cos.c k_rem_pio2.c \
k_sin.c k_standard.c k_tan.c logb.c rndint.c s_asinh.c s_atan.c s_cbrt.c \
e_scalb.c e_sinh.c e_sqrt.c k_cos.c k_rem_pio2.c \
k_sin.c k_standard.c k_tan.c s_asinh.c s_atan.c s_cbrt.c \
s_ceil.c s_copysign.c s_cos.c s_erf.c s_expm1.c s_fabs.c s_finite.c s_floor.c \
s_frexp.c s_ilogb.c s_ldexp.c s_lib_version.c s_log1p.c s_logb.c s_matherr.c \
s_modf.c s_nextafter.c s_rint.c s_scalbn.c s_signgam.c s_significand.c s_sin.c \
s_tan.c s_tanh.c scalb.c sign.c w_acos.c w_acosh.c w_asin.c w_atan2.c w_atanh.c \
s_tan.c s_tanh.c w_acos.c w_acosh.c w_asin.c w_atan2.c w_atanh.c \
w_cabs.c w_cosh.c w_drem.c w_exp.c w_fmod.c w_gamma.c w_gamma_r.c w_hypot.c \
w_j0.c w_j1.c w_jn.c w_lgamma.c w_lgamma_r.c w_log.c w_log10.c w_pow.c \
w_remainder.c w_scalb.c w_sinh.c w_sqrt.c libm.c
@ -31,7 +31,7 @@ w_cosh.c w_exp.c w_fmod.c w_log.c w_log10.c w_pow.c w_sinh.c w_sqrt.c libm.c
endif
libmath_a_SOURCES += fpP.h mathP.h
libmath_a_SOURCES += mathP.h
if CONFIG_KBUILD
rtai_math.ko: @RTAI_KBUILD_ENV@

179
base/math/ceilfloor.c
View file

@ -1,179 +0,0 @@
#if defined(__ppc__)
/*******************************************************************************
* *
* File ceilfloor.c, *
* Function ceil(x) and floor(x), *
* Implementation of ceil and floor for the PowerPC. *
* *
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on November 1991, *
* *
* based on math.h, library code for Macintoshes with a 68881/68882 *
* by Jim Thomas. *
* *
* W A R N I N G: This routine expects a 64 bit double model. *
* *
* December 03 1992: first rs6000 port. *
* July 14 1993: comment changes and addition of #pragma fenv_access. *
* May 06 1997: port of the ibm/taligent ceil and floor routines. *
* April 11 2001: first port to os x using gcc. *
* June 13 2001: replaced __setflm with in-line assembly *
* *
*******************************************************************************/
#if !defined(__ppc__)
#define asm(x)
#endif
static const double twoTo52 = 4503599627370496.0;
static const unsigned long signMask = 0x80000000ul;
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;
/*******************************************************************************
* Functions needed for the computation. *
*******************************************************************************/
/*******************************************************************************
* Ceil(x) returns the smallest integer not less than x. *
*******************************************************************************/
double ceil ( double x )
{
DblInHex xInHex,OldEnvironment;
register double y;
register unsigned long int xhi;
register int target;
xInHex.dbl = x;
xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x|
target = ( xInHex.words.hi < signMask );
if ( xhi < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xhi < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case
return ( x );
else
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 1.0 );
else
return ( -0.0 );
}
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{
y = ( x + twoTo52 ) - twoTo52; // round at binary pt.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary pt.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
/*******************************************************************************
* Floor(x) returns the largest integer not greater than x. *
*******************************************************************************/
double floor ( double x )
{
DblInHex xInHex,OldEnvironment;
register double y;
register unsigned long int xhi;
register long int target;
xInHex.dbl = x;
xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x|
target = ( xInHex.words.hi < signMask );
if ( xhi < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xhi < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case
return ( x );
else
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 0.0 );
else
return ( -1.0 );
}
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{
y = ( x + twoTo52 ) - twoTo52; // round at binary pt.
if ( y > x )
return ( y - 1.0 );
else
return ( y );
}
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary pt.
if ( y > x )
return ( y - 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
#endif /* __ppc__ */

112
base/math/fpP.h
View file

@ -1,112 +0,0 @@
/*******************************************************************************
* *
* File fpP.h, *
* All pack 4 dependencies for the MathLib elems plus some defines used *
* throughout MathLib. *
* *
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on October 1991, *
* *
* W A R N I N G: This routine expects a 64 bit double model. *
* *
*******************************************************************************/
#define NoException 0
/*******************************************************************************
* Values of constants. *
*******************************************************************************/
//#define SgnMask 0x8000
#define dSgnMask 0x80000000
#define sSgnMask 0x7FFFFFFF
//#define ExpMask 0x7FFF
#define dExpMask 0x7FF00000
#define sExpMask 0xFF000000
/* according to rounding BIG & SMALL are: */
#define BIG 1.1e+300 /* used to deliver ±° or largest number, */
#define SMALL 1.1e-300 /* used to deliver ±0 or smallest number. */
#define InfExp 0x7FF
#define dMaxExp 0x7FF00000
#define MaxExpP1 1024
#define MaxExp 1023
#define DenormLimit -52
//#define ManMask 0x80000000
#define dManMask 0x00080000
//#define IsItDenorm 0x80000000
#define dIsItDenorm 0x00080000
//#define xIsItSNaN 0x40000000
#define dIsItSNaN 0x00080000
#define dHighMan 0x000FFFFF
#define dFirstBitSet 0x00080000
#define BIAS 0x3FF
//#define GetSign 0x8000
#define dGetSign 0x80000000
#define sGetSign 0x80000000
//#define Infinity(x) ( x.hex.exponent & ExpMask ) == ExpMask
#define dInfinity(x) ( x.hex.high & dExpMask ) == dExpMask
#define sInfinity(x) ( ( x.hexsgl << 1 ) & sExpMask ) == sExpMask
//#define Exponent(x) x.hex.exponent & ExpMask
#define dExponent(x) x.hex.high & dExpMask
#define sExponent(x) ( ( x.hexsgl << 1 ) & sExpMask )
#define sZero(x) ( x.hexsgl & sSgnMask ) == 0
//#define Sign(x) ( x.hex.exponent & SgnMask ) == SgnMask
/*******************************************************************************
* Types used in the auxiliary functions. *
*******************************************************************************/
typedef struct /* Hex representation of a double. */
{
#if defined(__BIG_ENDIAN__)
unsigned long int high;
unsigned long int low;
#else
unsigned long int low;
unsigned long int high;
#endif
} dHexParts;
typedef union
{
unsigned char byties[8];
double dbl;
} DblInHex;
//enum boolean { FALSE, TRUE };
/*******************************************************************************
* Macros to access long subfields of a double value. *
*******************************************************************************/
#define highpartd(x) *((long *) &x)
#define lowpartd(x) *((long *) &x + 1)
enum {
FP_SNAN = 0, /* signaling NaN
*/
FP_QNAN = 1, /* quiet NaN
*/
FP_INFINITE = 2, /* + or - infinity
*/
FP_ZERO = 3, /* + or - zero
*/
FP_NORMAL = 4, /* all normal numbers
*/
FP_SUBNORMAL = 5 /* denormal numbers
*/
};

239
base/math/fpmacros.c
View file

@ -1,239 +0,0 @@
/***********************************************************************
** File: fpmacros.c
**
** Contains: C source code for implementations of floating-point
** functions which involve float format numbers, as
** defined in header <fp.h>. In particular, this file
** contains implementations of functions
** __fpclassify(d,f), __isnormal(d,f), __isfinite(d,f),
** __isnan(d,f), and __signbit(d,f). This file targets
** PowerPC platforms.
**
** Written by: Robert A. Murley, Ali Sazegari
**
** Copyright: c 2001 by Apple Computer, Inc., all rights reserved
**
** Change History (most recent first):
**
** 07 Jul 01 ram First created from fpfloatfunc.c, fp.c,
** classify.c and sign.c in MathLib v3 Mac OS9.
**
***********************************************************************/
#include "fpP.h"
#define SIGN_MASK 0x80000000
#define NSIGN_MASK 0x7fffffff
#define FEXP_MASK 0x7f800000
#define FFRAC_MASK 0x007fffff
/***********************************************************************
long int __fpclassifyf(float x) returns the classification code of the
argument x, as defined in <fp.h>.
Exceptions: INVALID signaled if x is a signaling NaN; in this case,
the FP_QNAN code is returned.
Calls: none
***********************************************************************/
long int __fpclassifyf ( float x )
{
unsigned long int iexp;
union {
unsigned long int lval;
float fval;
} z;
z.fval = x;
iexp = z.lval & FEXP_MASK; /* isolate float exponent */
if (iexp == FEXP_MASK) { /* NaN or INF case */
if ((z.lval & 0x007fffff) == 0)
return (long int) FP_INFINITE;
else if ((z.lval & 0x00400000) != 0)
return (long int) FP_QNAN;
else
return (long int) FP_SNAN;
}
if (iexp != 0) /* normal float */
return (long int) FP_NORMAL;
if (x == 0.0)
return (long int) FP_ZERO; /* zero */
else
return (long int) FP_SUBNORMAL; /* must be subnormal */
}
/***********************************************************************
Function __fpclassify,
Implementation of classify of a double number for the PowerPC.
Exceptions: INVALID signaled if x is a signaling NaN; in this case,
the FP_QNAN code is returned.
Calls: none
***********************************************************************/
long int __fpclassify ( double arg )
{
register unsigned long int exponent;
union
{
dHexParts hex;
double dbl;
} x;
x.dbl = arg;
exponent = x.hex.high & dExpMask;
if ( exponent == dExpMask )
{
if ( ( ( x.hex.high & dHighMan ) | x.hex.low ) == 0 )
return (long int) FP_INFINITE;
else
return ( x.hex.high & 0x00080000 ) ? FP_QNAN : FP_SNAN;
}
else if ( exponent != 0)
return (long int) FP_NORMAL;
else {
if ( arg == 0.0 )
return (long int) FP_ZERO;
else
return (long int) FP_SUBNORMAL;
}
}
/***********************************************************************
long int __isnormalf(float x) returns nonzero if and only if x is a
normalized float number and zero otherwise.
Exceptions: INVALID is raised if x is a signaling NaN; in this case,
zero is returned.
Calls: none
***********************************************************************/
long int __isnormalf ( float x )
{
unsigned long int iexp;
union {
unsigned long int lval;
float fval;
} z;
z.fval = x;
iexp = z.lval & FEXP_MASK; /* isolate float exponent */
return ((iexp != FEXP_MASK) && (iexp != 0));
}
long int __isnorma ( double x )
{
return ( __fpclassify ( x ) == FP_NORMAL );
}
/***********************************************************************
long int __isfinitef(float x) returns nonzero if and only if x is a
finite (normal, subnormal, or zero) float number and zero otherwise.
Exceptions: INVALID is raised if x is a signaling NaN; in this case,
zero is returned.
Calls: none
***********************************************************************/
long int __isfinitef ( float x )
{
union {
unsigned long int lval;
float fval;
} z;
z.fval = x;
return ((z.lval & FEXP_MASK) != FEXP_MASK);
}
long int __isfinite ( double x )
{
return ( __fpclassify ( x ) >= FP_ZERO );
}
/***********************************************************************
long int __isnanf(float x) returns nonzero if and only if x is a
NaN and zero otherwise.
Exceptions: INVALID is raised if x is a signaling NaN; in this case,
nonzero is returned.
Calls: none
***********************************************************************/
long int __isnanf ( float x )
{
union {
unsigned long int lval;
float fval;
} z;
z.fval = x;
return (((z.lval&FEXP_MASK) == FEXP_MASK) && ((z.lval&FFRAC_MASK) != 0));
}
long int __isnan ( double x )
{
long int class = __fpclassify(x);
return ( ( class == FP_SNAN ) || ( class == FP_QNAN ) );
}
/***********************************************************************
long int __signbitf(float x) returns nonzero if and only if the sign
bit of x is set and zero otherwise.
Exceptions: INVALID is raised if x is a signaling NaN.
Calls: none
***********************************************************************/
long int __signbitf ( float x )
{
union {
unsigned long int lval;
float fval;
} z;
z.fval = x;
return ((z.lval & SIGN_MASK) != 0);
}
/***********************************************************************
Function sign of a double.
Implementation of sign bit for the PowerPC.
Calls: none
***********************************************************************/
long int __signbit ( double arg )
{
union
{
dHexParts hex;
double dbl;
} x;
long int sign;
x.dbl = arg;
sign = ( ( x.hex.high & dSgnMask ) == dSgnMask ) ? 1 : 0;
return sign;
}

73
base/math/frexpldexp.c
View file

@ -1,73 +0,0 @@
#if defined(__ppc__)
/*******************************************************************************
* *
* File frexpldexp.c, *
* Functions frexp(x) and ldexp(x), *
* Implementation of frexp and ldexp functions for the PowerPC. *
* *
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on January 1991, *
* *
* W A R N I N G: This routine expects a 64 bit double model. *
* *
* December03 1992: first rs6000 implementation. *
* October 05 1993: added special cases for NaN and ° in frexp. *
* May 27 1997: improved the performance of frexp by eliminating the *
* switch statement. *
* June 13 2001: (ram) rewrote frexp to eliminate calls to scalb and *
* logb. *
* *
*******************************************************************************/
#include <limits.h>
#include <math.h>
static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
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;
double ldexp ( double value, int exp )
{
if ( exp > SHRT_MAX )
exp = SHRT_MAX;
else if ( exp < -SHRT_MAX )
exp = -SHRT_MAX;
return scalb ( value, exp );
}
double frexp ( double value, int *eptr )
{
DblInHex argument;
unsigned long int valueHead;
argument.dbl = value;
valueHead = argument.words.hi & 0x7fffffffUL; // valueHead <- |x|
*eptr = 0;
if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 )
return value; // 0, inf, or NaN
if ( valueHead < 0x00100000 )
{ // denorm
argument.dbl = two54 * value;
valueHead = argument.words.hi &0x7fffffff;
*eptr = -54;
}
*eptr += ( valueHead >> 20 ) - 1022;
argument.words.hi = ( argument.words.hi & 0x800fffff ) | 0x3fe00000;
return argument.dbl;
}
#endif /* __ppc__ */

104
base/math/logb.c
View file

@ -1,104 +0,0 @@
#if defined(__ppc__)
/*******************************************************************************
* *
* File logb.c, *
* Functions logb. *
* Implementation of logb for the PowerPC. *
* *
* Copyright © 1991 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on June 1991, *
* *
* August 26 1991: removed CFront Version 1.1d17 warnings. *
* August 27 1991: no errors reported by the test suite. *
* November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and *
* + or - infinity to constants. *
* November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to *
* improve performance. *
* February 07 1992: changed bit operations to macros ( object size is *
* unchanged ). *
* September24 1992: took the "#include support.h" out. *
* December 03 1992: first rs/6000 port. *
* August 30 1992: set the divide by zero for the zero argument case. *
* October 05 1993: corrected the environment. *
* October 17 1994: replaced all environmental functions with __setflm. *
* May 28 1997: made speed improvements. *
* April 30 2001: forst mac os x port using gcc. *
* *
********************************************************************************
* The C math library offers a similar function called "frexp". It is *
* different in details from logb, but similar in spirit. This current *
* implementation of logb follows the recommendation in IEEE Standard 854 *
* which is different in its handling of denormalized numbers from the IEEE *
* Standard 754. *
*******************************************************************************/
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 twoTo52 = 4.50359962737049600e15; // 0x1p52
static const double klTod = 4503601774854144.0; // 0x1.000008p52
static const unsigned long int signMask = 0x80000000ul;
static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }};
/*******************************************************************************
********************************************************************************
* L O G B *
********************************************************************************
*******************************************************************************/
double logb ( double x )
{
DblInHex xInHex;
long int shiftedExp;
xInHex.dbl = x;
shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
if ( shiftedExp == 2047 )
{ // NaN or INF
if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
return x; // NaN or +INF return x
else
return -x; // -INF returns +INF
}
if ( shiftedExp != 0 ) // normal number
shiftedExp -= 1023; // unbias exponent
else if ( x == 0.0 )
{ // zero
xInHex.words.hi = 0x0UL; // return -infinity
return ( minusInf.dbl );
}
else
{ // subnormal number
xInHex.dbl *= twoTo52; // scale up
shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
shiftedExp -= 1075; // unbias exponent
}
if ( shiftedExp == 0 ) // zero result
return ( 0.0 );
else
{ // nonzero result
xInHex.dbl = klTod;
xInHex.words.lo += shiftedExp;
return ( xInHex.dbl - klTod );
}
}
#endif /* __ppc__ */

632
base/math/rndint.c
View file

@ -1,632 +0,0 @@
/*******************************************************************************
** File: rndint.c
**
** Contains: C source code for implementations of floating-point
** functions which round to integral value or format, as
** defined in header <fp.h>. In particular, this file
** contains implementations of functions rint, nearbyint,
** rinttol, round, roundtol, trunc, modf and modfl. This file
** targets PowerPC or Power platforms.
**
** Written by: A. Sazegari, Apple AltiVec Group
** Created originally by Jon Okada, Apple Numerics Group
**
** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
**
** Change History (most recent first):
**
** 13 Jul 01 ram replaced --setflm calls with inline assembly
** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm
** definition.
** 1. removed double_t, put in double for now.
** 2. removed iclass from nearbyint.
** 3. removed wrong comments intrunc.
** 4.
** 13 May 97 ali made performance improvements in rint, rinttol, roundtol
** and trunc by folding some of the taligent ideas into this
** implementation. nearbyint is faster than the one in taligent,
** rint is more elegant, but slower by %30 than the taligent one.
** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c
** 15 Sep 94 ali Major overhaul and performance improvements of all functions.
** 20 Jul 94 PAF New faster version
** 16 Jul 93 ali Added the modfl function.
** 18 Feb 93 ali Changed the return value of fenv functions
** feclearexcept and feraiseexcept to their new
** NCEG X3J11.1/93-001 definitions.
** 16 Dec 92 JPO Removed __itrunc implementation to a
** separate file.
** 15 Dec 92 JPO Added __itrunc implementation and modified
** rinttol to include conversion from double
** to long int format. Modified roundtol to
** call __itrunc.
** 10 Dec 92 JPO Added modf (double) implementation.
** 04 Dec 92 JPO First created.
**
*******************************************************************************/
#include <limits.h>
#include <math.h>
#if !defined(__ppc__)
#define asm(x)
#endif
#define SET_INVALID 0x01000000UL
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 unsigned long int signMask = 0x80000000ul;
static const double twoTo52 = 4503599627370496.0;
static const double doubleToLong = 4503603922337792.0; // 2^52
static const DblInHex Huge = {{ 0x7FF00000, 0x00000000 }};
static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
/*******************************************************************************
* *
* The function rint rounds its double argument to integral value *
* according to the current rounding direction and returns the result in *
* double format. This function signals inexact if an ordered return *
* value is not equal to the operand. *
* *
********************************************************************************
* *
* This function calls: fabs. *
* *
*******************************************************************************/
/*******************************************************************************
* First, an elegant implementation. *
********************************************************************************
*
*double rint ( double x )
* {
* double y;
*
* y = twoTo52.fval;
*
* if ( fabs ( x ) >= y ) // huge case is exact
* return x;
* if ( x < 0 ) y = -y; // negative case
* y = ( x + y ) - y; // force rounding
* if ( y == 0.0 ) // zero results mirror sign of x
* y = copysign ( y, x );
* return ( y );
* }
********************************************************************************
* Now a bit twidling version that is about %30 faster. *
*******************************************************************************/
#if defined(__ppc__)
double rint ( double x )
{
DblInHex argument;
register double y;
unsigned long int xHead;
register long int target;
argument.dbl = x;
xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
target = ( argument.words.hi < signMask ); // flags positive sign
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( target )
y = ( x + twoTo52 ) - twoTo52; // round at binary point
else
y = ( x - twoTo52 ) + twoTo52; // round at binary point
if ( y == 0.0 )
{ // fix sign of zero result
if ( target )
return ( 0.0 );
else
return ( -0.0 );
}
return y;
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt.
else
return ( ( x - twoTo52 ) + twoTo52 );
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
#endif /* __ppc__ */
/*******************************************************************************
* *
* The function nearbyint rounds its double argument to integral value *
* according to the current rounding direction and returns the result in *
* double format. This function does not signal inexact. *
* *
********************************************************************************
* *
* This function calls fabs and copysign. *
* *
*******************************************************************************/
double nearbyint ( double x )
{
double y;
#if defined(__ppc__)
double OldEnvironment;
#endif /* __ppc__ */
y = twoTo52;
asm ("mffs %0" : "=f" (OldEnvironment)); /* get the environement */
if ( fabs ( x ) >= y ) /* huge case is exact */
return x;
if ( x < 0 ) y = -y; /* negative case */
y = ( x + y ) - y; /* force rounding */
if ( y == 0.0 ) /* zero results mirror sign of x */
y = copysign ( y, x );
// restore old flags
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
return ( y );
}
/*******************************************************************************
* *
* The function rinttol converts its double argument to integral value *
* according to the current rounding direction and returns the result in *
* long int format. This conversion signals invalid if the argument is a *
* NaN or the rounded intermediate result is out of range of the *
* destination long int format, and it delivers an unspecified result in *
* this case. This function signals inexact if the rounded result is *
* within range of the long int format but unequal to the operand. *
* *
*******************************************************************************/
long int rinttol ( double x )
{
register double y;
DblInHex argument, OldEnvironment;
unsigned long int xHead;
register long int target;
argument.dbl = x;
target = ( argument.words.hi < signMask ); // flag positive sign
xHead = argument.words.hi & 0x7ffffffful; // high 32 bits of x
if ( target )
/*******************************************************************************
* Sign of x is positive. *
*******************************************************************************/
{
if ( xHead < 0x41dffffful )
{ // x is safely in long range
y = ( x + twoTo52 ) - twoTo52; // round at binary point
argument.dbl = y + doubleToLong; // force result into argument.words.lo
return ( ( long ) argument.words.lo );
}
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead > 0x41dffffful )
{ // x is safely out of long range
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX );
}
/*******************************************************************************
* x > 0.0 and may or may not be out of range of long. *
*******************************************************************************/
y = ( x + twoTo52 ) - twoTo52; // do rounding
if ( y > ( double ) LONG_MAX )
{ // out of range of long
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX );
}
argument.dbl = y + doubleToLong; // in range
return ( ( long ) argument.words.lo ); // return result & flags
}
/*******************************************************************************
* Sign of x is negative. *
*******************************************************************************/
if ( xHead < 0x41e00000ul )
{ // x is safely in long range
y = ( x - twoTo52 ) + twoTo52;
argument.dbl = y + doubleToLong;
return ( ( long ) argument.words.lo );
}
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead > 0x41e00000ul )
{ // x is safely out of long range
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN );
}
/*******************************************************************************
* x < 0.0 and may or may not be out of range of long. *
*******************************************************************************/
y = ( x - twoTo52 ) + twoTo52; // do rounding
if ( y < ( double ) LONG_MIN )
{ // out of range of long
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN );
}
argument.dbl = y + doubleToLong; // in range
return ( ( long ) argument.words.lo ); // return result & flags
}
/*******************************************************************************
* *
* The function round rounds its double argument to integral value *
* according to the "add half to the magnitude and truncate" rounding of *
* Pascal's Round function and FORTRAN's ANINT function and returns the *
* result in double format. This function signals inexact if an ordered *
* return value is not equal to the operand. *
* *
*******************************************************************************/
double round ( double x )
{
DblInHex argument, OldEnvironment;
register double y, z;
register unsigned long int xHead;
register long int target;
argument.dbl = x;
xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^52? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
if ( xHead < 0x3fe00000ul )
/*******************************************************************************
* Is |x| < 0.5? *
*******************************************************************************/
{
if ( ( xHead | argument.words.lo ) != 0ul )
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 0.0 );
else
return ( -0.0 );
}
/*******************************************************************************
* Is 0.5 ² |x| < 1.0? *
*******************************************************************************/
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target )
return ( 1.0 );
else
return ( -1.0 );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{ // positive x
y = ( x + twoTo52 ) - twoTo52; // round at binary point
if ( y == x ) // exact case
return ( x );
z = x + 0.5; // inexact case
y = ( z + twoTo52 ) - twoTo52; // round at binary point
if ( y > z )
return ( y - 1.0 );
else
return ( y );
}
/*******************************************************************************
* Is x < 0? *
*******************************************************************************/
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary point
if ( y == x )
return ( x );
z = x - 0.5;
y = ( z - twoTo52 ) + twoTo52; // round at binary point
if ( y < z )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
/*******************************************************************************
* *
* The function roundtol converts its double argument to integral format *
* according to the "add half to the magnitude and chop" rounding mode of *
* Pascal's Round function and FORTRAN's NINT function. This conversion *
* signals invalid if the argument is a NaN or the rounded intermediate *
* result is out of range of the destination long int format, and it *
* delivers an unspecified result in this case. This function signals *
* inexact if the rounded result is within range of the long int format but *
* unequal to the operand. *
* *
*******************************************************************************/
long int roundtol ( double x )
{
register double y, z;
DblInHex argument, OldEnvironment;
register unsigned long int xhi;
register long int target;
#if defined(__ppc__)
const DblInHex kTZ = {{ 0x0, 0x1 }};
const DblInHex kUP = {{ 0x0, 0x2 }};
#endif /* __ppc__ */
argument.dbl = x;
xhi = argument.words.hi & 0x7ffffffful; // high 32 bits of x
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xhi > 0x41e00000ul )
/*******************************************************************************
* Is x is out of long range or NaN? *
*******************************************************************************/
{
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
if ( target ) // pin result
return ( LONG_MAX );
else
return ( LONG_MIN );
}
if ( target )
/*******************************************************************************
* Is sign of x is "+"? *
*******************************************************************************/
{
if ( x < 2147483647.5 )
/*******************************************************************************
* x is in the range of a long. *
*******************************************************************************/
{
y = ( x + doubleToLong ) - doubleToLong; // round at binary point
if ( y != x )
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
z = x + 0.5; // truncate x + 0.5
argument.dbl = z + doubleToLong;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( ( long ) argument.words.lo );
}
argument.dbl = y + doubleToLong; // force result into argument.words.lo
return ( ( long ) argument.words.lo ); // return long result
}
/*******************************************************************************
* Rounded positive x is out of the range of a long. *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MAX ); // return pinned result
}
/*******************************************************************************
* x < 0.0 and may or may not be out of the range of a long. *
*******************************************************************************/
if ( x > -2147483648.5 )
/*******************************************************************************
* x is in the range of a long. *
*******************************************************************************/
{
y = ( x + doubleToLong ) - doubleToLong; // round at binary point
if ( y != x )
{ // inexact case
asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // save environment
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
z = x - 0.5; // truncate x - 0.5
argument.dbl = z + doubleToLong;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( ( long ) argument.words.lo );
}
argument.dbl = y + doubleToLong;
return ( ( long ) argument.words.lo ); // return long result
}
/*******************************************************************************
* Rounded negative x is out of the range of a long. *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= SET_INVALID;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
return ( LONG_MIN ); // return pinned result
}
/*******************************************************************************
* *
* The function trunc truncates its double argument to integral value *
* and returns the result in double format. This function signals *
* inexact if an ordered return value is not equal to the operand. *
* *
*******************************************************************************/
double trunc ( double x )
{
DblInHex argument,OldEnvironment;
register double y;
register unsigned long int xhi;
register long int target;
argument.dbl = x;
xhi = argument.words.hi & 0x7fffffffUL; // xhi <- high half of |x|
target = ( argument.words.hi < signMask ); // flag positive sign
if ( xhi < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^53? *
*******************************************************************************/
{
if ( xhi < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
if ( ( xhi | argument.words.lo ) != 0ul )
{ // raise deserved INEXACT
asm ("mffs %0" : "=f" (OldEnvironment.dbl));
OldEnvironment.words.lo |= 0x02000000ul;
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
}
if ( target ) // return properly signed zero
return ( 0.0 );
else
return ( -0.0 );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
if ( target )
{
y = ( x + twoTo52 ) - twoTo52; // round at binary point
if ( y > x )
return ( y - 1.0 );
else
return ( y );
}
else
{
y = ( x - twoTo52 ) + twoTo52; // round at binary point.
if ( y < x )
return ( y + 1.0 );
else
return ( y );
}
}
/*******************************************************************************
* Is |x| >= 2.0^52 or x is a NaN. *
*******************************************************************************/
return ( x );
}
/*******************************************************************************
* The modf family of functions separate a floating-point number into its *
* fractional and integral parts, returning the fractional part and writing *
* the integral part in floating-point format to the object pointed to by a *
* pointer argument. If the input argument is integral or infinite in *
* value, the return value is a zero with the sign of the input argument. *
* The modf family of functions raises no floating-point exceptions. older *
* implemenation set the INVALID flag due to signaling NaN input. *
* *
*******************************************************************************/
/*******************************************************************************
* modf is the double implementation. *
*******************************************************************************/
#if defined(__ppc__)
double modf ( double x, double *iptr )
{
register double OldEnvironment, xtrunc;
register unsigned long int xHead, signBit;
DblInHex argument;
argument.dbl = x;
xHead = argument.words.hi & 0x7ffffffful; // |x| high bit pattern
signBit = ( argument.words.hi & 0x80000000ul ); // isolate sign bit
if (xHead == 0x7ff81fe0)
signBit = signBit | 0;
if ( xHead < 0x43300000ul )
/*******************************************************************************
* Is |x| < 2.0^53? *
*******************************************************************************/
{
if ( xHead < 0x3ff00000ul )
/*******************************************************************************
* Is |x| < 1.0? *
*******************************************************************************/
{
argument.words.hi = signBit; // truncate to zero
argument.words.lo = 0ul;
*iptr = argument.dbl;
return ( x );
}
/*******************************************************************************
* Is 1.0 < |x| < 2.0^52? *
*******************************************************************************/
asm ("mffs %0" : "=f" (OldEnvironment)); // save environment
// round toward zero
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
if ( signBit == 0ul ) // truncate to integer
xtrunc = ( x + twoTo52 ) - twoTo52;
else
xtrunc = ( x - twoTo52 ) + twoTo52;
// restore caller's env
asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
*iptr = xtrunc; // store integral part
if ( x != xtrunc ) // nonzero fraction
return ( x - xtrunc );
else
{ // zero with x's sign
argument.words.hi = signBit;
argument.words.lo = 0ul;
return ( argument.dbl );
}
}
*iptr = x; // x is integral or NaN
if ( x != x ) // NaN is returned
return x;
else
{ // zero with x's sign
argument.words.hi = signBit;
argument.words.lo = 0ul;
return ( argument.dbl );
}
}
#endif /* __ppc__ */

87
base/math/scalb.c
View file

@ -1,87 +0,0 @@
#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__ */

58
base/math/sign.c
View file

@ -1,58 +0,0 @@
#if defined(__ppc__)
/*******************************************************************************
* *
* File sign.c, *
* Functions copysign and __signbitd. *
* For PowerPC based machines. *
* *
* Copyright © 1991, 2001 Apple Computer, Inc. All rights reserved. *
* *
* Written by Ali Sazegari, started on June 1991. *
* *
* August 26 1991: no CFront Version 1.1d17 warnings. *
* September 06 1991: passes the test suite with invalid raised on *
* signaling nans. sane rom code behaves the same. *
* September 24 1992: took the Ò#include support.hÓ out. *
* Dcember 02 1992: PowerPC port. *
* July 20 1994: __fabs added *
* July 21 1994: deleted unnecessary functions: neg, COPYSIGNnew, *
* and SIGNNUMnew. *
* April 11 2001: first port to os x using gcc. *
* removed fabs and deffered to gcc for direct *
* instruction generation. *
* *
*******************************************************************************/
#include "fpP.h"
/*******************************************************************************
* *
* Function copysign. *
* Implementation of copysign for the PowerPC. *
* *
********************************************************************************
* Note: The order of the operands in this function is reversed from that *
* suggested in the IEEE standard 754. *
*******************************************************************************/
double copysign ( double arg2, double arg1 )
{
union
{
dHexParts hex;
double dbl;
} x, y;
/*******************************************************************************
* No need to flush NaNs out. *
*******************************************************************************/
x.dbl = arg1;
y.dbl = arg2;
y.hex.high = y.hex.high & 0x7FFFFFFF;
y.hex.high = ( y.hex.high | ( x.hex.high & dSgnMask ) );
return y.dbl;
}
#endif /* __ppc__ */