bb3f4825c4
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.
104 lines
4.4 KiB
C
104 lines
4.4 KiB
C
|
|
/*
|
|
* IBM Accurate Mathematical Library
|
|
* written by International Business Machines Corp.
|
|
* Copyright (C) 2001 Free Software Foundation
|
|
*
|
|
* 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
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* 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
|
|
* GNU Lesser General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* along with this program; if not, write to the Free Software
|
|
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
|
|
*/
|
|
/****************************************************************************/
|
|
/* 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"
|
|
|
|
/****************************************************************************/
|
|
/* 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 */
|
|
/****************************************************************************/
|
|
|
|
double fastiroot(double);
|
|
|
|
void __mpsqrt(mp_no *x, mp_no *y, int p) {
|
|
#include "mpsqrt.h"
|
|
|
|
int i,m,ex,ey;
|
|
double dx,dy;
|
|
mp_no
|
|
mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
|
|
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
|
|
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
|
|
mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
|
|
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
|
|
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
|
|
mp_no mpxn,mpz,mpu,mpt1,mpt2;
|
|
|
|
/* Prepare multi-precision 1/2 and 3/2 */
|
|
mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD;
|
|
mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD;
|
|
|
|
ex=EX; 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);
|
|
|
|
m=mp[p];
|
|
for (i=0; i<m; i++) {
|
|
__mul(&mpu,&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;
|
|
|
|
return;
|
|
}
|
|
|
|
/***********************************************************/
|
|
/* Compute a double precision approximation for 1/sqrt(x) */
|
|
/* with the relative error bounded by 2**-51. */
|
|
/***********************************************************/
|
|
double fastiroot(double x) {
|
|
union {int i[2]; double d;} p,q;
|
|
double y,z, t;
|
|
int n;
|
|
static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
|
|
|
|
p.d = x;
|
|
p.i[HIGH_HALF] = (p.i[HIGH_HALF] & 0x3FFFFFFF ) | 0x3FE00000 ;
|
|
q.d = x;
|
|
y = p.d;
|
|
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 */
|
|
p.i[HIGH_HALF] -= n;
|
|
t = x*p.d;
|
|
return p.d*(1.5 - 0.5*p.d*t);
|
|
}
|