aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r--sysdeps/ieee754/cbrt.c120
-rw-r--r--sysdeps/ieee754/drem.c107
-rw-r--r--sysdeps/ieee754/logb.c48
-rw-r--r--sysdeps/ieee754/sqrt.c120
4 files changed, 0 insertions, 395 deletions
diff --git a/sysdeps/ieee754/cbrt.c b/sysdeps/ieee754/cbrt.c
deleted file mode 100644
index fe5fb95511..0000000000
--- a/sysdeps/ieee754/cbrt.c
+++ /dev/null
@@ -1,120 +0,0 @@
-/*
- * Copyright (c) 1985, 1993
- * The Regents of the University of California. All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- * notice, this list of conditions and the following disclaimer in the
- * documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- * must display the following acknowledgement:
- * This product includes software developed by the University of
- * California, Berkeley and its contributors.
- * 4. Neither the name of the University nor the names of its contributors
- * may be used to endorse or promote products derived from this software
- * without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- */
-
-#ifndef lint
-static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93";
-#endif /* not lint */
-
-#include <sys/cdefs.h>
-
-/* kahan's cube root (53 bits IEEE double precision)
- * for IEEE machines only
- * coded in C by K.C. Ng, 4/30/85
- *
- * Accuracy:
- * better than 0.667 ulps according to an error analysis. Maximum
- * error observed was 0.666 ulps in an 1,000,000 random arguments test.
- *
- * Warning: this code is semi machine dependent; the ordering of words in
- * a floating point number must be known in advance. I assume that the
- * long interger at the address of a floating point number will be the
- * leading 32 bits of that floating point number (i.e., sign, exponent,
- * and the 20 most significant bits).
- * On a National machine, it has different ordering; therefore, this code
- * must be compiled with flag -DNATIONAL.
- */
-#if !defined(vax)&&!defined(tahoe)
-
-static const unsigned long
- B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
- B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
-static const double
- C= 19./35.,
- D= -864./1225.,
- E= 99./70.,
- F= 45./28.,
- G= 5./14.;
-
-double cbrt(x)
-double x;
-{
- double r,s,t=0.0,w;
- unsigned long *px = (unsigned long *) &x,
- *pt = (unsigned long *) &t,
- mexp,sign;
-
-#ifdef national /* ordering of words in a floating points number */
- const int n0=1,n1=0;
-#else /* national */
- const int n0=0,n1=1;
-#endif /* national */
-
- mexp=px[n0]&0x7ff00000;
- if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */
- if(x==0.0) return(x); /* cbrt(0) is itself */
-
- sign=px[n0]&0x80000000; /* sign= sign(x) */
- px[n0] ^= sign; /* x=|x| */
-
-
- /* rough cbrt to 5 bits */
- if(mexp==0) /* subnormal number */
- {pt[n0]=0x43500000; /* set t= 2**54 */
- t*=x; pt[n0]=pt[n0]/3+B2;
- }
- else
- pt[n0]=px[n0]/3+B1;
-
-
- /* new cbrt to 23 bits, may be implemented in single precision */
- r=t*t/x;
- s=C+r*t;
- t*=G+F/(s+E+D/s);
-
- /* chopped to 20 bits and make it larger than cbrt(x) */
- pt[n1]=0; pt[n0]+=0x00000001;
-
-
- /* one step newton iteration to 53 bits with error less than 0.667 ulps */
- s=t*t; /* t*t is exact */
- r=x/s;
- w=t+t;
- r=(r-t)/(w+r); /* r-t is exact */
- t=t+t*r;
-
-
- /* retore the sign bit */
- pt[n0] |= sign;
- return(t);
-}
-#endif
diff --git a/sysdeps/ieee754/drem.c b/sysdeps/ieee754/drem.c
deleted file mode 100644
index cab3a04535..0000000000
--- a/sysdeps/ieee754/drem.c
+++ /dev/null
@@ -1,107 +0,0 @@
-/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
-This file is part of the GNU C Library.
-
-The GNU C Library is free software; you can redistribute it and/or
-modify it under the terms of the GNU Library General Public License as
-published by the Free Software Foundation; either version 2 of the
-License, or (at your option) any later version.
-
-The GNU C Library is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-Library General Public License for more details.
-
-You should have received a copy of the GNU Library General Public
-License along with the GNU C Library; see the file COPYING.LIB. If
-not, write to the Free Software Foundation, Inc., 675 Mass Ave,
-Cambridge, MA 02139, USA. */
-
-/*
- * Copyright (c) 1985 Regents of the University of California.
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms are permitted provided
- * that: (1) source distributions retain this entire copyright notice and
- * comment, and (2) distributions including binaries display the following
- * acknowledgement: ``This product includes software developed by the
- * University of California, Berkeley and its contributors'' in the
- * documentation or other materials provided with the distribution and in
- * all advertising materials mentioning features or use of this software.
- * Neither the name of the University nor the names of its contributors may
- * be used to endorse or promote products derived from this software without
- * specific prior written permission.
- * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- */
-
-#include <ansidecl.h>
-#include <math.h>
-#include <float.h>
-#include "ieee754.h"
-
-/* Return the remainder of X/Y. */
-double
-DEFUN(__drem, (x, y),
- double x AND double y)
-{
- union ieee754_double ux, uy;
-
- ux.d = x;
- uy.d = y;
-#define x ux.d
-#define y uy.d
-
- uy.ieee.negative = 0;
-
- if (!__finite (x) || y == 0.0)
- return NAN;
- else if (__isnan (y))
- return y;
- else if (__isinf (y))
- return x;
- else if (uy.ieee.exponent <= 1)
- {
- /* Subnormal (or almost subnormal) Y value. */
- double b = __scalb (1.0, 54);
- y *= b;
- x = __drem (x, y);
- x *= b;
- return __drem (x, y) / b;
- }
- else if (y >= 1.7e308 / 2)
- {
- y /= 2;
- x /= 2;
- return __drem (x, y) * 2;
- }
- else
- {
- union ieee754_double a;
- double b;
- unsigned int negative = ux.ieee.negative;
- a.d = y + y;
- b = y / 2;
- ux.ieee.negative = 0;
- while (x > a.d)
- {
- unsigned short int k = ux.ieee.exponent - a.ieee.exponent;
- union ieee754_double tmp;
- tmp.d = a.d;
- tmp.ieee.exponent += k;
- if (x < tmp.d)
- --tmp.ieee.exponent;
- x -= tmp.d;
- }
- if (x > b)
- {
- x -= y;
- if (x >= b)
- x -= y;
- }
- ux.ieee.negative ^= negative;
- return x;
- }
-}
-
-weak_alias (__drem, drem)
diff --git a/sysdeps/ieee754/logb.c b/sysdeps/ieee754/logb.c
deleted file mode 100644
index 918de9ebad..0000000000
--- a/sysdeps/ieee754/logb.c
+++ /dev/null
@@ -1,48 +0,0 @@
-/* Copyright (C) 1992, 1995 Free Software Foundation, Inc.
-This file is part of the GNU C Library.
-
-The GNU C Library is free software; you can redistribute it and/or
-modify it under the terms of the GNU Library General Public License as
-published by the Free Software Foundation; either version 2 of the
-License, or (at your option) any later version.
-
-The GNU C Library is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-Library General Public License for more details.
-
-You should have received a copy of the GNU Library General Public
-License along with the GNU C Library; see the file COPYING.LIB. If
-not, write to the Free Software Foundation, Inc., 675 Mass Ave,
-Cambridge, MA 02139, USA. */
-
-#include <ansidecl.h>
-#include <math.h>
-#include <float.h>
-#include "ieee754.h"
-
-/* Return the base 2 signed integral exponent of X. */
-double
-DEFUN(__logb, (x), double x)
-{
- union ieee754_double u;
-
- if (__isnan (x))
- return x;
- else if (__isinf (x))
- return HUGE_VAL;
- else if (x == 0.0)
- return - HUGE_VAL;
-
- u.d = x;
-
- if (u.ieee.exponent == 0)
- /* A denormalized number.
- Multiplying by 2 ** DBL_MANT_DIG normalizes it;
- we then subtract the DBL_MANT_DIG we added to the exponent. */
- return (__logb (x * ldexp (1.0, DBL_MANT_DIG)) - DBL_MANT_DIG);
-
- return (int) u.ieee.exponent - (DBL_MAX_EXP - 1);
-}
-
-weak_alias (__logb, logb)
diff --git a/sysdeps/ieee754/sqrt.c b/sysdeps/ieee754/sqrt.c
deleted file mode 100644
index 7e350e0d91..0000000000
--- a/sysdeps/ieee754/sqrt.c
+++ /dev/null
@@ -1,120 +0,0 @@
-/*
- * Copyright (c) 1985 Regents of the University of California.
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms are permitted provided
- * that: (1) source distributions retain this entire copyright notice and
- * comment, and (2) distributions including binaries display the following
- * acknowledgement: ``This product includes software developed by the
- * University of California, Berkeley and its contributors'' in the
- * documentation or other materials provided with the distribution and in
- * all advertising materials mentioning features or use of this software.
- * Neither the name of the University nor the names of its contributors may
- * be used to endorse or promote products derived from this software without
- * specific prior written permission.
- * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
- * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
- */
-
-#include <ansidecl.h>
-#include <errno.h>
-#include <math.h>
-
-/* Return the square root of X. */
-double
-DEFUN (sqrt, (x), double x)
-{
- double q, s, b, r, t;
- CONST double zero = 0.0;
- int m, n, i;
-
- /* sqrt (NaN) is NaN; sqrt (+-0) is +-0. */
- if (__isnan (x) || x == zero)
- return x;
-
- if (x < zero)
- return zero / zero;
-
- /* sqrt (Inf) is Inf. */
- if (__isinf (x))
- return x;
-
- /* Scale X to [1,4). */
- n = __logb (x);
- x = __scalb (x, -n);
- m = __logb (x);
- if (m != 0)
- /* Subnormal number. */
- x = __scalb (x, -m);
-
- m += n;
- n = m / 2;
-
- if ((n + n) != m)
- {
- x *= 2;
- --m;
- n = m / 2;
- }
-
- /* Generate sqrt (X) bit by bit (accumulating in Q). */
- q = 1.0;
- s = 4.0;
- x -= 1.0;
- r = 1;
- for (i = 1; i <= 51; i++)
- {
- t = s + 1;
- x *= 4;
- r /= 2;
- if (t <= x)
- {
- s = t + t + 2, x -= t;
- q += r;
- }
- else
- s *= 2;
- }
-
- /* Generate the last bit and determine the final rounding. */
- r /= 2;
- x *= 4;
- if (x == zero)
- goto end;
- (void) (100 + r); /* Trigger inexact flag. */
- if (s < x)
- {
- q += r;
- x -= s;
- s += 2;
- s *= 2;
- x *= 4;
- t = (x - s) - 5;
- b = 1.0 + 3 * r / 4;
- if (b == 1.0)
- goto end; /* B == 1: Round to zero. */
- b = 1.0 + r / 4;
- if (b > 1.0)
- t = 1; /* B > 1: Round to +Inf. */
- if (t >= 0)
- q += r;
- } /* Else round to nearest. */
- else
- {
- s *= 2;
- x *= 4;
- t = (x - s) - 1;
- b = 1.0 + 3 * r / 4;
- if (b == 1.0)
- goto end;
- b = 1.0 + r / 4;
- if (b > 1.0)
- t = 1;
- if (t >= 0)
- q += r;
- }
-
-end:
- return __scalb (q, n);
-}