Commit 1d92226b by Jakub Jelinek Committed by Jakub Jelinek

re PR fortran/47642 (real(kind=16) - libquadmath - segfault on amd64 FreeBSD)

	PR fortran/47642
	* libquadmath.texi (quadmath_snprintf): Document.
	(quadmath_flt128tostr): Remove.
	* Makefile.am (libquadmath_la_SOURCES): Add printf/*.c.
	Remove  quadmath_io.c, gdtoa/gdtoa.c, gdtoa/g__fmt.c,
	gdtoa/g_Qfmt.c, gdtoa/dmisc.c and gdtoa/ulp.c.
	* quadmath.h (quadmath_snprintf): New prototype.
	(quadmath_flt128tostr): Remove.
	* quadmath_weak.h (quadmath_snprintf): Add.
	(quadmath_flt128tostr): Remove.
	* configure.ac: New AC_CHECK_HEADERS headers: langinfo.h, wchar.h,
	wctype.h, limits.h, ctype.h, printf.h, errno.h.
	(AC_USE_SYSTEM_EXTENSIONS): Add.
	(HAVE_HIDDEN_VISIBILITY, HAVE_PRINTF_HOOKS,
	USE_LOCALE_SUPPORT, USE_I18N_NUMBER_H): New checks.
	* quadmath.map (QUADMATH_1.0): Add quadmath_snprintf.  Remove
	quadmath_flt128tostr.
	* printf/printf_fphex.c: New file.
	* printf/_itowa.h: New file.
	* printf/mul_n.c: New file.
	* printf/quadmath-printf.h: New file.
	* printf/submul_1.c: New file.
	* printf/quadmath-printf.c: New file.
	* printf/gmp-impl.h: New file.
	* printf/lshift.c: New file.
	* printf/fpioconst.h: New file.
	* printf/add_n.c: New file.
	* printf/cmp.c: New file.
	* printf/sub_n.c: New file.
	* printf/mul.c: New file.
	* printf/divrem.c: New file.
	* printf/addmul_1.c: New file.
	* printf/printf_fp.c: New file.
	* printf/_itoa.h: New file.
	* printf/fpioconst.c: New file.
	* printf/_i18n_number.h: New file.
	* printf/flt1282mpn.c: New file.
	* printf/rshift.c: New file.
	* printf/mul_1.c: New file.
	* quadmath_io.c: Removed.
	* gdtoa/gdtoa.c: Removed.
	* gdtoa/g__fmt.c: Removed.
	* gdtoa/g_Qfmt.c: Removed.
	* gdtoa/dmisc.c: Removed.
	* gdtoa/ulp.c: Removed.
	* config.h.in: Regenerated.
	* configure: Regenerated.
	* Makefile.in: Regenerated.

	* io/write_float.def (DTOAQ): Use quadmath_snprintf instead of
	quadmath_flt128tostr.
	* io/transfer128.c (tmp2): Initialize to quadmath_snprintf instead
	of quadmath_flt128tostr.

From-SVN: r170135
parent 944f4bb3
2011-02-14 Jakub Jelinek <jakub@redhat.com>
PR fortran/47642
* io/write_float.def (DTOAQ): Use quadmath_snprintf instead of
quadmath_flt128tostr.
* io/transfer128.c (tmp2): Initialize to quadmath_snprintf instead
of quadmath_flt128tostr.
2011-02-13 Ralf Wildenhues <Ralf.Wildenhues@gmx.de> 2011-02-13 Ralf Wildenhues <Ralf.Wildenhues@gmx.de>
* Makefile.in: Regenerate. * Makefile.in: Regenerate.
......
/* Copyright (C) 2010 /* Copyright (C) 2010, 2011
Free Software Foundation, Inc. Free Software Foundation, Inc.
This file is part of the GNU Fortran runtime library (libgfortran). This file is part of the GNU Fortran runtime library (libgfortran).
...@@ -62,12 +62,12 @@ export_proto(transfer_complex128_write); ...@@ -62,12 +62,12 @@ export_proto(transfer_complex128_write);
/* Make sure that libquadmath is pulled in. The functions strtoflt128 /* Make sure that libquadmath is pulled in. The functions strtoflt128
and quadmath_flt128tostr are weakly referrenced in convert_real and and quadmath_snprintf are weakly referrenced in convert_real and
write_float; the pointer assignment with USED attribute make sure write_float; the pointer assignment with USED attribute make sure
that there is a non-weakref dependence if the quadmath functions that there is a non-weakref dependence if the quadmath functions
are used. That avoids segfault when libquadmath is statically linked. */ are used. That avoids segfault when libquadmath is statically linked. */
static void __attribute__((used)) *tmp1 = strtoflt128; static void __attribute__((used)) *tmp1 = strtoflt128;
static void __attribute__((used)) *tmp2 = quadmath_flt128tostr; static void __attribute__((used)) *tmp2 = quadmath_snprintf;
void void
transfer_real128 (st_parameter_dt *dtp, void *p, int kind) transfer_real128 (st_parameter_dt *dtp, void *p, int kind)
......
/* Copyright (C) 2007, 2008, 2009, 2010 Free Software Foundation, Inc. /* Copyright (C) 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
Contributed by Andy Vaught Contributed by Andy Vaught
Write float code factoring to this file by Jerry DeLisle Write float code factoring to this file by Jerry DeLisle
F2003 I/O support contributed by Jerry DeLisle F2003 I/O support contributed by Jerry DeLisle
...@@ -1000,7 +1000,9 @@ sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \ ...@@ -1000,7 +1000,9 @@ sprintf (buffer, "%+-#" STR(MIN_FIELD_WIDTH) ".*" \
#if defined(GFC_REAL_16_IS_FLOAT128) #if defined(GFC_REAL_16_IS_FLOAT128)
#define DTOAQ \ #define DTOAQ \
__qmath_(quadmath_flt128tostr) (buffer, size, ndigits - 1, tmp); __qmath_(quadmath_snprintf) (buffer, sizeof buffer, \
"%+-#" STR(MIN_FIELD_WIDTH) ".*" \
"Qe", ndigits - 1, tmp);
#endif #endif
#define WRITE_FLOAT(x,y)\ #define WRITE_FLOAT(x,y)\
......
2011-02-14 Jakub Jelinek <jakub@redhat.com>
PR fortran/47642
* libquadmath.texi (quadmath_snprintf): Document.
(quadmath_flt128tostr): Remove.
* Makefile.am (libquadmath_la_SOURCES): Add printf/*.c.
Remove quadmath_io.c, gdtoa/gdtoa.c, gdtoa/g__fmt.c,
gdtoa/g_Qfmt.c, gdtoa/dmisc.c and gdtoa/ulp.c.
* quadmath.h (quadmath_snprintf): New prototype.
(quadmath_flt128tostr): Remove.
* quadmath_weak.h (quadmath_snprintf): Add.
(quadmath_flt128tostr): Remove.
* configure.ac: New AC_CHECK_HEADERS headers: langinfo.h, wchar.h,
wctype.h, limits.h, ctype.h, printf.h, errno.h.
(AC_USE_SYSTEM_EXTENSIONS): Add.
(HAVE_HIDDEN_VISIBILITY, HAVE_PRINTF_HOOKS,
USE_LOCALE_SUPPORT, USE_I18N_NUMBER_H): New checks.
* quadmath.map (QUADMATH_1.0): Add quadmath_snprintf. Remove
quadmath_flt128tostr.
* printf/printf_fphex.c: New file.
* printf/_itowa.h: New file.
* printf/mul_n.c: New file.
* printf/quadmath-printf.h: New file.
* printf/submul_1.c: New file.
* printf/quadmath-printf.c: New file.
* printf/gmp-impl.h: New file.
* printf/lshift.c: New file.
* printf/fpioconst.h: New file.
* printf/add_n.c: New file.
* printf/cmp.c: New file.
* printf/sub_n.c: New file.
* printf/mul.c: New file.
* printf/divrem.c: New file.
* printf/addmul_1.c: New file.
* printf/printf_fp.c: New file.
* printf/_itoa.h: New file.
* printf/fpioconst.c: New file.
* printf/_i18n_number.h: New file.
* printf/flt1282mpn.c: New file.
* printf/rshift.c: New file.
* printf/mul_1.c: New file.
* quadmath_io.c: Removed.
* gdtoa/gdtoa.c: Removed.
* gdtoa/g__fmt.c: Removed.
* gdtoa/g_Qfmt.c: Removed.
* gdtoa/dmisc.c: Removed.
* gdtoa/ulp.c: Removed.
* config.h.in: Regenerated.
* configure: Regenerated.
* Makefile.in: Regenerated.
2011-02-13 Ralf Wildenhues <Ralf.Wildenhues@gmx.de> 2011-02-13 Ralf Wildenhues <Ralf.Wildenhues@gmx.de>
* Makefile.in: Regenerate. * Makefile.in: Regenerate.
......
...@@ -45,10 +45,9 @@ libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include ...@@ -45,10 +45,9 @@ libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
libquadmath_la_SOURCES = \ libquadmath_la_SOURCES = \
gdtoa/arith.h gdtoa/gdtoa_fltrnds.h gdtoa/gd_qnan.h gdtoa/gdtoaimp.h \ gdtoa/arith.h gdtoa/gdtoa_fltrnds.h gdtoa/gd_qnan.h gdtoa/gdtoaimp.h \
gdtoa/gdtoa.h quadmath-imp.h \ gdtoa/gdtoa.h quadmath-imp.h \
gdtoa/dmisc.c gdtoa/gdtoa.c gdtoa/hd_init.c gdtoa/smisc.c gdtoa/sum.c \ gdtoa/hd_init.c gdtoa/smisc.c gdtoa/sum.c \
gdtoa/g_Qfmt.c gdtoa/gethex.c gdtoa/hexnan.c gdtoa/strtodg.c \ gdtoa/gethex.c gdtoa/hexnan.c gdtoa/strtodg.c \
gdtoa/ulp.c gdtoa/g__fmt.c gdtoa/gmisc.c gdtoa/misc.c gdtoa/strtopQ.c \ gdtoa/gmisc.c gdtoa/misc.c gdtoa/strtopQ.c \
quadmath_io.c \
math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \ math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \
math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \ math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \ math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
...@@ -65,7 +64,11 @@ libquadmath_la_SOURCES = \ ...@@ -65,7 +64,11 @@ libquadmath_la_SOURCES = \
math/cacoshq.c math/cacosq.c math/casinhq.c math/casinq.c \ math/cacoshq.c math/cacosq.c math/casinhq.c math/casinq.c \
math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \ math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \
math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \ math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \
math/llrintq.c math/log2q.c math/lrintq.c math/nearbyintq.c math/remquoq.c math/llrintq.c math/log2q.c math/lrintq.c math/nearbyintq.c math/remquoq.c \
printf/addmul_1.c printf/add_n.c printf/cmp.c printf/divrem.c \
printf/flt1282mpn.c printf/fpioconst.c printf/lshift.c printf/mul_1.c \
printf/mul_n.c printf/mul.c printf/printf_fphex.c printf/printf_fp.c \
printf/quadmath-printf.c printf/rshift.c printf/submul_1.c printf/sub_n.c
# Work around what appears to be a GNU make bug handling MAKEFLAGS # Work around what appears to be a GNU make bug handling MAKEFLAGS
......
...@@ -3,9 +3,15 @@ ...@@ -3,9 +3,15 @@
/* libm includes cbrtl */ /* libm includes cbrtl */
#undef HAVE_CBRTL #undef HAVE_CBRTL
/* Define to 1 if you have the <ctype.h> header file. */
#undef HAVE_CTYPE_H
/* Define to 1 if you have the <dlfcn.h> header file. */ /* Define to 1 if you have the <dlfcn.h> header file. */
#undef HAVE_DLFCN_H #undef HAVE_DLFCN_H
/* Define to 1 if you have the <errno.h> header file. */
#undef HAVE_ERRNO_H
/* libm includes feholdexcept */ /* libm includes feholdexcept */
#undef HAVE_FEHOLDEXCEPT #undef HAVE_FEHOLDEXCEPT
...@@ -24,12 +30,27 @@ ...@@ -24,12 +30,27 @@
/* libm includes feupdateenv */ /* libm includes feupdateenv */
#undef HAVE_FEUPDATEENV #undef HAVE_FEUPDATEENV
/* __attribute__((visibility ("hidden"))) supported */
#undef HAVE_HIDDEN_VISIBILITY
/* Define to 1 if you have the <inttypes.h> header file. */ /* Define to 1 if you have the <inttypes.h> header file. */
#undef HAVE_INTTYPES_H #undef HAVE_INTTYPES_H
/* Define to 1 if you have the <langinfo.h> header file. */
#undef HAVE_LANGINFO_H
/* Define to 1 if you have the <limits.h> header file. */
#undef HAVE_LIMITS_H
/* Define to 1 if you have the <memory.h> header file. */ /* Define to 1 if you have the <memory.h> header file. */
#undef HAVE_MEMORY_H #undef HAVE_MEMORY_H
/* Define to 1 if you have the <printf.h> header file. */
#undef HAVE_PRINTF_H
/* GNU C Library stype printf hooks supported */
#undef HAVE_PRINTF_HOOKS
/* libm includes sqrtl */ /* libm includes sqrtl */
#undef HAVE_SQRTL #undef HAVE_SQRTL
...@@ -54,6 +75,12 @@ ...@@ -54,6 +75,12 @@
/* Define to 1 if you have the <unistd.h> header file. */ /* Define to 1 if you have the <unistd.h> header file. */
#undef HAVE_UNISTD_H #undef HAVE_UNISTD_H
/* Define to 1 if you have the <wchar.h> header file. */
#undef HAVE_WCHAR_H
/* Define to 1 if you have the <wctype.h> header file. */
#undef HAVE_WCTYPE_H
/* Define to the sub-directory in which libtool stores uninstalled libraries. /* Define to the sub-directory in which libtool stores uninstalled libraries.
*/ */
#undef LT_OBJDIR #undef LT_OBJDIR
...@@ -85,5 +112,43 @@ ...@@ -85,5 +112,43 @@
/* Define to 1 if you have the ANSI C header files. */ /* Define to 1 if you have the ANSI C header files. */
#undef STDC_HEADERS #undef STDC_HEADERS
/* whether i18n number rewriting can be supported */
#undef USE_I18N_NUMBER_H
/* whether nl_langinfo is sufficiently supported */
#undef USE_LOCALE_SUPPORT
/* Enable extensions on AIX 3, Interix. */
#ifndef _ALL_SOURCE
# undef _ALL_SOURCE
#endif
/* Enable GNU extensions on systems that have them. */
#ifndef _GNU_SOURCE
# undef _GNU_SOURCE
#endif
/* Enable threading extensions on Solaris. */
#ifndef _POSIX_PTHREAD_SEMANTICS
# undef _POSIX_PTHREAD_SEMANTICS
#endif
/* Enable extensions on HP NonStop. */
#ifndef _TANDEM_SOURCE
# undef _TANDEM_SOURCE
#endif
/* Enable general extensions on Solaris. */
#ifndef __EXTENSIONS__
# undef __EXTENSIONS__
#endif
/* Version number of package */ /* Version number of package */
#undef VERSION #undef VERSION
/* Define to 1 if on MINIX. */
#undef _MINIX
/* Define to 2 if the system does not provide POSIX.1 features except with
this defined. */
#undef _POSIX_1_SOURCE
/* Define to 1 if you need to in order for `stat' and other things to work. */
#undef _POSIX_SOURCE
This source diff could not be displayed because it is too large. You can view the blob instead.
...@@ -42,6 +42,8 @@ AC_MSG_RESULT($version_specific_libs) ...@@ -42,6 +42,8 @@ AC_MSG_RESULT($version_specific_libs)
GCC_NO_EXECUTABLES GCC_NO_EXECUTABLES
AC_USE_SYSTEM_EXTENSIONS
# See if makeinfo has been installed and is modern enough # See if makeinfo has been installed and is modern enough
# that we can use it. # that we can use it.
ACX_CHECK_PROG_VER([MAKEINFO], [makeinfo], [--version], ACX_CHECK_PROG_VER([MAKEINFO], [makeinfo], [--version],
...@@ -110,7 +112,7 @@ esac ...@@ -110,7 +112,7 @@ esac
AC_SUBST(toolexecdir) AC_SUBST(toolexecdir)
AC_SUBST(toolexeclibdir) AC_SUBST(toolexeclibdir)
AC_CHECK_HEADERS(fenv.h) AC_CHECK_HEADERS(fenv.h langinfo.h wchar.h wctype.h limits.h ctype.h printf.h errno.h)
# If available, sqrtl and cbrtl speed up the calculation - # If available, sqrtl and cbrtl speed up the calculation -
# but they are not required # but they are not required
...@@ -146,6 +148,16 @@ else ...@@ -146,6 +148,16 @@ else
fi fi
fi fi
# Check for hidden visibility (copied from libssp).
AC_MSG_CHECKING([whether hidden visibility is supported])
AC_TRY_COMPILE([
void __attribute__((visibility ("hidden"))) bar (void) {}],,
[quadmath_hidden=yes],[quadmath_hidden=no])
AC_MSG_RESULT($quadmath_hidden)
if test x$quadmath_hidden = xyes; then
AC_DEFINE([HAVE_HIDDEN_VISIBILITY],[1],[__attribute__((visibility ("hidden"))) supported])
fi
# Check for symbol versioning (copied from libssp). # Check for symbol versioning (copied from libssp).
AC_MSG_CHECKING([whether symbol versioning is supported]) AC_MSG_CHECKING([whether symbol versioning is supported])
if test x$gcc_no_link = xyes; then if test x$gcc_no_link = xyes; then
...@@ -213,6 +225,75 @@ AC_CACHE_CHECK([whether __float128 is supported], [libquad_cv_have_float128], ...@@ -213,6 +225,75 @@ AC_CACHE_CHECK([whether __float128 is supported], [libquad_cv_have_float128],
])]) ])])
AM_CONDITIONAL(BUILD_LIBQUADMATH, [test "x$libquad_cv_have_float128" = xyes]) AM_CONDITIONAL(BUILD_LIBQUADMATH, [test "x$libquad_cv_have_float128" = xyes])
# Check for printf hook support.
AC_MSG_CHECKING([whether printf hooks are supported])
AC_TRY_COMPILE([
#include <printf.h>
#include <stdarg.h>
#include <stdlib.h>
extern void flt128_va (void *, va_list *);
extern int flt128_ais (const struct printf_info *, size_t, int *, int *);
extern int flt128_printf_fp (FILE *, const struct printf_info *, const void *const *);
],[
int pa_flt128 = register_printf_type (flt128_va);
int mod_Q = register_printf_modifier (L"Q");
int res = register_printf_specifier ('f', flt128_printf_fp, flt128_ais);
],
[quadmath_printf_hooks=yes],[quadmath_printf_hooks=no])
AC_MSG_RESULT($quadmath_printf_hooks)
if test x$quadmath_printf_hooks = xyes; then
AC_DEFINE([HAVE_PRINTF_HOOKS],[1],[GNU C Library stype printf hooks supported])
fi
# Check for whether locale support for quadmath_snprintf or Q printf hooks
# should be provided.
AC_MSG_CHECKING([whether locale support for quadmath_snprintf should be added])
AC_TRY_COMPILE([#include <langinfo.h>],[
const char *s;
s = nl_langinfo (DECIMAL_POINT);
s = nl_langinfo (MON_DECIMAL_POINT);
s = nl_langinfo (GROUPING);
s = nl_langinfo (MON_GROUPING);
s = nl_langinfo (THOUSANDS_SEP);
s = nl_langinfo (MON_THOUSANDS_SEP);
s = nl_langinfo (_NL_NUMERIC_DECIMAL_POINT_WC);
s = nl_langinfo (_NL_MONETARY_DECIMAL_POINT_WC);
s = nl_langinfo (_NL_NUMERIC_THOUSANDS_SEP_WC);
s = nl_langinfo (_NL_MONETARY_THOUSANDS_SEP_WC);
s = nl_langinfo (_NL_CTYPE_MB_CUR_MAX);
(void) s;
],
[quadmath_use_locale_support=yes],[quadmath_use_locale_support=no])
AC_MSG_RESULT($quadmath_use_locale_support)
if test x$quadmath_use_locale_support = xyes; then
AC_DEFINE([USE_LOCALE_SUPPORT],[1],[whether nl_langinfo is sufficiently supported])
fi
# Check for whether i18n number rewriting support for quadmath_snprintf
# or Q printf hooks should be provided.
AC_MSG_CHECKING([whether i18n number rewriting support for quadmath_snprintf should be added])
AC_TRY_COMPILE([#include <langinfo.h>
#include <limits.h>
#include <string.h>
#include <wchar.h>
#include <wctype.h>],[
const char *s;
char decimal[MB_LEN_MAX];
wctrans_t map = wctrans ("to_outpunct");
wint_t wdecimal = towctrans (L'.', map);
mbstate_t state;
memset (&state, '\0', sizeof (state));
wcrtomb (decimal, wdecimal, &state);
s = nl_langinfo (_NL_CTYPE_OUTDIGIT0_MB);
s = nl_langinfo (_NL_CTYPE_OUTDIGIT0_WC);
(void) s;
],
[quadmath_use_i18n_number_h=yes],[quadmath_use_i18n_number_h=no])
AC_MSG_RESULT($quadmath_use_i18n_number_h)
if test x$quadmath_use_i18n_number_h = xyes; then
AC_DEFINE([USE_I18N_NUMBER_H],[1],[whether i18n number rewriting can be supported])
fi
AC_CACHE_SAVE AC_CACHE_SAVE
if test ${multilib} = yes; then if test ${multilib} = yes; then
......
/****************************************************************
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);
*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);
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;
}
/****************************************************************
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
char*
#ifdef KR_headers
g_Qfmt(buf, V, ndig, bufsize) char *buf; char *V; int ndig; size_t bufsize;
#else
g_Qfmt(char *buf, void *V, int ndig, size_t bufsize)
#endif
{
static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0 };
char *b, *s, *se;
ULong bits[4], *L, sign;
int decpt, ex, i, mode;
#ifdef Honor_FLT_ROUNDS
#include "gdtoa_fltrnds.h"
#else
#define fpi &fpi0
#endif
if (ndig < 0)
ndig = 0;
if (bufsize < ndig + 10)
return 0;
L = (ULong*)V;
sign = L[_0] & 0x80000000L;
bits[3] = L[_0] & 0xffff;
bits[2] = L[_1];
bits[1] = L[_2];
bits[0] = L[_3];
b = buf;
if ( (ex = (L[_0] & 0x7fff0000L) >> 16) !=0) {
if (ex == 0x7fff) {
/* Infinity or NaN */
if (bits[0] | bits[1] | bits[2] | bits[3])
b = strcpy(b, "NaN");
else {
b = buf;
if (sign)
*b++ = '-';
b = strcpy(b, "Infinity");
}
return b;
}
i = STRTOG_Normal;
bits[3] |= 0x10000;
}
else if (bits[0] | bits[1] | bits[2] | bits[3]) {
i = STRTOG_Denormal;
ex = 1;
}
else {
#ifndef IGNORE_ZERO_SIGN
if (sign)
*b++ = '-';
#endif
*b++ = '0';
*b = 0;
return b;
}
ex -= 0x3fff + 112;
mode = 2;
if (ndig <= 0) {
if (bufsize < 48)
return 0;
mode = 0;
}
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
// FXC -- modifications
return g__fmt(buf, s, se, decpt, sign, bufsize);
// FXC -- end of modifications
}
/****************************************************************
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
char *
#ifdef KR_headers
g__fmt(b, s, se, decpt, sign, blen) char *b; char *s; char *se; int decpt; ULong sign; size_t blen;
#else
g__fmt(char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
#endif
{
int i, j, k;
char *be, *s0;
size_t len;
#ifdef USE_LOCALE
#ifdef NO_LOCALE_CACHE
char *decimalpoint = localeconv()->decimal_point;
size_t dlen = strlen(decimalpoint);
#else
char *decimalpoint;
static char *decimalpoint_cache;
static size_t dlen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
dlen = strlen(s0);
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
}
decimalpoint = s0;
#endif
#else
#define dlen 0
#endif
s0 = s;
len = (se-s) + dlen + 6; /* 6 = sign + e+dd + trailing null */
if (blen < len)
goto ret0;
be = b + blen - 1;
if (sign)
*b++ = '-';
if (decpt <= -4 || decpt > se - s + 5) {
*b++ = *s++;
if (*s) {
#ifdef USE_LOCALE
while((*b = *decimalpoint++))
++b;
#else
*b++ = '.';
#endif
while((*b = *s++) !=0)
b++;
}
*b++ = 'e';
/* sprintf(b, "%+.2d", decpt - 1); */
if (--decpt < 0) {
*b++ = '-';
decpt = -decpt;
}
else
*b++ = '+';
for(j = 2, k = 10; 10*k <= decpt; j++, k *= 10){}
for(;;) {
i = decpt / k;
if (b >= be)
goto ret0;
*b++ = i + '0';
if (--j <= 0)
break;
decpt -= i*k;
decpt *= 10;
}
*b = 0;
}
else if (decpt <= 0) {
#ifdef USE_LOCALE
while((*b = *decimalpoint++))
++b;
#else
*b++ = '.';
#endif
if (be < b - decpt + (se - s))
goto ret0;
for(; decpt < 0; decpt++)
*b++ = '0';
while((*b = *s++) != 0)
b++;
}
else {
while((*b = *s++) != 0) {
b++;
if (--decpt == 0 && *s) {
#ifdef USE_LOCALE
while(*b = *decimalpoint++)
++b;
#else
*b++ = '.';
#endif
}
}
if (b + decpt > be) {
ret0:
b = 0;
goto ret;
}
for(; decpt > 0; decpt--)
*b++ = '0';
*b = 0;
}
ret:
freedtoa(s0);
return b;
}
/****************************************************************
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);
}
...@@ -246,7 +246,7 @@ The following mathematical functions are available: ...@@ -246,7 +246,7 @@ The following mathematical functions are available:
@menu @menu
* @code{strtoflt128}: strtoflt128, Convert from string * @code{strtoflt128}: strtoflt128, Convert from string
* @code{quadmath_flt128tostr}: quadmath_flt128tostr, Convert to string * @code{quadmath_snprintf}: quadmath_snprintf, Convert to string
@end menu @end menu
...@@ -289,43 +289,73 @@ int main () ...@@ -289,43 +289,73 @@ int main ()
@end table @end table
@node quadmath_flt128tostr @node quadmath_snprintf
@section @code{quadmath_flt128tostr} --- Convert to string @section @code{quadmath_snprintf} --- Convert to string
The function @code{quadmath_flt128tostr} converts a @code{__float128} floating-point The function @code{quadmath_snprintf} converts a @code{__float128} floating-point
number into a string. number into a string. It is a specialized alternative to @code{snprintf}, where
the format string is restricted to a single conversion specifier with @code{Q}
modifier and conversion specifier @code{e}, @code{E}, @code{f}, @code{F}, @code{g},
@code{G}, @code{a} or @code{A}, with no extra characters before or after the
conversion specifier. The @code{%m$} or @code{*m$} style must not be used in
the format.
@table @asis @table @asis
@item Syntax @item Syntax
@code{void quadmath_flt128tostr (char *s, size_t size, size_t n, __float128 x)} @code{int quadmath_snprintf (char *s, size_t size, const char *format, ...)}
@item @emph{Arguments}: @item @emph{Arguments}:
@multitable @columnfractions .15 .70 @multitable @columnfractions .15 .70
@item @var{s} @tab output string @item @var{s} @tab output string
@item @var{size} @tab byte size of the string, including tailing NUL @item @var{size} @tab byte size of the string, including tailing NUL
@item @var{n} @tab number of digits after the decimal point @item @var{format} @tab conversion specifier string
@item @var{x} @tab the number to be converted
@end multitable @end multitable
@item Example @item Example
@smallexample @smallexample
#include <quadmath.h> #include <quadmath.h>
#include <stdlib.h>
#include <stdio.h>
int main () int main ()
@{ @{
__float128 r; __float128 r;
char str[200]; int prec = 20;
int width = 46;
char buf[128];
r = 2.0q; r = 2.0q;
r = sqrtq(r); r = sqrtq (r);
quadmath_flt128tostr (str, sizeof (str), 20, r); int n = quadmath_snprintf (buf, sizeof buf, "%+-#*.20Qe", width, r);
printf("%s\n", str); if ((size_t) n < sizeof buf)
/* Prints: +1.41421356237309504880e+00 */ printf ("%s\n", buf);
/* Prints: +1.41421356237309504880e+00 */
quadmath_snprintf (buf, sizeof buf, "%Qa", r);
if ((size_t) n < sizeof buf)
printf ("%s\n", buf);
/* Prints: 0x1.6a09e667f3bcc908b2fb1366ea96p+0 */
n = quadmath_snprintf (NULL, 0, "%+-#46.*Qe", prec, r);
if (n > -1)
@{
char *str = malloc (n + 1);
if (str)
@{
quadmath_snprintf (str, n + 1, "%+-#46.*Qe", prec, r);
printf ("%s\n", str);
/* Prints: +1.41421356237309504880e+00 */
@}
free (str);
@}
return 0; return 0;
@} @}
@end smallexample @end smallexample
@end table @end table
On some targets when supported by the C library hooks are installed
for @code{printf} family of functions, so that @code{printf ("%Qe", 1.2Q);}
etc.@: works too.
@c --------------------------------------------------------------------- @c ---------------------------------------------------------------------
@c GNU Free Documentation License @c GNU Free Documentation License
......
/* Copyright (C) 2000, 2004, 2008 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@gnu.org>, 2000.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
/* Look up the value of the next multibyte character and return its numerical
value if it is one of the digits known in the locale. If *DECIDED is
-1 this means it is not yet decided which form it is and we have to
search through all available digits. Otherwise we know which script
the digits are from. */
static inline char *
outdigit_value (char *s, int n)
{
const char *outdigit;
size_t dlen;
assert (0 <= n && n <= 9);
outdigit = nl_langinfo (_NL_CTYPE_OUTDIGIT0_MB + n);
dlen = strlen (outdigit);
s -= dlen;
while (dlen-- > 0)
s[dlen] = outdigit[dlen];
return s;
}
/* Look up the value of the next multibyte character and return its numerical
value if it is one of the digits known in the locale. If *DECIDED is
-1 this means it is not yet decided which form it is and we have to
search through all available digits. Otherwise we know which script
the digits are from. */
static inline wchar_t
outdigitwc_value (int n)
{
assert (0 <= n && n <= 9);
return nl_langinfo_wc (_NL_CTYPE_OUTDIGIT0_WC + n);
}
static char *
_i18n_number_rewrite (char *w, char *rear_ptr, char *end)
{
char decimal[MB_LEN_MAX + 1];
char thousands[MB_LEN_MAX + 1];
/* "to_outpunct" is a map from ASCII decimal point and thousands-sep
to their equivalent in locale. This is defined for locales which
use extra decimal point and thousands-sep. */
wctrans_t map = wctrans ("to_outpunct");
wint_t wdecimal = towctrans (L_('.'), map);
wint_t wthousands = towctrans (L_(','), map);
if (__builtin_expect (map != NULL, 0))
{
mbstate_t state;
memset (&state, '\0', sizeof (state));
size_t n = wcrtomb (decimal, wdecimal, &state);
if (n == (size_t) -1)
memcpy (decimal, ".", 2);
else
decimal[n] = '\0';
memset (&state, '\0', sizeof (state));
n = wcrtomb (thousands, wthousands, &state);
if (n == (size_t) -1)
memcpy (thousands, ",", 2);
else
thousands[n] = '\0';
}
/* Copy existing string so that nothing gets overwritten. */
char *src;
int use_alloca = (rear_ptr - w) < 4096;
if (__builtin_expect (use_alloca, 1))
src = (char *) alloca (rear_ptr - w);
else
{
src = (char *) malloc (rear_ptr - w);
if (src == NULL)
/* If we cannot allocate the memory don't rewrite the string.
It is better than nothing. */
return w;
}
memcpy (src, w, rear_ptr - w);
char *s = src + (rear_ptr - w);
w = end;
/* Process all characters in the string. */
while (--s >= src)
{
if (*s >= '0' && *s <= '9')
{
if (sizeof (char) == 1)
w = (char *) outdigit_value ((char *) w, *s - '0');
else
*--w = (char) outdigitwc_value (*s - '0');
}
else if (__builtin_expect (map == NULL, 1) || (*s != '.' && *s != ','))
*--w = *s;
else
{
if (sizeof (char) == 1)
{
const char *outpunct = *s == '.' ? decimal : thousands;
size_t dlen = strlen (outpunct);
w -= dlen;
while (dlen-- > 0)
w[dlen] = outpunct[dlen];
}
else
*--w = *s == '.' ? (char) wdecimal : (char) wthousands;
}
}
if (! use_alloca)
free (src);
return w;
}
/* Internal function for converting integers to ASCII.
Copyright (C) 1994-1999,2002,2003,2007 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
#ifndef _ITOA_H
#define _ITOA_H
/* Convert VALUE into ASCII in base BASE (2..16).
Write backwards starting the character just before BUFLIM.
Return the address of the first (left-to-right) character in the number.
Use upper case letters iff UPPER_CASE is nonzero. */
static const char _itoa_lower_digits[16] = "0123456789abcdef";
static const char _itoa_upper_digits[16] = "0123456789ABCDEF";
static inline char * __attribute__ ((unused, always_inline))
_itoa_word (unsigned long value, char *buflim,
unsigned int base, int upper_case)
{
const char *digits = (upper_case ? _itoa_upper_digits : _itoa_lower_digits);
switch (base)
{
# define SPECIAL(Base) \
case Base: \
do \
*--buflim = digits[value % Base]; \
while ((value /= Base) != 0); \
break
SPECIAL (10);
SPECIAL (16);
SPECIAL (8);
default:
do
*--buflim = digits[value % base];
while ((value /= base) != 0);
}
return buflim;
}
static inline char * __attribute__ ((unused, always_inline))
_itoa (uint64_t value, char *buflim,
unsigned int base, int upper_case)
{
const char *digits = (upper_case ? _itoa_upper_digits : _itoa_lower_digits);
switch (base)
{
SPECIAL (10);
SPECIAL (16);
SPECIAL (8);
default:
do
*--buflim = digits[value % base];
while ((value /= base) != 0);
}
return buflim;
}
# undef SPECIAL
#endif /* itoa.h */
/* Internal function for converting integers to ASCII.
Copyright (C) 1994,95,96,97,98,99,2002,2003 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
#ifndef _ITOWA_H
#define _ITOWA_H 1
/* Convert VALUE into ASCII in base BASE (2..16).
Write backwards starting the character just before BUFLIM.
Return the address of the first (left-to-right) character in the number.
Use upper case letters iff UPPER_CASE is nonzero. */
static const wchar_t _itowa_lower_digits[16] = L_("0123456789abcdef");
static const wchar_t _itowa_upper_digits[16] = L_("0123456789ABCDEF");
static inline wchar_t *
__attribute__ ((unused, always_inline))
_itowa_word (unsigned long value, wchar_t *buflim,
unsigned int base, int upper_case)
{
const wchar_t *digits = (upper_case
? _itowa_upper_digits : _itowa_lower_digits);
wchar_t *bp = buflim;
switch (base)
{
#define SPECIAL(Base) \
case Base: \
do \
*--bp = digits[value % Base]; \
while ((value /= Base) != 0); \
break
SPECIAL (10);
SPECIAL (16);
SPECIAL (8);
default:
do
*--bp = digits[value % base];
while ((value /= base) != 0);
}
return bp;
}
static inline wchar_t *
__attribute__ ((unused, always_inline))
_itowa (uint64_t value, wchar_t *buflim,
unsigned int base, int upper_case)
{
const wchar_t *digits = (upper_case
? _itowa_upper_digits : _itowa_lower_digits);
wchar_t *bp = buflim;
switch (base)
{
SPECIAL (10);
SPECIAL (16);
SPECIAL (8);
default:
do
*--bp = digits[value % base];
while ((value /= base) != 0);
}
return bp;
}
#undef SPECIAL
#endif /* itowa.h */
/* mpn_add_n -- Add two limb vectors of equal, non-zero length.
Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
mp_limb_t
#if __STDC__
mpn_add_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
#else
mpn_add_n (res_ptr, s1_ptr, s2_ptr, size)
register mp_ptr res_ptr;
register mp_srcptr s1_ptr;
register mp_srcptr s2_ptr;
mp_size_t size;
#endif
{
register mp_limb_t x, y, cy;
register mp_size_t j;
/* The loop counter and index J goes from -SIZE to -1. This way
the loop becomes faster. */
j = -size;
/* Offset the base pointers to compensate for the negative indices. */
s1_ptr -= j;
s2_ptr -= j;
res_ptr -= j;
cy = 0;
do
{
y = s2_ptr[j];
x = s1_ptr[j];
y += cy; /* add previous carry to one addend */
cy = (y < cy); /* get out carry from that addition */
y = x + y; /* add other addend */
cy = (y < x) + cy; /* get out carry from that add, combine */
res_ptr[j] = y;
}
while (++j != 0);
return cy;
}
/* mpn_addmul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
by S2_LIMB, add the S1_SIZE least significant limbs of the product to the
limb vector pointed to by RES_PTR. Return the most significant limb of
the product, adjusted for carry-out from the addition.
Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
mp_limb_t
mpn_addmul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
register mp_ptr res_ptr;
register mp_srcptr s1_ptr;
mp_size_t s1_size;
register mp_limb_t s2_limb;
{
register mp_limb_t cy_limb;
register mp_size_t j;
register mp_limb_t prod_high, prod_low;
register mp_limb_t x;
/* The loop counter and index J goes from -SIZE to -1. This way
the loop becomes faster. */
j = -s1_size;
/* Offset the base pointers to compensate for the negative indices. */
res_ptr -= j;
s1_ptr -= j;
cy_limb = 0;
do
{
umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
prod_low += cy_limb;
cy_limb = (prod_low < cy_limb) + prod_high;
x = res_ptr[j];
prod_low = x + prod_low;
cy_limb += (prod_low < x);
res_ptr[j] = prod_low;
}
while (++j != 0);
return cy_limb;
}
/* mpn_cmp -- Compare two low-level natural-number integers.
Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Compare OP1_PTR/OP1_SIZE with OP2_PTR/OP2_SIZE.
There are no restrictions on the relative sizes of
the two arguments.
Return 1 if OP1 > OP2, 0 if they are equal, and -1 if OP1 < OP2. */
int
#if __STDC__
mpn_cmp (mp_srcptr op1_ptr, mp_srcptr op2_ptr, mp_size_t size)
#else
mpn_cmp (op1_ptr, op2_ptr, size)
mp_srcptr op1_ptr;
mp_srcptr op2_ptr;
mp_size_t size;
#endif
{
mp_size_t i;
mp_limb_t op1_word, op2_word;
for (i = size - 1; i >= 0; i--)
{
op1_word = op1_ptr[i];
op2_word = op2_ptr[i];
if (op1_word != op2_word)
goto diff;
}
return 0;
diff:
/* This can *not* be simplified to
op2_word - op2_word
since that expression might give signed overflow. */
return (op1_word > op2_word) ? 1 : -1;
}
/* mpn_divrem -- Divide natural numbers, producing both remainder and
quotient.
Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
the NSIZE-DSIZE least significant quotient limbs at QP
and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
non-zero, generate that many fraction bits and append them after the
other quotient limbs.
Return the most significant limb of the quotient, this is always 0 or 1.
Preconditions:
0. NSIZE >= DSIZE.
1. The most significant bit of the divisor must be set.
2. QP must either not overlap with the input operands at all, or
QP + DSIZE >= NP must hold true. (This means that it's
possible to put the quotient in the high part of NUM, right after the
remainder in NUM.
3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. */
mp_limb_t
#if __STDC__
mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
mp_ptr np, mp_size_t nsize,
mp_srcptr dp, mp_size_t dsize)
#else
mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
mp_ptr qp;
mp_size_t qextra_limbs;
mp_ptr np;
mp_size_t nsize;
mp_srcptr dp;
mp_size_t dsize;
#endif
{
mp_limb_t most_significant_q_limb = 0;
switch (dsize)
{
case 0:
/* We are asked to divide by zero, so go ahead and do it! (To make
the compiler not remove this statement, return the value.) */
return 1 / dsize;
case 1:
{
mp_size_t i;
mp_limb_t n1;
mp_limb_t d;
d = dp[0];
n1 = np[nsize - 1];
if (n1 >= d)
{
n1 -= d;
most_significant_q_limb = 1;
}
qp += qextra_limbs;
for (i = nsize - 2; i >= 0; i--)
udiv_qrnnd (qp[i], n1, n1, np[i], d);
qp -= qextra_limbs;
for (i = qextra_limbs - 1; i >= 0; i--)
udiv_qrnnd (qp[i], n1, n1, 0, d);
np[0] = n1;
}
break;
case 2:
{
mp_size_t i;
mp_limb_t n1, n0, n2;
mp_limb_t d1, d0;
np += nsize - 2;
d1 = dp[1];
d0 = dp[0];
n1 = np[1];
n0 = np[0];
if (n1 >= d1 && (n1 > d1 || n0 >= d0))
{
sub_ddmmss (n1, n0, n1, n0, d1, d0);
most_significant_q_limb = 1;
}
for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
{
mp_limb_t q;
mp_limb_t r;
if (i >= qextra_limbs)
np--;
else
np[0] = 0;
if (n1 == d1)
{
/* Q should be either 111..111 or 111..110. Need special
treatment of this rare case as normal division would
give overflow. */
q = ~(mp_limb_t) 0;
r = n0 + d1;
if (r < d1) /* Carry in the addition? */
{
add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
qp[i] = q;
continue;
}
n1 = d0 - (d0 != 0);
n0 = -d0;
}
else
{
udiv_qrnnd (q, r, n1, n0, d1);
umul_ppmm (n1, n0, d0, q);
}
n2 = np[0];
q_test:
if (n1 > r || (n1 == r && n0 > n2))
{
/* The estimated Q was too large. */
q--;
sub_ddmmss (n1, n0, n1, n0, 0, d0);
r += d1;
if (r >= d1) /* If not carry, test Q again. */
goto q_test;
}
qp[i] = q;
sub_ddmmss (n1, n0, r, n2, n1, n0);
}
np[1] = n1;
np[0] = n0;
}
break;
default:
{
mp_size_t i;
mp_limb_t dX, d1, n0;
np += nsize - dsize;
dX = dp[dsize - 1];
d1 = dp[dsize - 2];
n0 = np[dsize - 1];
if (n0 >= dX)
{
if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
{
mpn_sub_n (np, np, dp, dsize);
n0 = np[dsize - 1];
most_significant_q_limb = 1;
}
}
for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
{
mp_limb_t q;
mp_limb_t n1, n2;
mp_limb_t cy_limb;
if (i >= qextra_limbs)
{
np--;
n2 = np[dsize];
}
else
{
n2 = np[dsize - 1];
MPN_COPY_DECR (np + 1, np, dsize);
np[0] = 0;
}
if (n0 == dX)
/* This might over-estimate q, but it's probably not worth
the extra code here to find out. */
q = ~(mp_limb_t) 0;
else
{
mp_limb_t r;
udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
umul_ppmm (n1, n0, d1, q);
while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
{
q--;
r += dX;
if (r < dX) /* I.e. "carry in previous addition?" */
break;
n1 -= n0 < d1;
n0 -= d1;
}
}
/* Possible optimization: We already have (q * n0) and (1 * n1)
after the calculation of q. Taking advantage of that, we
could make this loop make two iterations less. */
cy_limb = mpn_submul_1 (np, dp, dsize, q);
if (n2 != cy_limb)
{
mpn_add_n (np, np, dp, dsize);
q--;
}
qp[i] = q;
n0 = np[dsize - 1];
}
}
}
return most_significant_q_limb;
}
/* Copyright (C) 1995,1996,1997,1998,1999,2002,2003
Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
#include <float.h>
#include <math.h>
#include <stdlib.h>
#include "gmp-impl.h"
/* Convert a `__float128' in IEEE854 quad-precision format to a
multi-precision integer representing the significand scaled up by its
number of bits (113 for long double) and an integral power of two
(MPN frexpl). */
mp_size_t
mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size,
int *expt, int *is_neg,
__float128 value)
{
ieee854_float128 u;
u.value = value;
*is_neg = u.ieee.negative;
*expt = (int) u.ieee.exponent - IEEE854_FLOAT128_BIAS;
#if BITS_PER_MP_LIMB == 32
res_ptr[0] = u.ieee.mant_low; /* Low-order 32 bits of fraction. */
res_ptr[1] = (u.ieee.mant_low >> 32);
res_ptr[2] = u.ieee.mant_high;
res_ptr[3] = (u.ieee.mant_high >> 32); /* High-order 32 bits. */
#define N 4
#elif BITS_PER_MP_LIMB == 64
res_ptr[0] = u.ieee.mant_low;
res_ptr[1] = u.ieee.mant_high;
#define N 2
#else
#error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
#endif
/* The format does not fill the last limb. There are some zeros. */
#define NUM_LEADING_ZEROS (BITS_PER_MP_LIMB \
- (FLT128_MANT_DIG - ((N - 1) * BITS_PER_MP_LIMB)))
if (u.ieee.exponent == 0)
{
/* A biased exponent of zero is a special case.
Either it is a zero or it is a denormal number. */
if (res_ptr[0] == 0 && res_ptr[1] == 0
&& res_ptr[N - 2] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=4. */
/* It's zero. */
*expt = 0;
else
{
/* It is a denormal number, meaning it has no implicit leading
one bit, and its exponent is in fact the format minimum. */
int cnt;
#if N == 2
if (res_ptr[N - 1] != 0)
{
count_leading_zeros (cnt, res_ptr[N - 1]);
cnt -= NUM_LEADING_ZEROS;
res_ptr[N - 1] = res_ptr[N - 1] << cnt
| (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
res_ptr[0] <<= cnt;
*expt = FLT128_MIN_EXP - 1 - cnt;
}
else
{
count_leading_zeros (cnt, res_ptr[0]);
if (cnt >= NUM_LEADING_ZEROS)
{
res_ptr[N - 1] = res_ptr[0] << (cnt - NUM_LEADING_ZEROS);
res_ptr[0] = 0;
}
else
{
res_ptr[N - 1] = res_ptr[0] >> (NUM_LEADING_ZEROS - cnt);
res_ptr[0] <<= BITS_PER_MP_LIMB - (NUM_LEADING_ZEROS - cnt);
}
*expt = FLT128_MIN_EXP - 1
- (BITS_PER_MP_LIMB - NUM_LEADING_ZEROS) - cnt;
}
#else
int j, k, l;
for (j = N - 1; j > 0; j--)
if (res_ptr[j] != 0)
break;
count_leading_zeros (cnt, res_ptr[j]);
cnt -= NUM_LEADING_ZEROS;
l = N - 1 - j;
if (cnt < 0)
{
cnt += BITS_PER_MP_LIMB;
l--;
}
if (!cnt)
for (k = N - 1; k >= l; k--)
res_ptr[k] = res_ptr[k-l];
else
{
for (k = N - 1; k > l; k--)
res_ptr[k] = res_ptr[k-l] << cnt
| res_ptr[k-l-1] >> (BITS_PER_MP_LIMB - cnt);
res_ptr[k--] = res_ptr[0] << cnt;
}
for (; k >= 0; k--)
res_ptr[k] = 0;
*expt = FLT128_MIN_EXP - 1 - l * BITS_PER_MP_LIMB - cnt;
#endif
}
}
else
/* Add the implicit leading one bit for a normalized number. */
res_ptr[N - 1] |= (mp_limb_t) 1 << (FLT128_MANT_DIG - 1
- ((N - 1) * BITS_PER_MP_LIMB));
return N;
}
/* Header file for constants used in floating point <-> decimal conversions.
Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003
Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
#ifndef _FPIOCONST_H
#define _FPIOCONST_H
#include <float.h>
#include <math.h>
/* These values are used by __printf_fp, where they are noncritical (if the
value is not large enough, it will just be slower); and by
strtof/strtod/strtold, where it is critical (it's used for overflow
detection).
XXX These should be defined in <float.h>. For the time being, we have the
IEEE754 values here. */
#define FLT128_MAX_10_EXP_LOG 12 /* = floor(log_2(FLT128_MAX_10_EXP)) */
/* The array with the number representation. */
#define __tens __quadmath_tens
extern const mp_limb_t __tens[] attribute_hidden;
/* Table of powers of ten. This is used by __printf_fp and by
strtof/strtod/strtold. */
struct mp_power
{
size_t arrayoff; /* Offset in `__tens'. */
mp_size_t arraysize; /* Size of the array. */
int p_expo; /* Exponent of the number 10^(2^i). */
int m_expo; /* Exponent of the number 10^-(2^i-1). */
};
#define _fpioconst_pow10 __quadmath_fpioconst_pow10
extern const struct mp_power _fpioconst_pow10[FLT128_MAX_10_EXP_LOG + 1]
attribute_hidden;
/* The constants in the array `_fpioconst_pow10' have an offset. */
#if BITS_PER_MP_LIMB == 32
# define _FPIO_CONST_OFFSET 2
#else
# define _FPIO_CONST_OFFSET 1
#endif
#endif /* fpioconst.h */
/* Include file for internal GNU MP types and definitions.
Copyright (C) 1991, 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include <stdlib.h>
#include "quadmath-imp.h"
#undef alloca
#define alloca __builtin_alloca
#define ABS(x) (x >= 0 ? x : -x)
#ifndef MIN
#define MIN(l,o) ((l) < (o) ? (l) : (o))
#endif
#ifndef MAX
#define MAX(h,i) ((h) > (i) ? (h) : (i))
#endif
#define BITS_PER_MP_LIMB (__SIZEOF_LONG__ * __CHAR_BIT__)
#define BYTES_PER_MP_LIMB (BITS_PER_MP_LIMB / __CHAR_BIT__)
typedef unsigned long int mp_limb_t;
typedef long int mp_limb_signed_t;
typedef mp_limb_t * mp_ptr;
typedef const mp_limb_t * mp_srcptr;
typedef long int mp_size_t;
typedef long int mp_exp_t;
/* Define stuff for longlong.h. */
typedef unsigned int UQItype __attribute__ ((mode (QI)));
typedef int SItype __attribute__ ((mode (SI)));
typedef unsigned int USItype __attribute__ ((mode (SI)));
typedef int DItype __attribute__ ((mode (DI)));
typedef unsigned int UDItype __attribute__ ((mode (DI)));
typedef mp_limb_t UWtype;
typedef unsigned int UHWtype;
#define W_TYPE_SIZE BITS_PER_MP_LIMB
#ifdef HAVE_HIDDEN_VISIBILITY
#define attribute_hidden __attribute__((__visibility__ ("hidden")))
#else
#define attribute_hidden
#endif
#include "../../gcc/longlong.h"
/* Copy NLIMBS *limbs* from SRC to DST. */
#define MPN_COPY_INCR(DST, SRC, NLIMBS) \
do { \
mp_size_t __i; \
for (__i = 0; __i < (NLIMBS); __i++) \
(DST)[__i] = (SRC)[__i]; \
} while (0)
#define MPN_COPY_DECR(DST, SRC, NLIMBS) \
do { \
mp_size_t __i; \
for (__i = (NLIMBS) - 1; __i >= 0; __i--) \
(DST)[__i] = (SRC)[__i]; \
} while (0)
#define MPN_COPY MPN_COPY_INCR
/* Zero NLIMBS *limbs* AT DST. */
#define MPN_ZERO(DST, NLIMBS) \
do { \
mp_size_t __i; \
for (__i = 0; __i < (NLIMBS); __i++) \
(DST)[__i] = 0; \
} while (0)
#define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
do { \
if ((size) < KARATSUBA_THRESHOLD) \
impn_mul_n_basecase (prodp, up, vp, size); \
else \
impn_mul_n (prodp, up, vp, size, tspace); \
} while (0);
#define __MPN(x) __quadmath_mpn_##x
/* Internal mpn calls */
#define impn_mul_n_basecase __MPN(impn_mul_n_basecase)
#define impn_mul_n __MPN(impn_mul_n)
/* Prototypes for internal mpn calls. */
void impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp,
mp_size_t size) attribute_hidden;
void impn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size,
mp_ptr tspace) attribute_hidden;
#define mpn_add_n __MPN(add_n)
#define mpn_addmul_1 __MPN(addmul_1)
#define mpn_cmp __MPN(cmp)
#define mpn_divrem __MPN(divrem)
#define mpn_lshift __MPN(lshift)
#define mpn_mul __MPN(mul)
#define mpn_mul_1 __MPN(mul_1)
#define mpn_rshift __MPN(rshift)
#define mpn_sub_n __MPN(sub_n)
#define mpn_submul_1 __MPN(submul_1)
mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t)
attribute_hidden;
mp_limb_t mpn_addmul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
attribute_hidden;
int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t) attribute_hidden;
mp_limb_t mpn_divrem (mp_ptr, mp_size_t, mp_ptr, mp_size_t, mp_srcptr,
mp_size_t) attribute_hidden;
mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int)
attribute_hidden;
mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t)
attribute_hidden;
mp_limb_t mpn_mul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
attribute_hidden;
mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int)
attribute_hidden;
mp_limb_t mpn_sub_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t)
attribute_hidden;
mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t)
attribute_hidden;
#define mpn_extract_flt128 __MPN(extract_flt128)
mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, int *expt,
int *is_neg, __float128 value) attribute_hidden;
#define mpn_divmod(qp,np,nsize,dp,dsize) mpn_divrem (qp,0,np,nsize,dp,dsize)
static inline mp_limb_t
mpn_add_1 (register mp_ptr res_ptr,
register mp_srcptr s1_ptr,
register mp_size_t s1_size,
register mp_limb_t s2_limb)
{
register mp_limb_t x;
x = *s1_ptr++;
s2_limb = x + s2_limb;
*res_ptr++ = s2_limb;
if (s2_limb < x)
{
while (--s1_size != 0)
{
x = *s1_ptr++ + 1;
*res_ptr++ = x;
if (x != 0)
goto fin;
}
return 1;
}
fin:
if (res_ptr != s1_ptr)
{
mp_size_t i;
for (i = 0; i < s1_size - 1; i++)
res_ptr[i] = s1_ptr[i];
}
return 0;
}
/* mpn_lshift -- Shift left low level.
Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the left
and store the USIZE least significant digits of the result at WP.
Return the bits shifted out from the most significant digit.
Argument constraints:
1. 0 < CNT < BITS_PER_MP_LIMB
2. If the result is to be written over the input, WP must be >= UP.
*/
mp_limb_t
#if __STDC__
mpn_lshift (register mp_ptr wp,
register mp_srcptr up, mp_size_t usize,
register unsigned int cnt)
#else
mpn_lshift (wp, up, usize, cnt)
register mp_ptr wp;
register mp_srcptr up;
mp_size_t usize;
register unsigned int cnt;
#endif
{
register mp_limb_t high_limb, low_limb;
register unsigned sh_1, sh_2;
register mp_size_t i;
mp_limb_t retval;
#ifdef DEBUG
if (usize == 0 || cnt == 0)
abort ();
#endif
sh_1 = cnt;
#if 0
if (sh_1 == 0)
{
if (wp != up)
{
/* Copy from high end to low end, to allow specified input/output
overlapping. */
for (i = usize - 1; i >= 0; i--)
wp[i] = up[i];
}
return 0;
}
#endif
wp += 1;
sh_2 = BITS_PER_MP_LIMB - sh_1;
i = usize - 1;
low_limb = up[i];
retval = low_limb >> sh_2;
high_limb = low_limb;
while (--i >= 0)
{
low_limb = up[i];
wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
high_limb = low_limb;
}
wp[i] = high_limb << sh_1;
return retval;
}
/* mpn_mul -- Multiply two natural numbers.
Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
and v (pointed to by VP, with VSIZE limbs), and store the result at
PRODP. USIZE + VSIZE limbs are always stored, but if the input
operands are normalized. Return the most significant limb of the
result.
NOTE: The space pointed to by PRODP is overwritten before finished
with U and V, so overlap is an error.
Argument constraints:
1. USIZE >= VSIZE.
2. PRODP != UP and PRODP != VP, i.e. the destination
must be distinct from the multiplier and the multiplicand. */
/* If KARATSUBA_THRESHOLD is not already defined, define it to a
value which is good on most machines. */
#ifndef KARATSUBA_THRESHOLD
#define KARATSUBA_THRESHOLD 32
#endif
mp_limb_t
#if __STDC__
mpn_mul (mp_ptr prodp,
mp_srcptr up, mp_size_t usize,
mp_srcptr vp, mp_size_t vsize)
#else
mpn_mul (prodp, up, usize, vp, vsize)
mp_ptr prodp;
mp_srcptr up;
mp_size_t usize;
mp_srcptr vp;
mp_size_t vsize;
#endif
{
mp_ptr prod_endp = prodp + usize + vsize - 1;
mp_limb_t cy;
mp_ptr tspace;
if (vsize < KARATSUBA_THRESHOLD)
{
/* Handle simple cases with traditional multiplication.
This is the most critical code of the entire function. All
multiplies rely on this, both small and huge. Small ones arrive
here immediately. Huge ones arrive here as this is the base case
for Karatsuba's recursive algorithm below. */
mp_size_t i;
mp_limb_t cy_limb;
mp_limb_t v_limb;
if (vsize == 0)
return 0;
/* Multiply by the first limb in V separately, as the result can be
stored (not added) to PROD. We also avoid a loop for zeroing. */
v_limb = vp[0];
if (v_limb <= 1)
{
if (v_limb == 1)
MPN_COPY (prodp, up, usize);
else
MPN_ZERO (prodp, usize);
cy_limb = 0;
}
else
cy_limb = mpn_mul_1 (prodp, up, usize, v_limb);
prodp[usize] = cy_limb;
prodp++;
/* For each iteration in the outer loop, multiply one limb from
U with one limb from V, and add it to PROD. */
for (i = 1; i < vsize; i++)
{
v_limb = vp[i];
if (v_limb <= 1)
{
cy_limb = 0;
if (v_limb == 1)
cy_limb = mpn_add_n (prodp, prodp, up, usize);
}
else
cy_limb = mpn_addmul_1 (prodp, up, usize, v_limb);
prodp[usize] = cy_limb;
prodp++;
}
return cy_limb;
}
tspace = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
MPN_MUL_N_RECURSE (prodp, up, vp, vsize, tspace);
prodp += vsize;
up += vsize;
usize -= vsize;
if (usize >= vsize)
{
mp_ptr tp = (mp_ptr) alloca (2 * vsize * BYTES_PER_MP_LIMB);
do
{
MPN_MUL_N_RECURSE (tp, up, vp, vsize, tspace);
cy = mpn_add_n (prodp, prodp, tp, vsize);
mpn_add_1 (prodp + vsize, tp + vsize, vsize, cy);
prodp += vsize;
up += vsize;
usize -= vsize;
}
while (usize >= vsize);
}
/* True: usize < vsize. */
/* Make life simple: Recurse. */
if (usize != 0)
{
mpn_mul (tspace, vp, vsize, up, usize);
cy = mpn_add_n (prodp, prodp, tspace, vsize);
mpn_add_1 (prodp + vsize, tspace + vsize, usize, cy);
}
return *prod_endp;
}
/* mpn_mul_1 -- Multiply a limb vector with a single limb and
store the product in a second limb vector.
Copyright (C) 1991, 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
mp_limb_t
mpn_mul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
register mp_ptr res_ptr;
register mp_srcptr s1_ptr;
mp_size_t s1_size;
register mp_limb_t s2_limb;
{
register mp_limb_t cy_limb;
register mp_size_t j;
register mp_limb_t prod_high, prod_low;
/* The loop counter and index J goes from -S1_SIZE to -1. This way
the loop becomes faster. */
j = -s1_size;
/* Offset the base pointers to compensate for the negative indices. */
s1_ptr -= j;
res_ptr -= j;
cy_limb = 0;
do
{
umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
prod_low += cy_limb;
cy_limb = (prod_low < cy_limb) + prod_high;
res_ptr[j] = prod_low;
}
while (++j != 0);
return cy_limb;
}
/* mpn_mul_n -- Multiply two natural numbers of length n.
Copyright (C) 1991, 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
both with SIZE limbs, and store the result at PRODP. 2 * SIZE limbs are
always stored. Return the most significant limb.
Argument constraints:
1. PRODP != UP and PRODP != VP, i.e. the destination
must be distinct from the multiplier and the multiplicand. */
/* If KARATSUBA_THRESHOLD is not already defined, define it to a
value which is good on most machines. */
#ifndef KARATSUBA_THRESHOLD
#define KARATSUBA_THRESHOLD 32
#endif
/* The code can't handle KARATSUBA_THRESHOLD smaller than 2. */
#if KARATSUBA_THRESHOLD < 2
#undef KARATSUBA_THRESHOLD
#define KARATSUBA_THRESHOLD 2
#endif
/* Handle simple cases with traditional multiplication.
This is the most critical code of multiplication. All multiplies rely
on this, both small and huge. Small ones arrive here immediately. Huge
ones arrive here as this is the base case for Karatsuba's recursive
algorithm below. */
void
#if __STDC__
impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
#else
impn_mul_n_basecase (prodp, up, vp, size)
mp_ptr prodp;
mp_srcptr up;
mp_srcptr vp;
mp_size_t size;
#endif
{
mp_size_t i;
mp_limb_t cy_limb;
mp_limb_t v_limb;
/* Multiply by the first limb in V separately, as the result can be
stored (not added) to PROD. We also avoid a loop for zeroing. */
v_limb = vp[0];
if (v_limb <= 1)
{
if (v_limb == 1)
MPN_COPY (prodp, up, size);
else
MPN_ZERO (prodp, size);
cy_limb = 0;
}
else
cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
prodp[size] = cy_limb;
prodp++;
/* For each iteration in the outer loop, multiply one limb from
U with one limb from V, and add it to PROD. */
for (i = 1; i < size; i++)
{
v_limb = vp[i];
if (v_limb <= 1)
{
cy_limb = 0;
if (v_limb == 1)
cy_limb = mpn_add_n (prodp, prodp, up, size);
}
else
cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
prodp[size] = cy_limb;
prodp++;
}
}
void
#if __STDC__
impn_mul_n (mp_ptr prodp,
mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
#else
impn_mul_n (prodp, up, vp, size, tspace)
mp_ptr prodp;
mp_srcptr up;
mp_srcptr vp;
mp_size_t size;
mp_ptr tspace;
#endif
{
if ((size & 1) != 0)
{
/* The size is odd, the code code below doesn't handle that.
Multiply the least significant (size - 1) limbs with a recursive
call, and handle the most significant limb of S1 and S2
separately. */
/* A slightly faster way to do this would be to make the Karatsuba
code below behave as if the size were even, and let it check for
odd size in the end. I.e., in essence move this code to the end.
Doing so would save us a recursive call, and potentially make the
stack grow a lot less. */
mp_size_t esize = size - 1; /* even size */
mp_limb_t cy_limb;
MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
cy_limb = mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
prodp[esize + esize] = cy_limb;
cy_limb = mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
prodp[esize + size] = cy_limb;
}
else
{
/* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
Split U in two pieces, U1 and U0, such that
U = U0 + U1*(B**n),
and V in V1 and V0, such that
V = V0 + V1*(B**n).
UV is then computed recursively using the identity
2n n n n
UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
1 1 1 0 0 1 0 0
Where B = 2**BITS_PER_MP_LIMB. */
mp_size_t hsize = size >> 1;
mp_limb_t cy;
int negflg;
/*** Product H. ________________ ________________
|_____U1 x V1____||____U0 x V0_____| */
/* Put result in upper part of PROD and pass low part of TSPACE
as new TSPACE. */
MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
/*** Product M. ________________
|_(U1-U0)(V0-V1)_| */
if (mpn_cmp (up + hsize, up, hsize) >= 0)
{
mpn_sub_n (prodp, up + hsize, up, hsize);
negflg = 0;
}
else
{
mpn_sub_n (prodp, up, up + hsize, hsize);
negflg = 1;
}
if (mpn_cmp (vp + hsize, vp, hsize) >= 0)
{
mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
negflg ^= 1;
}
else
{
mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
/* No change of NEGFLG. */
}
/* Read temporary operands from low part of PROD.
Put result in low part of TSPACE using upper part of TSPACE
as new TSPACE. */
MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
/*** Add/copy product H. */
MPN_COPY (prodp + hsize, prodp + size, hsize);
cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
/*** Add product M (if NEGFLG M is a negative number). */
if (negflg)
cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
else
cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
/*** Product L. ________________ ________________
|________________||____U0 x V0_____| */
/* Read temporary operands from low part of PROD.
Put result in low part of TSPACE using upper part of TSPACE
as new TSPACE. */
MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
/*** Add/copy Product L (twice). */
cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
if (cy)
mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
MPN_COPY (prodp, tspace, hsize);
cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
if (cy)
mpn_add_1 (prodp + size, prodp + size, size, 1);
}
}
/* GCC Quad-Precision Math Library
Copyright (C) 2011 Free Software Foundation, Inc.
Written by Jakub Jelinek <jakub@redhat.com>
This file is part of the libquadmath library.
Libquadmath is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
Libquadmath is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with libquadmath; see the file COPYING.LIB. If
not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
Boston, MA 02110-1301, USA. */
#include <stdlib.h>
#include <stdio.h>
#ifdef HAVE_LIMITS_H
#include <limits.h>
#endif
#ifdef HAVE_LANGINFO_H
#include <langinfo.h>
#endif
#ifdef HAVE_CTYPE_H
#include <ctype.h>
#endif
#ifdef HAVE_WCHAR_H
#include <wchar.h>
#endif
#ifdef HAVE_WCTYPE_H
#include <wctype.h>
#endif
#ifdef HAVE_PRINTF_HOOKS
#include <printf.h>
#endif
#include "quadmath-imp.h"
#include "gmp-impl.h"
#ifdef HAVE_WCHAR_H
#define L_(x) L##x
#else
#define L_(x) x
#undef wchar_t
#undef wint_t
#undef putwc
#undef WEOF
#define wchar_t char
#define wint_t int
#define putwc(c,f) putc(c,f)
#define WEOF EOF
#endif
#ifndef HAVE_CTYPE_H
/* Won't work for EBCDIC. */
#undef isupper
#undef isdigit
#undef tolower
#define isupper(x) ((x) >= 'A' && (x) <= 'Z')
#define isdigit(x) ((x) >= '0' && (x) <= '9')
#define tolower(x) (isupper (x) ? (x) - 'A' + 'a' : (x))
#endif
#ifndef CHAR_MAX
#ifdef __CHAR_UNSIGNED__
#define CHAR_MAX (2 * __SCHAR_MAX__ + 1)
#else
#define CHAR_MAX __SCHAR_MAX__
#endif
#endif
#ifndef HAVE_PRINTF_HOOKS
#define printf_info __quadmath_printf_info
struct printf_info
{
int prec; /* Precision. */
int width; /* Width. */
wchar_t spec; /* Format letter. */
unsigned int is_long_double:1;/* L flag. */
unsigned int is_short:1; /* h flag. */
unsigned int is_long:1; /* l flag. */
unsigned int alt:1; /* # flag. */
unsigned int space:1; /* Space flag. */
unsigned int left:1; /* - flag. */
unsigned int showsign:1; /* + flag. */
unsigned int group:1; /* ' flag. */
unsigned int extra:1; /* For special use. */
unsigned int is_char:1; /* hh flag. */
unsigned int wide:1; /* Nonzero for wide character streams. */
unsigned int i18n:1; /* I flag. */
unsigned short int user; /* Bits for user-installed modifiers. */
wchar_t pad; /* Padding character. */
};
#endif
struct __quadmath_printf_file
{
FILE *fp;
char *str;
size_t size;
size_t len;
int file_p;
};
int
__quadmath_printf_fp (struct __quadmath_printf_file *fp,
const struct printf_info *info,
const void *const *args) attribute_hidden;
int
__quadmath_printf_fphex (struct __quadmath_printf_file *fp,
const struct printf_info *info,
const void *const *args) attribute_hidden;
size_t __quadmath_do_pad (struct __quadmath_printf_file *fp, int wide,
int c, size_t n) attribute_hidden;
static inline __attribute__((__unused__)) size_t
__quadmath_do_put (struct __quadmath_printf_file *fp, int wide,
const char *s, size_t n)
{
size_t len;
if (fp->file_p)
{
if (wide)
{
size_t cnt;
const wchar_t *ls = (const wchar_t *) s;
for (cnt = 0; cnt < n; cnt++)
if (putwc (ls[cnt], fp->fp) == WEOF)
break;
return cnt;
}
return fwrite (s, 1, n, fp->fp);
}
len = MIN (fp->size, n);
memcpy (fp->str, s, len);
fp->str += len;
fp->size -= len;
fp->len += n;
return n;
}
static inline __attribute__((__unused__)) int
__quadmath_do_putc (struct __quadmath_printf_file *fp, int wide,
wchar_t c)
{
if (fp->file_p)
return wide ? (int) putwc (c, fp->fp) : putc (c, fp->fp);
if (fp->size)
{
*(fp->str++) = c;
fp->size--;
}
fp->len++;
return (unsigned char) c;
}
#define PUT(f, s, n) __quadmath_do_put (f, wide, s, n)
#define PAD(f, c, n) __quadmath_do_pad (f, wide, c, n)
#define PUTC(c, f) __quadmath_do_putc (f, wide, c)
#define nl_langinfo_wc(x) \
({ union { const char *mb; wchar_t wc; } u; u.mb = nl_langinfo (x); u.wc; })
/* mpn_rshift -- Shift right a low-level natural-number integer.
Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
and store the USIZE least significant limbs of the result at WP.
The bits shifted out to the right are returned.
Argument constraints:
1. 0 < CNT < BITS_PER_MP_LIMB
2. If the result is to be written over the input, WP must be <= UP.
*/
mp_limb_t
#if __STDC__
mpn_rshift (register mp_ptr wp,
register mp_srcptr up, mp_size_t usize,
register unsigned int cnt)
#else
mpn_rshift (wp, up, usize, cnt)
register mp_ptr wp;
register mp_srcptr up;
mp_size_t usize;
register unsigned int cnt;
#endif
{
register mp_limb_t high_limb, low_limb;
register unsigned sh_1, sh_2;
register mp_size_t i;
mp_limb_t retval;
#ifdef DEBUG
if (usize == 0 || cnt == 0)
abort ();
#endif
sh_1 = cnt;
#if 0
if (sh_1 == 0)
{
if (wp != up)
{
/* Copy from low end to high end, to allow specified input/output
overlapping. */
for (i = 0; i < usize; i++)
wp[i] = up[i];
}
return usize;
}
#endif
wp -= 1;
sh_2 = BITS_PER_MP_LIMB - sh_1;
high_limb = up[0];
retval = high_limb << sh_2;
low_limb = high_limb;
for (i = 1; i < usize; i++)
{
high_limb = up[i];
wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
low_limb = high_limb;
}
wp[i] = low_limb >> sh_1;
return retval;
}
/* mpn_sub_n -- Subtract two limb vectors of equal, non-zero length.
Copyright (C) 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
mp_limb_t
#if __STDC__
mpn_sub_n (mp_ptr res_ptr, mp_srcptr s1_ptr, mp_srcptr s2_ptr, mp_size_t size)
#else
mpn_sub_n (res_ptr, s1_ptr, s2_ptr, size)
register mp_ptr res_ptr;
register mp_srcptr s1_ptr;
register mp_srcptr s2_ptr;
mp_size_t size;
#endif
{
register mp_limb_t x, y, cy;
register mp_size_t j;
/* The loop counter and index J goes from -SIZE to -1. This way
the loop becomes faster. */
j = -size;
/* Offset the base pointers to compensate for the negative indices. */
s1_ptr -= j;
s2_ptr -= j;
res_ptr -= j;
cy = 0;
do
{
y = s2_ptr[j];
x = s1_ptr[j];
y += cy; /* add previous carry to subtrahend */
cy = (y < cy); /* get out carry from that addition */
y = x - y; /* main subtract */
cy = (y > x) + cy; /* get out carry from the subtract, combine */
res_ptr[j] = y;
}
while (++j != 0);
return cy;
}
/* mpn_submul_1 -- multiply the S1_SIZE long limb vector pointed to by S1_PTR
by S2_LIMB, subtract the S1_SIZE least significant limbs of the product
from the limb vector pointed to by RES_PTR. Return the most significant
limb of the product, adjusted for carry-out from the subtraction.
Copyright (C) 1992, 1993, 1994, 1996, 2005 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp-impl.h"
mp_limb_t
mpn_submul_1 (res_ptr, s1_ptr, s1_size, s2_limb)
register mp_ptr res_ptr;
register mp_srcptr s1_ptr;
mp_size_t s1_size;
register mp_limb_t s2_limb;
{
register mp_limb_t cy_limb;
register mp_size_t j;
register mp_limb_t prod_high, prod_low;
register mp_limb_t x;
/* The loop counter and index J goes from -SIZE to -1. This way
the loop becomes faster. */
j = -s1_size;
/* Offset the base pointers to compensate for the negative indices. */
res_ptr -= j;
s1_ptr -= j;
cy_limb = 0;
do
{
umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
prod_low += cy_limb;
cy_limb = (prod_low < cy_limb) + prod_high;
x = res_ptr[j];
prod_low = x - prod_low;
cy_limb += (prod_low > x);
res_ptr[j] = prod_low;
}
while (++j != 0);
return cy_limb;
}
...@@ -132,8 +132,8 @@ extern __complex128 ctanhq (__complex128) __quadmath_throw; ...@@ -132,8 +132,8 @@ extern __complex128 ctanhq (__complex128) __quadmath_throw;
/* Prototypes for string <-> __float128 conversion functions */ /* Prototypes for string <-> __float128 conversion functions */
extern __float128 strtoflt128 (const char *, char **) __quadmath_throw; extern __float128 strtoflt128 (const char *, char **) __quadmath_throw;
extern void quadmath_flt128tostr (char *, size_t, size_t, __float128) extern int quadmath_snprintf (char *str, size_t size,
__quadmath_throw; const char *format, ...) __quadmath_throw;
/* Macros */ /* Macros */
......
...@@ -92,7 +92,7 @@ QUADMATH_1.0 { ...@@ -92,7 +92,7 @@ QUADMATH_1.0 {
ctanhq; ctanhq;
strtoflt128; strtoflt128;
quadmath_flt128tostr; quadmath_snprintf;
local: local:
*; *;
}; };
/* GCC Quad-Precision Math Library
Copyright (C) 2010, 2011 Free Software Foundation, Inc.
Written by Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
This file is part of the libquadmath library.
Libquadmath is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
Libquadmath is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with libquadmath; see the file COPYING.LIB. If
not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
Boston, MA 02110-1301, USA. */
#include "quadmath.h"
#include <stdio.h>
#include <string.h>
#define ABS(x) ((x) >= 0 ? (x) : -(x))
static void
format (char * res, const __float128 x, size_t n)
{
char buffer[1024];
char *p;
memset (buffer, 0, sizeof(buffer));
g_Qfmt (buffer, &x, n + 1, sizeof(buffer) - 3);
p = buffer + (*buffer == '-' ? 1 : 0);
/* The sign is the easiest part. */
res[0] = (signbitq (x) ? '-' : '+');
if (*p == '.')
{
/* We have a number smaller than 1, without exponent. */
int exp = 0;
char *c;
for (c = p+1; *c == '0'; c++)
exp++;
/* We move the string "exp" characters left. */
size_t l = strlen (p+1+exp);
memcpy (res + 2, p + 1 + exp, l);
memset (res + 2 + l, '0', n - l + 1);
sprintf (res + n + 3, "e-%02d", exp + 1);
res[1] = res[2];
res[2] = '.';
return;
}
/* Now, do we already have an exponent. */
char *c;
for (c = p; *c && *c != 'e'; c++)
;
if (*c)
{
int exp = strtol (c + 1, NULL, 10);
size_t l = c - p;
memcpy (res + 1, p, l);
if (l <= n + 1)
memset (res + 1 + l, '0', (int) n - l + 2);
sprintf (res + n + 3, "e%c%02d", exp >= 0 ? '+' : '-', ABS(exp));
return;
}
else
{
/* If we have no exponent, normalize and add the exponent. */
for (c = p; *c && *c != '.'; c++)
;
res[1] = *p;
res[2] = '.';
size_t l = c - p;
memcpy (res + 3, p + 1, l);
size_t l2 = strlen (c + 1);
memcpy (res + 2 + l, c + 1, l2);
memset (res + 2 + l + l2, '0', n - (l + l2) + 1);
sprintf (res + n + 3, "e+%02d", l - 1);
return;
}
}
void
quadmath_flt128tostr (char *s, size_t size, size_t n, __float128 x)
{
char buffer[1024];
memset (buffer, 0, sizeof(buffer));
format (buffer, x, n);
memcpy (s, buffer, size);
}
...@@ -132,6 +132,6 @@ __qmath3 (ctanhq) ...@@ -132,6 +132,6 @@ __qmath3 (ctanhq)
/* Prototypes for string <-> flt128 conversion functions. */ /* Prototypes for string <-> flt128 conversion functions. */
__qmath3 (strtoflt128) __qmath3 (strtoflt128)
__qmath3 (quadmath_flt128tostr) __qmath3 (quadmath_snprintf)
#endif #endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment