From 554fd8c5195424bdbcabf5de30fdc183aba391bd Mon Sep 17 00:00:00 2001 From: upstream source tree Date: Sun, 15 Mar 2015 20:14:05 -0400 Subject: obtained gcc-4.6.4.tar.bz2 from upstream website; verified gcc-4.6.4.tar.bz2.sig; imported gcc-4.6.4 source tree from verified upstream tarball. downloading a git-generated archive based on the 'upstream' tag should provide you with a source tree that is binary identical to the one extracted from the above tarball. if you have obtained the source via the command 'git clone', however, do note that line-endings of files in your working directory might differ from line-endings of the respective files in the upstream repository. --- libquadmath/printf/_i18n_number.h | 139 ++++ libquadmath/printf/_itoa.h | 77 ++ libquadmath/printf/_itowa.h | 83 +++ libquadmath/printf/add_n.c | 62 ++ libquadmath/printf/addmul_1.c | 64 ++ libquadmath/printf/cmp.c | 56 ++ libquadmath/printf/divrem.c | 244 +++++++ libquadmath/printf/flt1282mpn.c | 136 ++++ libquadmath/printf/fpioconst.c | 463 ++++++++++++ libquadmath/printf/fpioconst.h | 64 ++ libquadmath/printf/gmp-impl.h | 181 +++++ libquadmath/printf/lshift.c | 87 +++ libquadmath/printf/mul.c | 148 ++++ libquadmath/printf/mul_1.c | 58 ++ libquadmath/printf/mul_n.c | 220 ++++++ libquadmath/printf/printf_fp.c | 1306 ++++++++++++++++++++++++++++++++++ libquadmath/printf/printf_fphex.c | 463 ++++++++++++ libquadmath/printf/quadmath-printf.c | 422 +++++++++++ libquadmath/printf/quadmath-printf.h | 186 +++++ libquadmath/printf/rshift.c | 88 +++ libquadmath/printf/sub_n.c | 62 ++ libquadmath/printf/submul_1.c | 64 ++ 22 files changed, 4673 insertions(+) create mode 100644 libquadmath/printf/_i18n_number.h create mode 100644 libquadmath/printf/_itoa.h create mode 100644 libquadmath/printf/_itowa.h create mode 100644 libquadmath/printf/add_n.c create mode 100644 libquadmath/printf/addmul_1.c create mode 100644 libquadmath/printf/cmp.c create mode 100644 libquadmath/printf/divrem.c create mode 100644 libquadmath/printf/flt1282mpn.c create mode 100644 libquadmath/printf/fpioconst.c create mode 100644 libquadmath/printf/fpioconst.h create mode 100644 libquadmath/printf/gmp-impl.h create mode 100644 libquadmath/printf/lshift.c create mode 100644 libquadmath/printf/mul.c create mode 100644 libquadmath/printf/mul_1.c create mode 100644 libquadmath/printf/mul_n.c create mode 100644 libquadmath/printf/printf_fp.c create mode 100644 libquadmath/printf/printf_fphex.c create mode 100644 libquadmath/printf/quadmath-printf.c create mode 100644 libquadmath/printf/quadmath-printf.h create mode 100644 libquadmath/printf/rshift.c create mode 100644 libquadmath/printf/sub_n.c create mode 100644 libquadmath/printf/submul_1.c (limited to 'libquadmath/printf') diff --git a/libquadmath/printf/_i18n_number.h b/libquadmath/printf/_i18n_number.h new file mode 100644 index 000000000..80065d88b --- /dev/null +++ b/libquadmath/printf/_i18n_number.h @@ -0,0 +1,139 @@ +/* Copyright (C) 2000, 2004, 2008 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 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; +} diff --git a/libquadmath/printf/_itoa.h b/libquadmath/printf/_itoa.h new file mode 100644 index 000000000..a0cd2b05c --- /dev/null +++ b/libquadmath/printf/_itoa.h @@ -0,0 +1,77 @@ +/* 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 */ diff --git a/libquadmath/printf/_itowa.h b/libquadmath/printf/_itowa.h new file mode 100644 index 000000000..4717b5c65 --- /dev/null +++ b/libquadmath/printf/_itowa.h @@ -0,0 +1,83 @@ +/* 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 */ diff --git a/libquadmath/printf/add_n.c b/libquadmath/printf/add_n.c new file mode 100644 index 000000000..749cf3ee8 --- /dev/null +++ b/libquadmath/printf/add_n.c @@ -0,0 +1,62 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/addmul_1.c b/libquadmath/printf/addmul_1.c new file mode 100644 index 000000000..f527f9848 --- /dev/null +++ b/libquadmath/printf/addmul_1.c @@ -0,0 +1,64 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/cmp.c b/libquadmath/printf/cmp.c new file mode 100644 index 000000000..a4be43e2a --- /dev/null +++ b/libquadmath/printf/cmp.c @@ -0,0 +1,56 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/divrem.c b/libquadmath/printf/divrem.c new file mode 100644 index 000000000..944d1a0c0 --- /dev/null +++ b/libquadmath/printf/divrem.c @@ -0,0 +1,244 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/flt1282mpn.c b/libquadmath/printf/flt1282mpn.c new file mode 100644 index 000000000..0105314ef --- /dev/null +++ b/libquadmath/printf/flt1282mpn.c @@ -0,0 +1,136 @@ +/* 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 +#include +#include +#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; +} diff --git a/libquadmath/printf/fpioconst.c b/libquadmath/printf/fpioconst.c new file mode 100644 index 000000000..8c67e6f90 --- /dev/null +++ b/libquadmath/printf/fpioconst.c @@ -0,0 +1,463 @@ +/* Table of MP integer constants 10^(2^i), used for floating point <-> decimal. + 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 +#include "gmp-impl.h" /* This defines BITS_PER_MP_LIMB. */ +#include "fpioconst.h" + +/* First page : 32-bit limbs + Second page : 64-bit limbs + Last page : table of pointers + */ + +#if BITS_PER_MP_LIMB == 32 + +/* Table with constants of 10^(2^i), i=0..12 for 32-bit limbs. */ + +const mp_limb_t __tens[] = +{ +#define TENS_P0_IDX 0 +#define TENS_P0_SIZE 3 + [TENS_P0_IDX] = 0x00000000, 0x00000000, 0x0000000a, + +#define TENS_P1_IDX (TENS_P0_IDX + TENS_P0_SIZE) +#define TENS_P1_SIZE 3 + [TENS_P1_IDX] = 0x00000000, 0x00000000, 0x00000064, + +#define TENS_P2_IDX (TENS_P1_IDX + TENS_P1_SIZE) +#define TENS_P2_SIZE 3 + [TENS_P2_IDX] = 0x00000000, 0x00000000, 0x00002710, + +#define TENS_P3_IDX (TENS_P2_IDX + TENS_P2_SIZE) +#define TENS_P3_SIZE 3 + [TENS_P3_IDX] = 0x00000000, 0x00000000, 0x05f5e100, + +#define TENS_P4_IDX (TENS_P3_IDX + TENS_P3_SIZE) +#define TENS_P4_SIZE 4 + [TENS_P4_IDX] = 0x00000000, 0x00000000, 0x6fc10000, 0x002386f2, + +#define TENS_P5_IDX (TENS_P4_IDX + TENS_P4_SIZE) +#define TENS_P5_SIZE 6 + [TENS_P5_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x85acef81, 0x2d6d415b, + 0x000004ee, + +#define TENS_P6_IDX (TENS_P5_IDX + TENS_P5_SIZE) +#define TENS_P6_SIZE 9 + [TENS_P6_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0xbf6a1f01, + 0x6e38ed64, 0xdaa797ed, 0xe93ff9f4, 0x00184f03, + +#define TENS_P7_IDX (TENS_P6_IDX + TENS_P6_SIZE) +#define TENS_P7_SIZE 16 + [TENS_P7_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x2e953e01, 0x03df9909, 0x0f1538fd, 0x2374e42f, 0xd3cff5ec, + 0xc404dc08, 0xbccdb0da, 0xa6337f19, 0xe91f2603, 0x0000024e, + +#define TENS_P8_IDX (TENS_P7_IDX + TENS_P7_SIZE) +#define TENS_P8_SIZE 29 + [TENS_P8_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x982e7c01, + 0xbed3875b, 0xd8d99f72, 0x12152f87, 0x6bde50c6, 0xcf4a6e70, 0xd595d80f, + 0x26b2716e, 0xadc666b0, 0x1d153624, 0x3c42d35a, 0x63ff540e, 0xcc5573c0, + 0x65f9ef17, 0x55bc28f2, 0x80dcc7f7, 0xf46eeddc, 0x5fdcefce, 0x000553f7, + +#ifndef __NO_LONG_DOUBLE_MATH +# define TENS_P9_IDX (TENS_P8_IDX + TENS_P8_SIZE) +# define TENS_P9_SIZE 56 + [TENS_P9_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0xfc6cf801, 0x77f27267, 0x8f9546dc, 0x5d96976f, 0xb83a8a97, + 0xc31e1ad9, 0x46c40513, 0x94e65747, 0xc88976c1, 0x4475b579, 0x28f8733b, + 0xaa1da1bf, 0x703ed321, 0x1e25cfea, 0xb21a2f22, 0xbc51fb2e, 0x96e14f5d, + 0xbfa3edac, 0x329c57ae, 0xe7fc7153, 0xc3fc0695, 0x85a91924, 0xf95f635e, + 0xb2908ee0, 0x93abade4, 0x1366732a, 0x9449775c, 0x69be5b0e, 0x7343afac, + 0xb099bc81, 0x45a71d46, 0xa2699748, 0x8cb07303, 0x8a0b1f13, 0x8cab8a97, + 0xc1d238d9, 0x633415d4, 0x0000001c, + +# define TENS_P10_IDX (TENS_P9_IDX + TENS_P9_SIZE) +# define TENS_P10_SIZE 109 + [TENS_P10_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x2919f001, + 0xf55b2b72, 0x6e7c215b, 0x1ec29f86, 0x991c4e87, 0x15c51a88, 0x140ac535, + 0x4c7d1e1a, 0xcc2cd819, 0x0ed1440e, 0x896634ee, 0x7de16cfb, 0x1e43f61f, + 0x9fce837d, 0x231d2b9c, 0x233e55c7, 0x65dc60d7, 0xf451218b, 0x1c5cd134, + 0xc9635986, 0x922bbb9f, 0xa7e89431, 0x9f9f2a07, 0x62be695a, 0x8e1042c4, + 0x045b7a74, 0x1abe1de3, 0x8ad822a5, 0xba34c411, 0xd814b505, 0xbf3fdeb3, + 0x8fc51a16, 0xb1b896bc, 0xf56deeec, 0x31fb6bfd, 0xb6f4654b, 0x101a3616, + 0x6b7595fb, 0xdc1a47fe, 0x80d98089, 0x80bda5a5, 0x9a202882, 0x31eb0f66, + 0xfc8f1f90, 0x976a3310, 0xe26a7b7e, 0xdf68368a, 0x3ce3a0b8, 0x8e4262ce, + 0x75a351a2, 0x6cb0b6c9, 0x44597583, 0x31b5653f, 0xc356e38a, 0x35faaba6, + 0x0190fba0, 0x9fc4ed52, 0x88bc491b, 0x1640114a, 0x005b8041, 0xf4f3235e, + 0x1e8d4649, 0x36a8de06, 0x73c55349, 0xa7e6bd2a, 0xc1a6970c, 0x47187094, + 0xd2db49ef, 0x926c3f5b, 0xae6209d4, 0x2d433949, 0x34f4a3c6, 0xd4305d94, + 0xd9d61a05, 0x00000325, + +# define TENS_P11_IDX (TENS_P10_IDX + TENS_P10_SIZE) +# define TENS_P11_SIZE 215 + [TENS_P11_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x1333e001, 0xe3096865, 0xb27d4d3f, 0x49e28dcf, 0xec2e4721, + 0xee87e354, 0xb6067584, 0x368b8abb, 0xa5e5a191, 0x2ed56d55, 0xfd827773, + 0xea50d142, 0x51b78db2, 0x98342c9e, 0xc850dabc, 0x866ed6f1, 0x19342c12, + 0x92794987, 0xd2f869c2, 0x66912e4a, 0x71c7fd8f, 0x57a7842d, 0x235552eb, + 0xfb7fedcc, 0xf3861ce0, 0x38209ce1, 0x9713b449, 0x34c10134, 0x8c6c54de, + 0xa7a8289c, 0x2dbb6643, 0xe3cb64f3, 0x8074ff01, 0xe3892ee9, 0x10c17f94, + 0xa8f16f92, 0xa8281ed6, 0x967abbb3, 0x5a151440, 0x9952fbed, 0x13b41e44, + 0xafe609c3, 0xa2bca416, 0xf111821f, 0xfb1264b4, 0x91bac974, 0xd6c7d6ab, + 0x8e48ff35, 0x4419bd43, 0xc4a65665, 0x685e5510, 0x33554c36, 0xab498697, + 0x0dbd21fe, 0x3cfe491d, 0x982da466, 0xcbea4ca7, 0x9e110c7b, 0x79c56b8a, + 0x5fc5a047, 0x84d80e2e, 0x1aa9f444, 0x730f203c, 0x6a57b1ab, 0xd752f7a6, + 0x87a7dc62, 0x944545ff, 0x40660460, 0x77c1a42f, 0xc9ac375d, 0xe866d7ef, + 0x744695f0, 0x81428c85, 0xa1fc6b96, 0xd7917c7b, 0x7bf03c19, 0x5b33eb41, + 0x5715f791, 0x8f6cae5f, 0xdb0708fd, 0xb125ac8e, 0x785ce6b7, 0x56c6815b, + 0x6f46eadb, 0x4eeebeee, 0x195355d8, 0xa244de3c, 0x9d7389c0, 0x53761abd, + 0xcf99d019, 0xde9ec24b, 0x0d76ce39, 0x70beb181, 0x2e55ecee, 0xd5f86079, + 0xf56d9d4b, 0xfb8886fb, 0x13ef5a83, 0x408f43c5, 0x3f3389a4, 0xfad37943, + 0x58ccf45c, 0xf82df846, 0x415c7f3e, 0x2915e818, 0x8b3d5cf4, 0x6a445f27, + 0xf8dbb57a, 0xca8f0070, 0x8ad803ec, 0xb2e87c34, 0x038f9245, 0xbedd8a6c, + 0xc7c9dee0, 0x0eac7d56, 0x2ad3fa14, 0xe0de0840, 0xf775677c, 0xf1bd0ad5, + 0x92be221e, 0x87fa1fb9, 0xce9d04a4, 0xd2c36fa9, 0x3f6f7024, 0xb028af62, + 0x907855ee, 0xd83e49d6, 0x4efac5dc, 0xe7151aab, 0x77cd8c6b, 0x0a753b7d, + 0x0af908b4, 0x8c983623, 0xe50f3027, 0x94222771, 0x1d08e2d6, 0xf7e928e6, + 0xf2ee5ca6, 0x1b61b93c, 0x11eb962b, 0x9648b21c, 0xce2bcba1, 0x34f77154, + 0x7bbebe30, 0xe526a319, 0x8ce329ac, 0xde4a74d2, 0xb5dc53d5, 0x0009e8b3, + +# define TENS_P12_IDX (TENS_P11_IDX + TENS_P11_SIZE) +# define TENS_P12_SIZE 428 + [TENS_P12_IDX] = 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, + 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x2a67c001, + 0xd4724e8d, 0x8efe7ae7, 0xf89a1e90, 0xef084117, 0x54e05154, 0x13b1bb51, + 0x506be829, 0xfb29b172, 0xe599574e, 0xf0da6146, 0x806c0ed3, 0xb86ae5be, + 0x45155e93, 0xc0591cc2, 0x7e1e7c34, 0x7c4823da, 0x1d1f4cce, 0x9b8ba1e8, + 0xd6bfdf75, 0xe341be10, 0xc2dfae78, 0x016b67b2, 0x0f237f1a, 0x3dbeabcd, + 0xaf6a2574, 0xcab3e6d7, 0x142e0e80, 0x61959127, 0x2c234811, 0x87009701, + 0xcb4bf982, 0xf8169c84, 0x88052f8c, 0x68dde6d4, 0xbc131761, 0xff0b0905, + 0x54ab9c41, 0x7613b224, 0x1a1c304e, 0x3bfe167b, 0x441c2d47, 0x4f6cea9c, + 0x78f06181, 0xeb659fb8, 0x30c7ae41, 0x947e0d0e, 0xa1ebcad7, 0xd97d9556, + 0x2130504d, 0x1a8309cb, 0xf2acd507, 0x3f8ec72a, 0xfd82373a, 0x95a842bc, + 0x280f4d32, 0xf3618ac0, 0x811a4f04, 0x6dc3a5b4, 0xd3967a1b, 0x15b8c898, + 0xdcfe388f, 0x454eb2a0, 0x8738b909, 0x10c4e996, 0x2bd9cc11, 0x3297cd0c, + 0x655fec30, 0xae0725b1, 0xf4090ee8, 0x037d19ee, 0x398c6fed, 0x3b9af26b, + 0xc994a450, 0xb5341743, 0x75a697b2, 0xac50b9c1, 0x3ccb5b92, 0xffe06205, + 0xa8329761, 0xdfea5242, 0xeb83cadb, 0xe79dadf7, 0x3c20ee69, 0x1e0a6817, + 0x7021b97a, 0x743074fa, 0x176ca776, 0x77fb8af6, 0xeca19beb, 0x92baf1de, + 0xaf63b712, 0xde35c88b, 0xa4eb8f8c, 0xe137d5e9, 0x40b464a0, 0x87d1cde8, + 0x42923bbd, 0xcd8f62ff, 0x2e2690f3, 0x095edc16, 0x59c89f1b, 0x1fa8fd5d, + 0x5138753d, 0x390a2b29, 0x80152f18, 0x2dd8d925, 0xf984d83e, 0x7a872e74, + 0xc19e1faf, 0xed4d542d, 0xecf9b5d0, 0x9462ea75, 0xc53c0adf, 0x0caea134, + 0x37a2d439, 0xc8fa2e8a, 0x2181327e, 0x6e7bb827, 0x2d240820, 0x50be10e0, + 0x5893d4b8, 0xab312bb9, 0x1f2b2322, 0x440b3f25, 0xbf627ede, 0x72dac789, + 0xb608b895, 0x78787e2a, 0x86deb3f0, 0x6fee7aab, 0xbb9373f4, 0x27ecf57b, + 0xf7d8b57e, 0xfca26a9f, 0x3d04e8d2, 0xc9df13cb, 0x3172826a, 0xcd9e8d7c, + 0xa8fcd8e0, 0xb2c39497, 0x307641d9, 0x1cc939c1, 0x2608c4cf, 0xb6d1c7bf, + 0x3d326a7e, 0xeeaf19e6, 0x8e13e25f, 0xee63302b, 0x2dfe6d97, 0x25971d58, + 0xe41d3cc4, 0x0a80627c, 0xab8db59a, 0x9eea37c8, 0xe90afb77, 0x90ca19cf, + 0x9ee3352c, 0x3613c850, 0xfe78d682, 0x788f6e50, 0x5b060904, 0xb71bd1a4, + 0x3fecb534, 0xb32c450c, 0x20c33857, 0xa6e9cfda, 0x0239f4ce, 0x48497187, + 0xa19adb95, 0xb492ed8a, 0x95aca6a8, 0x4dcd6cd9, 0xcf1b2350, 0xfbe8b12a, + 0x1a67778c, 0x38eb3acc, 0xc32da383, 0xfb126ab1, 0xa03f40a8, 0xed5bf546, + 0xe9ce4724, 0x4c4a74fd, 0x73a130d8, 0xd9960e2d, 0xa2ebd6c1, 0x94ab6feb, + 0x6f233b7c, 0x49126080, 0x8e7b9a73, 0x4b8c9091, 0xd298f999, 0x35e836b5, + 0xa96ddeff, 0x96119b31, 0x6b0dd9bc, 0xc6cc3f8d, 0x282566fb, 0x72b882e7, + 0xd6769f3b, 0xa674343d, 0x00fc509b, 0xdcbf7789, 0xd6266a3f, 0xae9641fd, + 0x4e89541b, 0x11953407, 0x53400d03, 0x8e0dd75a, 0xe5b53345, 0x108f19ad, + 0x108b89bc, 0x41a4c954, 0xe03b2b63, 0x437b3d7f, 0x97aced8e, 0xcbd66670, + 0x2c5508c2, 0x650ebc69, 0x5c4f2ef0, 0x904ff6bf, 0x9985a2df, 0x9faddd9e, + 0x5ed8d239, 0x25585832, 0xe3e51cb9, 0x0ff4f1d4, 0x56c02d9a, 0x8c4ef804, + 0xc1a08a13, 0x13fd01c8, 0xe6d27671, 0xa7c234f4, 0x9d0176cc, 0xd0d73df2, + 0x4d8bfa89, 0x544f10cd, 0x2b17e0b2, 0xb70a5c7d, 0xfd86fe49, 0xdf373f41, + 0x214495bb, 0x84e857fd, 0x00d313d5, 0x0496fcbe, 0xa4ba4744, 0xe8cac982, + 0xaec29e6e, 0x87ec7038, 0x7000a519, 0xaeee333b, 0xff66e42c, 0x8afd6b25, + 0x03b4f63b, 0xbd7991dc, 0x5ab8d9c7, 0x2ed4684e, 0x48741a6c, 0xaf06940d, + 0x2fdc6349, 0xb03d7ecd, 0xe974996f, 0xac7867f9, 0x52ec8721, 0xbcdd9d4a, + 0x8edd2d00, 0x3557de06, 0x41c759f8, 0x3956d4b9, 0xa75409f2, 0x123cd8a1, + 0xb6100fab, 0x3e7b21e2, 0x2e8d623b, 0x92959da2, 0xbca35f77, 0x200c03a5, + 0x35fcb457, 0x1bb6c6e4, 0xf74eb928, 0x3d5d0b54, 0x87cc1d21, 0x4964046f, + 0x18ae4240, 0xd868b275, 0x8bd2b496, 0x1c5563f4, 0xc234d8f5, 0xf868e970, + 0xf9151fff, 0xae7be4a2, 0x271133ee, 0xbb0fd922, 0x25254932, 0xa60a9fc0, + 0x104bcd64, 0x30290145, 0x00000062 +#endif /* !__NO_LONG_DOUBLE_MATH */ +}; + +#elif BITS_PER_MP_LIMB == 64 + +/* Table with constants of 10^(2^i), i=0..12 for 64-bit limbs. */ + +const mp_limb_t __tens[] = +{ +#define TENS_P0_IDX 0 +#define TENS_P0_SIZE 2 + [TENS_P0_IDX] = 0x0000000000000000ull, 0x000000000000000aull, + +#define TENS_P1_IDX (TENS_P0_IDX + TENS_P0_SIZE) +#define TENS_P1_SIZE 2 + [TENS_P1_IDX] = 0x0000000000000000ull, 0x0000000000000064ull, + +#define TENS_P2_IDX (TENS_P1_IDX + TENS_P1_SIZE) +#define TENS_P2_SIZE 2 + [TENS_P2_IDX] = 0x0000000000000000ull, 0x0000000000002710ull, + +#define TENS_P3_IDX (TENS_P2_IDX + TENS_P2_SIZE) +#define TENS_P3_SIZE 2 + [TENS_P3_IDX] = 0x0000000000000000ull, 0x0000000005f5e100ull, + +#define TENS_P4_IDX (TENS_P3_IDX + TENS_P3_SIZE) +#define TENS_P4_SIZE 2 + [TENS_P4_IDX] = 0x0000000000000000ull, 0x002386f26fc10000ull, + +#define TENS_P5_IDX (TENS_P4_IDX + TENS_P4_SIZE) +#define TENS_P5_SIZE 3 + [TENS_P5_IDX] = 0x0000000000000000ull, 0x85acef8100000000ull, + 0x000004ee2d6d415bull, + +#define TENS_P6_IDX (TENS_P5_IDX + TENS_P5_SIZE) +#define TENS_P6_SIZE 5 + [TENS_P6_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x6e38ed64bf6a1f01ull, 0xe93ff9f4daa797edull, 0x0000000000184f03ull, + +#define TENS_P7_IDX (TENS_P6_IDX + TENS_P6_SIZE) +#define TENS_P7_SIZE 8 + [TENS_P7_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x03df99092e953e01ull, 0x2374e42f0f1538fdull, + 0xc404dc08d3cff5ecull, 0xa6337f19bccdb0daull, 0x0000024ee91f2603ull, + +#define TENS_P8_IDX (TENS_P7_IDX + TENS_P7_SIZE) +#define TENS_P8_SIZE 15 + [TENS_P8_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0xbed3875b982e7c01ull, 0x12152f87d8d99f72ull, 0xcf4a6e706bde50c6ull, + 0x26b2716ed595d80full, 0x1d153624adc666b0ull, 0x63ff540e3c42d35aull, + 0x65f9ef17cc5573c0ull, 0x80dcc7f755bc28f2ull, 0x5fdcefcef46eeddcull, + 0x00000000000553f7ull, +#if FLT128_MAX_EXP > 1024 +# define TENS_P9_IDX (TENS_P8_IDX + TENS_P8_SIZE) +# define TENS_P9_SIZE 28 + [TENS_P9_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x77f27267fc6cf801ull, 0x5d96976f8f9546dcull, + 0xc31e1ad9b83a8a97ull, 0x94e6574746c40513ull, 0x4475b579c88976c1ull, + 0xaa1da1bf28f8733bull, 0x1e25cfea703ed321ull, 0xbc51fb2eb21a2f22ull, + 0xbfa3edac96e14f5dull, 0xe7fc7153329c57aeull, 0x85a91924c3fc0695ull, + 0xb2908ee0f95f635eull, 0x1366732a93abade4ull, 0x69be5b0e9449775cull, + 0xb099bc817343afacull, 0xa269974845a71d46ull, 0x8a0b1f138cb07303ull, + 0xc1d238d98cab8a97ull, 0x0000001c633415d4ull, + +# define TENS_P10_IDX (TENS_P9_IDX + TENS_P9_SIZE) +# define TENS_P10_SIZE 55 + [TENS_P10_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0xf55b2b722919f001ull, 0x1ec29f866e7c215bull, 0x15c51a88991c4e87ull, + 0x4c7d1e1a140ac535ull, 0x0ed1440ecc2cd819ull, 0x7de16cfb896634eeull, + 0x9fce837d1e43f61full, 0x233e55c7231d2b9cull, 0xf451218b65dc60d7ull, + 0xc96359861c5cd134ull, 0xa7e89431922bbb9full, 0x62be695a9f9f2a07ull, + 0x045b7a748e1042c4ull, 0x8ad822a51abe1de3ull, 0xd814b505ba34c411ull, + 0x8fc51a16bf3fdeb3ull, 0xf56deeecb1b896bcull, 0xb6f4654b31fb6bfdull, + 0x6b7595fb101a3616ull, 0x80d98089dc1a47feull, 0x9a20288280bda5a5ull, + 0xfc8f1f9031eb0f66ull, 0xe26a7b7e976a3310ull, 0x3ce3a0b8df68368aull, + 0x75a351a28e4262ceull, 0x445975836cb0b6c9ull, 0xc356e38a31b5653full, + 0x0190fba035faaba6ull, 0x88bc491b9fc4ed52ull, 0x005b80411640114aull, + 0x1e8d4649f4f3235eull, 0x73c5534936a8de06ull, 0xc1a6970ca7e6bd2aull, + 0xd2db49ef47187094ull, 0xae6209d4926c3f5bull, 0x34f4a3c62d433949ull, + 0xd9d61a05d4305d94ull, 0x0000000000000325ull, + +# define TENS_P11_IDX (TENS_P10_IDX + TENS_P10_SIZE) +# define TENS_P11_SIZE 108 + [TENS_P11_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0xe30968651333e001ull, 0x49e28dcfb27d4d3full, + 0xee87e354ec2e4721ull, 0x368b8abbb6067584ull, 0x2ed56d55a5e5a191ull, + 0xea50d142fd827773ull, 0x98342c9e51b78db2ull, 0x866ed6f1c850dabcull, + 0x9279498719342c12ull, 0x66912e4ad2f869c2ull, 0x57a7842d71c7fd8full, + 0xfb7fedcc235552ebull, 0x38209ce1f3861ce0ull, 0x34c101349713b449ull, + 0xa7a8289c8c6c54deull, 0xe3cb64f32dbb6643ull, 0xe3892ee98074ff01ull, + 0xa8f16f9210c17f94ull, 0x967abbb3a8281ed6ull, 0x9952fbed5a151440ull, + 0xafe609c313b41e44ull, 0xf111821fa2bca416ull, 0x91bac974fb1264b4ull, + 0x8e48ff35d6c7d6abull, 0xc4a656654419bd43ull, 0x33554c36685e5510ull, + 0x0dbd21feab498697ull, 0x982da4663cfe491dull, 0x9e110c7bcbea4ca7ull, + 0x5fc5a04779c56b8aull, 0x1aa9f44484d80e2eull, 0x6a57b1ab730f203cull, + 0x87a7dc62d752f7a6ull, 0x40660460944545ffull, 0xc9ac375d77c1a42full, + 0x744695f0e866d7efull, 0xa1fc6b9681428c85ull, 0x7bf03c19d7917c7bull, + 0x5715f7915b33eb41ull, 0xdb0708fd8f6cae5full, 0x785ce6b7b125ac8eull, + 0x6f46eadb56c6815bull, 0x195355d84eeebeeeull, 0x9d7389c0a244de3cull, + 0xcf99d01953761abdull, 0x0d76ce39de9ec24bull, 0x2e55ecee70beb181ull, + 0xf56d9d4bd5f86079ull, 0x13ef5a83fb8886fbull, 0x3f3389a4408f43c5ull, + 0x58ccf45cfad37943ull, 0x415c7f3ef82df846ull, 0x8b3d5cf42915e818ull, + 0xf8dbb57a6a445f27ull, 0x8ad803ecca8f0070ull, 0x038f9245b2e87c34ull, + 0xc7c9dee0bedd8a6cull, 0x2ad3fa140eac7d56ull, 0xf775677ce0de0840ull, + 0x92be221ef1bd0ad5ull, 0xce9d04a487fa1fb9ull, 0x3f6f7024d2c36fa9ull, + 0x907855eeb028af62ull, 0x4efac5dcd83e49d6ull, 0x77cd8c6be7151aabull, + 0x0af908b40a753b7dull, 0xe50f30278c983623ull, 0x1d08e2d694222771ull, + 0xf2ee5ca6f7e928e6ull, 0x11eb962b1b61b93cull, 0xce2bcba19648b21cull, + 0x7bbebe3034f77154ull, 0x8ce329ace526a319ull, 0xb5dc53d5de4a74d2ull, + 0x000000000009e8b3ull, + +# define TENS_P12_IDX (TENS_P11_IDX + TENS_P11_SIZE) +# define TENS_P12_SIZE 214 + [TENS_P12_IDX] = 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0x0000000000000000ull, 0x0000000000000000ull, 0x0000000000000000ull, + 0xd4724e8d2a67c001ull, 0xf89a1e908efe7ae7ull, 0x54e05154ef084117ull, + 0x506be82913b1bb51ull, 0xe599574efb29b172ull, 0x806c0ed3f0da6146ull, + 0x45155e93b86ae5beull, 0x7e1e7c34c0591cc2ull, 0x1d1f4cce7c4823daull, + 0xd6bfdf759b8ba1e8ull, 0xc2dfae78e341be10ull, 0x0f237f1a016b67b2ull, + 0xaf6a25743dbeabcdull, 0x142e0e80cab3e6d7ull, 0x2c23481161959127ull, + 0xcb4bf98287009701ull, 0x88052f8cf8169c84ull, 0xbc13176168dde6d4ull, + 0x54ab9c41ff0b0905ull, 0x1a1c304e7613b224ull, 0x441c2d473bfe167bull, + 0x78f061814f6cea9cull, 0x30c7ae41eb659fb8ull, 0xa1ebcad7947e0d0eull, + 0x2130504dd97d9556ull, 0xf2acd5071a8309cbull, 0xfd82373a3f8ec72aull, + 0x280f4d3295a842bcull, 0x811a4f04f3618ac0ull, 0xd3967a1b6dc3a5b4ull, + 0xdcfe388f15b8c898ull, 0x8738b909454eb2a0ull, 0x2bd9cc1110c4e996ull, + 0x655fec303297cd0cull, 0xf4090ee8ae0725b1ull, 0x398c6fed037d19eeull, + 0xc994a4503b9af26bull, 0x75a697b2b5341743ull, 0x3ccb5b92ac50b9c1ull, + 0xa8329761ffe06205ull, 0xeb83cadbdfea5242ull, 0x3c20ee69e79dadf7ull, + 0x7021b97a1e0a6817ull, 0x176ca776743074faull, 0xeca19beb77fb8af6ull, + 0xaf63b71292baf1deull, 0xa4eb8f8cde35c88bull, 0x40b464a0e137d5e9ull, + 0x42923bbd87d1cde8ull, 0x2e2690f3cd8f62ffull, 0x59c89f1b095edc16ull, + 0x5138753d1fa8fd5dull, 0x80152f18390a2b29ull, 0xf984d83e2dd8d925ull, + 0xc19e1faf7a872e74ull, 0xecf9b5d0ed4d542dull, 0xc53c0adf9462ea75ull, + 0x37a2d4390caea134ull, 0x2181327ec8fa2e8aull, 0x2d2408206e7bb827ull, + 0x5893d4b850be10e0ull, 0x1f2b2322ab312bb9ull, 0xbf627ede440b3f25ull, + 0xb608b89572dac789ull, 0x86deb3f078787e2aull, 0xbb9373f46fee7aabull, + 0xf7d8b57e27ecf57bull, 0x3d04e8d2fca26a9full, 0x3172826ac9df13cbull, + 0xa8fcd8e0cd9e8d7cull, 0x307641d9b2c39497ull, 0x2608c4cf1cc939c1ull, + 0x3d326a7eb6d1c7bfull, 0x8e13e25feeaf19e6ull, 0x2dfe6d97ee63302bull, + 0xe41d3cc425971d58ull, 0xab8db59a0a80627cull, 0xe90afb779eea37c8ull, + 0x9ee3352c90ca19cfull, 0xfe78d6823613c850ull, 0x5b060904788f6e50ull, + 0x3fecb534b71bd1a4ull, 0x20c33857b32c450cull, 0x0239f4cea6e9cfdaull, + 0xa19adb9548497187ull, 0x95aca6a8b492ed8aull, 0xcf1b23504dcd6cd9ull, + 0x1a67778cfbe8b12aull, 0xc32da38338eb3accull, 0xa03f40a8fb126ab1ull, + 0xe9ce4724ed5bf546ull, 0x73a130d84c4a74fdull, 0xa2ebd6c1d9960e2dull, + 0x6f233b7c94ab6febull, 0x8e7b9a7349126080ull, 0xd298f9994b8c9091ull, + 0xa96ddeff35e836b5ull, 0x6b0dd9bc96119b31ull, 0x282566fbc6cc3f8dull, + 0xd6769f3b72b882e7ull, 0x00fc509ba674343dull, 0xd6266a3fdcbf7789ull, + 0x4e89541bae9641fdull, 0x53400d0311953407ull, 0xe5b533458e0dd75aull, + 0x108b89bc108f19adull, 0xe03b2b6341a4c954ull, 0x97aced8e437b3d7full, + 0x2c5508c2cbd66670ull, 0x5c4f2ef0650ebc69ull, 0x9985a2df904ff6bfull, + 0x5ed8d2399faddd9eull, 0xe3e51cb925585832ull, 0x56c02d9a0ff4f1d4ull, + 0xc1a08a138c4ef804ull, 0xe6d2767113fd01c8ull, 0x9d0176cca7c234f4ull, + 0x4d8bfa89d0d73df2ull, 0x2b17e0b2544f10cdull, 0xfd86fe49b70a5c7dull, + 0x214495bbdf373f41ull, 0x00d313d584e857fdull, 0xa4ba47440496fcbeull, + 0xaec29e6ee8cac982ull, 0x7000a51987ec7038ull, 0xff66e42caeee333bull, + 0x03b4f63b8afd6b25ull, 0x5ab8d9c7bd7991dcull, 0x48741a6c2ed4684eull, + 0x2fdc6349af06940dull, 0xe974996fb03d7ecdull, 0x52ec8721ac7867f9ull, + 0x8edd2d00bcdd9d4aull, 0x41c759f83557de06ull, 0xa75409f23956d4b9ull, + 0xb6100fab123cd8a1ull, 0x2e8d623b3e7b21e2ull, 0xbca35f7792959da2ull, + 0x35fcb457200c03a5ull, 0xf74eb9281bb6c6e4ull, 0x87cc1d213d5d0b54ull, + 0x18ae42404964046full, 0x8bd2b496d868b275ull, 0xc234d8f51c5563f4ull, + 0xf9151ffff868e970ull, 0x271133eeae7be4a2ull, 0x25254932bb0fd922ull, + 0x104bcd64a60a9fc0ull, 0x0000006230290145ull +#endif +}; + +#else +# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for" +#endif + +/* Each of array variable above defines one mpn integer which is a power of 10. + This table points to those variables, indexed by the exponent. */ + +const struct mp_power _fpioconst_pow10[FLT128_MAX_10_EXP_LOG + 1] = +{ + { TENS_P0_IDX, TENS_P0_SIZE, 4, }, + { TENS_P1_IDX, TENS_P1_SIZE, 7, 4 }, + { TENS_P2_IDX, TENS_P2_SIZE, 14, 10 }, + { TENS_P3_IDX, TENS_P3_SIZE, 27, 24 }, + { TENS_P4_IDX, TENS_P4_SIZE, 54, 50 }, + { TENS_P5_IDX, TENS_P5_SIZE, 107, 103 }, + { TENS_P6_IDX, TENS_P6_SIZE, 213, 210 }, + { TENS_P7_IDX, TENS_P7_SIZE, 426, 422 }, + { TENS_P8_IDX, TENS_P8_SIZE, 851, 848 }, +#if FLT128_MAX_EXP > 1024 + { TENS_P9_IDX, TENS_P9_SIZE, 1701, 1698 }, + { TENS_P10_IDX, TENS_P10_SIZE, 3402, 3399 }, + { TENS_P11_IDX, TENS_P11_SIZE, 6804, 6800 }, + { TENS_P12_IDX, TENS_P12_SIZE, 13607, 13604 } +#endif +}; + +#if LAST_POW10 > _LAST_POW10 +# error "Need to expand 10^(2^i) table for i up to" LAST_POW10 +#endif diff --git a/libquadmath/printf/fpioconst.h b/libquadmath/printf/fpioconst.h new file mode 100644 index 000000000..418755547 --- /dev/null +++ b/libquadmath/printf/fpioconst.h @@ -0,0 +1,64 @@ +/* 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 +#include + + +/* 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 . 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 */ diff --git a/libquadmath/printf/gmp-impl.h b/libquadmath/printf/gmp-impl.h new file mode 100644 index 000000000..ca49e1966 --- /dev/null +++ b/libquadmath/printf/gmp-impl.h @@ -0,0 +1,181 @@ +/* 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 +#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_construct_float128 __MPN(construct_float128) +__float128 mpn_construct_float128 (mp_srcptr frac_ptr, int expt, int sign) + 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; +} diff --git a/libquadmath/printf/lshift.c b/libquadmath/printf/lshift.c new file mode 100644 index 000000000..58aa8d464 --- /dev/null +++ b/libquadmath/printf/lshift.c @@ -0,0 +1,87 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/mul.c b/libquadmath/printf/mul.c new file mode 100644 index 000000000..d31fa36fa --- /dev/null +++ b/libquadmath/printf/mul.c @@ -0,0 +1,148 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/mul_1.c b/libquadmath/printf/mul_1.c new file mode 100644 index 000000000..48a273f07 --- /dev/null +++ b/libquadmath/printf/mul_1.c @@ -0,0 +1,58 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/mul_n.c b/libquadmath/printf/mul_n.c new file mode 100644 index 000000000..c4bc1bed2 --- /dev/null +++ b/libquadmath/printf/mul_n.c @@ -0,0 +1,220 @@ +/* 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 +#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); + } +} diff --git a/libquadmath/printf/printf_fp.c b/libquadmath/printf/printf_fp.c new file mode 100644 index 000000000..eb663726d --- /dev/null +++ b/libquadmath/printf/printf_fp.c @@ -0,0 +1,1306 @@ +/* Floating point output for `printf'. + Copyright (C) 1995-2003, 2006-2008, 2011 Free Software Foundation, Inc. + + This file is part of the GNU C Library. + Written by Ulrich Drepper , 1995. + + 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 +#include +#include +#include +#include +#include +#define NDEBUG +#include +#ifdef HAVE_ERRNO_H +#include +#endif +#include +#include +#include "quadmath-printf.h" +#include "fpioconst.h" + +#ifdef USE_I18N_NUMBER_H +#include "_i18n_number.h" +#endif + + +/* Macros for doing the actual output. */ + +#define outchar(ch) \ + do \ + { \ + register const int outc = (ch); \ + if (PUTC (outc, fp) == EOF) \ + { \ + if (buffer_malloced) \ + free (wbuffer); \ + return -1; \ + } \ + ++done; \ + } while (0) + +#define PRINT(ptr, wptr, len) \ + do \ + { \ + register size_t outlen = (len); \ + if (len > 20) \ + { \ + if (PUT (fp, wide ? (const char *) wptr : ptr, outlen) != outlen) \ + { \ + if (buffer_malloced) \ + free (wbuffer); \ + return -1; \ + } \ + ptr += outlen; \ + done += outlen; \ + } \ + else \ + { \ + if (wide) \ + while (outlen-- > 0) \ + outchar (*wptr++); \ + else \ + while (outlen-- > 0) \ + outchar (*ptr++); \ + } \ + } while (0) + +#define PADN(ch, len) \ + do \ + { \ + if (PAD (fp, ch, len) != len) \ + { \ + if (buffer_malloced) \ + free (wbuffer); \ + return -1; \ + } \ + done += len; \ + } \ + while (0) + + +/* We use the GNU MP library to handle large numbers. + + An MP variable occupies a varying number of entries in its array. We keep + track of this number for efficiency reasons. Otherwise we would always + have to process the whole array. */ +#define MPN_VAR(name) mp_limb_t *name; mp_size_t name##size + +#define MPN_ASSIGN(dst,src) \ + memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t)) +#define MPN_GE(u,v) \ + (u##size > v##size || (u##size == v##size && mpn_cmp (u, v, u##size) >= 0)) + +extern mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, + int *expt, int *is_neg, + __float128 value) attribute_hidden; +static unsigned int guess_grouping (unsigned int intdig_max, + const char *grouping); + + +static wchar_t *group_number (wchar_t *buf, wchar_t *bufend, + unsigned int intdig_no, const char *grouping, + wchar_t thousands_sep, int ngroups); + + +int +__quadmath_printf_fp (struct __quadmath_printf_file *fp, + const struct printf_info *info, + const void *const *args) +{ + /* The floating-point value to output. */ + __float128 fpnum; + + /* Locale-dependent representation of decimal point. */ + const char *decimal; + wchar_t decimalwc; + + /* Locale-dependent thousands separator and grouping specification. */ + const char *thousands_sep = NULL; + wchar_t thousands_sepwc = L_('\0'); + const char *grouping; + + /* "NaN" or "Inf" for the special cases. */ + const char *special = NULL; + const wchar_t *wspecial = NULL; + + /* We need just a few limbs for the input before shifting to the right + position. */ + mp_limb_t fp_input[(FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB]; + /* We need to shift the contents of fp_input by this amount of bits. */ + int to_shift = 0; + + /* The fraction of the floting-point value in question */ + MPN_VAR(frac); + /* and the exponent. */ + int exponent; + /* Sign of the exponent. */ + int expsign = 0; + /* Sign of float number. */ + int is_neg = 0; + + /* Scaling factor. */ + MPN_VAR(scale); + + /* Temporary bignum value. */ + MPN_VAR(tmp); + + /* Digit which is result of last hack_digit() call. */ + wchar_t digit; + + /* The type of output format that will be used: 'e'/'E' or 'f'. */ + int type; + + /* Counter for number of written characters. */ + int done = 0; + + /* General helper (carry limb). */ + mp_limb_t cy; + + /* Nonzero if this is output on a wide character stream. */ + int wide = info->wide; + + /* Buffer in which we produce the output. */ + wchar_t *wbuffer = NULL; + /* Flag whether wbuffer is malloc'ed or not. */ + int buffer_malloced = 0; + + auto wchar_t hack_digit (void); + + wchar_t hack_digit (void) + { + mp_limb_t hi; + + if (expsign != 0 && type == 'f' && exponent-- > 0) + hi = 0; + else if (scalesize == 0) + { + hi = frac[fracsize - 1]; + frac[fracsize - 1] = mpn_mul_1 (frac, frac, fracsize - 1, 10); + } + else + { + if (fracsize < scalesize) + hi = 0; + else + { + hi = mpn_divmod (tmp, frac, fracsize, scale, scalesize); + tmp[fracsize - scalesize] = hi; + hi = tmp[0]; + + fracsize = scalesize; + while (fracsize != 0 && frac[fracsize - 1] == 0) + --fracsize; + if (fracsize == 0) + { + /* We're not prepared for an mpn variable with zero + limbs. */ + fracsize = 1; + return L_('0') + hi; + } + } + + mp_limb_t _cy = mpn_mul_1 (frac, frac, fracsize, 10); + if (_cy != 0) + frac[fracsize++] = _cy; + } + + return L_('0') + hi; + } + + /* Figure out the decimal point character. */ +#ifdef USE_NL_LANGINFO + if (info->extra == 0) + decimal = nl_langinfo (DECIMAL_POINT); + else + { + decimal = nl_langinfo (MON_DECIMAL_POINT); + if (*decimal == '\0') + decimal = nl_langinfo (DECIMAL_POINT); + } + /* The decimal point character must never be zero. */ + assert (*decimal != '\0'); +#elif defined USE_LOCALECONV + const struct lconv *lc = localeconv (); + if (info->extra == 0) + decimal = lc->decimal_point; + else + { + decimal = lc->mon_decimal_point; + if (decimal == NULL || *decimal == '\0') + decimal = lc->decimal_point; + } + if (decimal == NULL || *decimal == '\0') + decimal = "."; +#else + decimal = "."; +#endif +#ifdef USE_NL_LANGINFO_WC + if (info->extra == 0) + decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); + else + { + decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC); + if (decimalwc == L_('\0')) + decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); + } + /* The decimal point character must never be zero. */ + assert (decimalwc != L_('\0')); +#else + decimalwc = L_('.'); +#endif + +#if defined USE_NL_LANGINFO && defined USE_NL_LANGINFO_WC + if (info->group) + { + if (info->extra == 0) + grouping = nl_langinfo (GROUPING); + else + grouping = nl_langinfo (MON_GROUPING); + + if (*grouping <= 0 || *grouping == CHAR_MAX) + grouping = NULL; + else + { + /* Figure out the thousands separator character. */ + if (wide) + { + if (info->extra == 0) + thousands_sepwc = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC); + else + thousands_sepwc = nl_langinfo_wc (_NL_MONETARY_THOUSANDS_SEP_WC); + + if (thousands_sepwc == L_('\0')) + grouping = NULL; + } + else + { + if (info->extra == 0) + thousands_sep = nl_langinfo (THOUSANDS_SEP); + else + thousands_sep = nl_langinfo (MON_THOUSANDS_SEP); + if (*thousands_sep == '\0') + grouping = NULL; + } + } + } + else +#elif defined USE_NL_LANGINFO + if (info->group && !wide) + { + if (info->extra == 0) + grouping = nl_langinfo (GROUPING); + else + grouping = nl_langinfo (MON_GROUPING); + + if (*grouping <= 0 || *grouping == CHAR_MAX) + grouping = NULL; + else + { + /* Figure out the thousands separator character. */ + if (info->extra == 0) + thousands_sep = nl_langinfo (THOUSANDS_SEP); + else + thousands_sep = nl_langinfo (MON_THOUSANDS_SEP); + + if (*thousands_sep == '\0') + grouping = NULL; + } + } + else +#elif defined USE_LOCALECONV + if (info->group && !wide) + { + if (info->extra == 0) + grouping = lc->grouping; + else + grouping = lc->mon_grouping; + + if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX) + grouping = NULL; + else + { + /* Figure out the thousands separator character. */ + if (info->extra == 0) + thousands_sep = lc->thousands_sep; + else + thousands_sep = lc->mon_thousands_sep; + + if (thousands_sep == NULL || *thousands_sep == '\0') + grouping = NULL; + } + } + else +#endif + grouping = NULL; + if (grouping != NULL && !wide) + /* If we are printing multibyte characters and there is a + multibyte representation for the thousands separator, + we must ensure the wide character thousands separator + is available, even if it is fake. */ + thousands_sepwc = (wchar_t) 0xfffffffe; + + /* Fetch the argument value. */ + { + fpnum = **(const __float128 **) args[0]; + + /* Check for special values: not a number or infinity. */ + if (isnanq (fpnum)) + { + ieee854_float128 u = { .value = fpnum }; + is_neg = u.ieee.negative != 0; + if (isupper (info->spec)) + { + special = "NAN"; + wspecial = L_("NAN"); + } + else + { + special = "nan"; + wspecial = L_("nan"); + } + } + else if (isinfq (fpnum)) + { + is_neg = fpnum < 0; + if (isupper (info->spec)) + { + special = "INF"; + wspecial = L_("INF"); + } + else + { + special = "inf"; + wspecial = L_("inf"); + } + } + else + { + fracsize = mpn_extract_flt128 (fp_input, + (sizeof (fp_input) / + sizeof (fp_input[0])), + &exponent, &is_neg, fpnum); + to_shift = 1 + fracsize * BITS_PER_MP_LIMB - FLT128_MANT_DIG; + } + } + + if (special) + { + int width = info->width; + + if (is_neg || info->showsign || info->space) + --width; + width -= 3; + + if (!info->left && width > 0) + PADN (' ', width); + + if (is_neg) + outchar ('-'); + else if (info->showsign) + outchar ('+'); + else if (info->space) + outchar (' '); + + PRINT (special, wspecial, 3); + + if (info->left && width > 0) + PADN (' ', width); + + return done; + } + + + /* We need three multiprecision variables. Now that we have the exponent + of the number we can allocate the needed memory. It would be more + efficient to use variables of the fixed maximum size but because this + would be really big it could lead to memory problems. */ + { + mp_size_t bignum_size = ((ABS (exponent) + BITS_PER_MP_LIMB - 1) + / BITS_PER_MP_LIMB + + (FLT128_MANT_DIG / BITS_PER_MP_LIMB > 2 ? 8 : 4)) + * sizeof (mp_limb_t); + frac = (mp_limb_t *) alloca (bignum_size); + tmp = (mp_limb_t *) alloca (bignum_size); + scale = (mp_limb_t *) alloca (bignum_size); + } + + /* We now have to distinguish between numbers with positive and negative + exponents because the method used for the one is not applicable/efficient + for the other. */ + scalesize = 0; + if (exponent > 2) + { + /* |FP| >= 8.0. */ + int scaleexpo = 0; + int explog = FLT128_MAX_10_EXP_LOG; + int exp10 = 0; + const struct mp_power *powers = &_fpioconst_pow10[explog + 1]; + int cnt_h, cnt_l, i; + + if ((exponent + to_shift) % BITS_PER_MP_LIMB == 0) + { + MPN_COPY_DECR (frac + (exponent + to_shift) / BITS_PER_MP_LIMB, + fp_input, fracsize); + fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB; + } + else + { + cy = mpn_lshift (frac + (exponent + to_shift) / BITS_PER_MP_LIMB, + fp_input, fracsize, + (exponent + to_shift) % BITS_PER_MP_LIMB); + fracsize += (exponent + to_shift) / BITS_PER_MP_LIMB; + if (cy) + frac[fracsize++] = cy; + } + MPN_ZERO (frac, (exponent + to_shift) / BITS_PER_MP_LIMB); + + assert (powers > &_fpioconst_pow10[0]); + do + { + --powers; + + /* The number of the product of two binary numbers with n and m + bits respectively has m+n or m+n-1 bits. */ + if (exponent >= scaleexpo + powers->p_expo - 1) + { + if (scalesize == 0) + { + if (FLT128_MANT_DIG > _FPIO_CONST_OFFSET * BITS_PER_MP_LIMB) + { +#define _FPIO_CONST_SHIFT \ + (((FLT128_MANT_DIG + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) \ + - _FPIO_CONST_OFFSET) + /* 64bit const offset is not enough for + IEEE quad long double. */ + tmpsize = powers->arraysize + _FPIO_CONST_SHIFT; + memcpy (tmp + _FPIO_CONST_SHIFT, + &__tens[powers->arrayoff], + tmpsize * sizeof (mp_limb_t)); + MPN_ZERO (tmp, _FPIO_CONST_SHIFT); + /* Adjust exponent, as scaleexpo will be this much + bigger too. */ + exponent += _FPIO_CONST_SHIFT * BITS_PER_MP_LIMB; + } + else + { + tmpsize = powers->arraysize; + memcpy (tmp, &__tens[powers->arrayoff], + tmpsize * sizeof (mp_limb_t)); + } + } + else + { + cy = mpn_mul (tmp, scale, scalesize, + &__tens[powers->arrayoff + + _FPIO_CONST_OFFSET], + powers->arraysize - _FPIO_CONST_OFFSET); + tmpsize = scalesize + powers->arraysize - _FPIO_CONST_OFFSET; + if (cy == 0) + --tmpsize; + } + + if (MPN_GE (frac, tmp)) + { + int cnt; + MPN_ASSIGN (scale, tmp); + count_leading_zeros (cnt, scale[scalesize - 1]); + scaleexpo = (scalesize - 2) * BITS_PER_MP_LIMB - cnt - 1; + exp10 |= 1 << explog; + } + } + --explog; + } + while (powers > &_fpioconst_pow10[0]); + exponent = exp10; + + /* Optimize number representations. We want to represent the numbers + with the lowest number of bytes possible without losing any + bytes. Also the highest bit in the scaling factor has to be set + (this is a requirement of the MPN division routines). */ + if (scalesize > 0) + { + /* Determine minimum number of zero bits at the end of + both numbers. */ + for (i = 0; scale[i] == 0 && frac[i] == 0; i++) + ; + + /* Determine number of bits the scaling factor is misplaced. */ + count_leading_zeros (cnt_h, scale[scalesize - 1]); + + if (cnt_h == 0) + { + /* The highest bit of the scaling factor is already set. So + we only have to remove the trailing empty limbs. */ + if (i > 0) + { + MPN_COPY_INCR (scale, scale + i, scalesize - i); + scalesize -= i; + MPN_COPY_INCR (frac, frac + i, fracsize - i); + fracsize -= i; + } + } + else + { + if (scale[i] != 0) + { + count_trailing_zeros (cnt_l, scale[i]); + if (frac[i] != 0) + { + int cnt_l2; + count_trailing_zeros (cnt_l2, frac[i]); + if (cnt_l2 < cnt_l) + cnt_l = cnt_l2; + } + } + else + count_trailing_zeros (cnt_l, frac[i]); + + /* Now shift the numbers to their optimal position. */ + if (i == 0 && BITS_PER_MP_LIMB - cnt_h > cnt_l) + { + /* We cannot save any memory. So just roll both numbers + so that the scaling factor has its highest bit set. */ + + (void) mpn_lshift (scale, scale, scalesize, cnt_h); + cy = mpn_lshift (frac, frac, fracsize, cnt_h); + if (cy != 0) + frac[fracsize++] = cy; + } + else if (BITS_PER_MP_LIMB - cnt_h <= cnt_l) + { + /* We can save memory by removing the trailing zero limbs + and by packing the non-zero limbs which gain another + free one. */ + + (void) mpn_rshift (scale, scale + i, scalesize - i, + BITS_PER_MP_LIMB - cnt_h); + scalesize -= i + 1; + (void) mpn_rshift (frac, frac + i, fracsize - i, + BITS_PER_MP_LIMB - cnt_h); + fracsize -= frac[fracsize - i - 1] == 0 ? i + 1 : i; + } + else + { + /* We can only save the memory of the limbs which are zero. + The non-zero parts occupy the same number of limbs. */ + + (void) mpn_rshift (scale, scale + (i - 1), + scalesize - (i - 1), + BITS_PER_MP_LIMB - cnt_h); + scalesize -= i; + (void) mpn_rshift (frac, frac + (i - 1), + fracsize - (i - 1), + BITS_PER_MP_LIMB - cnt_h); + fracsize -= frac[fracsize - (i - 1) - 1] == 0 ? i : i - 1; + } + } + } + } + else if (exponent < 0) + { + /* |FP| < 1.0. */ + int exp10 = 0; + int explog = FLT128_MAX_10_EXP_LOG; + const struct mp_power *powers = &_fpioconst_pow10[explog + 1]; + + /* Now shift the input value to its right place. */ + cy = mpn_lshift (frac, fp_input, fracsize, to_shift); + frac[fracsize++] = cy; + assert (cy == 1 || (frac[fracsize - 2] == 0 && frac[0] == 0)); + + expsign = 1; + exponent = -exponent; + + assert (powers != &_fpioconst_pow10[0]); + do + { + --powers; + + if (exponent >= powers->m_expo) + { + int i, incr, cnt_h, cnt_l; + mp_limb_t topval[2]; + + /* The mpn_mul function expects the first argument to be + bigger than the second. */ + if (fracsize < powers->arraysize - _FPIO_CONST_OFFSET) + cy = mpn_mul (tmp, &__tens[powers->arrayoff + + _FPIO_CONST_OFFSET], + powers->arraysize - _FPIO_CONST_OFFSET, + frac, fracsize); + else + cy = mpn_mul (tmp, frac, fracsize, + &__tens[powers->arrayoff + _FPIO_CONST_OFFSET], + powers->arraysize - _FPIO_CONST_OFFSET); + tmpsize = fracsize + powers->arraysize - _FPIO_CONST_OFFSET; + if (cy == 0) + --tmpsize; + + count_leading_zeros (cnt_h, tmp[tmpsize - 1]); + incr = (tmpsize - fracsize) * BITS_PER_MP_LIMB + + BITS_PER_MP_LIMB - 1 - cnt_h; + + assert (incr <= powers->p_expo); + + /* If we increased the exponent by exactly 3 we have to test + for overflow. This is done by comparing with 10 shifted + to the right position. */ + if (incr == exponent + 3) + { + if (cnt_h <= BITS_PER_MP_LIMB - 4) + { + topval[0] = 0; + topval[1] + = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4 - cnt_h); + } + else + { + topval[0] = ((mp_limb_t) 10) << (BITS_PER_MP_LIMB - 4); + topval[1] = 0; + (void) mpn_lshift (topval, topval, 2, + BITS_PER_MP_LIMB - cnt_h); + } + } + + /* We have to be careful when multiplying the last factor. + If the result is greater than 1.0 be have to test it + against 10.0. If it is greater or equal to 10.0 the + multiplication was not valid. This is because we cannot + determine the number of bits in the result in advance. */ + if (incr < exponent + 3 + || (incr == exponent + 3 && + (tmp[tmpsize - 1] < topval[1] + || (tmp[tmpsize - 1] == topval[1] + && tmp[tmpsize - 2] < topval[0])))) + { + /* The factor is right. Adapt binary and decimal + exponents. */ + exponent -= incr; + exp10 |= 1 << explog; + + /* If this factor yields a number greater or equal to + 1.0, we must not shift the non-fractional digits down. */ + if (exponent < 0) + cnt_h += -exponent; + + /* Now we optimize the number representation. */ + for (i = 0; tmp[i] == 0; ++i); + if (cnt_h == BITS_PER_MP_LIMB - 1) + { + MPN_COPY (frac, tmp + i, tmpsize - i); + fracsize = tmpsize - i; + } + else + { + count_trailing_zeros (cnt_l, tmp[i]); + + /* Now shift the numbers to their optimal position. */ + if (i == 0 && BITS_PER_MP_LIMB - 1 - cnt_h > cnt_l) + { + /* We cannot save any memory. Just roll the + number so that the leading digit is in a + separate limb. */ + + cy = mpn_lshift (frac, tmp, tmpsize, cnt_h + 1); + fracsize = tmpsize + 1; + frac[fracsize - 1] = cy; + } + else if (BITS_PER_MP_LIMB - 1 - cnt_h <= cnt_l) + { + (void) mpn_rshift (frac, tmp + i, tmpsize - i, + BITS_PER_MP_LIMB - 1 - cnt_h); + fracsize = tmpsize - i; + } + else + { + /* We can only save the memory of the limbs which + are zero. The non-zero parts occupy the same + number of limbs. */ + + (void) mpn_rshift (frac, tmp + (i - 1), + tmpsize - (i - 1), + BITS_PER_MP_LIMB - 1 - cnt_h); + fracsize = tmpsize - (i - 1); + } + } + } + } + --explog; + } + while (powers != &_fpioconst_pow10[1] && exponent > 0); + /* All factors but 10^-1 are tested now. */ + if (exponent > 0) + { + int cnt_l; + + cy = mpn_mul_1 (tmp, frac, fracsize, 10); + tmpsize = fracsize; + assert (cy == 0 || tmp[tmpsize - 1] < 20); + + count_trailing_zeros (cnt_l, tmp[0]); + if (cnt_l < MIN (4, exponent)) + { + cy = mpn_lshift (frac, tmp, tmpsize, + BITS_PER_MP_LIMB - MIN (4, exponent)); + if (cy != 0) + frac[tmpsize++] = cy; + } + else + (void) mpn_rshift (frac, tmp, tmpsize, MIN (4, exponent)); + fracsize = tmpsize; + exp10 |= 1; + assert (frac[fracsize - 1] < 10); + } + exponent = exp10; + } + else + { + /* This is a special case. We don't need a factor because the + numbers are in the range of 1.0 <= |fp| < 8.0. We simply + shift it to the right place and divide it by 1.0 to get the + leading digit. (Of course this division is not really made.) */ + assert (0 <= exponent && exponent < 3 && + exponent + to_shift < BITS_PER_MP_LIMB); + + /* Now shift the input value to its right place. */ + cy = mpn_lshift (frac, fp_input, fracsize, (exponent + to_shift)); + frac[fracsize++] = cy; + exponent = 0; + } + + { + int width = info->width; + wchar_t *wstartp, *wcp; + size_t chars_needed; + int expscale; + int intdig_max, intdig_no = 0; + int fracdig_min; + int fracdig_max; + int dig_max; + int significant; + int ngroups = 0; + char spec = tolower (info->spec); + + if (spec == 'e') + { + type = info->spec; + intdig_max = 1; + fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec; + chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4; + /* d . ddd e +- ddd */ + dig_max = __INT_MAX__; /* Unlimited. */ + significant = 1; /* Does not matter here. */ + } + else if (spec == 'f') + { + type = 'f'; + fracdig_min = fracdig_max = info->prec < 0 ? 6 : info->prec; + dig_max = __INT_MAX__; /* Unlimited. */ + significant = 1; /* Does not matter here. */ + if (expsign == 0) + { + intdig_max = exponent + 1; + /* This can be really big! */ /* XXX Maybe malloc if too big? */ + chars_needed = (size_t) exponent + 1 + 1 + (size_t) fracdig_max; + } + else + { + intdig_max = 1; + chars_needed = 1 + 1 + (size_t) fracdig_max; + } + } + else + { + dig_max = info->prec < 0 ? 6 : (info->prec == 0 ? 1 : info->prec); + if ((expsign == 0 && exponent >= dig_max) + || (expsign != 0 && exponent > 4)) + { + if ('g' - 'G' == 'e' - 'E') + type = 'E' + (info->spec - 'G'); + else + type = isupper (info->spec) ? 'E' : 'e'; + fracdig_max = dig_max - 1; + intdig_max = 1; + chars_needed = 1 + 1 + (size_t) fracdig_max + 1 + 1 + 4; + } + else + { + type = 'f'; + intdig_max = expsign == 0 ? exponent + 1 : 0; + fracdig_max = dig_max - intdig_max; + /* We need space for the significant digits and perhaps + for leading zeros when < 1.0. The number of leading + zeros can be as many as would be required for + exponential notation with a negative two-digit + exponent, which is 4. */ + chars_needed = (size_t) dig_max + 1 + 4; + } + fracdig_min = info->alt ? fracdig_max : 0; + significant = 0; /* We count significant digits. */ + } + + if (grouping) + { + /* Guess the number of groups we will make, and thus how + many spaces we need for separator characters. */ + ngroups = guess_grouping (intdig_max, grouping); + /* Allocate one more character in case rounding increases the + number of groups. */ + chars_needed += ngroups + 1; + } + + /* Allocate buffer for output. We need two more because while rounding + it is possible that we need two more characters in front of all the + other output. If the amount of memory we have to allocate is too + large use `malloc' instead of `alloca'. */ + if (__builtin_expect (chars_needed >= (size_t) -1 / sizeof (wchar_t) - 2 + || chars_needed < fracdig_max, 0)) + { + /* Some overflow occurred. */ +#if defined HAVE_ERRNO_H && defined ERANGE + errno = ERANGE; +#endif + return -1; + } + size_t wbuffer_to_alloc = (2 + chars_needed) * sizeof (wchar_t); + buffer_malloced = wbuffer_to_alloc >= 4096; + if (__builtin_expect (buffer_malloced, 0)) + { + wbuffer = (wchar_t *) malloc (wbuffer_to_alloc); + if (wbuffer == NULL) + /* Signal an error to the caller. */ + return -1; + } + else + wbuffer = (wchar_t *) alloca (wbuffer_to_alloc); + wcp = wstartp = wbuffer + 2; /* Let room for rounding. */ + + /* Do the real work: put digits in allocated buffer. */ + if (expsign == 0 || type != 'f') + { + assert (expsign == 0 || intdig_max == 1); + while (intdig_no < intdig_max) + { + ++intdig_no; + *wcp++ = hack_digit (); + } + significant = 1; + if (info->alt + || fracdig_min > 0 + || (fracdig_max > 0 && (fracsize > 1 || frac[0] != 0))) + *wcp++ = decimalwc; + } + else + { + /* |fp| < 1.0 and the selected type is 'f', so put "0." + in the buffer. */ + *wcp++ = L_('0'); + --exponent; + *wcp++ = decimalwc; + } + + /* Generate the needed number of fractional digits. */ + int fracdig_no = 0; + int added_zeros = 0; + while (fracdig_no < fracdig_min + added_zeros + || (fracdig_no < fracdig_max && (fracsize > 1 || frac[0] != 0))) + { + ++fracdig_no; + *wcp = hack_digit (); + if (*wcp++ != L_('0')) + significant = 1; + else if (significant == 0) + { + ++fracdig_max; + if (fracdig_min > 0) + ++added_zeros; + } + } + + /* Do rounding. */ + digit = hack_digit (); + if (digit > L_('4')) + { + wchar_t *wtp = wcp; + + if (digit == L_('5') + && ((*(wcp - 1) != decimalwc && (*(wcp - 1) & 1) == 0) + || ((*(wcp - 1) == decimalwc && (*(wcp - 2) & 1) == 0)))) + { + /* This is the critical case. */ + if (fracsize == 1 && frac[0] == 0) + /* Rest of the number is zero -> round to even. + (IEEE 754-1985 4.1 says this is the default rounding.) */ + goto do_expo; + else if (scalesize == 0) + { + /* Here we have to see whether all limbs are zero since no + normalization happened. */ + size_t lcnt = fracsize; + while (lcnt >= 1 && frac[lcnt - 1] == 0) + --lcnt; + if (lcnt == 0) + /* Rest of the number is zero -> round to even. + (IEEE 754-1985 4.1 says this is the default rounding.) */ + goto do_expo; + } + } + + if (fracdig_no > 0) + { + /* Process fractional digits. Terminate if not rounded or + radix character is reached. */ + int removed = 0; + while (*--wtp != decimalwc && *wtp == L_('9')) + { + *wtp = L_('0'); + ++removed; + } + if (removed == fracdig_min && added_zeros > 0) + --added_zeros; + if (*wtp != decimalwc) + /* Round up. */ + (*wtp)++; + else if (__builtin_expect (spec == 'g' && type == 'f' && info->alt + && wtp == wstartp + 1 + && wstartp[0] == L_('0'), + 0)) + /* This is a special case: the rounded number is 1.0, + the format is 'g' or 'G', and the alternative format + is selected. This means the result must be "1.". */ + --added_zeros; + } + + if (fracdig_no == 0 || *wtp == decimalwc) + { + /* Round the integer digits. */ + if (*(wtp - 1) == decimalwc) + --wtp; + + while (--wtp >= wstartp && *wtp == L_('9')) + *wtp = L_('0'); + + if (wtp >= wstartp) + /* Round up. */ + (*wtp)++; + else + /* It is more critical. All digits were 9's. */ + { + if (type != 'f') + { + *wstartp = '1'; + exponent += expsign == 0 ? 1 : -1; + + /* The above exponent adjustment could lead to 1.0e-00, + e.g. for 0.999999999. Make sure exponent 0 always + uses + sign. */ + if (exponent == 0) + expsign = 0; + } + else if (intdig_no == dig_max) + { + /* This is the case where for type %g the number fits + really in the range for %f output but after rounding + the number of digits is too big. */ + *--wstartp = decimalwc; + *--wstartp = L_('1'); + + if (info->alt || fracdig_no > 0) + { + /* Overwrite the old radix character. */ + wstartp[intdig_no + 2] = L_('0'); + ++fracdig_no; + } + + fracdig_no += intdig_no; + intdig_no = 1; + fracdig_max = intdig_max - intdig_no; + ++exponent; + /* Now we must print the exponent. */ + type = isupper (info->spec) ? 'E' : 'e'; + } + else + { + /* We can simply add another another digit before the + radix. */ + *--wstartp = L_('1'); + ++intdig_no; + } + + /* While rounding the number of digits can change. + If the number now exceeds the limits remove some + fractional digits. */ + if (intdig_no + fracdig_no > dig_max) + { + wcp -= intdig_no + fracdig_no - dig_max; + fracdig_no -= intdig_no + fracdig_no - dig_max; + } + } + } + } + + do_expo: + /* Now remove unnecessary '0' at the end of the string. */ + while (fracdig_no > fracdig_min + added_zeros && *(wcp - 1) == L_('0')) + { + --wcp; + --fracdig_no; + } + /* If we eliminate all fractional digits we perhaps also can remove + the radix character. */ + if (fracdig_no == 0 && !info->alt && *(wcp - 1) == decimalwc) + --wcp; + + if (grouping) + { + /* Rounding might have changed the number of groups. We allocated + enough memory but we need here the correct number of groups. */ + if (intdig_no != intdig_max) + ngroups = guess_grouping (intdig_no, grouping); + + /* Add in separator characters, overwriting the same buffer. */ + wcp = group_number (wstartp, wcp, intdig_no, grouping, thousands_sepwc, + ngroups); + } + + /* Write the exponent if it is needed. */ + if (type != 'f') + { + if (__builtin_expect (expsign != 0 && exponent == 4 && spec == 'g', 0)) + { + /* This is another special case. The exponent of the number is + really smaller than -4, which requires the 'e'/'E' format. + But after rounding the number has an exponent of -4. */ + assert (wcp >= wstartp + 1); + assert (wstartp[0] == L_('1')); + memcpy (wstartp, L_("0.0001"), 6 * sizeof (wchar_t)); + wstartp[1] = decimalwc; + if (wcp >= wstartp + 2) + { + size_t cnt; + for (cnt = 0; cnt < wcp - (wstartp + 2); cnt++) + wstartp[6 + cnt] = L_('0'); + wcp += 4; + } + else + wcp += 5; + } + else + { + *wcp++ = (wchar_t) type; + *wcp++ = expsign ? L_('-') : L_('+'); + + /* Find the magnitude of the exponent. */ + expscale = 10; + while (expscale <= exponent) + expscale *= 10; + + if (exponent < 10) + /* Exponent always has at least two digits. */ + *wcp++ = L_('0'); + else + do + { + expscale /= 10; + *wcp++ = L_('0') + (exponent / expscale); + exponent %= expscale; + } + while (expscale > 10); + *wcp++ = L_('0') + exponent; + } + } + + /* Compute number of characters which must be filled with the padding + character. */ + if (is_neg || info->showsign || info->space) + --width; + width -= wcp - wstartp; + + if (!info->left && info->pad != '0' && width > 0) + PADN (info->pad, width); + + if (is_neg) + outchar ('-'); + else if (info->showsign) + outchar ('+'); + else if (info->space) + outchar (' '); + + if (!info->left && info->pad == '0' && width > 0) + PADN ('0', width); + + { + char *buffer = NULL; + char *buffer_end __attribute__((__unused__)) = NULL; + char *cp = NULL; + char *tmpptr; + + if (! wide) + { + /* Create the single byte string. */ + size_t decimal_len; + size_t thousands_sep_len; + wchar_t *copywc; +#ifdef USE_I18N_NUMBER_H + size_t factor = (info->i18n + ? nl_langinfo_wc (_NL_CTYPE_MB_CUR_MAX) + : 1); +#else + size_t factor = 1; +#endif + + decimal_len = strlen (decimal); + + if (thousands_sep == NULL) + thousands_sep_len = 0; + else + thousands_sep_len = strlen (thousands_sep); + + size_t nbuffer = (2 + chars_needed * factor + decimal_len + + ngroups * thousands_sep_len); + if (__builtin_expect (buffer_malloced, 0)) + { + buffer = (char *) malloc (nbuffer); + if (buffer == NULL) + { + /* Signal an error to the caller. */ + free (wbuffer); + return -1; + } + } + else + buffer = (char *) alloca (nbuffer); + buffer_end = buffer + nbuffer; + + /* Now copy the wide character string. Since the character + (except for the decimal point and thousands separator) must + be coming from the ASCII range we can esily convert the + string without mapping tables. */ + for (cp = buffer, copywc = wstartp; copywc < wcp; ++copywc) + if (*copywc == decimalwc) + memcpy (cp, decimal, decimal_len), cp += decimal_len; + else if (*copywc == thousands_sepwc) + memcpy (cp, thousands_sep, thousands_sep_len), cp += thousands_sep_len; + else + *cp++ = (char) *copywc; + } + + tmpptr = buffer; +#if USE_I18N_NUMBER_H + if (__builtin_expect (info->i18n, 0)) + { + tmpptr = _i18n_number_rewrite (tmpptr, cp, buffer_end); + cp = buffer_end; + assert ((uintptr_t) buffer <= (uintptr_t) tmpptr); + assert ((uintptr_t) tmpptr < (uintptr_t) buffer_end); + } +#endif + + PRINT (tmpptr, wstartp, wide ? wcp - wstartp : cp - tmpptr); + + /* Free the memory if necessary. */ + if (__builtin_expect (buffer_malloced, 0)) + { + free (buffer); + free (wbuffer); + } + } + + if (info->left && width > 0) + PADN (info->pad, width); + } + return done; +} + +/* Return the number of extra grouping characters that will be inserted + into a number with INTDIG_MAX integer digits. */ + +static unsigned int +guess_grouping (unsigned int intdig_max, const char *grouping) +{ + unsigned int groups; + + /* We treat all negative values like CHAR_MAX. */ + + if (*grouping == CHAR_MAX || *grouping <= 0) + /* No grouping should be done. */ + return 0; + + groups = 0; + while (intdig_max > (unsigned int) *grouping) + { + ++groups; + intdig_max -= *grouping++; + + if (*grouping == 0) + { + /* Same grouping repeats. */ + groups += (intdig_max - 1) / grouping[-1]; + break; + } + else if (*grouping == CHAR_MAX || *grouping <= 0) + /* No more grouping should be done. */ + break; + } + + return groups; +} + +/* Group the INTDIG_NO integer digits of the number in [BUF,BUFEND). + There is guaranteed enough space past BUFEND to extend it. + Return the new end of buffer. */ + +static wchar_t * +group_number (wchar_t *buf, wchar_t *bufend, unsigned int intdig_no, + const char *grouping, wchar_t thousands_sep, int ngroups) +{ + wchar_t *p; + + if (ngroups == 0) + return bufend; + + /* Move the fractional part down. */ + memmove (buf + intdig_no + ngroups, buf + intdig_no, + (bufend - (buf + intdig_no)) * sizeof (wchar_t)); + + p = buf + intdig_no + ngroups - 1; + do + { + unsigned int len = *grouping++; + do + *p-- = buf[--intdig_no]; + while (--len > 0); + *p-- = thousands_sep; + + if (*grouping == 0) + /* Same grouping repeats. */ + --grouping; + else if (*grouping == CHAR_MAX || *grouping <= 0) + /* No more grouping should be done. */ + break; + } while (intdig_no > (unsigned int) *grouping); + + /* Copy the remaining ungrouped digits. */ + do + *p-- = buf[--intdig_no]; + while (p > buf); + + return bufend + ngroups; +} diff --git a/libquadmath/printf/printf_fphex.c b/libquadmath/printf/printf_fphex.c new file mode 100644 index 000000000..941e93307 --- /dev/null +++ b/libquadmath/printf/printf_fphex.c @@ -0,0 +1,463 @@ +/* Print floating point number in hexadecimal notation according to ISO C99. + Copyright (C) 1997-2002,2004,2006 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + 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 +#include +#include +#include +#include +#define NDEBUG +#include +#include "quadmath-printf.h" +#include "_itoa.h" +#include "_itowa.h" + + +/* Macros for doing the actual output. */ + +#define outchar(ch) \ + do \ + { \ + register const int outc = (ch); \ + if (PUTC (outc, fp) == EOF) \ + return -1; \ + ++done; \ + } while (0) + +#define PRINT(ptr, wptr, len) \ + do \ + { \ + register size_t outlen = (len); \ + if (wide) \ + while (outlen-- > 0) \ + outchar (*wptr++); \ + else \ + while (outlen-- > 0) \ + outchar (*ptr++); \ + } while (0) + +#define PADN(ch, len) \ + do \ + { \ + if (PAD (fp, ch, len) != len) \ + return -1; \ + done += len; \ + } \ + while (0) + + + +int +__quadmath_printf_fphex (struct __quadmath_printf_file *fp, + const struct printf_info *info, + const void *const *args) +{ + /* The floating-point value to output. */ + ieee854_float128 fpnum; + + /* Locale-dependent representation of decimal point. */ + const char *decimal; + wchar_t decimalwc; + + /* "NaN" or "Inf" for the special cases. */ + const char *special = NULL; + const wchar_t *wspecial = NULL; + + /* Buffer for the generated number string for the mantissa. The + maximal size for the mantissa is 128 bits. */ + char numbuf[32]; + char *numstr; + char *numend; + wchar_t wnumbuf[32]; + wchar_t *wnumstr; + wchar_t *wnumend; + int negative; + + /* The maximal exponent of two in decimal notation has 5 digits. */ + char expbuf[5]; + char *expstr; + wchar_t wexpbuf[5]; + wchar_t *wexpstr; + int expnegative; + int exponent; + + /* Non-zero is mantissa is zero. */ + int zero_mantissa; + + /* The leading digit before the decimal point. */ + char leading; + + /* Precision. */ + int precision = info->prec; + + /* Width. */ + int width = info->width; + + /* Number of characters written. */ + int done = 0; + + /* Nonzero if this is output on a wide character stream. */ + int wide = info->wide; + + /* Figure out the decimal point character. */ +#ifdef USE_NL_LANGINFO + if (info->extra == 0) + decimal = nl_langinfo (DECIMAL_POINT); + else + { + decimal = nl_langinfo (MON_DECIMAL_POINT); + if (*decimal == '\0') + decimal = nl_langinfo (DECIMAL_POINT); + } + /* The decimal point character must never be zero. */ + assert (*decimal != '\0'); +#elif defined USE_LOCALECONV + const struct lconv *lc = localeconv (); + if (info->extra == 0) + decimal = lc->decimal_point; + else + { + decimal = lc->mon_decimal_point; + if (decimal == NULL || *decimal == '\0') + decimal = lc->decimal_point; + } + if (decimal == NULL || *decimal == '\0') + decimal = "."; +#else + decimal = "."; +#endif +#ifdef USE_NL_LANGINFO_WC + if (info->extra == 0) + decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); + else + { + decimalwc = nl_langinfo_wc (_NL_MONETARY_DECIMAL_POINT_WC); + if (decimalwc == L_('\0')) + decimalwc = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); + } + /* The decimal point character must never be zero. */ + assert (decimalwc != L_('\0')); +#else + decimalwc = L_('.'); +#endif + + /* Fetch the argument value. */ + { + fpnum.value = **(const __float128 **) args[0]; + + /* Check for special values: not a number or infinity. */ + if (isnanq (fpnum.value)) + { + negative = fpnum.ieee.negative != 0; + if (isupper (info->spec)) + { + special = "NAN"; + wspecial = L_("NAN"); + } + else + { + special = "nan"; + wspecial = L_("nan"); + } + } + else + { + if (isinfq (fpnum.value)) + { + if (isupper (info->spec)) + { + special = "INF"; + wspecial = L_("INF"); + } + else + { + special = "inf"; + wspecial = L_("inf"); + } + } + + negative = signbitq (fpnum.value); + } + } + + if (special) + { + int width = info->width; + + if (negative || info->showsign || info->space) + --width; + width -= 3; + + if (!info->left && width > 0) + PADN (' ', width); + + if (negative) + outchar ('-'); + else if (info->showsign) + outchar ('+'); + else if (info->space) + outchar (' '); + + PRINT (special, wspecial, 3); + + if (info->left && width > 0) + PADN (' ', width); + + return done; + } + + { + /* We have 112 bits of mantissa plus one implicit digit. Since + 112 bits are representable without rest using hexadecimal + digits we use only the implicit digits for the number before + the decimal point. */ + uint64_t num0, num1; + + assert (sizeof (long double) == 16); + + num0 = fpnum.ieee.mant_high; + num1 = fpnum.ieee.mant_low; + + zero_mantissa = (num0|num1) == 0; + + if (sizeof (unsigned long int) > 6) + { + numstr = _itoa_word (num1, numbuf + sizeof numbuf, 16, + info->spec == 'A'); + wnumstr = _itowa_word (num1, + wnumbuf + sizeof (wnumbuf) / sizeof (wchar_t), + 16, info->spec == 'A'); + } + else + { + numstr = _itoa (num1, numbuf + sizeof numbuf, 16, + info->spec == 'A'); + wnumstr = _itowa (num1, + wnumbuf + sizeof (wnumbuf) / sizeof (wchar_t), + 16, info->spec == 'A'); + } + + while (numstr > numbuf + (sizeof numbuf - 64 / 4)) + { + *--numstr = '0'; + *--wnumstr = L_('0'); + } + + if (sizeof (unsigned long int) > 6) + { + numstr = _itoa_word (num0, numstr, 16, info->spec == 'A'); + wnumstr = _itowa_word (num0, wnumstr, 16, info->spec == 'A'); + } + else + { + numstr = _itoa (num0, numstr, 16, info->spec == 'A'); + wnumstr = _itowa (num0, wnumstr, 16, info->spec == 'A'); + } + + /* Fill with zeroes. */ + while (numstr > numbuf + (sizeof numbuf - 112 / 4)) + { + *--numstr = '0'; + *--wnumstr = L_('0'); + } + + leading = fpnum.ieee.exponent == 0 ? '0' : '1'; + + exponent = fpnum.ieee.exponent; + + if (exponent == 0) + { + if (zero_mantissa) + expnegative = 0; + else + { + /* This is a denormalized number. */ + expnegative = 1; + exponent = IEEE854_FLOAT128_BIAS - 1; + } + } + else if (exponent >= IEEE854_FLOAT128_BIAS) + { + expnegative = 0; + exponent -= IEEE854_FLOAT128_BIAS; + } + else + { + expnegative = 1; + exponent = -(exponent - IEEE854_FLOAT128_BIAS); + } + } + + /* Look for trailing zeroes. */ + if (! zero_mantissa) + { + wnumend = &wnumbuf[sizeof wnumbuf / sizeof wnumbuf[0]]; + numend = &numbuf[sizeof numbuf / sizeof numbuf[0]]; + while (wnumend[-1] == L_('0')) + { + --wnumend; + --numend; + } + + if (precision == -1) + precision = numend - numstr; + else if (precision < numend - numstr + && (numstr[precision] > '8' + || (('A' < '0' || 'a' < '0') + && numstr[precision] < '0') + || (numstr[precision] == '8' + && (precision + 1 < numend - numstr + /* Round to even. */ + || (precision > 0 + && ((numstr[precision - 1] & 1) + ^ (isdigit (numstr[precision - 1]) == 0))) + || (precision == 0 + && ((leading & 1) + ^ (isdigit (leading) == 0))))))) + { + /* Round up. */ + int cnt = precision; + while (--cnt >= 0) + { + char ch = numstr[cnt]; + /* We assume that the digits and the letters are ordered + like in ASCII. This is true for the rest of GNU, too. */ + if (ch == '9') + { + wnumstr[cnt] = (wchar_t) info->spec; + numstr[cnt] = info->spec; /* This is tricky, + think about it! */ + break; + } + else if (tolower (ch) < 'f') + { + ++numstr[cnt]; + ++wnumstr[cnt]; + break; + } + else + { + numstr[cnt] = '0'; + wnumstr[cnt] = L_('0'); + } + } + if (cnt < 0) + { + /* The mantissa so far was fff...f Now increment the + leading digit. Here it is again possible that we + get an overflow. */ + if (leading == '9') + leading = info->spec; + else if (tolower (leading) < 'f') + ++leading; + else + { + leading = '1'; + if (expnegative) + { + exponent -= 4; + if (exponent <= 0) + { + exponent = -exponent; + expnegative = 0; + } + } + else + exponent += 4; + } + } + } + } + else + { + if (precision == -1) + precision = 0; + numend = numstr; + wnumend = wnumstr; + } + + /* Now we can compute the exponent string. */ + expstr = _itoa_word (exponent, expbuf + sizeof expbuf, 10, 0); + wexpstr = _itowa_word (exponent, + wexpbuf + sizeof wexpbuf / sizeof (wchar_t), 10, 0); + + /* Now we have all information to compute the size. */ + width -= ((negative || info->showsign || info->space) + /* Sign. */ + + 2 + 1 + 0 + precision + 1 + 1 + /* 0x h . hhh P ExpoSign. */ + + ((expbuf + sizeof expbuf) - expstr)); + /* Exponent. */ + + /* Count the decimal point. + A special case when the mantissa or the precision is zero and the `#' + is not given. In this case we must not print the decimal point. */ + if (precision > 0 || info->alt) + width -= wide ? 1 : strlen (decimal); + + if (!info->left && info->pad != '0' && width > 0) + PADN (' ', width); + + if (negative) + outchar ('-'); + else if (info->showsign) + outchar ('+'); + else if (info->space) + outchar (' '); + + outchar ('0'); + if ('X' - 'A' == 'x' - 'a') + outchar (info->spec + ('x' - 'a')); + else + outchar (info->spec == 'A' ? 'X' : 'x'); + + if (!info->left && info->pad == '0' && width > 0) + PADN ('0', width); + + outchar (leading); + + if (precision > 0 || info->alt) + { + const wchar_t *wtmp = &decimalwc; + PRINT (decimal, wtmp, wide ? 1 : strlen (decimal)); + } + + if (precision > 0) + { + ssize_t tofill = precision - (numend - numstr); + PRINT (numstr, wnumstr, MIN (numend - numstr, precision)); + if (tofill > 0) + PADN ('0', tofill); + } + + if ('P' - 'A' == 'p' - 'a') + outchar (info->spec + ('p' - 'a')); + else + outchar (info->spec == 'A' ? 'P' : 'p'); + + outchar (expnegative ? '-' : '+'); + + PRINT (expstr, wexpstr, (expbuf + sizeof expbuf) - expstr); + + if (info->left && info->pad != '0' && width > 0) + PADN (info->pad, width); + + return done; +} diff --git a/libquadmath/printf/quadmath-printf.c b/libquadmath/printf/quadmath-printf.c new file mode 100644 index 000000000..b70f432cc --- /dev/null +++ b/libquadmath/printf/quadmath-printf.c @@ -0,0 +1,422 @@ +/* GCC Quad-Precision Math Library + Copyright (C) 2011 Free Software Foundation, Inc. + Written by Jakub Jelinek + +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 +#include +#include +#include +#include "quadmath-printf.h" + +/* Read a simple integer from a string and update the string pointer. + It is assumed that the first character is a digit. */ +static unsigned int +read_int (const char **pstr) +{ + unsigned int retval = (unsigned char) **pstr - '0'; + + while (isdigit ((unsigned char) *++(*pstr))) + { + retval *= 10; + retval += (unsigned char) **pstr - '0'; + } + + return retval; +} + +#define PADSIZE 16 +static char const blanks[PADSIZE] = +{' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '}; +static char const zeroes[PADSIZE] = +{'0','0','0','0','0','0','0','0','0','0','0','0','0','0','0','0'}; +static wchar_t const wblanks[PADSIZE] = +{ + L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), + L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' '), L_(' ') +}; +static wchar_t const wzeroes[PADSIZE] = +{ + L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), + L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0'), L_('0') +}; + +attribute_hidden size_t +__quadmath_do_pad (struct __quadmath_printf_file *fp, int wide, int c, + size_t n) +{ + ssize_t i; + char padbuf[PADSIZE]; + wchar_t wpadbuf[PADSIZE]; + const char *padstr; + size_t w, written = 0; + if (wide) + { + if (c == ' ') + padstr = (const char *) wblanks; + else if (c == '0') + padstr = (const char *) wzeroes; + else + { + padstr = (const char *) wpadbuf; + for (i = 0; i < PADSIZE; i++) + wpadbuf[i] = c; + } + } + else + { + if (c == ' ') + padstr = blanks; + else if (c == '0') + padstr = zeroes; + else + { + padstr = (const char *) padbuf; + for (i = 0; i < PADSIZE; i++) + padbuf[i] = c; + } + } + for (i = n; i >= PADSIZE; i -= PADSIZE) + { + w = PUT (fp, (char *) padstr, PADSIZE); + written += w; + if (w != PADSIZE) + return written; + } + if (i > 0) + { + w = PUT (fp, (char *) padstr, i); + written += w; + } + return written; +} + +/* This is a stripped down version of snprintf, which just handles + a single %eEfFgGaA format entry with Q modifier. % has to be + the first character of the format string, no $ can be used. */ +int +quadmath_snprintf (char *str, size_t size, const char *format, ...) +{ + struct printf_info info; + va_list ap; + __float128 fpnum, *fpnum_addr = &fpnum, **fpnum_addr2 = &fpnum_addr; + struct __quadmath_printf_file qfp; + + if (*format++ != '%') + return -1; + + /* Clear information structure. */ + memset (&info, '\0', sizeof info); + /* info.alt = 0; + info.space = 0; + info.left = 0; + info.showsign = 0; + info.group = 0; + info.i18n = 0; + info.extra = 0; */ + info.pad = ' '; + /* info.wide = 0; */ + + /* Check for spec modifiers. */ + do + { + switch (*format) + { + case ' ': + /* Output a space in place of a sign, when there is no sign. */ + info.space = 1; + continue; + case '+': + /* Always output + or - for numbers. */ + info.showsign = 1; + continue; + case '-': + /* Left-justify things. */ + info.left = 1; + continue; + case '#': + /* Use the "alternate form": + Hex has 0x or 0X, FP always has a decimal point. */ + info.alt = 1; + continue; + case '0': + /* Pad with 0s. */ + info.pad = '0'; + continue; + case '\'': + /* Show grouping in numbers if the locale information + indicates any. */ + info.group = 1; + continue; + case 'I': + /* Use the internationalized form of the output. Currently + means to use the `outdigits' of the current locale. */ + info.i18n = 1; + continue; + default: + break; + } + break; + } + while (*++format); + + if (info.left) + info.pad = ' '; + + va_start (ap, format); + + /* Get the field width. */ + /* info.width = 0; */ + if (*format == '*') + { + /* The field width is given in an argument. + A negative field width indicates left justification. */ + ++format; + info.width = va_arg (ap, int); + } + else if (isdigit (*format)) + /* Constant width specification. */ + info.width = read_int (&format); + + /* Get the precision. */ + /* -1 means none given; 0 means explicit 0. */ + info.prec = -1; + if (*format == '.') + { + ++format; + if (*format == '*') + { + /* The precision is given in an argument. */ + ++format; + + info.prec = va_arg (ap, int); + } + else if (isdigit (*format)) + info.prec = read_int (&format); + else + /* "%.?" is treated like "%.0?". */ + info.prec = 0; + } + + /* Check for type modifiers. */ + /* info.is_long_double = 0; + info.is_short = 0; + info.is_long = 0; + info.is_char = 0; + info.user = 0; */ + + /* We require Q modifier. */ + if (*format++ != 'Q') + { + va_end (ap); + return -1; + } + + /* Get the format specification. */ + info.spec = (wchar_t) *format++; + if (info.spec == L_('\0') || *format != '\0') + { + va_end (ap); + return -1; + } + + switch (info.spec) + { + case L_('e'): + case L_('E'): + case L_('f'): + case L_('F'): + case L_('g'): + case L_('G'): + case L_('a'): + case L_('A'): + break; + default: + va_end (ap); + return -1; + } + + fpnum = va_arg (ap, __float128); + va_end (ap); + + qfp.fp = NULL; + qfp.str = str; + qfp.size = size ? size - 1 : 0; + qfp.len = 0; + qfp.file_p = 0; + + if (info.spec == L_('a') || info.spec == L_('A')) + __quadmath_printf_fphex (&qfp, &info, (const void *const *)&fpnum_addr2); + else + __quadmath_printf_fp (&qfp, &info, (const void *const *)&fpnum_addr2); + + if (size) + *qfp.str = '\0'; + + return qfp.len; +} + +#ifdef HAVE_PRINTF_HOOKS +static int pa_flt128; +int mod_Q attribute_hidden; + +static void +flt128_va (void *mem, va_list *ap) +{ + __float128 d = va_arg (*ap, __float128); + memcpy (mem, &d, sizeof (d)); +} + +static int +flt128_ais (const struct printf_info *info, size_t n __attribute__ ((unused)), + int *argtype, int *size) +{ + if (info->user & mod_Q) + { + argtype[0] = pa_flt128; + size[0] = sizeof (__float128); + return 1; + } +#if __GLIBC__ < 2 || (__GLIBC__ == 2 && __GLIBC_MINOR__ <= 13) + /* Workaround bug in glibc printf hook handling. */ + size[0] = -1; + switch (info->spec) + { + case L_('i'): + case L_('d'): + case L_('u'): + case L_('o'): + case L_('X'): + case L_('x'): +#if __LONG_MAX__ != __LONG_LONG_MAX__ + if (info->is_long_double) + argtype[0] = PA_INT|PA_FLAG_LONG_LONG; + else +#endif + if (info->is_long) + argtype[0] = PA_INT|PA_FLAG_LONG; + else if (info->is_short) + argtype[0] = PA_INT|PA_FLAG_SHORT; + else if (info->is_char) + argtype[0] = PA_CHAR; + else + argtype[0] = PA_INT; + return 1; + case L_('e'): + case L_('E'): + case L_('f'): + case L_('F'): + case L_('g'): + case L_('G'): + case L_('a'): + case L_('A'): + if (info->is_long_double) + argtype[0] = PA_DOUBLE|PA_FLAG_LONG_DOUBLE; + else + argtype[0] = PA_DOUBLE; + return 1; + case L_('c'): + argtype[0] = PA_CHAR; + return 1; + case L_('C'): + argtype[0] = PA_WCHAR; + return 1; + case L_('s'): + argtype[0] = PA_STRING; + return 1; + case L_('S'): + argtype[0] = PA_WSTRING; + return 1; + case L_('p'): + argtype[0] = PA_POINTER; + return 1; + case L_('n'): + argtype[0] = PA_INT|PA_FLAG_PTR; + return 1; + + case L_('m'): + default: + /* An unknown spec will consume no args. */ + return 0; + } +#endif + return -1; +} + +static int +flt128_printf_fp (FILE *fp, const struct printf_info *info, + const void *const *args) +{ + struct __quadmath_printf_file qpf + = { .fp = fp, .str = NULL, .size = 0, .len = 0, .file_p = 1 }; + + if ((info->user & mod_Q) == 0) + return -2; + + return __quadmath_printf_fp (&qpf, info, args); +} + +static int +flt128_printf_fphex (FILE *fp, const struct printf_info *info, + const void *const *args) +{ + struct __quadmath_printf_file qpf + = { .fp = fp, .str = NULL, .size = 0, .len = 0, .file_p = 1 }; + + if ((info->user & mod_Q) == 0) + return -2; + + return __quadmath_printf_fphex (&qpf, info, args); +} + +__attribute__((constructor)) static void +register_printf_flt128 (void) +{ + pa_flt128 = register_printf_type (flt128_va); + if (pa_flt128 == -1) + return; + mod_Q = register_printf_modifier (L_("Q")); + if (mod_Q == -1) + return; + register_printf_specifier ('f', flt128_printf_fp, flt128_ais); + register_printf_specifier ('F', flt128_printf_fp, flt128_ais); + register_printf_specifier ('e', flt128_printf_fp, flt128_ais); + register_printf_specifier ('E', flt128_printf_fp, flt128_ais); + register_printf_specifier ('g', flt128_printf_fp, flt128_ais); + register_printf_specifier ('G', flt128_printf_fp, flt128_ais); + register_printf_specifier ('a', flt128_printf_fphex, flt128_ais); + register_printf_specifier ('A', flt128_printf_fphex, flt128_ais); +} + +__attribute__((destructor)) static void +unregister_printf_flt128 (void) +{ + /* No way to unregister printf type and modifier currently, + and only one printf specifier can be registered right now. */ + if (pa_flt128 == -1 || mod_Q == -1) + return; + register_printf_specifier ('f', NULL, NULL); + register_printf_specifier ('F', NULL, NULL); + register_printf_specifier ('e', NULL, NULL); + register_printf_specifier ('E', NULL, NULL); + register_printf_specifier ('g', NULL, NULL); + register_printf_specifier ('G', NULL, NULL); + register_printf_specifier ('a', NULL, NULL); + register_printf_specifier ('A', NULL, NULL); +} +#endif diff --git a/libquadmath/printf/quadmath-printf.h b/libquadmath/printf/quadmath-printf.h new file mode 100644 index 000000000..32ebcec92 --- /dev/null +++ b/libquadmath/printf/quadmath-printf.h @@ -0,0 +1,186 @@ +/* GCC Quad-Precision Math Library + Copyright (C) 2011 Free Software Foundation, Inc. + Written by Jakub Jelinek + +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 +#include +#ifdef HAVE_LIMITS_H +#include +#endif +#ifdef HAVE_LANGINFO_H +#include +#endif +#ifdef HAVE_CTYPE_H +#include +#endif +#ifdef HAVE_WCHAR_H +#include +#endif +#ifdef HAVE_WCTYPE_H +#include +#endif +#ifdef HAVE_PRINTF_HOOKS +#include +#endif +#ifdef HAVE_LOCALE_H +#include +#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 isxdigit +#undef tolower +#define isupper(x) \ + ({__typeof(x) __is_x = (x); __is_x >= 'A' && __is_x <= 'Z'; }) +#define isdigit(x) \ + ({__typeof(x) __is_x = (x); __is_x >= '0' && __is_x <= '9'; }) +#define isxdigit(x) \ + ({__typeof(x) __is_x = (x); \ + (__is_x >= '0' && __is_x <= '9') \ + || ((x) >= 'A' && (x) <= 'F') \ + || ((x) >= 'a' && (x) <= 'f'); }) +#define tolower(x) \ + ({__typeof(x) __is_x = (x); \ + (__is_x >= 'A' && __is_x <= 'Z') ? __is_x - 'A' + 'a' : __is_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; }) + +#undef _itoa +#define _itoa __quadmath_itoa + +#undef NAN +#define NAN __builtin_nanf ("") diff --git a/libquadmath/printf/rshift.c b/libquadmath/printf/rshift.c new file mode 100644 index 000000000..17fd914d9 --- /dev/null +++ b/libquadmath/printf/rshift.c @@ -0,0 +1,88 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/sub_n.c b/libquadmath/printf/sub_n.c new file mode 100644 index 000000000..92e718731 --- /dev/null +++ b/libquadmath/printf/sub_n.c @@ -0,0 +1,62 @@ +/* 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 +#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; +} diff --git a/libquadmath/printf/submul_1.c b/libquadmath/printf/submul_1.c new file mode 100644 index 000000000..31903c628 --- /dev/null +++ b/libquadmath/printf/submul_1.c @@ -0,0 +1,64 @@ +/* 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 +#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; +} -- cgit v1.2.3