2001-03-12 00:04:52 +00:00
|
|
|
/*
|
|
|
|
* IBM Accurate Mathematical Library
|
2002-07-06 06:36:39 +00:00
|
|
|
* written by International Business Machines Corp.
|
2014-01-01 11:03:15 +00:00
|
|
|
* Copyright (C) 2001-2014 Free Software Foundation, Inc.
|
2001-03-12 00:04:52 +00:00
|
|
|
*
|
|
|
|
* This program is free software; you can redistribute it and/or modify
|
|
|
|
* it under the terms of the GNU Lesser General Public License as published by
|
2002-08-26 22:40:48 +00:00
|
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
2001-03-12 00:04:52 +00:00
|
|
|
* (at your option) any later version.
|
2001-03-12 07:57:09 +00:00
|
|
|
*
|
2001-03-12 00:04:52 +00:00
|
|
|
* This program 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
|
2002-08-20 21:51:55 +00:00
|
|
|
* GNU Lesser General Public License for more details.
|
2001-03-12 00:04:52 +00:00
|
|
|
*
|
|
|
|
* You should have received a copy of the GNU Lesser General Public License
|
2012-02-09 23:18:22 +00:00
|
|
|
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
2001-03-12 00:04:52 +00:00
|
|
|
*/
|
|
|
|
/****************************************************************************/
|
|
|
|
/* MODULE_NAME:mpsqrt.c */
|
|
|
|
/* */
|
|
|
|
/* FUNCTION:mpsqrt */
|
|
|
|
/* fastiroot */
|
|
|
|
/* */
|
|
|
|
/* FILES NEEDED:endian.h mpa.h mpsqrt.h */
|
|
|
|
/* mpa.c */
|
|
|
|
/* Multi-Precision square root function subroutine for precision p >= 4. */
|
|
|
|
/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
|
|
|
|
/* */
|
|
|
|
/****************************************************************************/
|
|
|
|
#include "endian.h"
|
|
|
|
#include "mpa.h"
|
|
|
|
|
2011-10-25 04:56:33 +00:00
|
|
|
#ifndef SECTION
|
|
|
|
# define SECTION
|
|
|
|
#endif
|
|
|
|
|
|
|
|
#include "mpsqrt.h"
|
|
|
|
|
2001-03-12 00:04:52 +00:00
|
|
|
/****************************************************************************/
|
|
|
|
/* Multi-Precision square root function subroutine for precision p >= 4. */
|
|
|
|
/* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */
|
|
|
|
/* Routine receives two pointers to Multi Precision numbers: */
|
|
|
|
/* x (left argument) and y (next argument). Routine also receives precision */
|
|
|
|
/* p as integer. Routine computes sqrt(*x) and stores result in *y */
|
|
|
|
/****************************************************************************/
|
|
|
|
|
2013-02-27 05:58:20 +00:00
|
|
|
static double fastiroot (double);
|
2001-03-12 00:04:52 +00:00
|
|
|
|
2011-10-25 04:56:33 +00:00
|
|
|
void
|
|
|
|
SECTION
|
2013-02-27 05:58:20 +00:00
|
|
|
__mpsqrt (mp_no *x, mp_no *y, int p)
|
|
|
|
{
|
|
|
|
int i, m, ey;
|
|
|
|
double dx, dy;
|
2013-03-21 07:37:44 +00:00
|
|
|
static const mp_no mphalf = {0, {1.0, HALFRAD}};
|
|
|
|
static const mp_no mp3halfs = {1, {1.0, 1.0, HALFRAD}};
|
2013-02-27 05:58:20 +00:00
|
|
|
mp_no mpxn, mpz, mpu, mpt1, mpt2;
|
2001-03-12 00:04:52 +00:00
|
|
|
|
2013-02-27 05:58:20 +00:00
|
|
|
ey = EX / 2;
|
|
|
|
__cpy (x, &mpxn, p);
|
|
|
|
mpxn.e -= (ey + ey);
|
|
|
|
__mp_dbl (&mpxn, &dx, p);
|
|
|
|
dy = fastiroot (dx);
|
|
|
|
__dbl_mp (dy, &mpu, p);
|
|
|
|
__mul (&mpxn, &mphalf, &mpz, p);
|
2001-03-12 00:04:52 +00:00
|
|
|
|
2013-02-27 05:58:20 +00:00
|
|
|
m = __mpsqrt_mp[p];
|
|
|
|
for (i = 0; i < m; i++)
|
|
|
|
{
|
|
|
|
__sqr (&mpu, &mpt1, p);
|
|
|
|
__mul (&mpt1, &mpz, &mpt2, p);
|
|
|
|
__sub (&mp3halfs, &mpt2, &mpt1, p);
|
|
|
|
__mul (&mpu, &mpt1, &mpt2, p);
|
|
|
|
__cpy (&mpt2, &mpu, p);
|
|
|
|
}
|
|
|
|
__mul (&mpxn, &mpu, y, p);
|
|
|
|
EY += ey;
|
2001-03-12 00:04:52 +00:00
|
|
|
}
|
|
|
|
|
|
|
|
/***********************************************************/
|
|
|
|
/* Compute a double precision approximation for 1/sqrt(x) */
|
|
|
|
/* with the relative error bounded by 2**-51. */
|
|
|
|
/***********************************************************/
|
2011-10-25 04:56:33 +00:00
|
|
|
static double
|
|
|
|
SECTION
|
2013-02-27 05:58:20 +00:00
|
|
|
fastiroot (double x)
|
|
|
|
{
|
|
|
|
union
|
|
|
|
{
|
|
|
|
int i[2];
|
|
|
|
double d;
|
|
|
|
} p, q;
|
|
|
|
double y, z, t;
|
Update.
2003-11-28 Ulrich Drepper <drepper@redhat.com>
* sysdeps/x86_64/fpu/libm-test-ulps: Add some more minor changes
to compensate other setup.
2003-11-27 Andreas Jaeger <aj@suse.de>
* sysdeps/x86_64/fpu/libm-test-ulps: Add ulps for new atan2 test.
* math/libm-test.inc (atan2_test): Add test that run infinitly.
Reported by "Willus" <etc231etc231@willus.com>.
2003-11-27 Michael Matz <matz@suse.de>
* sysdeps/ieee754/dbl-64/mpsqrt.c (fastiroot): Fix 64-bit problem
with wrong types.
2003-11-28 Jakub Jelinek <jakub@redhat.com>
* posix/regexec.c (acquire_init_state_context): Make inline.
Add always_inline attribute.
(check_matching): Add BE macro. Move if (cur_state->has_backref)
into if (dfa->nbackref).
(sift_states_backward): Fix comment.
(transit_state): Add BE macro. Move if (next_state->has_backref)
into if (dfa->nbackref && next_state). Don't check for next_state
!= NULL twice.
* posix/regcomp.c (peek_token): Use opr.ctx_type instead of opr.idx
for ANCHOR.
(parse_expression): Only call init_word_char if word context will be
needed.
* posix/bug-regex11.c (tests): Add new tests.
* posix/tst-regex.c: Include getopt.h.
(timing): New variable.
(main): Set timing to 1 if --timing argument is present.
Add 2 new tests.
(run_test, run_test_backwards): Handle timing.
2003-11-27 Jakub Jelinek <jakub@redhat.com>
* posix/regex_internal.h (re_string_t): Remove mbs_case field.
Add offsets, valid_raw_len, raw_len, raw_stop, mbs_allocated and
offsets_needed fields. Change icase, is_utf8 and map_notascii
type from int bitfield to unsigned char.
(MBS_ALLOCATED, MBS_CASE_ALLOCATED): Remove.
(build_wcs_upper_buffer): Change prototype to return int.
(re_string_peek_byte_case, re_string_fetch_byte_case): Remove
defines, add prototypes.
* posix/regex_internal.c (re_string_allocate): Don't initialize
stop here. Don't initialize mbs_case. Set valid_raw_len.
Use mbs_allocated instead of MBS_* macros.
(re_string_construct): Don't initialize stop and valid_len here.
Don't initialize mbs_case. Use mbs_allocated instead of MBS_*
macros. Reallocate buffers if build_wcs_upper_buffer converted
too few bytes. Set valid_len to bufs_len only for single byte
no translation and set in that case valid_raw_len as well.
(re_string_realloc_buffers): Reallocate offsets if not NULL.
Use mbs_allocated instead of MBS_ALLOCATED. Don't reallocate
mbs_case.
(re_string_construct_common): Initialize raw_len, mbs_allocated,
stop and raw_stop.
(build_wcs_buffer): Apply pstr->trans before mbrtowc instead of
after it. Set valid_raw_len. Don't set mbs_case.
(build_wcs_upper_buffer): Return REG_NOERROR or REG_ESPACE.
Only use the fast path if !pstr->offsets_needed. Apply pstr->trans
before mbrtowc instead of after it. If upper case character
uses different number of bytes than lower case, goto to the
slow path. Don't call towupper unnecessarily twice. Set
valid_raw_len as well. Handle in the slow path the case if
lower and upper case use different number of characters.
Don't set mbs_case.
(re_string_skip_chars): Use valid_raw_len instead of valid_len.
(build_upper_buffer): Don't set mbs_case. Add BE macro. Set
valid_raw_len.
(re_string_translate_buffer): Set mbs instead of mbs_case. Set
valid_raw_len.
(re_string_reconstruct): Use raw_len/raw_stop to initialize
len/stop. Clear valid_raw_len and offsets_needed when clearing
valid_len. Use mbs_allocated instead of MBS_* macros.
Check original offset against valid_raw_len instead of valid_len.
Remove mbs_case handling. Adjust valid_raw_len together with
valid_len. If is_utf8 and looking for tip context, apply
pstr->trans first. If buffers start with partial multi-byte
character, initialize mbs array as well if mbs_allocated.
Check return value of build_wcs_upper_buffer.
(re_string_peek_byte_case): New function.
(re_string_fetch_byte_case): New function.
(re_string_destruct): Use mbs_allocated instead of MBS_ALLOCATED.
Don't free mbs_case. Free offsets.
* posix/regcomp.c (init_dfa): Only check if charset name is UTF-8
if mb_cur_max == 6.
* posix/regexec.c (re_search_internal): Initialize input.raw_stop
as well. Use valid_raw_len instead of valid_len when looking
through fastmap. Adjust registers through input.offsets.
(extend_buffers): Allow build_wcs_upper_buffer to fail.
* posix/bug-regex18.c (tests): Enable #ifdefed out tests. Add new
tests.
2003-11-29 06:13:09 +00:00
|
|
|
int n;
|
2013-02-27 05:58:20 +00:00
|
|
|
static const double c0 = 0.99674, c1 = -0.53380;
|
|
|
|
static const double c2 = 0.45472, c3 = -0.21553;
|
2001-03-12 07:57:09 +00:00
|
|
|
|
2001-03-12 00:04:52 +00:00
|
|
|
p.d = x;
|
2013-02-27 05:58:20 +00:00
|
|
|
p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF) | 0x3FE00000;
|
2001-03-12 00:04:52 +00:00
|
|
|
q.d = x;
|
|
|
|
y = p.d;
|
2013-02-27 05:58:20 +00:00
|
|
|
z = y - 1.0;
|
|
|
|
n = (q.i[HIGH_HALF] - p.i[HIGH_HALF]) >> 1;
|
|
|
|
z = ((c3 * z + c2) * z + c1) * z + c0; /* 2**-7 */
|
|
|
|
z = z * (1.5 - 0.5 * y * z * z); /* 2**-14 */
|
|
|
|
p.d = z * (1.5 - 0.5 * y * z * z); /* 2**-28 */
|
2001-03-12 00:04:52 +00:00
|
|
|
p.i[HIGH_HALF] -= n;
|
2013-02-27 05:58:20 +00:00
|
|
|
t = x * p.d;
|
|
|
|
return p.d * (1.5 - 0.5 * p.d * t);
|
2001-03-12 00:04:52 +00:00
|
|
|
}
|