diff options
-rw-r--r-- | extract.c | 92 |
1 files changed, 47 insertions, 45 deletions
@@ -1,61 +1,63 @@ -#include <stdio.h> -#include <math.h> -#include <stdlib.h> -#include <strings.h> +/* mpfr_extract -- bit-extraction function for the binary splitting algorithm + +Copyright (C) 2000 Free Software Foundation. + +This file is part of the MPFR Library. + +The MPFR Library 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. + +The MPFR 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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + #include "gmp.h" #include "gmp-impl.h" #include "mpfr.h" -#include "assert.h" +/* given 0 <= |p| < 1, this function extracts limbs of p and puts them in y. + It is mainly designed for the "binary splitting" algorithm together with + generic.c. -int + More precisely, if B = 2^BITS_PER_MP_LIMB: + - for i=0, y = floor(p * B) + - for i>0, y = (p * B^(2^i)) mod B^(2^(i-1)) + */ +void #if __STDC__ -mpfr_extract(mpz_ptr y, mpfr_srcptr p, int i) +mpfr_extract (mpz_ptr y, mpfr_srcptr p, unsigned int i) #else -mpfr_extract(y,p,i) -mpz_ptr y; -mpfr_srcptr p; -int i; +mpfr_extract (y, p, i) + mpz_ptr y; + mpfr_srcptr p; + unsigned int i; #endif { - int two_i=1 << i; + int two_i = 1 << i; int two_i_2 = i ? two_i / 2 : 1; - int size; - int j; + mp_size_t size_p = MPFR_ABSSIZE(p); - /* TODO: gestion de l'infini. */ + /* as 0 <= |p| < 1, we don't have to care with infinities, NaN, ... */ - if (MPFR_ABSSIZE(p) < two_i){ - int j; - y->_mp_alloc=two_i_2 ; - y->_mp_size=two_i_2 ; - PTR(y) = (*_mp_allocate_func)(two_i_2 * sizeof(mp_limb_t)); + _mpz_realloc (y, two_i_2); + if (size_p < two_i) { MPN_ZERO (PTR(y), two_i_2); - assert(y->_mp_d!=NULL); - for(j= 0; j < MPFR_ABSSIZE(p) - two_i_2 ; j++){ - y->_mp_d[j + two_i - MPFR_ABSSIZE(p)] = p->_mp_d[j]; - } - } else - { - PTR(y) = (*_mp_allocate_func)(two_i_2 * sizeof(mp_limb_t)); - memcpy(y -> _mp_d, p->_mp_d+MPFR_ABSSIZE(p) - two_i, two_i_2 * sizeof(mp_limb_t)); - y->_mp_alloc=two_i_2 ; - y->_mp_size=two_i_2 ; - } - - size = MPFR_ABSSIZE(y); - for (j = 0; j < size; j++) - { - if (y->_mp_d[j]) - { - if MPFR_ISNEG(p) - mpz_neg(y,y); - return 0; - } - } - y->_mp_size=0; - - return 0; + if (size_p >= two_i_2) + MPN_COPY (PTR(y) + two_i - size_p, MPFR_MANT(p), size_p - two_i_2); + } + else + MPN_COPY (PTR(y), MPFR_MANT(p) + size_p - two_i, two_i_2); + + MPN_NORMALIZE (PTR(y), two_i_2); + SIZ(y) = (MPFR_ISNEG(p)) ? -two_i_2 : two_i_2; } |