generic-library/bionic
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

am bfd89274: am 4abaa576: Merge "Switch to gdtoa."

Browse Source 
* commit 'bfd8927447b76bca6393c1e7c735a1333aa60786':
  Switch to gdtoa.
This commit is contained in:
Elliott Hughes 2014-04-17 23:34:44 +00:00 committed by Android Git Automerger
commit fb7b42df1b
33 changed files with 7907 additions and 2855 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
libc/Android.mk
View File

@ -85,7 +85,6 @@ libc_common_src_files := \
stdlib/getenv.c \
stdlib/putenv.c \
stdlib/setenv.c \
stdlib/strtod.c \
unistd/syslog.c \
# Fortify implementations of libc functions.
@ -208,7 +207,6 @@ libc_bionic_src_files := \
bionic/strerror.cpp \
bionic/strerror_r.cpp \
bionic/strsignal.cpp \
bionic/strtof.cpp \
bionic/strtold.cpp \
bionic/stubs.cpp \
bionic/symlink.cpp \
@ -311,6 +309,33 @@ libc_upstream_netbsd_src_files := \
upstream-netbsd/lib/libc/thread-stub/__isthreaded.c \
upstream-netbsd/lib/libc/unistd/killpg.c \
libc_upstream_openbsd_gdtoa_src_files := \
upstream-openbsd/gdtoa_support.cpp \
upstream-openbsd/lib/libc/gdtoa/dmisc.c \
upstream-openbsd/lib/libc/gdtoa/dtoa.c \
upstream-openbsd/lib/libc/gdtoa/gdtoa.c \
upstream-openbsd/lib/libc/gdtoa/gethex.c \
upstream-openbsd/lib/libc/gdtoa/gmisc.c \
upstream-openbsd/lib/libc/gdtoa/hd_init.c \
upstream-openbsd/lib/libc/gdtoa/hdtoa.c \
upstream-openbsd/lib/libc/gdtoa/hexnan.c \
upstream-openbsd/lib/libc/gdtoa/ldtoa.c \
upstream-openbsd/lib/libc/gdtoa/misc.c \
upstream-openbsd/lib/libc/gdtoa/smisc.c \
upstream-openbsd/lib/libc/gdtoa/strtod.c \
upstream-openbsd/lib/libc/gdtoa/strtodg.c \
upstream-openbsd/lib/libc/gdtoa/strtof.c \
upstream-openbsd/lib/libc/gdtoa/strtord.c \
upstream-openbsd/lib/libc/gdtoa/sum.c \
upstream-openbsd/lib/libc/gdtoa/ulp.c \
libc_upstream_openbsd_gdtoa_src_files_32 := \
$(libc_upstream_openbsd_gdtoa_src_files) \
libc_upstream_openbsd_gdtoa_src_files_64 := \
$(libc_upstream_openbsd_gdtoa_src_files) \
upstream-openbsd/lib/libc/gdtoa/strtorQ.c \
libc_upstream_openbsd_src_files := \
upstream-openbsd/lib/libc/gen/alarm.c \
upstream-openbsd/lib/libc/gen/exec.c \
@ -614,6 +639,34 @@ $(eval $(call patch-up-arch-specific-flags,LOCAL_CFLAGS,libc_common_cflags))
include $(BUILD_STATIC_LIBRARY)
# ========================================================
# libc_gdtoa.a - upstream OpenBSD C library gdtoa code
# ========================================================
#
# These files are built with the openbsd-compat.h header file
# automatically included.
include $(CLEAR_VARS)
LOCAL_SRC_FILES_32 := $(libc_upstream_openbsd_gdtoa_src_files_32)
LOCAL_SRC_FILES_64 := $(libc_upstream_openbsd_gdtoa_src_files_64)
LOCAL_CFLAGS := \
$(libc_common_cflags) \
-I$(LOCAL_PATH)/upstream-openbsd \
-I$(LOCAL_PATH)/upstream-openbsd/lib/libc/include \
-include upstream-openbsd/openbsd-compat.h \
LOCAL_CONLYFLAGS := $(libc_common_conlyflags)
LOCAL_CPPFLAGS := $(libc_common_cppflags)
LOCAL_C_INCLUDES := $(libc_common_c_includes)
LOCAL_MODULE := libc_gdtoa
LOCAL_ADDITIONAL_DEPENDENCIES := $(libc_common_additional_dependencies)
LOCAL_SYSTEM_SHARED_LIBRARIES :=
$(eval $(call patch-up-arch-specific-flags,LOCAL_CFLAGS,libc_common_cflags))
include $(BUILD_STATIC_LIBRARY)
# ========================================================
# libc_bionic.a - home-grown C library code
# ========================================================
@ -670,12 +723,14 @@ LOCAL_CFLAGS := $(libc_common_cflags)
LOCAL_CONLYFLAGS := $(libc_common_conlyflags)
LOCAL_CPPFLAGS := $(libc_common_cppflags)
LOCAL_C_INCLUDES := $(libc_common_c_includes)
LOCAL_ARM_MODE := arm # Work around arm linker bug http://b/14090368.
LOCAL_MODULE := libc_common
LOCAL_ADDITIONAL_DEPENDENCIES := $(libc_common_additional_dependencies)
LOCAL_WHOLE_STATIC_LIBRARIES := \
libc_bionic \
libc_dns \
libc_freebsd \
libc_gdtoa \
libc_netbsd \
libc_openbsd \
libc_stack_protector \
@ -773,7 +828,7 @@ LOCAL_SRC_FILES := \
bionic/pthread_debug.cpp \
bionic/libc_init_dynamic.cpp \
LOCAL_MODULE:= libc
LOCAL_MODULE := libc
LOCAL_ADDITIONAL_DEPENDENCIES := $(libc_common_additional_dependencies)
LOCAL_REQUIRED_MODULES := tzdata
@ -786,7 +841,7 @@ LOCAL_STRIP_MODULE := keep_symbols
# ensures that symbols that are pulled into those new libraries from libgcc.a are not declared
# external; if that were the case, then libc would not pull those symbols from libgcc.a as it
# should, instead relying on the external symbols from the dependent libraries. That would
# create an "cloaked" dependency on libgcc.a in libc though the libraries, which is not what
# create a "cloaked" dependency on libgcc.a in libc though the libraries, which is not what
# you wanted!
LOCAL_SHARED_LIBRARIES := libdl
@ -834,7 +889,7 @@ LOCAL_SRC_FILES := \
bionic/malloc_debug_leak.cpp \
bionic/malloc_debug_check.cpp \
LOCAL_MODULE:= libc_malloc_debug_leak
LOCAL_MODULE := libc_malloc_debug_leak
LOCAL_ADDITIONAL_DEPENDENCIES := $(libc_common_additional_dependencies)
LOCAL_SHARED_LIBRARIES := libc libdl
@ -865,7 +920,7 @@ LOCAL_C_INCLUDES := $(libc_common_c_includes)
LOCAL_SRC_FILES := \
bionic/malloc_debug_qemu.cpp
LOCAL_MODULE:= libc_malloc_debug_qemu
LOCAL_MODULE := libc_malloc_debug_qemu
LOCAL_ADDITIONAL_DEPENDENCIES := $(libc_common_additional_dependencies)
LOCAL_SHARED_LIBRARIES := libc libdl

110
libc/arch-arm/include/machine/ieee.h
View File

@ -45,49 +45,6 @@
* @(#)ieee.h 8.1 (Berkeley) 6/11/93
*/
/*
* ieee.h defines the machine-dependent layout of the machine's IEEE
* floating point.
*/
/*
* Define the number of bits in each fraction and exponent.
*
* k k+1
* Note that 1.0 x 2 == 0.1 x 2 and that denorms are represented
*
* (-exp_bias+1)
* as fractions that look like 0.fffff x 2 . This means that
*
* -126
* the number 0.10000 x 2 , for instance, is the same as the normalized
*
* -127 -128
* float 1.0 x 2 . Thus, to represent 2 , we need one leading zero
*
* -129
* in the fraction; to represent 2 , we need two, and so on. This
*
* (-exp_bias-fracbits+1)
* implies that the smallest denormalized number is 2
*
* for whichever format we are talking about: for single precision, for
*
* -126 -149
* instance, we get .00000000000000000000001 x 2 , or 1.0 x 2 , and
*
* -149 == -127 - 23 + 1.
*/
/*
* The ARM has two sets of FP data formats. The FPA supports 32-bit, 64-bit
* and 96-bit IEEE formats, with the words in big-endian order. VFP supports
* 32-bin and 64-bit IEEE formats with the words in the CPU's native byte
* order.
*
* The FPA also has two packed decimal formats, but we ignore them here.
*/
#define SNG_EXPBITS 8
#define SNG_FRACBITS 23
@ -96,97 +53,30 @@
#define DBL_FRACLBITS 32
#define DBL_FRACBITS 52
#ifndef __VFP_FP__
#define E80_EXPBITS 15
#define E80_FRACHBITS 31
#define E80_FRACLBITS 32
#define E80_FRACBITS 64
#define EXT_EXPBITS 15
#define EXT_FRACHBITS 16
#define EXT_FRACHMBITS 32
#define EXT_FRACLMBITS 32
#define EXT_FRACLBITS 32
#define EXT_FRACBITS 112
#endif
struct ieee_single {
u_int sng_frac:23;
u_int sng_exp:8;
u_int sng_sign:1;
};
#ifdef __VFP_FP__
struct ieee_double {
u_int dbl_fracl;
u_int dbl_frach:20;
u_int dbl_exp:11;
u_int dbl_sign:1;
};
#else /* !__VFP_FP__ */
struct ieee_double {
u_int dbl_frach:20;
u_int dbl_exp:11;
u_int dbl_sign:1;
u_int dbl_fracl;
};
union ieee_double_u {
double dblu_d;
struct ieee_double dblu_dbl;
};
struct ieee_e80 {
u_int e80_exp:15;
u_int e80_zero:16;
u_int e80_sign:1;
u_int e80_frach:31;
u_int e80_j:1;
u_int e80_fracl;
};
struct ieee_ext {
u_int ext_frach:16;
u_int ext_exp:15;
u_int ext_sign:1;
u_int ext_frachm;
u_int ext_fraclm;
u_int ext_fracl;
};
#endif /* !__VFP_FP__ */
/*
* Floats whose exponent is in [1..INFNAN) (of whatever type) are
* `normal'. Floats whose exponent is INFNAN are either Inf or NaN.
* Floats whose exponent is zero are either zero (iff all fraction
* bits are zero) or subnormal values.
*
* A NaN is a `signalling NaN' if its QUIETNAN bit is clear in its
* high fraction; if the bit is set, it is a `quiet NaN'.
*/
#define SNG_EXP_INFNAN 255
#define DBL_EXP_INFNAN 2047
#ifndef __VFP_FP__
#define E80_EXP_INFNAN 32767
#define EXT_EXP_INFNAN 32767
#endif /* !__VFP_FP__ */
#if 0
#define SNG_QUIETNAN (1 << 22)
#define DBL_QUIETNAN (1 << 19)
#ifndef __VFP_FP__
#define E80_QUIETNAN (1 << 15)
#define EXT_QUIETNAN (1 << 15)
#endif /* !__VFP_FP__ */
#endif
/*
* Exponent biases.
*/
#define SNG_EXP_BIAS 127
#define DBL_EXP_BIAS 1023
#ifndef __VFP_FP__
#define E80_EXP_BIAS 16383
#define EXT_EXP_BIAS 16383
#endif /* !__VFP_FP__ */

33
libc/arch-arm64/include/machine/ieee.h
View File

@ -45,11 +45,6 @@
* @(#)ieee.h 8.1 (Berkeley) 6/11/93
*/
/*
* ieee.h defines the machine-dependent layout of the machine's IEEE
* floating point.
*/
#define SNG_EXPBITS 8
#define SNG_FRACBITS 23
@ -58,6 +53,22 @@
#define DBL_FRACLBITS 32
#define DBL_FRACBITS 52
#define EXT_EXPBITS 15
#define EXT_FRACHBITS 16
#define EXT_FRACHMBITS 32
#define EXT_FRACLMBITS 32
#define EXT_FRACLBITS 32
#define EXT_FRACBITS 112
#define EXT_IMPLICIT_NBIT
#define EXT_TO_ARRAY32(p, a) do { \
(a)[0] = (uint32_t)(p)->ext_fracl; \
(a)[1] = (uint32_t)(p)->ext_fraclm; \
(a)[2] = (uint32_t)(p)->ext_frachm; \
(a)[3] = (uint32_t)(p)->ext_frach; \
} while(0)
struct ieee_single {
u_int sng_frac:23;
u_int sng_exp:8;
@ -71,14 +82,20 @@ struct ieee_double {
u_int dbl_sign:1;
};
struct ieee_ext {
u_int ext_fracl;
u_int ext_fraclm;
u_int ext_frachm;
u_int ext_frach:16;
u_int ext_exp:15;
u_int ext_sign:1;
};
/*
* Floats whose exponent is in [1..INFNAN) (of whatever type) are
* `normal'. Floats whose exponent is INFNAN are either Inf or NaN.
* Floats whose exponent is zero are either zero (iff all fraction
* bits are zero) or subnormal values.
*
* A NaN is a `signalling NaN' if its QUIETNAN bit is clear in its
* high fraction; if the bit is set, it is a `quiet NaN'.
*/
#define SNG_EXP_INFNAN 255
#define DBL_EXP_INFNAN 2047

34
libc/bionic/strtof.cpp
View File

@ -1,34 +0,0 @@
/*
* Copyright (C) 2014 The Android Open Source Project
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* * Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* * 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.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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
* COPYRIGHT OWNER 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.
*/
#include <stdlib.h>
float strtof(const char* s, char** end_ptr) {
// TODO: upgrade to a non-hack implementation.
return strtod(s, end_ptr);
}

11
libc/bionic/strtold.cpp
View File

@ -28,7 +28,16 @@
#include <stdlib.h>
extern "C" int __strtorQ(const char*, char**, int, void*);
long double strtold(const char* s, char** end_ptr) {
// TODO: this is fine for LP32 where double == long double, but is broken on LP64.
#if __LP64__
long double result;
// TODO: use the current rounding mode?
__strtorQ(s, end_ptr, 1 /* FPI_Round_near */, &result);
return result;
#else
// This is fine for LP32 where long double is just double.
return strtod(s, end_ptr);
#endif
}

6
libc/stdio/vfprintf.c
View File

@ -146,7 +146,7 @@ __sbprintf(FILE *fp, const char *fmt, va_list ap)
#define DEFPREC 6
static char *cvt(double, int, int, char *, int *, int, int *);
extern void freedtoa(char *);
extern void __freedtoa(char *);
static int exponent(char *, int, int);
#else /* no FLOATING_POINT */
#define BUF 40
@ -550,7 +550,7 @@ reswitch: switch (ch) {
break;
}
if (dtoaresult != NULL) freedtoa(dtoaresult);
if (dtoaresult != NULL) __freedtoa(dtoaresult);
flags |= FPT;
dtoaresult = cp = cvt(_double, prec, flags, &softsign,
&expt, ch, &ndig);
@ -846,7 +846,7 @@ done:
error:
#ifdef FLOATING_POINT
if (dtoaresult != NULL) {
freedtoa(dtoaresult);
__freedtoa(dtoaresult);
}
#endif
if (argtable != NULL && argtable != statargtable) {

2689
libc/stdlib/strtod.c
View File

File diff suppressed because it is too large Load Diff

28
libc/upstream-openbsd/arith.h Normal file
View File

@ -0,0 +1,28 @@
/*
* Copyright (C) 2014 The Android Open Source Project
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#define IEEE_8087
#define Arith_Kind_ASL 1
#if __LP64__
#define Long int
#define Intcast (int)(long)
#define Double_Align
#define X64_bit_pointers
#endif
#define INFNAN_CHECK
#define MULTIPLE_THREADS

48
libc/upstream-openbsd/gd_qnan.h Normal file
View File

@ -0,0 +1,48 @@
/*
* Copyright (C) 2014 The Android Open Source Project
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#if __arm__
#define f_QNAN 0xffffffff
#define d_QNAN0 0xffffffff
#define d_QNAN1 0xffffffff
#elif __mips__
#define f_QNAN 0x7fbfffff
#define d_QNAN0 0x7ff7ffff
#define d_QNAN1 0xffffffff
#else
#define f_QNAN 0xffc00000
#define d_QNAN0 0x00000000
#define d_QNAN1 0xfff80000
#endif
/* long double. */
#if __LP64__
#define ld_QNAN0 0x7fff8000
#define ld_QNAN1 0x00000000
#define ld_QNAN2 0x00000000
#define ld_QNAN3 0x00000000
#else
/* sizeof(long double) == sizeof(double), so we shouldn't be trying to use these constants. */
#endif

19
libc/upstream-openbsd/gdtoa_support.cpp Normal file
View File

@ -0,0 +1,19 @@
/*
* Copyright (C) 2014 The Android Open Source Project
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include <pthread.h>
__LIBC_HIDDEN__ pthread_mutex_t __dtoa_locks[] = { PTHREAD_MUTEX_INITIALIZER, PTHREAD_MUTEX_INITIALIZER };

224
libc/upstream-openbsd/lib/libc/gdtoa/dmisc.c Normal file
View File

@ -0,0 +1,224 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
#ifndef MULTIPLE_THREADS
char *dtoa_result;
#endif
char *
#ifdef KR_headers
rv_alloc(i) int i;
#else
rv_alloc(int i)
#endif
{
int j, k, *r;
j = sizeof(ULong);
for(k = 0;
sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
j <<= 1)
k++;
r = (int*)Balloc(k);
if (r == NULL)
return (
#ifndef MULTIPLE_THREADS
dtoa_result =
#endif
NULL);
*r = k;
return
#ifndef MULTIPLE_THREADS
dtoa_result =
#endif
(char *)(r+1);
}
char *
#ifdef KR_headers
nrv_alloc(s, rve, n) char *s, **rve; int n;
#else
nrv_alloc(char *s, char **rve, int n)
#endif
{
char *rv, *t;
t = rv = rv_alloc(n);
if (t == NULL)
return (NULL);
while((*t = *s++) !=0)
t++;
if (rve)
*rve = t;
return rv;
}
/* freedtoa(s) must be used to free values s returned by dtoa
* when MULTIPLE_THREADS is #defined. It should be used in all cases,
* but for consistency with earlier versions of dtoa, it is optional
* when MULTIPLE_THREADS is not defined.
*/
void
#ifdef KR_headers
freedtoa(s) char *s;
#else
freedtoa(char *s)
#endif
{
Bigint *b = (Bigint *)((int *)s - 1);
b->maxwds = 1 << (b->k = *(int*)b);
Bfree(b);
#ifndef MULTIPLE_THREADS
if (s == dtoa_result)
dtoa_result = 0;
#endif
}
int
quorem
#ifdef KR_headers
(b, S) Bigint *b, *S;
#else
(Bigint *b, Bigint *S)
#endif
{
int n;
ULong *bx, *bxe, q, *sx, *sxe;
#ifdef ULLong
ULLong borrow, carry, y, ys;
#else
ULong borrow, carry, y, ys;
#ifdef Pack_32
ULong si, z, zs;
#endif
#endif
n = S->wds;
#ifdef DEBUG
/*debug*/ if (b->wds > n)
/*debug*/ Bug("oversize b in quorem");
#endif
if (b->wds < n)
return 0;
sx = S->x;
sxe = sx + --n;
bx = b->x;
bxe = bx + n;
q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
#ifdef DEBUG
/*debug*/ if (q > 9)
/*debug*/ Bug("oversized quotient in quorem");
#endif
if (q) {
borrow = 0;
carry = 0;
do {
#ifdef ULLong
ys = *sx++ * (ULLong)q + carry;
carry = ys >> 32;
y = *bx - (ys & 0xffffffffUL) - borrow;
borrow = y >> 32 & 1UL;
*bx++ = y & 0xffffffffUL;
#else
#ifdef Pack_32
si = *sx++;
ys = (si & 0xffff) * q + carry;
zs = (si >> 16) * q + (ys >> 16);
carry = zs >> 16;
y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
z = (*bx >> 16) - (zs & 0xffff) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(bx, z, y);
#else
ys = *sx++ * q + carry;
carry = ys >> 16;
y = *bx - (ys & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
*bx++ = y & 0xffff;
#endif
#endif
}
while(sx <= sxe);
if (!*bxe) {
bx = b->x;
while(--bxe > bx && !*bxe)
--n;
b->wds = n;
}
}
if (cmp(b, S) >= 0) {
q++;
borrow = 0;
carry = 0;
bx = b->x;
sx = S->x;
do {
#ifdef ULLong
ys = *sx++ + carry;
carry = ys >> 32;
y = *bx - (ys & 0xffffffffUL) - borrow;
borrow = y >> 32 & 1UL;
*bx++ = y & 0xffffffffUL;
#else
#ifdef Pack_32
si = *sx++;
ys = (si & 0xffff) + carry;
zs = (si >> 16) + (ys >> 16);
carry = zs >> 16;
y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
z = (*bx >> 16) - (zs & 0xffff) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(bx, z, y);
#else
ys = *sx++ + carry;
carry = ys >> 16;
y = *bx - (ys & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
*bx++ = y & 0xffff;
#endif
#endif
}
while(sx <= sxe);
bx = b->x;
bxe = bx + n;
if (!*bxe) {
while(--bxe > bx && !*bxe)
--n;
b->wds = n;
}
}
return q;
}

839
libc/upstream-openbsd/lib/libc/gdtoa/dtoa.c Normal file
View File

@ -0,0 +1,839 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
*
* Inspired by "How to Print Floating-Point Numbers Accurately" by
* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
*
* Modifications:
* 1. Rather than iterating, we use a simple numeric overestimate
* to determine k = floor(log10(d)). We scale relevant
* quantities using O(log2(k)) rather than O(k) multiplications.
* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
* try to generate digits strictly left to right. Instead, we
* compute with fewer bits and propagate the carry if necessary
* when rounding the final digit up. This is often faster.
* 3. Under the assumption that input will be rounded nearest,
* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
* That is, we allow equality in stopping tests when the
* round-nearest rule will give the same floating-point value
* as would satisfaction of the stopping test with strict
* inequality.
* 4. We remove common factors of powers of 2 from relevant
* quantities.
* 5. When converting floating-point integers less than 1e16,
* we use floating-point arithmetic rather than resorting
* to multiple-precision integers.
* 6. When asked to produce fewer than 15 digits, we first try
* to get by with floating-point arithmetic; we resort to
* multiple-precision integer arithmetic only if we cannot
* guarantee that the floating-point calculation has given
* the correctly rounded result. For k requested digits and
* "uniformly" distributed input, the probability is
* something like 10^(k-15) that we must resort to the Long
* calculation.
*/
#ifdef Honor_FLT_ROUNDS
#undef Check_FLT_ROUNDS
#define Check_FLT_ROUNDS
#else
#define Rounding Flt_Rounds
#endif
char *
dtoa
#ifdef KR_headers
(d0, mode, ndigits, decpt, sign, rve)
double d0; int mode, ndigits, *decpt, *sign; char **rve;
#else
(double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
#endif
{
/* Arguments ndigits, decpt, sign are similar to those
of ecvt and fcvt; trailing zeros are suppressed from
the returned string. If not null, *rve is set to point
to the end of the return value. If d is +-Infinity or NaN,
then *decpt is set to 9999.
mode:
0 ==> shortest string that yields d when read in
and rounded to nearest.
1 ==> like 0, but with Steele & White stopping rule;
e.g. with IEEE P754 arithmetic , mode 0 gives
1e23 whereas mode 1 gives 9.999999999999999e22.
2 ==> max(1,ndigits) significant digits. This gives a
return value similar to that of ecvt, except
that trailing zeros are suppressed.
3 ==> through ndigits past the decimal point. This
gives a return value similar to that from fcvt,
except that trailing zeros are suppressed, and
ndigits can be negative.
4,5 ==> similar to 2 and 3, respectively, but (in
round-nearest mode) with the tests of mode 0 to
possibly return a shorter string that rounds to d.
With IEEE arithmetic and compilation with
-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
as modes 2 and 3 when FLT_ROUNDS != 1.
6-9 ==> Debugging modes similar to mode - 4: don't try
fast floating-point estimate (if applicable).
Values of mode other than 0-9 are treated as mode 0.
Sufficient space is allocated to the return value
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
spec_case, try_quick;
Long L;
#ifndef Sudden_Underflow
int denorm;
ULong x;
#endif
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
U d, d2, eps;
double ds;
char *s, *s0;
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
#ifdef Honor_FLT_ROUNDS /*{*/
int Rounding;
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
Rounding = Flt_Rounds;
#else /*}{*/
Rounding = 1;
switch(fegetround()) {
case FE_TOWARDZERO: Rounding = 0; break;
case FE_UPWARD: Rounding = 2; break;
case FE_DOWNWARD: Rounding = 3;
}
#endif /*}}*/
#endif /*}*/
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
freedtoa(dtoa_result);
dtoa_result = 0;
}
#endif
d.d = d0;
if (word0(&d) & Sign_bit) {
/* set sign for everything, including 0's and NaNs */
*sign = 1;
word0(&d) &= ~Sign_bit; /* clear sign bit */
}
else
*sign = 0;
#if defined(IEEE_Arith) + defined(VAX)
#ifdef IEEE_Arith
if ((word0(&d) & Exp_mask) == Exp_mask)
#else
if (word0(&d) == 0x8000)
#endif
{
/* Infinity or NaN */
*decpt = 9999;
#ifdef IEEE_Arith
if (!word1(&d) && !(word0(&d) & 0xfffff))
return nrv_alloc("Infinity", rve, 8);
#endif
return nrv_alloc("NaN", rve, 3);
}
#endif
#ifdef IBM
dval(&d) += 0; /* normalize */
#endif
if (!dval(&d)) {
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
#ifdef SET_INEXACT
try_quick = oldinexact = get_inexact();
inexact = 1;
#endif
#ifdef Honor_FLT_ROUNDS
if (Rounding >= 2) {
if (*sign)
Rounding = Rounding == 2 ? 0 : 2;
else
if (Rounding != 2)
Rounding = 0;
}
#endif
b = d2b(dval(&d), &be, &bbits);
if (b == NULL)
return (NULL);
#ifdef Sudden_Underflow
i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
#else
if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
#endif
dval(&d2) = dval(&d);
word0(&d2) &= Frac_mask1;
word0(&d2) |= Exp_11;
#ifdef IBM
if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
dval(&d2) /= 1 << j;
#endif
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
*
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
*
* We want k to be too large rather than too small.
* The error in the first-order Taylor series approximation
* is in our favor, so we just round up the constant enough
* to compensate for any error in the multiplication of
* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
* adding 1e-13 to the constant term more than suffices.
* Hence we adjust the constant term to 0.1760912590558.
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
i -= Bias;
#ifdef IBM
i <<= 2;
i += j;
#endif
#ifndef Sudden_Underflow
denorm = 0;
}
else {
/* d is denormalized */
i = bbits + be + (Bias + (P-1) - 1);
x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
: word1(&d) << (32 - i);
dval(&d2) = x;
word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
i -= (Bias + (P-1) - 1) + 1;
denorm = 1;
}
#endif
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
k = (int)ds;
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
if (k >= 0 && k <= Ten_pmax) {
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
j = bbits - i - 1;
if (j >= 0) {
b2 = 0;
s2 = j;
}
else {
b2 = -j;
s2 = 0;
}
if (k >= 0) {
b5 = 0;
s5 = k;
s2 += k;
}
else {
b2 -= k;
b5 = -k;
s5 = 0;
}
if (mode < 0 || mode > 9)
mode = 0;
#ifndef SET_INEXACT
#ifdef Check_FLT_ROUNDS
try_quick = Rounding == 1;
#else
try_quick = 1;
#endif
#endif /*SET_INEXACT*/
if (mode > 5) {
mode -= 4;
try_quick = 0;
}
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
i = 18;
ndigits = 0;
break;
case 2:
leftright = 0;
/* no break */
case 4:
if (ndigits <= 0)
ndigits = 1;
ilim = ilim1 = i = ndigits;
break;
case 3:
leftright = 0;
/* no break */
case 5:
i = ndigits + k + 1;
ilim = i;
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
s = s0 = rv_alloc(i);
if (s == NULL)
return (NULL);
#ifdef Honor_FLT_ROUNDS
if (mode > 1 && Rounding != 1)
leftright = 0;
#endif
if (ilim >= 0 && ilim <= Quick_max && try_quick) {
/* Try to get by with floating-point arithmetic. */
i = 0;
dval(&d2) = dval(&d);
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
if (k > 0) {
ds = tens[k&0xf];
j = k >> 4;
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
for(; j; j >>= 1, i++)
if (j & 1) {
ieps++;
ds *= bigtens[i];
}
dval(&d) /= ds;
}
else if (( j1 = -k )!=0) {
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(&d) *= bigtens[i];
}
}
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(&d) *= 10.;
ieps++;
}
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
#ifndef No_leftright
if (leftright) {
/* Use Steele & White method of only
* generating digits needed.
*/
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = dval(&d);
dval(&d) -= L;
*s++ = '0' + (int)L;
if (dval(&d) < dval(&eps))
goto ret1;
if (1. - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d));
if (!(dval(&d) -= L))
ilim = i;
*s++ = '0' + (int)L;
if (i == ilim) {
if (dval(&d) > 0.5 + dval(&eps))
goto bump_up;
else if (dval(&d) < 0.5 - dval(&eps)) {
while(*--s == '0');
s++;
goto ret1;
}
break;
}
}
#ifndef No_leftright
}
#endif
fast_failed:
s = s0;
dval(&d) = dval(&d2);
k = k0;
ilim = ilim0;
}
/* Do we have a "small" integer? */
if (be >= 0 && k <= Int_max) {
/* Yes. */
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d) / ds);
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(&d) < 0) {
L--;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (!dval(&d)) {
#ifdef SET_INEXACT
inexact = 0;
#endif
break;
}
if (i == ilim) {
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
switch(Rounding) {
case 0: goto ret1;
case 2: goto bump_up;
}
#endif
dval(&d) += dval(&d);
#ifdef ROUND_BIASED
if (dval(&d) >= ds)
#else
if (dval(&d) > ds || (dval(&d) == ds && L & 1))
#endif
{
bump_up:
while(*--s == '9')
if (s == s0) {
k++;
*s = '0';
break;
}
++*s++;
}
break;
}
}
goto ret1;
}
m2 = b2;
m5 = b5;
mhi = mlo = 0;
if (leftright) {
i =
#ifndef Sudden_Underflow
denorm ? be + (Bias + (P-1) - 1 + 1) :
#endif
#ifdef IBM
1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
#else
1 + P - bbits;
#endif
b2 += i;
s2 += i;
mhi = i2b(1);
if (mhi == NULL)
return (NULL);
}
if (m2 > 0 && s2 > 0) {
i = m2 < s2 ? m2 : s2;
b2 -= i;
m2 -= i;
s2 -= i;
}
if (b5 > 0) {
if (leftright) {
if (m5 > 0) {
mhi = pow5mult(mhi, m5);
if (mhi == NULL)
return (NULL);
b1 = mult(mhi, b);
if (b1 == NULL)
return (NULL);
Bfree(b);
b = b1;
}
if (( j = b5 - m5 )!=0) {
b = pow5mult(b, j);
if (b == NULL)
return (NULL);
}
}
else {
b = pow5mult(b, b5);
if (b == NULL)
return (NULL);
}
}
S = i2b(1);
if (S == NULL)
return (NULL);
if (s5 > 0) {
S = pow5mult(S, s5);
if (S == NULL)
return (NULL);
}
/* Check for special case that d is a normalized power of 2. */
spec_case = 0;
if ((mode < 2 || leftright)
#ifdef Honor_FLT_ROUNDS
&& Rounding == 1
#endif
) {
if (!word1(&d) && !(word0(&d) & Bndry_mask)
#ifndef Sudden_Underflow
&& word0(&d) & (Exp_mask & ~Exp_msk1)
#endif
) {
/* The special case */
b2 += Log2P;
s2 += Log2P;
spec_case = 1;
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
*
* Perhaps we should just compute leading 28 bits of S once
* and for all and pass them and a shift to quorem, so it
* can do shifts and ors to compute the numerator for q.
*/
#ifdef Pack_32
if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
i = 32 - i;
#else
if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
i = 16 - i;
#endif
if (i > 4) {
i -= 4;
b2 += i;
m2 += i;
s2 += i;
}
else if (i < 4) {
i += 28;
b2 += i;
m2 += i;
s2 += i;
}
if (b2 > 0) {
b = lshift(b, b2);
if (b == NULL)
return (NULL);
}
if (s2 > 0) {
S = lshift(S, s2);
if (S == NULL)
return (NULL);
}
if (k_check) {
if (cmp(b,S) < 0) {
k--;
b = multadd(b, 10, 0); /* we botched the k estimate */
if (b == NULL)
return (NULL);
if (leftright) {
mhi = multadd(mhi, 10, 0);
if (mhi == NULL)
return (NULL);
}
ilim = ilim1;
}
}
if (ilim <= 0 && (mode == 3 || mode == 5)) {
S = multadd(S,5,0);
if (S == NULL)
return (NULL);
if (ilim < 0 || cmp(b,S) <= 0) {
/* no digits, fcvt style */
no_digits:
k = -1 - ndigits;
goto ret;
}
one_digit:
*s++ = '1';
k++;
goto ret;
}
if (leftright) {
if (m2 > 0) {
mhi = lshift(mhi, m2);
if (mhi == NULL)
return (NULL);
}
/* Compute mlo -- check for special case
* that d is a normalized power of 2.
*/
mlo = mhi;
if (spec_case) {
mhi = Balloc(mhi->k);
if (mhi == NULL)
return (NULL);
Bcopy(mhi, mlo);
mhi = lshift(mhi, Log2P);
if (mhi == NULL)
return (NULL);
}
for(i = 1;;i++) {
dig = quorem(b,S) + '0';
/* Do we yet have the shortest decimal string
* that will round to d?
*/
j = cmp(b, mlo);
delta = diff(S, mhi);
if (delta == NULL)
return (NULL);
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
#ifdef Honor_FLT_ROUNDS
&& Rounding >= 1
#endif
) {
if (dig == '9')
goto round_9_up;
if (j > 0)
dig++;
#ifdef SET_INEXACT
else if (!b->x[0] && b->wds <= 1)
inexact = 0;
#endif
*s++ = dig;
goto ret;
}
#endif
if (j < 0 || (j == 0 && mode != 1
#ifndef ROUND_BIASED
&& !(word1(&d) & 1)
#endif
)) {
if (!b->x[0] && b->wds <= 1) {
#ifdef SET_INEXACT
inexact = 0;
#endif
goto accept_dig;
}
#ifdef Honor_FLT_ROUNDS
if (mode > 1)
switch(Rounding) {
case 0: goto accept_dig;
case 2: goto keep_dig;
}
#endif /*Honor_FLT_ROUNDS*/
if (j1 > 0) {
b = lshift(b, 1);
if (b == NULL)
return (NULL);
j1 = cmp(b, S);
#ifdef ROUND_BIASED
if (j1 >= 0 /*)*/
#else
if ((j1 > 0 || (j1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
}
accept_dig:
*s++ = dig;
goto ret;
}
if (j1 > 0) {
#ifdef Honor_FLT_ROUNDS
if (!Rounding)
goto accept_dig;
#endif
if (dig == '9') { /* possible if i == 1 */
round_9_up:
*s++ = '9';
goto roundoff;
}
*s++ = dig + 1;
goto ret;
}
#ifdef Honor_FLT_ROUNDS
keep_dig:
#endif
*s++ = dig;
if (i == ilim)
break;
b = multadd(b, 10, 0);
if (b == NULL)
return (NULL);
if (mlo == mhi) {
mlo = mhi = multadd(mhi, 10, 0);
if (mlo == NULL)
return (NULL);
}
else {
mlo = multadd(mlo, 10, 0);
if (mlo == NULL)
return (NULL);
mhi = multadd(mhi, 10, 0);
if (mhi == NULL)
return (NULL);
}
}
}
else
for(i = 1;; i++) {
*s++ = dig = quorem(b,S) + '0';
if (!b->x[0] && b->wds <= 1) {
#ifdef SET_INEXACT
inexact = 0;
#endif
goto ret;
}
if (i >= ilim)
break;
b = multadd(b, 10, 0);
if (b == NULL)
return (NULL);
}
/* Round off last digit */
#ifdef Honor_FLT_ROUNDS
switch(Rounding) {
case 0: goto trimzeros;
case 2: goto roundoff;
}
#endif
b = lshift(b, 1);
if (b == NULL)
return (NULL);
j = cmp(b, S);
#ifdef ROUND_BIASED
if (j >= 0)
#else
if (j > 0 || (j == 0 && dig & 1))
#endif
{
roundoff:
while(*--s == '9')
if (s == s0) {
k++;
*s++ = '1';
goto ret;
}
++*s++;
}
else {
#ifdef Honor_FLT_ROUNDS
trimzeros:
#endif
while(*--s == '0');
s++;
}
ret:
Bfree(S);
if (mhi) {
if (mlo && mlo != mhi)
Bfree(mlo);
Bfree(mhi);
}
ret1:
#ifdef SET_INEXACT
if (inexact) {
if (!oldinexact) {
word0(&d) = Exp_1 + (70 << Exp_shift);
word1(&d) = 0;
dval(&d) += 1.;
}
}
else if (!oldinexact)
clear_inexact();
#endif
Bfree(b);
*s = 0;
*decpt = k + 1;
if (rve)
*rve = s;
return s0;
}

829
libc/upstream-openbsd/lib/libc/gdtoa/gdtoa.c Normal file
View File

@ -0,0 +1,829 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
static Bigint *
#ifdef KR_headers
bitstob(bits, nbits, bbits) ULong *bits; int nbits; int *bbits;
#else
bitstob(ULong *bits, int nbits, int *bbits)
#endif
{
int i, k;
Bigint *b;
ULong *be, *x, *x0;
i = ULbits;
k = 0;
while(i < nbits) {
i <<= 1;
k++;
}
#ifndef Pack_32
if (!k)
k = 1;
#endif
b = Balloc(k);
if (b == NULL)
return (NULL);
be = bits + ((nbits - 1) >> kshift);
x = x0 = b->x;
do {
*x++ = *bits & ALL_ON;
#ifdef Pack_16
*x++ = (*bits >> 16) & ALL_ON;
#endif
} while(++bits <= be);
i = x - x0;
while(!x0[--i])
if (!i) {
b->wds = 0;
*bbits = 0;
goto ret;
}
b->wds = i + 1;
*bbits = i*ULbits + 32 - hi0bits(b->x[i]);
ret:
return b;
}
/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
*
* Inspired by "How to Print Floating-Point Numbers Accurately" by
* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
*
* Modifications:
* 1. Rather than iterating, we use a simple numeric overestimate
* to determine k = floor(log10(d)). We scale relevant
* quantities using O(log2(k)) rather than O(k) multiplications.
* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
* try to generate digits strictly left to right. Instead, we
* compute with fewer bits and propagate the carry if necessary
* when rounding the final digit up. This is often faster.
* 3. Under the assumption that input will be rounded nearest,
* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
* That is, we allow equality in stopping tests when the
* round-nearest rule will give the same floating-point value
* as would satisfaction of the stopping test with strict
* inequality.
* 4. We remove common factors of powers of 2 from relevant
* quantities.
* 5. When converting floating-point integers less than 1e16,
* we use floating-point arithmetic rather than resorting
* to multiple-precision integers.
* 6. When asked to produce fewer than 15 digits, we first try
* to get by with floating-point arithmetic; we resort to
* multiple-precision integer arithmetic only if we cannot
* guarantee that the floating-point calculation has given
* the correctly rounded result. For k requested digits and
* "uniformly" distributed input, the probability is
* something like 10^(k-15) that we must resort to the Long
* calculation.
*/
char *
gdtoa
#ifdef KR_headers
(fpi, be, bits, kindp, mode, ndigits, decpt, rve)
FPI *fpi; int be; ULong *bits;
int *kindp, mode, ndigits, *decpt; char **rve;
#else
(FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, int *decpt, char **rve)
#endif
{
/* Arguments ndigits and decpt are similar to the second and third
arguments of ecvt and fcvt; trailing zeros are suppressed from
the returned string. If not null, *rve is set to point
to the end of the return value. If d is +-Infinity or NaN,
then *decpt is set to 9999.
be = exponent: value = (integer represented by bits) * (2 to the power of be).
mode:
0 ==> shortest string that yields d when read in
and rounded to nearest.
1 ==> like 0, but with Steele & White stopping rule;
e.g. with IEEE P754 arithmetic , mode 0 gives
1e23 whereas mode 1 gives 9.999999999999999e22.
2 ==> max(1,ndigits) significant digits. This gives a
return value similar to that of ecvt, except
that trailing zeros are suppressed.
3 ==> through ndigits past the decimal point. This
gives a return value similar to that from fcvt,
except that trailing zeros are suppressed, and
ndigits can be negative.
4-9 should give the same return values as 2-3, i.e.,
4 <= mode <= 9 ==> same return as mode
2 + (mode & 1). These modes are mainly for
debugging; often they run slower but sometimes
faster than modes 2-3.
4,5,8,9 ==> left-to-right digit generation.
6-9 ==> don't try fast floating-point estimate
(if applicable).
Values of mode other than 0-9 are treated as mode 0.
Sufficient space is allocated to the return value
to hold the suppressed trailing zeros.
*/
int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
int j, j1, k, k0, k_check, kind, leftright, m2, m5, nbits;
int rdir, s2, s5, spec_case, try_quick;
Long L;
Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
double d2, ds;
char *s, *s0;
U d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
freedtoa(dtoa_result);
dtoa_result = 0;
}
#endif
inex = 0;
kind = *kindp &= ~STRTOG_Inexact;
switch(kind & STRTOG_Retmask) {
case STRTOG_Zero:
goto ret_zero;
case STRTOG_Normal:
case STRTOG_Denormal:
break;
case STRTOG_Infinite:
*decpt = -32768;
return nrv_alloc("Infinity", rve, 8);
case STRTOG_NaN:
*decpt = -32768;
return nrv_alloc("NaN", rve, 3);
default:
return 0;
}
b = bitstob(bits, nbits = fpi->nbits, &bbits);
if (b == NULL)
return (NULL);
be0 = be;
if ( (i = trailz(b)) !=0) {
rshift(b, i);
be += i;
bbits -= i;
}
if (!b->wds) {
Bfree(b);
ret_zero:
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
dval(&d) = b2d(b, &i);
i = be + bbits - 1;
word0(&d) &= Frac_mask1;
word0(&d) |= Exp_11;
#ifdef IBM
if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
dval(&d) /= 1 << j;
#endif
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
*
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
*
* We want k to be too large rather than too small.
* The error in the first-order Taylor series approximation
* is in our favor, so we just round up the constant enough
* to compensate for any error in the multiplication of
* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
* adding 1e-13 to the constant term more than suffices.
* Hence we adjust the constant term to 0.1760912590558.
* (We could get a more accurate k by invoking log10,
* but this is probably not worthwhile.)
*/
#ifdef IBM
i <<= 2;
i += j;
#endif
ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
/* correct assumption about exponent range */
if ((j = i) < 0)
j = -j;
if ((j -= 1077) > 0)
ds += j * 7e-17;
k = (int)ds;
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
#ifdef IBM
j = be + bbits - 1;
if ( (j1 = j & 3) !=0)
dval(&d) *= 1 << j1;
word0(&d) += j << Exp_shift - 2 & Exp_mask;
#else
word0(&d) += (be + bbits - 1) << Exp_shift;
#endif
if (k >= 0 && k <= Ten_pmax) {
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
j = bbits - i - 1;
if (j >= 0) {
b2 = 0;
s2 = j;
}
else {
b2 = -j;
s2 = 0;
}
if (k >= 0) {
b5 = 0;
s5 = k;
s2 += k;
}
else {
b2 -= k;
b5 = -k;
s5 = 0;
}
if (mode < 0 || mode > 9)
mode = 0;
try_quick = 1;
if (mode > 5) {
mode -= 4;
try_quick = 0;
}
else if (i >= -4 - Emin || i < Emin)
try_quick = 0;
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
i = (int)(nbits * .30103) + 3;
ndigits = 0;
break;
case 2:
leftright = 0;
/* no break */
case 4:
if (ndigits <= 0)
ndigits = 1;
ilim = ilim1 = i = ndigits;
break;
case 3:
leftright = 0;
/* no break */
case 5:
i = ndigits + k + 1;
ilim = i;
ilim1 = i - 1;
if (i <= 0)
i = 1;
}
s = s0 = rv_alloc(i);
if (s == NULL)
return (NULL);
if ( (rdir = fpi->rounding - 1) !=0) {
if (rdir < 0)
rdir = 2;
if (kind & STRTOG_Neg)
rdir = 3 - rdir;
}
/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
#ifndef IMPRECISE_INEXACT
&& k == 0
#endif
) {
/* Try to get by with floating-point arithmetic. */
i = 0;
d2 = dval(&d);
#ifdef IBM
if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
dval(&d) /= 1 << j;
#endif
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
if (k > 0) {
ds = tens[k&0xf];
j = k >> 4;
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
for(; j; j >>= 1, i++)
if (j & 1) {
ieps++;
ds *= bigtens[i];
}
}
else {
ds = 1.;
if ( (j1 = -k) !=0) {
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(&d) *= bigtens[i];
}
}
}
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(&d) *= 10.;
ieps++;
}
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
#ifndef No_leftright
if (leftright) {
/* Use Steele & White method of only
* generating digits needed.
*/
dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = (Long)(dval(&d)/ds);
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (dval(&d) < dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
if (ds - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
if ( (L = (Long)(dval(&d)/ds)) !=0)
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (i == ilim) {
ds *= 0.5;
if (dval(&d) > ds + dval(&eps))
goto bump_up;
else if (dval(&d) < ds - dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto clear_trailing0;
}
break;
}
}
#ifndef No_leftright
}
#endif
fast_failed:
s = s0;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
/* Do we have a "small" integer? */
if (be >= 0 && k <= Int_max) {
/* Yes. */
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(&d) *= 10.) {
L = dval(&d) / ds;
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(&d) < 0) {
L--;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (dval(&d) == 0.)
break;
if (i == ilim) {
if (rdir) {
if (rdir == 1)
goto bump_up;
inex = STRTOG_Inexlo;
goto ret1;
}
dval(&d) += dval(&d);
#ifdef ROUND_BIASED
if (dval(&d) >= ds)
#else
if (dval(&d) > ds || (dval(&d) == ds && L & 1))
#endif
{
bump_up:
inex = STRTOG_Inexhi;
while(*--s == '9')
if (s == s0) {
k++;
*s = '0';
break;
}
++*s++;
}
else {
inex = STRTOG_Inexlo;
clear_trailing0:
while(*--s == '0'){}
++s;
}
break;
}
}
goto ret1;
}
m2 = b2;
m5 = b5;
mhi = mlo = 0;
if (leftright) {
i = nbits - bbits;
if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
/* denormal */
i = be - fpi->emin + 1;
if (mode >= 2 && ilim > 0 && ilim < i)
goto small_ilim;
}
else if (mode >= 2) {
small_ilim:
j = ilim - 1;
if (m5 >= j)
m5 -= j;
else {
s5 += j -= m5;
b5 += j;
m5 = 0;
}
if ((i = ilim) < 0) {
m2 -= i;
i = 0;
}
}
b2 += i;
s2 += i;
mhi = i2b(1);
if (mhi == NULL)
return (NULL);
}
if (m2 > 0 && s2 > 0) {
i = m2 < s2 ? m2 : s2;
b2 -= i;
m2 -= i;
s2 -= i;
}
if (b5 > 0) {
if (leftright) {
if (m5 > 0) {
mhi = pow5mult(mhi, m5);
if (mhi == NULL)
return (NULL);
b1 = mult(mhi, b);
if (b1 == NULL)
return (NULL);
Bfree(b);
b = b1;
}
if ( (j = b5 - m5) !=0) {
b = pow5mult(b, j);
if (b == NULL)
return (NULL);
}
}
else {
b = pow5mult(b, b5);
if (b == NULL)
return (NULL);
}
}
S = i2b(1);
if (S == NULL)
return (NULL);
if (s5 > 0) {
S = pow5mult(S, s5);
if (S == NULL)
return (NULL);
}
/* Check for special case that d is a normalized power of 2. */
spec_case = 0;
if (mode < 2) {
if (bbits == 1 && be0 > fpi->emin + 1) {
/* The special case */
b2++;
s2++;
spec_case = 1;
}
}
/* Arrange for convenient computation of quotients:
* shift left if necessary so divisor has 4 leading 0 bits.
*
* Perhaps we should just compute leading 28 bits of S once
* and for all and pass them and a shift to quorem, so it
* can do shifts and ors to compute the numerator for q.
*/
i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
m2 += i;
if ((b2 += i) > 0) {
b = lshift(b, b2);
if (b == NULL)
return (NULL);
}
if ((s2 += i) > 0) {
S = lshift(S, s2);
if (S == NULL)
return (NULL);
}
if (k_check) {
if (cmp(b,S) < 0) {
k--;
b = multadd(b, 10, 0); /* we botched the k estimate */
if (b == NULL)
return (NULL);
if (leftright) {
mhi = multadd(mhi, 10, 0);
if (mhi == NULL)
return (NULL);
}
ilim = ilim1;
}
}
if (ilim <= 0 && mode > 2) {
S = multadd(S,5,0);
if (S == NULL)
return (NULL);
if (ilim < 0 || cmp(b,S) <= 0) {
/* no digits, fcvt style */
no_digits:
k = -1 - ndigits;
inex = STRTOG_Inexlo;
goto ret;
}
one_digit:
inex = STRTOG_Inexhi;
*s++ = '1';
k++;
goto ret;
}
if (leftright) {
if (m2 > 0) {
mhi = lshift(mhi, m2);
if (mhi == NULL)
return (NULL);
}
/* Compute mlo -- check for special case
* that d is a normalized power of 2.
*/
mlo = mhi;
if (spec_case) {
mhi = Balloc(mhi->k);
if (mhi == NULL)
return (NULL);
Bcopy(mhi, mlo);
mhi = lshift(mhi, 1);
if (mhi == NULL)
return (NULL);
}
for(i = 1;;i++) {
dig = quorem(b,S) + '0';
/* Do we yet have the shortest decimal string
* that will round to d?
*/
j = cmp(b, mlo);
delta = diff(S, mhi);
if (delta == NULL)
return (NULL);
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (j1 == 0 && !mode && !(bits[0] & 1) && !rdir) {
if (dig == '9')
goto round_9_up;
if (j <= 0) {
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
}
else {
dig++;
inex = STRTOG_Inexhi;
}
*s++ = dig;
goto ret;
}
#endif
if (j < 0 || (j == 0 && !mode
#ifndef ROUND_BIASED
&& !(bits[0] & 1)
#endif
)) {
if (rdir && (b->wds > 1 || b->x[0])) {
if (rdir == 2) {
inex = STRTOG_Inexlo;
goto accept;
}
while (cmp(S,mhi) > 0) {
*s++ = dig;
mhi1 = multadd(mhi, 10, 0);
if (mhi1 == NULL)
return (NULL);
if (mlo == mhi)
mlo = mhi1;
mhi = mhi1;
b = multadd(b, 10, 0);
if (b == NULL)
return (NULL);
dig = quorem(b,S) + '0';
}
if (dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
goto accept;
}
if (j1 > 0) {
b = lshift(b, 1);
if (b == NULL)
return (NULL);
j1 = cmp(b, S);
#ifdef ROUND_BIASED
if (j1 >= 0 /*)*/
#else
if ((j1 > 0 || (j1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
}
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
accept:
*s++ = dig;
goto ret;
}
if (j1 > 0 && rdir != 2) {
if (dig == '9') { /* possible if i == 1 */
round_9_up:
*s++ = '9';
inex = STRTOG_Inexhi;
goto roundoff;
}
inex = STRTOG_Inexhi;
*s++ = dig + 1;
goto ret;
}
*s++ = dig;
if (i == ilim)
break;
b = multadd(b, 10, 0);
if (b == NULL)
return (NULL);
if (mlo == mhi) {
mlo = mhi = multadd(mhi, 10, 0);
if (mlo == NULL)
return (NULL);
}
else {
mlo = multadd(mlo, 10, 0);
if (mlo == NULL)
return (NULL);
mhi = multadd(mhi, 10, 0);
if (mhi == NULL)
return (NULL);
}
}
}
else
for(i = 1;; i++) {
*s++ = dig = quorem(b,S) + '0';
if (i >= ilim)
break;
b = multadd(b, 10, 0);
if (b == NULL)
return (NULL);
}
/* Round off last digit */
if (rdir) {
if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
goto chopzeros;
goto roundoff;
}
b = lshift(b, 1);
if (b == NULL)
return (NULL);
j = cmp(b, S);
#ifdef ROUND_BIASED
if (j >= 0)
#else
if (j > 0 || (j == 0 && dig & 1))
#endif
{
roundoff:
inex = STRTOG_Inexhi;
while(*--s == '9')
if (s == s0) {
k++;
*s++ = '1';
goto ret;
}
++*s++;
}
else {
chopzeros:
if (b->wds > 1 || b->x[0])
inex = STRTOG_Inexlo;
while(*--s == '0'){}
++s;
}
ret:
Bfree(S);
if (mhi) {
if (mlo && mlo != mhi)
Bfree(mlo);
Bfree(mhi);
}
ret1:
Bfree(b);
*s = 0;
*decpt = k + 1;
if (rve)
*rve = s;
*kindp |= inex;
return s0;
}

155
libc/upstream-openbsd/lib/libc/gdtoa/gdtoa.h Normal file
View File

@ -0,0 +1,155 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#ifndef GDTOA_H_INCLUDED
#define GDTOA_H_INCLUDED
#include "arith.h"
#include <stddef.h> /* for size_t */
#ifndef Long
#define Long int
#endif
#ifndef ULong
typedef unsigned Long ULong;
#endif
#ifndef UShort
typedef unsigned short UShort;
#endif
#ifndef ANSI
#ifdef KR_headers
#define ANSI(x) ()
#define Void /*nothing*/
#else
#define ANSI(x) x
#define Void void
#endif
#endif /* ANSI */
#ifndef CONST
#ifdef KR_headers
#define CONST /* blank */
#else
#define CONST const
#endif
#endif /* CONST */
enum { /* return values from strtodg */
STRTOG_Zero = 0x000,
STRTOG_Normal = 0x001,
STRTOG_Denormal = 0x002,
STRTOG_Infinite = 0x003,
STRTOG_NaN = 0x004,
STRTOG_NaNbits = 0x005,
STRTOG_NoNumber = 0x006,
STRTOG_NoMemory = 0x007,
STRTOG_Retmask = 0x00f,
/* The following may be or-ed into one of the above values. */
STRTOG_Inexlo = 0x010, /* returned result rounded toward zero */
STRTOG_Inexhi = 0x020, /* returned result rounded away from zero */
STRTOG_Inexact = 0x030,
STRTOG_Underflow= 0x040,
STRTOG_Overflow = 0x080,
STRTOG_Neg = 0x100 /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
};
typedef struct
FPI {
int nbits;
int emin;
int emax;
int rounding;
int sudden_underflow;
} FPI;
enum { /* FPI.rounding values: same as FLT_ROUNDS */
FPI_Round_zero = 0,
FPI_Round_near = 1,
FPI_Round_up = 2,
FPI_Round_down = 3
};
#ifdef __cplusplus
extern "C" {
#endif
extern char* __dtoa ANSI((double d, int mode, int ndigits, int *decpt,
int *sign, char **rve));
extern char* __gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp,
int mode, int ndigits, int *decpt, char **rve));
extern void __freedtoa ANSI((char*));
extern float strtof ANSI((CONST char *, char **));
extern double strtod ANSI((CONST char *, char **));
extern int __strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*));
extern char* __g_ddfmt ANSI((char*, double*, int, size_t));
extern char* __g_dfmt ANSI((char*, double*, int, size_t));
extern char* __g_ffmt ANSI((char*, float*, int, size_t));
extern char* __g_Qfmt ANSI((char*, void*, int, size_t));
extern char* __g_xfmt ANSI((char*, void*, int, size_t));
extern char* __g_xLfmt ANSI((char*, void*, int, size_t));
extern int __strtoId ANSI((CONST char*, char**, double*, double*));
extern int __strtoIdd ANSI((CONST char*, char**, double*, double*));
extern int __strtoIf ANSI((CONST char*, char**, float*, float*));
extern int __strtoIQ ANSI((CONST char*, char**, void*, void*));
extern int __strtoIx ANSI((CONST char*, char**, void*, void*));
extern int __strtoIxL ANSI((CONST char*, char**, void*, void*));
extern int __strtord ANSI((CONST char*, char**, int, double*));
extern int __strtordd ANSI((CONST char*, char**, int, double*));
extern int __strtorf ANSI((CONST char*, char**, int, float*));
extern int __strtorQ ANSI((CONST char*, char**, int, void*));
extern int __strtorx ANSI((CONST char*, char**, int, void*));
extern int __strtorxL ANSI((CONST char*, char**, int, void*));
#if 1
extern int __strtodI ANSI((CONST char*, char**, double*));
extern int __strtopd ANSI((CONST char*, char**, double*));
extern int __strtopdd ANSI((CONST char*, char**, double*));
extern int __strtopf ANSI((CONST char*, char**, float*));
extern int __strtopQ ANSI((CONST char*, char**, void*));
extern int __strtopx ANSI((CONST char*, char**, void*));
extern int __strtopxL ANSI((CONST char*, char**, void*));
#else
#define __strtopd(s,se,x) strtord(s,se,1,x)
#define __strtopdd(s,se,x) strtordd(s,se,1,x)
#define __strtopf(s,se,x) strtorf(s,se,1,x)
#define __strtopQ(s,se,x) strtorQ(s,se,1,x)
#define __strtopx(s,se,x) strtorx(s,se,1,x)
#define __strtopxL(s,se,x) strtorxL(s,se,1,x)
#endif
#ifdef __cplusplus
}
#endif
#endif /* GDTOA_H_INCLUDED */

18
libc/upstream-openbsd/lib/libc/gdtoa/gdtoa_fltrnds.h Normal file
View File

@ -0,0 +1,18 @@
FPI *fpi, fpi1;
int Rounding;
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
Rounding = Flt_Rounds;
#else /*}{*/
Rounding = 1;
switch(fegetround()) {
case FE_TOWARDZERO: Rounding = 0; break;
case FE_UPWARD: Rounding = 2; break;
case FE_DOWNWARD: Rounding = 3;
}
#endif /*}}*/
fpi = &fpi0;
if (Rounding != 1) {
fpi1 = fpi0;
fpi = &fpi1;
fpi1.rounding = Rounding;
}

665
libc/upstream-openbsd/lib/libc/gdtoa/gdtoaimp.h Normal file
View File

@ -0,0 +1,665 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998-2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* This is a variation on dtoa.c that converts arbitary binary
floating-point formats to and from decimal notation. It uses
double-precision arithmetic internally, so there are still
various #ifdefs that adapt the calculations to the native
double-precision arithmetic (any of IEEE, VAX D_floating,
or IBM mainframe arithmetic).
Please send bug reports to David M. Gay (dmg at acm dot org,
with " at " changed at "@" and " dot " changed to ".").
*/
/* On a machine with IEEE extended-precision registers, it is
* necessary to specify double-precision (53-bit) rounding precision
* before invoking strtod or dtoa. If the machine uses (the equivalent
* of) Intel 80x87 arithmetic, the call
* _control87(PC_53, MCW_PC);
* does this with many compilers. Whether this or another call is
* appropriate depends on the compiler; for this to work, it may be
* necessary to #include "float.h" or another system-dependent header
* file.
*/
/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
*
* This strtod returns a nearest machine number to the input decimal
* string (or sets errno to ERANGE). With IEEE arithmetic, ties are
* broken by the IEEE round-even rule. Otherwise ties are broken by
* biased rounding (add half and chop).
*
* Inspired loosely by William D. Clinger's paper "How to Read Floating
* Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 112-126].
*
* Modifications:
*
* 1. We only require IEEE, IBM, or VAX double-precision
* arithmetic (not IEEE double-extended).
* 2. We get by with floating-point arithmetic in a case that
* Clinger missed -- when we're computing d * 10^n
* for a small integer d and the integer n is not too
* much larger than 22 (the maximum integer k for which
* we can represent 10^k exactly), we may be able to
* compute (d*10^k) * 10^(e-k) with just one roundoff.
* 3. Rather than a bit-at-a-time adjustment of the binary
* result in the hard case, we use floating-point
* arithmetic to determine the adjustment to within
* one bit; only in really hard cases do we need to
* compute a second residual.
* 4. Because of 3., we don't need a large table of powers of 10
* for ten-to-e (just some small tables, e.g. of 10^k
* for 0 <= k <= 22).
*/
/*
* #define IEEE_8087 for IEEE-arithmetic machines where the least
* significant byte has the lowest address.
* #define IEEE_MC68k for IEEE-arithmetic machines where the most
* significant byte has the lowest address.
* #define Long int on machines with 32-bit ints and 64-bit longs.
* #define Sudden_Underflow for IEEE-format machines without gradual
* underflow (i.e., that flush to zero on underflow).
* #define IBM for IBM mainframe-style floating-point arithmetic.
* #define VAX for VAX-style floating-point arithmetic (D_floating).
* #define No_leftright to omit left-right logic in fast floating-point
* computation of dtoa and gdtoa. This will cause modes 4 and 5 to be
* treated the same as modes 2 and 3 for some inputs.
* #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
* #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
* that use extended-precision instructions to compute rounded
* products and quotients) with IBM.
* #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
* that rounds toward +Infinity.
* #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
* rounding when the underlying floating-point arithmetic uses
* unbiased rounding. This prevent using ordinary floating-point
* arithmetic when the result could be computed with one rounding error.
* #define Inaccurate_Divide for IEEE-format with correctly rounded
* products but inaccurate quotients, e.g., for Intel i860.
* #define NO_LONG_LONG on machines that do not have a "long long"
* integer type (of >= 64 bits). On such machines, you can
* #define Just_16 to store 16 bits per 32-bit Long when doing
* high-precision integer arithmetic. Whether this speeds things
* up or slows things down depends on the machine and the number
* being converted. If long long is available and the name is
* something other than "long long", #define Llong to be the name,
* and if "unsigned Llong" does not work as an unsigned version of
* Llong, #define #ULLong to be the corresponding unsigned type.
* #define KR_headers for old-style C function headers.
* #define Bad_float_h if your system lacks a float.h or if it does not
* define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
* FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
* #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
* if memory is available and otherwise does something you deem
* appropriate. If MALLOC is undefined, malloc will be invoked
* directly -- and assumed always to succeed. Similarly, if you
* want something other than the system's free() to be called to
* recycle memory acquired from MALLOC, #define FREE to be the
* name of the alternate routine. (FREE or free is only called in
* pathological cases, e.g., in a gdtoa call after a gdtoa return in
* mode 3 with thousands of digits requested.)
* #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
* memory allocations from a private pool of memory when possible.
* When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
* unless #defined to be a different length. This default length
* suffices to get rid of MALLOC calls except for unusual cases,
* such as decimal-to-binary conversion of a very long string of
* digits. When converting IEEE double precision values, the
* longest string gdtoa can return is about 751 bytes long. For
* conversions by strtod of strings of 800 digits and all gdtoa
* conversions of IEEE doubles in single-threaded executions with
* 8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with
* 4-byte pointers, PRIVATE_MEM >= 7112 appears adequate.
* #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
* #defined automatically on IEEE systems. On such systems,
* when INFNAN_CHECK is #defined, strtod checks
* for Infinity and NaN (case insensitively).
* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
* strtodg also accepts (case insensitively) strings of the form
* NaN(x), where x is a string of hexadecimal digits (optionally
* preceded by 0x or 0X) and spaces; if there is only one string
* of hexadecimal digits, it is taken for the fraction bits of the
* resulting NaN; if there are two or more strings of hexadecimal
* digits, each string is assigned to the next available sequence
* of 32-bit words of fractions bits (starting with the most
* significant), right-aligned in each sequence.
* Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
* is consumed even when ... has the wrong form (in which case the
* "(...)" is consumed but ignored).
* #define MULTIPLE_THREADS if the system offers preemptively scheduled
* multiple threads. In this case, you must provide (or suitably
* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
* by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
* in pow5mult, ensures lazy evaluation of only one copy of high
* powers of 5; omitting this lock would introduce a small
* probability of wasting memory, but would otherwise be harmless.)
* You must also invoke freedtoa(s) to free the value s returned by
* dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
* #define IMPRECISE_INEXACT if you do not care about the setting of
* the STRTOG_Inexact bits in the special case of doing IEEE double
* precision conversions (which could also be done by the strtod in
* dtoa.c).
* #define NO_HEX_FP to disable recognition of C9x's hexadecimal
* floating-point constants.
* #define -DNO_ERRNO to suppress setting errno (in strtod.c and
* strtodg.c).
* #define NO_STRING_H to use private versions of memcpy.
* On some K&R systems, it may also be necessary to
* #define DECLARE_SIZE_T in this case.
* #define USE_LOCALE to use the current locale's decimal_point value.
*/
#ifndef GDTOAIMP_H_INCLUDED
#define GDTOAIMP_H_INCLUDED
#include "gdtoa.h"
#include "gd_qnan.h"
#ifdef Honor_FLT_ROUNDS
#include <fenv.h>
#endif
#ifdef DEBUG
#include "stdio.h"
#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
#endif
#include "stdlib.h"
#include "string.h"
#ifdef KR_headers
#define Char char
#else
#define Char void
#endif
#ifdef MALLOC
extern Char *MALLOC ANSI((size_t));
#else
#define MALLOC malloc
#endif
#undef IEEE_Arith
#undef Avoid_Underflow
#ifdef IEEE_MC68k
#define IEEE_Arith
#endif
#ifdef IEEE_8087
#define IEEE_Arith
#endif
#include "errno.h"
#ifdef Bad_float_h
#ifdef IEEE_Arith
#define DBL_DIG 15
#define DBL_MAX_10_EXP 308
#define DBL_MAX_EXP 1024
#define FLT_RADIX 2
#define DBL_MAX 1.7976931348623157e+308
#endif
#ifdef IBM
#define DBL_DIG 16
#define DBL_MAX_10_EXP 75
#define DBL_MAX_EXP 63
#define FLT_RADIX 16
#define DBL_MAX 7.2370055773322621e+75
#endif
#ifdef VAX
#define DBL_DIG 16
#define DBL_MAX_10_EXP 38
#define DBL_MAX_EXP 127
#define FLT_RADIX 2
#define DBL_MAX 1.7014118346046923e+38
#define n_bigtens 2
#endif
#ifndef LONG_MAX
#define LONG_MAX 2147483647
#endif
#else /* ifndef Bad_float_h */
#include "float.h"
#endif /* Bad_float_h */
#ifdef IEEE_Arith
#define Scale_Bit 0x10
#define n_bigtens 5
#endif
#ifdef IBM
#define n_bigtens 3
#endif
#ifdef VAX
#define n_bigtens 2
#endif
#ifndef __MATH_H__
#include "math.h"
#endif
#ifdef __cplusplus
extern "C" {
#endif
#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
#endif
typedef union { double d; ULong L[2]; } U;
#ifdef IEEE_8087
#define word0(x) (x)->L[1]
#define word1(x) (x)->L[0]
#else
#define word0(x) (x)->L[0]
#define word1(x) (x)->L[1]
#endif
#define dval(x) (x)->d
/* The following definition of Storeinc is appropriate for MIPS processors.
* An alternative that might be better on some machines is
* #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
*/
#if defined(IEEE_8087) + defined(VAX)
#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
((unsigned short *)a)[0] = (unsigned short)c, a++)
#else
#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
((unsigned short *)a)[1] = (unsigned short)c, a++)
#endif
/* #define P DBL_MANT_DIG */
/* Ten_pmax = floor(P*log(2)/log(5)) */
/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
#ifdef IEEE_Arith
#define Exp_shift 20
#define Exp_shift1 20
#define Exp_msk1 0x100000
#define Exp_msk11 0x100000
#define Exp_mask 0x7ff00000
#define P 53
#define Bias 1023
#define Emin (-1022)
#define Exp_1 0x3ff00000
#define Exp_11 0x3ff00000
#define Ebits 11
#define Frac_mask 0xfffff
#define Frac_mask1 0xfffff
#define Ten_pmax 22
#define Bletch 0x10
#define Bndry_mask 0xfffff
#define Bndry_mask1 0xfffff
#define LSB 1
#define Sign_bit 0x80000000
#define Log2P 1
#define Tiny0 0
#define Tiny1 1
#define Quick_max 14
#define Int_max 14
#ifndef Flt_Rounds
#ifdef FLT_ROUNDS
#define Flt_Rounds FLT_ROUNDS
#else
#define Flt_Rounds 1
#endif
#endif /*Flt_Rounds*/
#else /* ifndef IEEE_Arith */
#undef Sudden_Underflow
#define Sudden_Underflow
#ifdef IBM
#undef Flt_Rounds
#define Flt_Rounds 0
#define Exp_shift 24
#define Exp_shift1 24
#define Exp_msk1 0x1000000
#define Exp_msk11 0x1000000
#define Exp_mask 0x7f000000
#define P 14
#define Bias 65
#define Exp_1 0x41000000
#define Exp_11 0x41000000
#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
#define Frac_mask 0xffffff
#define Frac_mask1 0xffffff
#define Bletch 4
#define Ten_pmax 22
#define Bndry_mask 0xefffff
#define Bndry_mask1 0xffffff
#define LSB 1
#define Sign_bit 0x80000000
#define Log2P 4
#define Tiny0 0x100000
#define Tiny1 0
#define Quick_max 14
#define Int_max 15
#else /* VAX */
#undef Flt_Rounds
#define Flt_Rounds 1
#define Exp_shift 23
#define Exp_shift1 7
#define Exp_msk1 0x80
#define Exp_msk11 0x800000
#define Exp_mask 0x7f80
#define P 56
#define Bias 129
#define Emin (-127)
#define Exp_1 0x40800000
#define Exp_11 0x4080
#define Ebits 8
#define Frac_mask 0x7fffff
#define Frac_mask1 0xffff007f
#define Ten_pmax 24
#define Bletch 2
#define Bndry_mask 0xffff007f
#define Bndry_mask1 0xffff007f
#define LSB 0x10000
#define Sign_bit 0x8000
#define Log2P 1
#define Tiny0 0x80
#define Tiny1 0
#define Quick_max 15
#define Int_max 15
#endif /* IBM, VAX */
#endif /* IEEE_Arith */
#ifndef IEEE_Arith
#define ROUND_BIASED
#else
#ifdef ROUND_BIASED_without_Round_Up
#undef ROUND_BIASED
#define ROUND_BIASED
#endif
#endif
#ifdef RND_PRODQUOT
#define rounded_product(a,b) a = rnd_prod(a, b)
#define rounded_quotient(a,b) a = rnd_quot(a, b)
#ifdef KR_headers
extern double rnd_prod(), rnd_quot();
#else
extern double rnd_prod(double, double), rnd_quot(double, double);
#endif
#else
#define rounded_product(a,b) a *= b
#define rounded_quotient(a,b) a /= b
#endif
#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
#define Big1 0xffffffff
#undef Pack_16
#ifndef Pack_32
#define Pack_32
#endif
#ifdef NO_LONG_LONG
#undef ULLong
#ifdef Just_16
#undef Pack_32
#define Pack_16
/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
* This makes some inner loops simpler and sometimes saves work
* during multiplications, but it often seems to make things slightly
* slower. Hence the default is now to store 32 bits per Long.
*/
#endif
#else /* long long available */
#ifndef Llong
#define Llong long long
#endif
#ifndef ULLong
#define ULLong unsigned Llong
#endif
#endif /* NO_LONG_LONG */
#ifdef Pack_32
#define ULbits 32
#define kshift 5
#define kmask 31
#define ALL_ON 0xffffffff
#else
#define ULbits 16
#define kshift 4
#define kmask 15
#define ALL_ON 0xffff
#endif
#ifndef MULTIPLE_THREADS
#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
#define FREE_DTOA_LOCK(n) /*nothing*/
#else
#include "thread_private.h"
extern void *__dtoa_locks[];
#define ACQUIRE_DTOA_LOCK(n) _MUTEX_LOCK(&__dtoa_locks[n])
#define FREE_DTOA_LOCK(n) _MUTEX_UNLOCK(&__dtoa_locks[n])
#endif
#define Kmax 9
struct
Bigint {
struct Bigint *next;
int k, maxwds, sign, wds;
ULong x[1];
};
typedef struct Bigint Bigint;
#ifdef NO_STRING_H
#ifdef DECLARE_SIZE_T
typedef unsigned int size_t;
#endif
extern void memcpy_D2A ANSI((void*, const void*, size_t));
#define Bcopy(x,y) memcpy_D2A(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
#else /* !NO_STRING_H */
#define Bcopy(x,y) memcpy(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
#endif /* NO_STRING_H */
#define dtoa __dtoa
#define gdtoa __gdtoa
#define freedtoa __freedtoa
#define strtodg __strtodg
#define g_ddfmt __g_ddfmt
#define g_dfmt __g_dfmt
#define g_ffmt __g_ffmt
#define g_Qfmt __g_Qfmt
#define g_xfmt __g_xfmt
#define g_xLfmt __g_xLfmt
#define strtoId __strtoId
#define strtoIdd __strtoIdd
#define strtoIf __strtoIf
#define strtoIQ __strtoIQ
#define strtoIx __strtoIx
#define strtoIxL __strtoIxL
#define strtord __strtord
#define strtordd __strtordd
#define strtorf __strtorf
#define strtorQ __strtorQ
#define strtorx __strtorx
#define strtorxL __strtorxL
#define strtodI __strtodI
#define strtopd __strtopd
#define strtopdd __strtopdd
#define strtopf __strtopf
#define strtopQ __strtopQ
#define strtopx __strtopx
#define strtopxL __strtopxL
#define Balloc __Balloc_D2A
#define Bfree __Bfree_D2A
#define ULtoQ __ULtoQ_D2A
#define ULtof __ULtof_D2A
#define ULtod __ULtod_D2A
#define ULtodd __ULtodd_D2A
#define ULtox __ULtox_D2A
#define ULtoxL __ULtoxL_D2A
#define any_on __any_on_D2A
#define b2d __b2d_D2A
#define bigtens __bigtens_D2A
#define cmp __cmp_D2A
#define copybits __copybits_D2A
#define d2b __d2b_D2A
#define decrement __decrement_D2A
#define diff __diff_D2A
#define dtoa_result __dtoa_result_D2A
#define g__fmt __g__fmt_D2A
#define gethex __gethex_D2A
#define hexdig __hexdig_D2A
#define hexnan __hexnan_D2A
#define hi0bits(x) __hi0bits_D2A((ULong)(x))
#define hi0bits_D2A __hi0bits_D2A
#define i2b __i2b_D2A
#define increment __increment_D2A
#define lo0bits __lo0bits_D2A
#define lshift __lshift_D2A
#define match __match_D2A
#define mult __mult_D2A
#define multadd __multadd_D2A
#define nrv_alloc __nrv_alloc_D2A
#define pow5mult __pow5mult_D2A
#define quorem __quorem_D2A
#define ratio __ratio_D2A
#define rshift __rshift_D2A
#define rv_alloc __rv_alloc_D2A
#define s2b __s2b_D2A
#define set_ones __set_ones_D2A
#define strcp __strcp_D2A
#define strtoIg __strtoIg_D2A
#define sulp __sulp_D2A
#define sum __sum_D2A
#define tens __tens_D2A
#define tinytens __tinytens_D2A
#define tinytens __tinytens_D2A
#define trailz __trailz_D2A
#define ulp __ulp_D2A
extern char *dtoa_result;
extern CONST double bigtens[], tens[], tinytens[];
extern unsigned char hexdig[];
extern Bigint *Balloc ANSI((int));
extern void Bfree ANSI((Bigint*));
extern void ULtof ANSI((ULong*, ULong*, Long, int));
extern void ULtod ANSI((ULong*, ULong*, Long, int));
extern void ULtodd ANSI((ULong*, ULong*, Long, int));
extern void ULtoQ ANSI((ULong*, ULong*, Long, int));
extern void ULtox ANSI((UShort*, ULong*, Long, int));
extern void ULtoxL ANSI((ULong*, ULong*, Long, int));
extern ULong any_on ANSI((Bigint*, int));
extern double b2d ANSI((Bigint*, int*));
extern int cmp ANSI((Bigint*, Bigint*));
extern void copybits ANSI((ULong*, int, Bigint*));
extern Bigint *d2b ANSI((double, int*, int*));
extern void decrement ANSI((Bigint*));
extern Bigint *diff ANSI((Bigint*, Bigint*));
extern char *dtoa ANSI((double d, int mode, int ndigits,
int *decpt, int *sign, char **rve));
extern char *g__fmt ANSI((char*, char*, char*, int, ULong, size_t));
extern int gethex ANSI((CONST char**, FPI*, Long*, Bigint**, int));
extern void hexdig_init_D2A(Void);
extern int hexnan ANSI((CONST char**, FPI*, ULong*));
extern int hi0bits_D2A ANSI((ULong));
extern Bigint *i2b ANSI((int));
extern Bigint *increment ANSI((Bigint*));
extern int lo0bits ANSI((ULong*));
extern Bigint *lshift ANSI((Bigint*, int));
extern int match ANSI((CONST char**, char*));
extern Bigint *mult ANSI((Bigint*, Bigint*));
extern Bigint *multadd ANSI((Bigint*, int, int));
extern char *nrv_alloc ANSI((char*, char **, int));
extern Bigint *pow5mult ANSI((Bigint*, int));
extern int quorem ANSI((Bigint*, Bigint*));
extern double ratio ANSI((Bigint*, Bigint*));
extern void rshift ANSI((Bigint*, int));
extern char *rv_alloc ANSI((int));
extern Bigint *s2b ANSI((CONST char*, int, int, ULong, int));
extern Bigint *set_ones ANSI((Bigint*, int));
extern char *strcp ANSI((char*, const char*));
extern int strtoIg ANSI((CONST char*, char**, FPI*, Long*, Bigint**, int*));
extern double strtod ANSI((const char *s00, char **se));
extern Bigint *sum ANSI((Bigint*, Bigint*));
extern int trailz ANSI((Bigint*));
extern double ulp ANSI((U*));
#ifdef __cplusplus
}
#endif
/*
* NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c. Prior to
* 20050115, they used to be hard-wired here (to 0x7ff80000 and 0,
* respectively), but now are determined by compiling and running
* qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1.
* Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=...
* and -DNAN_WORD1=... values if necessary. This should still work.
* (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
*/
#ifdef IEEE_Arith
#ifndef NO_INFNAN_CHECK
#undef INFNAN_CHECK
#define INFNAN_CHECK
#endif
#ifdef IEEE_MC68k
#define _0 0
#define _1 1
#ifndef NAN_WORD0
#define NAN_WORD0 d_QNAN0
#endif
#ifndef NAN_WORD1
#define NAN_WORD1 d_QNAN1
#endif
#else
#define _0 1
#define _1 0
#ifndef NAN_WORD0
#define NAN_WORD0 d_QNAN1
#endif
#ifndef NAN_WORD1
#define NAN_WORD1 d_QNAN0
#endif
#endif
#else
#undef INFNAN_CHECK
#endif
#undef SI
#ifdef Sudden_Underflow
#define SI 1
#else
#define SI 0
#endif
#endif /* GDTOAIMP_H_INCLUDED */

360
libc/upstream-openbsd/lib/libc/gdtoa/gethex.c Normal file
View File

@ -0,0 +1,360 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
#ifdef USE_LOCALE
#include "locale.h"
#endif
int
#ifdef KR_headers
gethex(sp, fpi, exp, bp, sign)
CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign;
#else
gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
#endif
{
Bigint *b;
CONST unsigned char *decpt, *s0, *s, *s1;
int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
ULong L, lostbits, *x;
Long e, e1;
#ifdef USE_LOCALE
int i;
#ifdef NO_LOCALE_CACHE
const unsigned char *decimalpoint = (unsigned char*)localeconv()->decimal_point;
#else
const unsigned char *decimalpoint;
static unsigned char *decimalpoint_cache;
if (!(s0 = decimalpoint_cache)) {
s0 = (unsigned char*)localeconv()->decimal_point;
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strlcpy(decimalpoint_cache, s0, strlen(s0) + 1);
s0 = decimalpoint_cache;
}
}
decimalpoint = s0;
#endif
#endif
if (!hexdig['0'])
hexdig_init_D2A();
*bp = 0;
havedig = 0;
s0 = *(CONST unsigned char **)sp + 2;
while(s0[havedig] == '0')
havedig++;
s0 += havedig;
s = s0;
decpt = 0;
zret = 0;
e = 0;
if (hexdig[*s])
havedig++;
else {
zret = 1;
#ifdef USE_LOCALE
for(i = 0; decimalpoint[i]; ++i) {
if (s[i] != decimalpoint[i])
goto pcheck;
}
decpt = s += i;
#else
if (*s != '.')
goto pcheck;
decpt = ++s;
#endif
if (!hexdig[*s])
goto pcheck;
while(*s == '0')
s++;
if (hexdig[*s])
zret = 0;
havedig = 1;
s0 = s;
}
while(hexdig[*s])
s++;
#ifdef USE_LOCALE
if (*s == *decimalpoint && !decpt) {
for(i = 1; decimalpoint[i]; ++i) {
if (s[i] != decimalpoint[i])
goto pcheck;
}
decpt = s += i;
#else
if (*s == '.' && !decpt) {
decpt = ++s;
#endif
while(hexdig[*s])
s++;
}/*}*/
if (decpt)
e = -(((Long)(s-decpt)) << 2);
pcheck:
s1 = s;
big = esign = 0;
switch(*s) {
case 'p':
case 'P':
switch(*++s) {
case '-':
esign = 1;
/* no break */
case '+':
s++;
}
if ((n = hexdig[*s]) == 0 || n > 0x19) {
s = s1;
break;
}
e1 = n - 0x10;
while((n = hexdig[*++s]) !=0 && n <= 0x19) {
if (e1 & 0xf8000000)
big = 1;
e1 = 10*e1 + n - 0x10;
}
if (esign)
e1 = -e1;
e += e1;
}
*sp = (char*)s;
if (!havedig)
*sp = (char*)s0 - 1;
if (zret)
return STRTOG_Zero;
if (big) {
if (esign) {
switch(fpi->rounding) {
case FPI_Round_up:
if (sign)
break;
goto ret_tiny;
case FPI_Round_down:
if (!sign)
break;
goto ret_tiny;
}
goto retz;
ret_tiny:
b = Balloc(0);
if (b == NULL)
return (STRTOG_NoMemory);
b->wds = 1;
b->x[0] = 1;
goto dret;
}
switch(fpi->rounding) {
case FPI_Round_near:
goto ovfl1;
case FPI_Round_up:
if (!sign)
goto ovfl1;
goto ret_big;
case FPI_Round_down:
if (sign)
goto ovfl1;
goto ret_big;
}
ret_big:
nbits = fpi->nbits;
n0 = n = nbits >> kshift;
if (nbits & kmask)
++n;
for(j = n, k = 0; j >>= 1; ++k);
*bp = b = Balloc(k);
if (*bp == NULL)
return (STRTOG_NoMemory);
b->wds = n;
for(j = 0; j < n0; ++j)
b->x[j] = ALL_ON;
if (n > n0)
b->x[j] = ULbits >> (ULbits - (nbits & kmask));
*exp = fpi->emin;
return STRTOG_Normal | STRTOG_Inexlo;
}
n = s1 - s0 - 1;
for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
k++;
b = Balloc(k);
if (b == NULL)
return (STRTOG_NoMemory);
x = b->x;
n = 0;
L = 0;
#ifdef USE_LOCALE
for(i = 0; decimalpoint[i+1]; ++i);
#endif
while(s1 > s0) {
#ifdef USE_LOCALE
if (*--s1 == decimalpoint[i]) {
s1 -= i;
continue;
}
#else
if (*--s1 == '.')
continue;
#endif
if (n == ULbits) {
*x++ = L;
L = 0;
n = 0;
}
L |= (hexdig[*s1] & 0x0f) << n;
n += 4;
}
*x++ = L;
b->wds = n = x - b->x;
n = ULbits*n - hi0bits(L);
nbits = fpi->nbits;
lostbits = 0;
x = b->x;
if (n > nbits) {
n -= nbits;
if (any_on(b,n)) {
lostbits = 1;
k = n - 1;
if (x[k>>kshift] & 1 << (k & kmask)) {
lostbits = 2;
if (k > 0 && any_on(b,k))
lostbits = 3;
}
}
rshift(b, n);
e += n;
}
else if (n < nbits) {
n = nbits - n;
b = lshift(b, n);
if (b == NULL)
return (STRTOG_NoMemory);
e -= n;
x = b->x;
}
if (e > fpi->emax) {
ovfl:
Bfree(b);
ovfl1:
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
}
irv = STRTOG_Normal;
if (e < fpi->emin) {
irv = STRTOG_Denormal;
n = fpi->emin - e;
if (n >= nbits) {
switch (fpi->rounding) {
case FPI_Round_near:
if (n == nbits && (n < 2 || any_on(b,n-1)))
goto one_bit;
break;
case FPI_Round_up:
if (!sign)
goto one_bit;
break;
case FPI_Round_down:
if (sign) {
one_bit:
x[0] = b->wds = 1;
dret:
*bp = b;
*exp = fpi->emin;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Denormal | STRTOG_Inexhi
| STRTOG_Underflow;
}
}
Bfree(b);
retz:
#ifndef NO_ERRNO
errno = ERANGE;
#endif
return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
}
k = n - 1;
if (lostbits)
lostbits = 1;
else if (k > 0)
lostbits = any_on(b,k);
if (x[k>>kshift] & 1 << (k & kmask))
lostbits |= 2;
nbits -= n;
rshift(b,n);
e = fpi->emin;
}
if (lostbits) {
up = 0;
switch(fpi->rounding) {
case FPI_Round_zero:
break;
case FPI_Round_near:
if (lostbits & 2
&& (lostbits | x[0]) & 1)
up = 1;
break;
case FPI_Round_up:
up = 1 - sign;
break;
case FPI_Round_down:
up = sign;
}
if (up) {
k = b->wds;
b = increment(b);
if (b == NULL)
return (STRTOG_NoMemory);
x = b->x;
if (irv == STRTOG_Denormal) {
if (nbits == fpi->nbits - 1
&& x[nbits >> kshift] & 1 << (nbits & kmask))
irv = STRTOG_Normal;
}
else if (b->wds > k
|| ((n = nbits & kmask) !=0
&& hi0bits(x[k-1]) < 32-n)) {
rshift(b,1);
if (++e > fpi->emax)
goto ovfl;
}
irv |= STRTOG_Inexhi;
}
else
irv |= STRTOG_Inexlo;
}
*bp = b;
*exp = e;
return irv;
}

86
libc/upstream-openbsd/lib/libc/gdtoa/gmisc.c Normal file
View File

@ -0,0 +1,86 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
void
#ifdef KR_headers
rshift(b, k) Bigint *b; int k;
#else
rshift(Bigint *b, int k)
#endif
{
ULong *x, *x1, *xe, y;
int n;
x = x1 = b->x;
n = k >> kshift;
if (n < b->wds) {
xe = x + b->wds;
x += n;
if (k &= kmask) {
n = ULbits - k;
y = *x++ >> k;
while(x < xe) {
*x1++ = (y | (*x << n)) & ALL_ON;
y = *x++ >> k;
}
if ((*x1 = y) !=0)
x1++;
}
else
while(x < xe)
*x1++ = *x++;
}
if ((b->wds = x1 - b->x) == 0)
b->x[0] = 0;
}
int
#ifdef KR_headers
trailz(b) Bigint *b;
#else
trailz(Bigint *b)
#endif
{
ULong L, *x, *xe;
int n = 0;
x = b->x;
xe = x + b->wds;
for(n = 0; x < xe && !*x; x++)
n += ULbits;
if (x < xe) {
L = *x;
n += lo0bits(&L);
}
return n;
}

55
libc/upstream-openbsd/lib/libc/gdtoa/hd_init.c Normal file
View File

@ -0,0 +1,55 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
unsigned char hexdig[256];
static void
#ifdef KR_headers
htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
#else
htinit(unsigned char *h, unsigned char *s, int inc)
#endif
{
int i, j;
for(i = 0; (j = s[i]) !=0; i++)
h[j] = i + inc;
}
void
hexdig_init_D2A(Void)
{
#define USC (unsigned char *)
htinit(hexdig, USC "0123456789", 0x10);
htinit(hexdig, USC "abcdef", 0x10 + 10);
htinit(hexdig, USC "ABCDEF", 0x10 + 10);
}

332
libc/upstream-openbsd/lib/libc/gdtoa/hdtoa.c Normal file
View File

@ -0,0 +1,332 @@
/* $OpenBSD: hdtoa.c,v 1.2 2009/10/16 12:15:03 martynas Exp $ */
/*-
* Copyright (c) 2004, 2005 David Schultz <das@FreeBSD.ORG>
* 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.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
*/
#include <sys/types.h>
#include <machine/ieee.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include "gdtoaimp.h"
/* Strings values used by dtoa() */
#define INFSTR "Infinity"
#define NANSTR "NaN"
#define DBL_ADJ (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
#define LDBL_ADJ (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
/*
* Round up the given digit string. If the digit string is fff...f,
* this procedure sets it to 100...0 and returns 1 to indicate that
* the exponent needs to be bumped. Otherwise, 0 is returned.
*/
static int
roundup(char *s0, int ndigits)
{
char *s;
for (s = s0 + ndigits - 1; *s == 0xf; s--) {
if (s == s0) {
*s = 1;
return (1);
}
*s = 0;
}
++*s;
return (0);
}
/*
* Round the given digit string to ndigits digits according to the
* current rounding mode. Note that this could produce a string whose
* value is not representable in the corresponding floating-point
* type. The exponent pointed to by decpt is adjusted if necessary.
*/
static void
dorounding(char *s0, int ndigits, int sign, int *decpt)
{
int adjust = 0; /* do we need to adjust the exponent? */
switch (FLT_ROUNDS) {
case 0: /* toward zero */
default: /* implementation-defined */
break;
case 1: /* to nearest, halfway rounds to even */
if ((s0[ndigits] > 8) ||
(s0[ndigits] == 8 && s0[ndigits + 1] & 1))
adjust = roundup(s0, ndigits);
break;
case 2: /* toward +inf */
if (sign == 0)
adjust = roundup(s0, ndigits);
break;
case 3: /* toward -inf */
if (sign != 0)
adjust = roundup(s0, ndigits);
break;
}
if (adjust)
*decpt += 4;
}
/*
* This procedure converts a double-precision number in IEEE format
* into a string of hexadecimal digits and an exponent of 2. Its
* behavior is bug-for-bug compatible with dtoa() in mode 2, with the
* following exceptions:
*
* - An ndigits < 0 causes it to use as many digits as necessary to
* represent the number exactly.
* - The additional xdigs argument should point to either the string
* "0123456789ABCDEF" or the string "0123456789abcdef", depending on
* which case is desired.
* - This routine does not repeat dtoa's mistake of setting decpt
* to 9999 in the case of an infinity or NaN. INT_MAX is used
* for this purpose instead.
*
* Note that the C99 standard does not specify what the leading digit
* should be for non-zero numbers. For instance, 0x1.3p3 is the same
* as 0x2.6p2 is the same as 0x4.cp3. This implementation chooses the
* first digit so that subsequent digits are aligned on nibble
* boundaries (before rounding).
*
* Inputs: d, xdigs, ndigits
* Outputs: decpt, sign, rve
*/
char *
__hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
char **rve)
{
static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
struct ieee_double *p = (struct ieee_double *)&d;
char *s, *s0;
int bufsize;
*sign = p->dbl_sign;
switch (fpclassify(d)) {
case FP_NORMAL:
*decpt = p->dbl_exp - DBL_ADJ;
break;
case FP_ZERO:
*decpt = 1;
return (nrv_alloc("0", rve, 1));
case FP_SUBNORMAL:
d *= 0x1p514;
*decpt = p->dbl_exp - (514 + DBL_ADJ);
break;
case FP_INFINITE:
*decpt = INT_MAX;
return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
case FP_NAN:
*decpt = INT_MAX;
return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
default:
abort();
}
/* FP_NORMAL or FP_SUBNORMAL */
if (ndigits == 0) /* dtoa() compatibility */
ndigits = 1;
/*
* For simplicity, we generate all the digits even if the
* caller has requested fewer.
*/
bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
s0 = rv_alloc(bufsize);
if (s0 == NULL)
return (NULL);
/*
* We work from right to left, first adding any requested zero
* padding, then the least significant portion of the
* mantissa, followed by the most significant. The buffer is
* filled with the byte values 0x0 through 0xf, which are
* converted to xdigs[0x0] through xdigs[0xf] after the
* rounding phase.
*/
for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
*s = 0;
for (; s > s0 + sigfigs - (DBL_FRACLBITS / 4) - 1 && s > s0; s--) {
*s = p->dbl_fracl & 0xf;
p->dbl_fracl >>= 4;
}
for (; s > s0; s--) {
*s = p->dbl_frach & 0xf;
p->dbl_frach >>= 4;
}
/*
* At this point, we have snarfed all the bits in the
* mantissa, with the possible exception of the highest-order
* (partial) nibble, which is dealt with by the next
* statement. We also tack on the implicit normalization bit.
*/
*s = p->dbl_frach | (1U << ((DBL_MANT_DIG - 1) % 4));
/* If ndigits < 0, we are expected to auto-size the precision. */
if (ndigits < 0) {
for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
;
}
if (sigfigs > ndigits && s0[ndigits] != 0)
dorounding(s0, ndigits, p->dbl_sign, decpt);
s = s0 + ndigits;
if (rve != NULL)
*rve = s;
*s-- = '\0';
for (; s >= s0; s--)
*s = xdigs[(unsigned int)*s];
return (s0);
}
#if (LDBL_MANT_DIG > DBL_MANT_DIG)
/*
* This is the long double version of __hdtoa().
*/
char *
__hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
char **rve)
{
static const int sigfigs = (LDBL_MANT_DIG + 3) / 4;
struct ieee_ext *p = (struct ieee_ext *)&e;
char *s, *s0;
int bufsize;
*sign = p->ext_sign;
switch (fpclassify(e)) {
case FP_NORMAL:
*decpt = p->ext_exp - LDBL_ADJ;
break;
case FP_ZERO:
*decpt = 1;
return (nrv_alloc("0", rve, 1));
case FP_SUBNORMAL:
e *= 0x1p514L;
*decpt = p->ext_exp - (514 + LDBL_ADJ);
break;
case FP_INFINITE:
*decpt = INT_MAX;
return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
case FP_NAN:
*decpt = INT_MAX;
return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
default:
abort();
}
/* FP_NORMAL or FP_SUBNORMAL */
if (ndigits == 0) /* dtoa() compatibility */
ndigits = 1;
/*
* For simplicity, we generate all the digits even if the
* caller has requested fewer.
*/
bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
s0 = rv_alloc(bufsize);
if (s0 == NULL)
return (NULL);
/*
* We work from right to left, first adding any requested zero
* padding, then the least significant portion of the
* mantissa, followed by the most significant. The buffer is
* filled with the byte values 0x0 through 0xf, which are
* converted to xdigs[0x0] through xdigs[0xf] after the
* rounding phase.
*/
for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
*s = 0;
for (; s > s0 + sigfigs - (EXT_FRACLBITS / 4) - 1 && s > s0; s--) {
*s = p->ext_fracl & 0xf;
p->ext_fracl >>= 4;
}
#ifdef EXT_FRACHMBITS
for (; s > s0; s--) {
*s = p->ext_frachm & 0xf;
p->ext_frachm >>= 4;
}
#endif
#ifdef EXT_FRACLMBITS
for (; s > s0; s--) {
*s = p->ext_fraclm & 0xf;
p->ext_fraclm >>= 4;
}
#endif
for (; s > s0; s--) {
*s = p->ext_frach & 0xf;
p->ext_frach >>= 4;
}
/*
* At this point, we have snarfed all the bits in the
* mantissa, with the possible exception of the highest-order
* (partial) nibble, which is dealt with by the next
* statement. We also tack on the implicit normalization bit.
*/
*s = p->ext_frach | (1U << ((LDBL_MANT_DIG - 1) % 4));
/* If ndigits < 0, we are expected to auto-size the precision. */
if (ndigits < 0) {
for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
;
}
if (sigfigs > ndigits && s0[ndigits] != 0)
dorounding(s0, ndigits, p->ext_sign, decpt);
s = s0 + ndigits;
if (rve != NULL)
*rve = s;
*s-- = '\0';
for (; s >= s0; s--)
*s = xdigs[(unsigned int)*s];
return (s0);
}
#else /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
char *
__hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
char **rve)
{
return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
}
#endif /* (LDBL_MANT_DIG == DBL_MANT_DIG) */

150
libc/upstream-openbsd/lib/libc/gdtoa/hexnan.c Normal file
View File

@ -0,0 +1,150 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
static void
#ifdef KR_headers
L_shift(x, x1, i) ULong *x; ULong *x1; int i;
#else
L_shift(ULong *x, ULong *x1, int i)
#endif
{
int j;
i = 8 - i;
i <<= 2;
j = ULbits - i;
do {
*x |= x[1] << j;
x[1] >>= i;
} while(++x < x1);
}
int
#ifdef KR_headers
hexnan(sp, fpi, x0)
CONST char **sp; FPI *fpi; ULong *x0;
#else
hexnan( CONST char **sp, FPI *fpi, ULong *x0)
#endif
{
ULong c, h, *x, *x1, *xe;
CONST char *s;
int havedig, hd0, i, nbits;
if (!hexdig['0'])
hexdig_init_D2A();
nbits = fpi->nbits;
x = x0 + (nbits >> kshift);
if (nbits & kmask)
x++;
*--x = 0;
x1 = xe = x;
havedig = hd0 = i = 0;
s = *sp;
/* allow optional initial 0x or 0X */
while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(CONST unsigned char*)(s+3) > ' ')
s += 2;
while((c = *(CONST unsigned char*)++s)) {
if (!(h = hexdig[c])) {
if (c <= ' ') {
if (hd0 < havedig) {
if (x < x1 && i < 8)
L_shift(x, x1, i);
if (x <= x0) {
i = 8;
continue;
}
hd0 = havedig;
*--x = 0;
x1 = x;
i = 0;
}
while(*(CONST unsigned char*)(s+1) <= ' ')
++s;
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(CONST unsigned char*)(s+3) > ' ')
s += 2;
continue;
}
if (/*(*/ c == ')' && havedig) {
*sp = s + 1;
break;
}
#ifndef GDTOA_NON_PEDANTIC_NANCHECK
do {
if (/*(*/ c == ')') {
*sp = s + 1;
break;
}
} while((c = *++s));
#endif
return STRTOG_NaN;
}
havedig++;
if (++i > 8) {
if (x <= x0)
continue;
i = 1;
*--x = 0;
}
*x = (*x << 4) | (h & 0xf);
}
if (!havedig)
return STRTOG_NaN;
if (x < x1 && i < 8)
L_shift(x, x1, i);
if (x > x0) {
x1 = x0;
do *x1++ = *x++;
while(x <= xe);
do *x1++ = 0;
while(x1 <= xe);
}
else {
/* truncate high-order word if necessary */
if ( (i = nbits & (ULbits-1)) !=0)
*xe &= ((ULong)0xffffffff) >> (ULbits - i);
}
for(x1 = xe;; --x1) {
if (*x1 != 0)
break;
if (x1 == x0) {
*x1 = 1;
break;
}
}
return STRTOG_NaNbits;
}

124
libc/upstream-openbsd/lib/libc/gdtoa/ldtoa.c Normal file
View File

@ -0,0 +1,124 @@
/* $OpenBSD: ldtoa.c,v 1.1 2008/09/07 20:36:08 martynas Exp $ */
/*-
* Copyright (c) 2003 David Schultz <das@FreeBSD.ORG>
* 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.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
*/
#include <sys/types.h>
#ifndef __vax__
#include <machine/ieee.h>
#endif /* !__vax__ */
#include <float.h>
#include <inttypes.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
#include "gdtoaimp.h"
#if (LDBL_MANT_DIG > DBL_MANT_DIG)
/*
* ldtoa() is a wrapper for gdtoa() that makes it smell like dtoa(),
* except that the floating point argument is passed by reference.
* When dtoa() is passed a NaN or infinity, it sets expt to 9999.
* However, a long double could have a valid exponent of 9999, so we
* use INT_MAX in ldtoa() instead.
*/
char *
__ldtoa(long double *ld, int mode, int ndigits, int *decpt, int *sign,
char **rve)
{
FPI fpi = {
LDBL_MANT_DIG, /* nbits */
LDBL_MIN_EXP - LDBL_MANT_DIG, /* emin */
LDBL_MAX_EXP - LDBL_MANT_DIG, /* emax */
FLT_ROUNDS, /* rounding */
#ifdef Sudden_Underflow /* unused, but correct anyway */
1
#else
0
#endif
};
int be, kind;
char *ret;
struct ieee_ext *p = (struct ieee_ext *)ld;
uint32_t bits[(LDBL_MANT_DIG + 31) / 32];
void *vbits = bits;
/*
* gdtoa doesn't know anything about the sign of the number, so
* if the number is negative, we need to swap rounding modes of
* 2 (upwards) and 3 (downwards).
*/
*sign = p->ext_sign;
fpi.rounding ^= (fpi.rounding >> 1) & p->ext_sign;
be = p->ext_exp - (LDBL_MAX_EXP - 1) - (LDBL_MANT_DIG - 1);
EXT_TO_ARRAY32(p, bits);
switch (fpclassify(*ld)) {
case FP_NORMAL:
kind = STRTOG_Normal;
#ifdef EXT_IMPLICIT_NBIT
bits[LDBL_MANT_DIG / 32] |= 1 << ((LDBL_MANT_DIG - 1) % 32);
#endif /* EXT_IMPLICIT_NBIT */
break;
case FP_ZERO:
kind = STRTOG_Zero;
break;
case FP_SUBNORMAL:
kind = STRTOG_Denormal;
be++;
break;
case FP_INFINITE:
kind = STRTOG_Infinite;
break;
case FP_NAN:
kind = STRTOG_NaN;
break;
default:
abort();
}
ret = gdtoa(&fpi, be, vbits, &kind, mode, ndigits, decpt, rve);
if (*decpt == -32768)
*decpt = INT_MAX;
return ret;
}
#else /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
char *
__ldtoa(long double *ld, int mode, int ndigits, int *decpt, int *sign,
char **rve)
{
char *ret;
ret = dtoa((double)*ld, mode, ndigits, decpt, sign, rve);
if (*decpt == 9999)
*decpt = INT_MAX;
return ret;
}
#endif /* (LDBL_MANT_DIG == DBL_MANT_DIG) */

907
libc/upstream-openbsd/lib/libc/gdtoa/misc.c Normal file
View File

@ -0,0 +1,907 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
static Bigint *freelist[Kmax+1];
#ifndef Omit_Private_Memory
#ifndef PRIVATE_MEM
#define PRIVATE_MEM 2304
#endif
#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
#endif
Bigint *
Balloc
#ifdef KR_headers
(k) int k;
#else
(int k)
#endif
{
int x;
Bigint *rv;
#ifndef Omit_Private_Memory
unsigned int len;
#endif
ACQUIRE_DTOA_LOCK(0);
/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
/* but this case seems very unlikely. */
if (k <= Kmax && (rv = freelist[k]) !=0) {
freelist[k] = rv->next;
}
else {
x = 1 << k;
#ifdef Omit_Private_Memory
rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
if (rv == NULL)
return (NULL);
#else
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
/sizeof(double);
if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
rv = (Bigint*)pmem_next;
pmem_next += len;
}
else {
rv = (Bigint*)MALLOC(len*sizeof(double));
if (rv == NULL)
return (NULL);
}
#endif
rv->k = k;
rv->maxwds = x;
}
FREE_DTOA_LOCK(0);
rv->sign = rv->wds = 0;
return rv;
}
void
Bfree
#ifdef KR_headers
(v) Bigint *v;
#else
(Bigint *v)
#endif
{
if (v) {
if (v->k > Kmax)
#ifdef FREE
FREE((void*)v);
#else
free((void*)v);
#endif
else {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
}
}
}
int
lo0bits
#ifdef KR_headers
(y) ULong *y;
#else
(ULong *y)
#endif
{
int k;
ULong x = *y;
if (x & 7) {
if (x & 1)
return 0;
if (x & 2) {
*y = x >> 1;
return 1;
}
*y = x >> 2;
return 2;
}
k = 0;
if (!(x & 0xffff)) {
k = 16;
x >>= 16;
}
if (!(x & 0xff)) {
k += 8;
x >>= 8;
}
if (!(x & 0xf)) {
k += 4;
x >>= 4;
}
if (!(x & 0x3)) {
k += 2;
x >>= 2;
}
if (!(x & 1)) {
k++;
x >>= 1;
if (!x)
return 32;
}
*y = x;
return k;
}
Bigint *
multadd
#ifdef KR_headers
(b, m, a) Bigint *b; int m, a;
#else
(Bigint *b, int m, int a) /* multiply by m and add a */
#endif
{
int i, wds;
#ifdef ULLong
ULong *x;
ULLong carry, y;
#else
ULong carry, *x, y;
#ifdef Pack_32
ULong xi, z;
#endif
#endif
Bigint *b1;
wds = b->wds;
x = b->x;
i = 0;
carry = a;
do {
#ifdef ULLong
y = *x * (ULLong)m + carry;
carry = y >> 32;
*x++ = y & 0xffffffffUL;
#else
#ifdef Pack_32
xi = *x;
y = (xi & 0xffff) * m + carry;
z = (xi >> 16) * m + (y >> 16);
carry = z >> 16;
*x++ = (z << 16) + (y & 0xffff);
#else
y = *x * m + carry;
carry = y >> 16;
*x++ = y & 0xffff;
#endif
#endif
}
while(++i < wds);
if (carry) {
if (wds >= b->maxwds) {
b1 = Balloc(b->k+1);
if (b1 == NULL)
return (NULL);
Bcopy(b1, b);
Bfree(b);
b = b1;
}
b->x[wds++] = carry;
b->wds = wds;
}
return b;
}
int
hi0bits_D2A
#ifdef KR_headers
(x) ULong x;
#else
(ULong x)
#endif
{
int k = 0;
if (!(x & 0xffff0000)) {
k = 16;
x <<= 16;
}
if (!(x & 0xff000000)) {
k += 8;
x <<= 8;
}
if (!(x & 0xf0000000)) {
k += 4;
x <<= 4;
}
if (!(x & 0xc0000000)) {
k += 2;
x <<= 2;
}
if (!(x & 0x80000000)) {
k++;
if (!(x & 0x40000000))
return 32;
}
return k;
}
Bigint *
i2b
#ifdef KR_headers
(i) int i;
#else
(int i)
#endif
{
Bigint *b;
b = Balloc(1);
if (b == NULL)
return (NULL);
b->x[0] = i;
b->wds = 1;
return b;
}
Bigint *
mult
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
{
Bigint *c;
int k, wa, wb, wc;
ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
ULong y;
#ifdef ULLong
ULLong carry, z;
#else
ULong carry, z;
#ifdef Pack_32
ULong z2;
#endif
#endif
if (a->wds < b->wds) {
c = a;
a = b;
b = c;
}
k = a->k;
wa = a->wds;
wb = b->wds;
wc = wa + wb;
if (wc > a->maxwds)
k++;
c = Balloc(k);
if (c == NULL)
return (NULL);
for(x = c->x, xa = x + wc; x < xa; x++)
*x = 0;
xa = a->x;
xae = xa + wa;
xb = b->x;
xbe = xb + wb;
xc0 = c->x;
#ifdef ULLong
for(; xb < xbe; xc0++) {
if ( (y = *xb++) !=0) {
x = xa;
xc = xc0;
carry = 0;
do {
z = *x++ * (ULLong)y + *xc + carry;
carry = z >> 32;
*xc++ = z & 0xffffffffUL;
}
while(x < xae);
*xc = carry;
}
}
#else
#ifdef Pack_32
for(; xb < xbe; xb++, xc0++) {
if ( (y = *xb & 0xffff) !=0) {
x = xa;
xc = xc0;
carry = 0;
do {
z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
carry = z >> 16;
z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
carry = z2 >> 16;
Storeinc(xc, z2, z);
}
while(x < xae);
*xc = carry;
}
if ( (y = *xb >> 16) !=0) {
x = xa;
xc = xc0;
carry = 0;
z2 = *xc;
do {
z = (*x & 0xffff) * y + (*xc >> 16) + carry;
carry = z >> 16;
Storeinc(xc, z, z2);
z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
carry = z2 >> 16;
}
while(x < xae);
*xc = z2;
}
}
#else
for(; xb < xbe; xc0++) {
if ( (y = *xb++) !=0) {
x = xa;
xc = xc0;
carry = 0;
do {
z = *x++ * y + *xc + carry;
carry = z >> 16;
*xc++ = z & 0xffff;
}
while(x < xae);
*xc = carry;
}
}
#endif
#endif
for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
c->wds = wc;
return c;
}
static Bigint *p5s;
Bigint *
pow5mult
#ifdef KR_headers
(b, k) Bigint *b; int k;
#else
(Bigint *b, int k)
#endif
{
Bigint *b1, *p5, *p51;
int i;
static int p05[3] = { 5, 25, 125 };
if ( (i = k & 3) !=0) {
b = multadd(b, p05[i-1], 0);
if (b == NULL)
return (NULL);
}
if (!(k >>= 2))
return b;
if ((p5 = p5s) == 0) {
/* first time */
#ifdef MULTIPLE_THREADS
ACQUIRE_DTOA_LOCK(1);
if (!(p5 = p5s)) {
p5 = p5s = i2b(625);
if (p5 == NULL)
return (NULL);
p5->next = 0;
}
FREE_DTOA_LOCK(1);
#else
p5 = p5s = i2b(625);
if (p5 == NULL)
return (NULL);
p5->next = 0;
#endif
}
for(;;) {
if (k & 1) {
b1 = mult(b, p5);
if (b1 == NULL)
return (NULL);
Bfree(b);
b = b1;
}
if (!(k >>= 1))
break;
if ((p51 = p5->next) == 0) {
#ifdef MULTIPLE_THREADS
ACQUIRE_DTOA_LOCK(1);
if (!(p51 = p5->next)) {
p51 = p5->next = mult(p5,p5);
if (p51 == NULL)
return (NULL);
p51->next = 0;
}
FREE_DTOA_LOCK(1);
#else
p51 = p5->next = mult(p5,p5);
if (p51 == NULL)
return (NULL);
p51->next = 0;
#endif
}
p5 = p51;
}
return b;
}
Bigint *
lshift
#ifdef KR_headers
(b, k) Bigint *b; int k;
#else
(Bigint *b, int k)
#endif
{
int i, k1, n, n1;
Bigint *b1;
ULong *x, *x1, *xe, z;
n = k >> kshift;
k1 = b->k;
n1 = n + b->wds + 1;
for(i = b->maxwds; n1 > i; i <<= 1)
k1++;
b1 = Balloc(k1);
if (b1 == NULL)
return (NULL);
x1 = b1->x;
for(i = 0; i < n; i++)
*x1++ = 0;
x = b->x;
xe = x + b->wds;
if (k &= kmask) {
#ifdef Pack_32
k1 = 32 - k;
z = 0;
do {
*x1++ = *x << k | z;
z = *x++ >> k1;
}
while(x < xe);
if ((*x1 = z) !=0)
++n1;
#else
k1 = 16 - k;
z = 0;
do {
*x1++ = *x << k & 0xffff | z;
z = *x++ >> k1;
}
while(x < xe);
if (*x1 = z)
++n1;
#endif
}
else do
*x1++ = *x++;
while(x < xe);
b1->wds = n1 - 1;
Bfree(b);
return b1;
}
int
cmp
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
{
ULong *xa, *xa0, *xb, *xb0;
int i, j;
i = a->wds;
j = b->wds;
#ifdef DEBUG
if (i > 1 && !a->x[i-1])
Bug("cmp called with a->x[a->wds-1] == 0");
if (j > 1 && !b->x[j-1])
Bug("cmp called with b->x[b->wds-1] == 0");
#endif
if (i -= j)
return i;
xa0 = a->x;
xa = xa0 + j;
xb0 = b->x;
xb = xb0 + j;
for(;;) {
if (*--xa != *--xb)
return *xa < *xb ? -1 : 1;
if (xa <= xa0)
break;
}
return 0;
}
Bigint *
diff
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
{
Bigint *c;
int i, wa, wb;
ULong *xa, *xae, *xb, *xbe, *xc;
#ifdef ULLong
ULLong borrow, y;
#else
ULong borrow, y;
#ifdef Pack_32
ULong z;
#endif
#endif
i = cmp(a,b);
if (!i) {
c = Balloc(0);
if (c == NULL)
return (NULL);
c->wds = 1;
c->x[0] = 0;
return c;
}
if (i < 0) {
c = a;
a = b;
b = c;
i = 1;
}
else
i = 0;
c = Balloc(a->k);
if (c == NULL)
return (NULL);
c->sign = i;
wa = a->wds;
xa = a->x;
xae = xa + wa;
wb = b->wds;
xb = b->x;
xbe = xb + wb;
xc = c->x;
borrow = 0;
#ifdef ULLong
do {
y = (ULLong)*xa++ - *xb++ - borrow;
borrow = y >> 32 & 1UL;
*xc++ = y & 0xffffffffUL;
}
while(xb < xbe);
while(xa < xae) {
y = *xa++ - borrow;
borrow = y >> 32 & 1UL;
*xc++ = y & 0xffffffffUL;
}
#else
#ifdef Pack_32
do {
y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
while(xb < xbe);
while(xa < xae) {
y = (*xa & 0xffff) - borrow;
borrow = (y & 0x10000) >> 16;
z = (*xa++ >> 16) - borrow;
borrow = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
#else
do {
y = *xa++ - *xb++ - borrow;
borrow = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
while(xb < xbe);
while(xa < xae) {
y = *xa++ - borrow;
borrow = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
#endif
#endif
while(!*--xc)
wa--;
c->wds = wa;
return c;
}
double
b2d
#ifdef KR_headers
(a, e) Bigint *a; int *e;
#else
(Bigint *a, int *e)
#endif
{
ULong *xa, *xa0, w, y, z;
int k;
U d;
#ifdef VAX
ULong d0, d1;
#else
#define d0 word0(&d)
#define d1 word1(&d)
#endif
xa0 = a->x;
xa = xa0 + a->wds;
y = *--xa;
#ifdef DEBUG
if (!y) Bug("zero y in b2d");
#endif
k = hi0bits(y);
*e = 32 - k;
#ifdef Pack_32
if (k < Ebits) {
d0 = Exp_1 | y >> (Ebits - k);
w = xa > xa0 ? *--xa : 0;
d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
if (k -= Ebits) {
d0 = Exp_1 | y << k | z >> (32 - k);
y = xa > xa0 ? *--xa : 0;
d1 = z << k | y >> (32 - k);
}
else {
d0 = Exp_1 | y;
d1 = z;
}
#else
if (k < Ebits + 16) {
z = xa > xa0 ? *--xa : 0;
d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
w = xa > xa0 ? *--xa : 0;
y = xa > xa0 ? *--xa : 0;
d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
w = xa > xa0 ? *--xa : 0;
k -= Ebits + 16;
d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
y = xa > xa0 ? *--xa : 0;
d1 = w << k + 16 | y << k;
#endif
ret_d:
#ifdef VAX
word0(&d) = d0 >> 16 | d0 << 16;
word1(&d) = d1 >> 16 | d1 << 16;
#endif
return dval(&d);
}
#undef d0
#undef d1
Bigint *
d2b
#ifdef KR_headers
(dd, e, bits) double dd; int *e, *bits;
#else
(double dd, int *e, int *bits)
#endif
{
Bigint *b;
U d;
#ifndef Sudden_Underflow
int i;
#endif
int de, k;
ULong *x, y, z;
#ifdef VAX
ULong d0, d1;
#else
#define d0 word0(&d)
#define d1 word1(&d)
#endif
d.d = dd;
#ifdef VAX
d0 = word0(&d) >> 16 | word0(&d) << 16;
d1 = word1(&d) >> 16 | word1(&d) << 16;
#endif
#ifdef Pack_32
b = Balloc(1);
#else
b = Balloc(2);
#endif
if (b == NULL)
return (NULL);
x = b->x;
z = d0 & Frac_mask;
d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
#ifdef Sudden_Underflow
de = (int)(d0 >> Exp_shift);
#ifndef IBM
z |= Exp_msk11;
#endif
#else
if ( (de = (int)(d0 >> Exp_shift)) !=0)
z |= Exp_msk1;
#endif
#ifdef Pack_32
if ( (y = d1) !=0) {
if ( (k = lo0bits(&y)) !=0) {
x[0] = y | z << (32 - k);
z >>= k;
}
else
x[0] = y;
#ifndef Sudden_Underflow
i =
#endif
b->wds = (x[1] = z) !=0 ? 2 : 1;
}
else {
k = lo0bits(&z);
x[0] = z;
#ifndef Sudden_Underflow
i =
#endif
b->wds = 1;
k += 32;
}
#else
if ( (y = d1) !=0) {
if ( (k = lo0bits(&y)) !=0)
if (k >= 16) {
x[0] = y | z << 32 - k & 0xffff;
x[1] = z >> k - 16 & 0xffff;
x[2] = z >> k;
i = 2;
}
else {
x[0] = y & 0xffff;
x[1] = y >> 16 | z << 16 - k & 0xffff;
x[2] = z >> k & 0xffff;
x[3] = z >> k+16;
i = 3;
}
else {
x[0] = y & 0xffff;
x[1] = y >> 16;
x[2] = z & 0xffff;
x[3] = z >> 16;
i = 3;
}
}
else {
#ifdef DEBUG
if (!z)
Bug("Zero passed to d2b");
#endif
k = lo0bits(&z);
if (k >= 16) {
x[0] = z;
i = 0;
}
else {
x[0] = z & 0xffff;
x[1] = z >> 16;
i = 1;
}
k += 32;
}
while(!x[i])
--i;
b->wds = i + 1;
#endif
#ifndef Sudden_Underflow
if (de) {
#endif
#ifdef IBM
*e = (de - Bias - (P-1) << 2) + k;
*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
#else
*e = de - Bias - (P-1) + k;
*bits = P - k;
#endif
#ifndef Sudden_Underflow
}
else {
*e = de - Bias - (P-1) + 1 + k;
#ifdef Pack_32
*bits = 32*i - hi0bits(x[i-1]);
#else
*bits = (i+2)*16 - hi0bits(x[i]);
#endif
}
#endif
return b;
}
#undef d0
#undef d1
CONST double
#ifdef IEEE_Arith
bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
};
#else
#ifdef IBM
bigtens[] = { 1e16, 1e32, 1e64 };
CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
#else
bigtens[] = { 1e16, 1e32 };
CONST double tinytens[] = { 1e-16, 1e-32 };
#endif
#endif
CONST double
tens[] = {
1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1e20, 1e21, 1e22
#ifdef VAX
, 1e23, 1e24
#endif
};
char *
#ifdef KR_headers
strcp_D2A(a, b) char *a; char *b;
#else
strcp_D2A(char *a, CONST char *b)
#endif
{
while((*a = *b++))
a++;
return a;
}
#ifdef NO_STRING_H
Char *
#ifdef KR_headers
memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
#else
memcpy_D2A(void *a1, void *b1, size_t len)
#endif
{
char *a = (char*)a1, *ae = a + len;
char *b = (char*)b1, *a0 = a;
while(a < ae)
*a++ = *b++;
return a0;
}
#endif /* NO_STRING_H */

201
libc/upstream-openbsd/lib/libc/gdtoa/smisc.c Normal file
View File

@ -0,0 +1,201 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
Bigint *
s2b
#ifdef KR_headers
(s, nd0, nd, y9, dplen) CONST char *s; int dplen, nd0, nd; ULong y9;
#else
(CONST char *s, int nd0, int nd, ULong y9, int dplen)
#endif
{
Bigint *b;
int i, k;
Long x, y;
x = (nd + 8) / 9;
for(k = 0, y = 1; x > y; y <<= 1, k++) ;
#ifdef Pack_32
b = Balloc(k);
if (b == NULL)
return (NULL);
b->x[0] = y9;
b->wds = 1;
#else
b = Balloc(k+1);
if (b == NULL)
return (NULL);
b->x[0] = y9 & 0xffff;
b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
#endif
i = 9;
if (9 < nd0) {
s += 9;
do {
b = multadd(b, 10, *s++ - '0');
if (b == NULL)
return (NULL);
} while(++i < nd0);
s += dplen;
}
else
s += dplen + 9;
for(; i < nd; i++) {
b = multadd(b, 10, *s++ - '0');
if (b == NULL)
return (NULL);
}
return b;
}
double
ratio
#ifdef KR_headers
(a, b) Bigint *a, *b;
#else
(Bigint *a, Bigint *b)
#endif
{
U da, db;
int k, ka, kb;
dval(&da) = b2d(a, &ka);
dval(&db) = b2d(b, &kb);
k = ka - kb + ULbits*(a->wds - b->wds);
#ifdef IBM
if (k > 0) {
word0(&da) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(&da) *= 1 << k;
}
else {
k = -k;
word0(&db) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(&db) *= 1 << k;
}
#else
if (k > 0)
word0(&da) += k*Exp_msk1;
else {
k = -k;
word0(&db) += k*Exp_msk1;
}
#endif
return dval(&da) / dval(&db);
}
#ifdef INFNAN_CHECK
int
match
#ifdef KR_headers
(sp, t) char **sp, *t;
#else
(CONST char **sp, char *t)
#endif
{
int c, d;
CONST char *s = *sp;
while( (d = *t++) !=0) {
if ((c = *++s) >= 'A' && c <= 'Z')
c += 'a' - 'A';
if (c != d)
return 0;
}
*sp = s + 1;
return 1;
}
#endif /* INFNAN_CHECK */
void
#ifdef KR_headers
copybits(c, n, b) ULong *c; int n; Bigint *b;
#else
copybits(ULong *c, int n, Bigint *b)
#endif
{
ULong *ce, *x, *xe;
#ifdef Pack_16
int nw, nw1;
#endif
ce = c + ((n-1) >> kshift) + 1;
x = b->x;
#ifdef Pack_32
xe = x + b->wds;
while(x < xe)
*c++ = *x++;
#else
nw = b->wds;
nw1 = nw & 1;
for(xe = x + (nw - nw1); x < xe; x += 2)
Storeinc(c, x[1], x[0]);
if (nw1)
*c++ = *x;
#endif
while(c < ce)
*c++ = 0;
}
ULong
#ifdef KR_headers
any_on(b, k) Bigint *b; int k;
#else
any_on(Bigint *b, int k)
#endif
{
int n, nwds;
ULong *x, *x0, x1, x2;
x = b->x;
nwds = b->wds;
n = k >> kshift;
if (n > nwds)
n = nwds;
else if (n < nwds && (k &= kmask)) {
x1 = x2 = x[n];
x1 >>= k;
x1 <<= k;
if (x1 != x2)
return 1;
}
x0 = x;
x += n;
while(x > x0)
if (*--x)
return 1;
return 0;
}

1105
libc/upstream-openbsd/lib/libc/gdtoa/strtod.c Normal file
View File

File diff suppressed because it is too large Load Diff

1154
libc/upstream-openbsd/lib/libc/gdtoa/strtodg.c Normal file
View File

File diff suppressed because it is too large Load Diff

81
libc/upstream-openbsd/lib/libc/gdtoa/strtof.c Normal file
View File

@ -0,0 +1,81 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
float
#ifdef KR_headers
strtof(s, sp) CONST char *s; char **sp;
#else
strtof(CONST char *s, char **sp)
#endif
{
static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, SI };
ULong bits[1];
Long exp;
int k;
union { ULong L[1]; float f; } u;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
k = strtodg(s, sp, fpi, &exp, bits);
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
case STRTOG_Zero:
u.L[0] = 0;
break;
case STRTOG_Normal:
case STRTOG_NaNbits:
u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:
u.L[0] = bits[0];
break;
case STRTOG_NoMemory:
errno = ERANGE;
/* FALLTHROUGH */
case STRTOG_Infinite:
u.L[0] = 0x7f800000;
break;
case STRTOG_NaN:
u.L[0] = f_QNAN;
}
if (k & STRTOG_Neg)
u.L[0] |= 0x80000000L;
return u.f;
}

120
libc/upstream-openbsd/lib/libc/gdtoa/strtorQ.c Normal file
View File

@ -0,0 +1,120 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
#undef _0
#undef _1
/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
#ifdef IEEE_MC68k
#define _0 0
#define _1 1
#define _2 2
#define _3 3
#endif
#ifdef IEEE_8087
#define _0 3
#define _1 2
#define _2 1
#define _3 0
#endif
void
#ifdef KR_headers
ULtoQ(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
#else
ULtoQ(ULong *L, ULong *bits, Long exp, int k)
#endif
{
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
case STRTOG_Zero:
L[0] = L[1] = L[2] = L[3] = 0;
break;
case STRTOG_Normal:
case STRTOG_NaNbits:
L[_3] = bits[0];
L[_2] = bits[1];
L[_1] = bits[2];
L[_0] = (bits[3] & ~0x10000) | ((exp + 0x3fff + 112) << 16);
break;
case STRTOG_Denormal:
L[_3] = bits[0];
L[_2] = bits[1];
L[_1] = bits[2];
L[_0] = bits[3];
break;
case STRTOG_NoMemory:
errno = ERANGE;
/* FALLTHROUGH */
case STRTOG_Infinite:
L[_0] = 0x7fff0000;
L[_1] = L[_2] = L[_3] = 0;
break;
case STRTOG_NaN:
L[0] = ld_QNAN0;
L[1] = ld_QNAN1;
L[2] = ld_QNAN2;
L[3] = ld_QNAN3;
}
if (k & STRTOG_Neg)
L[_0] |= 0x80000000L;
}
int
#ifdef KR_headers
strtorQ(s, sp, rounding, L) CONST char *s; char **sp; int rounding; void *L;
#else
strtorQ(CONST char *s, char **sp, int rounding, void *L)
#endif
{
static FPI fpi0 = { 113, 1-16383-113+1, 32766-16383-113+1, 1, SI };
FPI *fpi, fpi1;
ULong bits[4];
Long exp;
int k;
fpi = &fpi0;
if (rounding != FPI_Round_near) {
fpi1 = fpi0;
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
ULtoQ((ULong*)L, bits, exp, k);
return k;
}

96
libc/upstream-openbsd/lib/libc/gdtoa/strtord.c Normal file
View File

@ -0,0 +1,96 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 2000 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
void
#ifdef KR_headers
ULtod(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
#else
ULtod(ULong *L, ULong *bits, Long exp, int k)
#endif
{
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
case STRTOG_Zero:
L[0] = L[1] = 0;
break;
case STRTOG_Denormal:
L[_1] = bits[0];
L[_0] = bits[1];
break;
case STRTOG_Normal:
case STRTOG_NaNbits:
L[_1] = bits[0];
L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
break;
case STRTOG_NoMemory:
errno = ERANGE;
/* FALLTHROUGH */
case STRTOG_Infinite:
L[_0] = 0x7ff00000;
L[_1] = 0;
break;
case STRTOG_NaN:
L[0] = d_QNAN0;
L[1] = d_QNAN1;
}
if (k & STRTOG_Neg)
L[_0] |= 0x80000000L;
}
int
#ifdef KR_headers
strtord(s, sp, rounding, d) CONST char *s; char **sp; int rounding; double *d;
#else
strtord(CONST char *s, char **sp, int rounding, double *d)
#endif
{
static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
FPI *fpi, fpi1;
ULong bits[2];
Long exp;
int k;
fpi = &fpi0;
if (rounding != FPI_Round_near) {
fpi1 = fpi0;
fpi1.rounding = rounding;
fpi = &fpi1;
}
k = strtodg(s, sp, fpi, &exp, bits);
ULtod((ULong*)d, bits, exp, k);
return k;
}

102
libc/upstream-openbsd/lib/libc/gdtoa/sum.c Normal file
View File

@ -0,0 +1,102 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
Bigint *
#ifdef KR_headers
sum(a, b) Bigint *a; Bigint *b;
#else
sum(Bigint *a, Bigint *b)
#endif
{
Bigint *c;
ULong carry, *xc, *xa, *xb, *xe, y;
#ifdef Pack_32
ULong z;
#endif
if (a->wds < b->wds) {
c = b; b = a; a = c;
}
c = Balloc(a->k);
if (c == NULL)
return (NULL);
c->wds = a->wds;
carry = 0;
xa = a->x;
xb = b->x;
xc = c->x;
xe = xc + b->wds;
#ifdef Pack_32
do {
y = (*xa & 0xffff) + (*xb & 0xffff) + carry;
carry = (y & 0x10000) >> 16;
z = (*xa++ >> 16) + (*xb++ >> 16) + carry;
carry = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
while(xc < xe);
xe += a->wds - b->wds;
while(xc < xe) {
y = (*xa & 0xffff) + carry;
carry = (y & 0x10000) >> 16;
z = (*xa++ >> 16) + carry;
carry = (z & 0x10000) >> 16;
Storeinc(xc, z, y);
}
#else
do {
y = *xa++ + *xb++ + carry;
carry = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
while(xc < xe);
xe += a->wds - b->wds;
while(xc < xe) {
y = *xa++ + carry;
carry = (y & 0x10000) >> 16;
*xc++ = y & 0xffff;
}
#endif
if (carry) {
if (c->wds == c->maxwds) {
b = Balloc(c->k + 1);
if (b == NULL)
return (NULL);
Bcopy(b, c);
Bfree(c);
c = b;
}
c->x[c->wds++] = 1;
}
return c;
}

70
libc/upstream-openbsd/lib/libc/gdtoa/ulp.c Normal file
View File

@ -0,0 +1,70 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 1998, 1999 by Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* Please send bug reports to David M. Gay (dmg at acm dot org,
* with " at " changed at "@" and " dot " changed to "."). */
#include "gdtoaimp.h"
double
ulp
#ifdef KR_headers
(x) U *x;
#else
(U *x)
#endif
{
Long L;
U a;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
if (L > 0) {
#endif
#ifdef IBM
L |= Exp_msk1 >> 4;
#endif
word0(&a) = L;
word1(&a) = 0;
#ifndef Sudden_Underflow
}
else {
L = -L >> Exp_shift;
if (L < Exp_shift) {
word0(&a) = 0x80000 >> L;
word1(&a) = 0;
}
else {
word0(&a) = 0;
L -= Exp_shift;
word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
}
}
#endif
return dval(&a);
}

26
libc/upstream-openbsd/thread_private.h Normal file
View File

@ -0,0 +1,26 @@
/*
* Copyright (C) 2014 The Android Open Source Project
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef _THREAD_PRIVATE_H_
#define _THREAD_PRIVATE_H_
#include <pthread.h>
/* Note that these aren't compatible with the usual OpenBSD ones which lazy-initialize! */
#define _MUTEX_LOCK(l) pthread_mutex_lock((pthread_mutex_t*) l)
#define _MUTEX_UNLOCK(l) pthread_mutex_unlock((pthread_mutex_t*) l)
#endif /* _THREAD_PRIVATE_H_ */

18
tests/stdio_test.cpp
View File

@ -306,19 +306,19 @@ TEST(stdio, snprintf_smoke) {
TEST(stdio, snprintf_f_special) {
char buf[BUFSIZ];
snprintf(buf, sizeof(buf), "%f", nanf(""));
EXPECT_STREQ("NaN", buf);
EXPECT_STRCASEEQ("NaN", buf);
snprintf(buf, sizeof(buf), "%f", HUGE_VALF);
EXPECT_STREQ("Inf", buf);
EXPECT_STRCASEEQ("Inf", buf);
}
TEST(stdio, snprintf_g_special) {
char buf[BUFSIZ];
snprintf(buf, sizeof(buf), "%g", nan(""));
EXPECT_STREQ("NaN", buf);
EXPECT_STRCASEEQ("NaN", buf);
snprintf(buf, sizeof(buf), "%g", HUGE_VAL);
EXPECT_STREQ("Inf", buf);
EXPECT_STRCASEEQ("Inf", buf);
}
TEST(stdio, snprintf_d_INT_MAX) {
@ -365,6 +365,16 @@ TEST(stdio, snprintf_lld_LLONG_MIN) {
EXPECT_STREQ("-9223372036854775808", buf);
}
TEST(stdio, snprintf_e) {
char buf[BUFSIZ];
snprintf(buf, sizeof(buf), "%e", 1.5);
EXPECT_STREQ("1.500000e+00", buf);
snprintf(buf, sizeof(buf), "%Le", 1.5l);
EXPECT_STREQ("1.500000e+00", buf);
}
TEST(stdio, popen) {
FILE* fp = popen("cat /proc/version", "r");
ASSERT_TRUE(fp != NULL);