hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

add fmpr_abs_bound_lt_2exp_si; some more simplification

Browse source
This commit is contained in:
Fredrik Johansson 2013-09-19 11:50:01 +01:00
parent 28159cd9a9
commit 991814c481
8 changed files with 173 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

7
doc/source/fmpr.rst
View file

@ -358,6 +358,13 @@ Comparisons
Sets *b* to the smallest integer such that `|x| < 2^b`.
If *x* is zero, infinity or NaN, the result is undefined.
.. function:: long fmpr_abs_bound_lt_2exp_si(const fmpr_t x)
Returns the smallest integer *b* such that `|x| < 2^b`, clamping
the result to lie between -*FMPR_PREC_EXACT* and *FMPR_PREC_EXACT*
inclusive. If *x* is zero, -*FMPR_PREC_EXACT* is returned,
and if *x* is infinity or NaN, *FMPR_PREC_EXACT* is returned.
Random number generation
-------------------------------------------------------------------------------

11
elefun/exp.c
View file

@ -36,16 +36,19 @@ elefun_exp_precomp(fmprb_t z, const fmprb_t x, long prec, int minus_one)
if (fmpr_cmp_2exp_si(fmprb_radref(x), -2) >= 0)
return 0;
if (fmpr_cmpabs_2exp_si(fmprb_midref(x), 128) >= 0)
/* magnitude clamped between +/- FMPR_PREC_EXACT */
mag = fmpr_abs_bound_lt_2exp_si(fmprb_midref(x));
/* too large */
if (mag >= 128)
return 0;
if (fmpr_cmpabs_2exp_si(fmprb_midref(x), -prec) < 0)
/* too small */
if (mag < -prec)
return 0;
fixed_wp = prec + 2 * FLINT_BIT_COUNT(prec);
mag = *fmpr_expref(fmprb_midref(x)) + fmpz_bits(fmpr_manref(fmprb_midref(x)));
if (mag > 0)
fixed_wp += mag; /* for argument reduction */
else if (minus_one)

31
fmpcb_poly/find_roots.c
View file

@ -25,36 +25,13 @@
#include "fmpcb_poly.h"
long _fmpr_get_exp(const fmpr_t x)
{
if (fmpr_is_zero(x))
{
return -FMPR_PREC_EXACT;
}
else if (fmpr_is_inf(x) || fmpr_is_nan(x))
{
return FMPR_PREC_EXACT;
}
else
{
long res;
fmpz_t t;
fmpz_init(t);
fmpz_set(t, fmpr_expref(x));
fmpz_add_ui(t, t, fmpz_bits(fmpr_manref(x)));
res = fmpz_get_si(t);
fmpz_clear(t);
return res;
}
}
long
_fmpcb_get_mid_mag(const fmpcb_t z)
{
long rm, im;
rm = _fmpr_get_exp(fmprb_midref(fmpcb_realref(z)));
im = _fmpr_get_exp(fmprb_midref(fmpcb_imagref(z)));
rm = fmpr_abs_bound_lt_2exp_si(fmprb_midref(fmpcb_realref(z)));
im = fmpr_abs_bound_lt_2exp_si(fmprb_midref(fmpcb_imagref(z)));
return FLINT_MAX(rm, im);
}
@ -64,8 +41,8 @@ _fmpcb_get_rad_mag(const fmpcb_t z)
{
long rm, im;
rm = _fmpr_get_exp(fmprb_radref(fmpcb_realref(z)));
im = _fmpr_get_exp(fmprb_radref(fmpcb_imagref(z)));
rm = fmpr_abs_bound_lt_2exp_si(fmprb_radref(fmpcb_realref(z)));
im = fmpr_abs_bound_lt_2exp_si(fmprb_radref(fmpcb_imagref(z)));
return FLINT_MAX(rm, im);
}

2
fmpr.h
View file

@ -612,6 +612,8 @@ fmpr_abs_bound_lt_2exp_fmpz(fmpz_t b, const fmpr_t x)
fmpz_add_ui(b, fmpr_expref(x), fmpz_bits(fmpr_manref(x)));
}
long fmpr_abs_bound_lt_2exp_si(const fmpr_t x);
static __inline__ void
fmpr_min(fmpr_t z, const fmpr_t a, const fmpr_t b)
{

60
fmpr/abs_bound_lt_2exp_si.c Normal file
View file

@ -0,0 +1,60 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "fmpr.h"
long
fmpr_abs_bound_lt_2exp_si(const fmpr_t x)
{
if (fmpr_is_special(x))
{
if (fmpr_is_zero(x))
return -FMPR_PREC_EXACT;
else
return FMPR_PREC_EXACT;
}
else
{
long res;
fmpz_t t;
fmpz_init(t);
fmpr_abs_bound_lt_2exp_fmpz(t, x);
if (fmpz_fits_si(t))
res = fmpz_get_si(t);
else
res = fmpz_sgn(t) < 0 ? -FMPR_PREC_EXACT : FMPR_PREC_EXACT;
fmpz_clear(t);
if (res < -FMPR_PREC_EXACT)
res = -FMPR_PREC_EXACT;
if (res > FMPR_PREC_EXACT)
res = FMPR_PREC_EXACT;
return res;
}
}

89
fmpr/test/t-abs_bound_lt_2exp_si.c Normal file
View file

@ -0,0 +1,89 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/
#include "fmpr.h"
int main()
{
long iter;
flint_rand_t state;
printf("abs_bound_lt_2exp_si....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
fmpr_t x;
fmpz_t b;
long c;
fmpr_init(x);
fmpz_init(b);
fmpr_randtest_special(x, state, 2 + n_randint(state, 1000), 100);
fmpr_abs_bound_lt_2exp_fmpz(b, x);
c = fmpr_abs_bound_lt_2exp_si(x);
if (c == -FMPR_PREC_EXACT)
{
if (!(fmpr_is_zero(x) || fmpz_cmp_si(b, -FMPR_PREC_EXACT) <= 0))
{
printf("FAIL (small/zero)\n\n");
abort();
}
}
else if (c == FMPR_PREC_EXACT)
{
if (!(fmpr_is_inf(x) || fmpr_is_nan(x) ||
fmpz_cmp_si(b, FMPR_PREC_EXACT) >= 0))
{
printf("FAIL (large/inf/nan)\n\n");
abort();
}
}
else
{
if (fmpz_cmp_si(b, c) != 0)
{
printf("FAIL (normal)\n\n");
printf("x = "); fmpr_print(x); printf("\n\n");
printf("b = "); fmpz_print(b); printf("\n\n");
printf("c = %ld\n\n", c);
abort();
}
}
fmpr_clear(x);
fmpz_clear(b);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

5
fmprb_poly/newton_refine_root.c
View file

@ -25,10 +25,9 @@
#include "fmprb_poly.h"
long _fmpr_mag(const fmpr_t c)
static __inline__ long _fmpr_mag(const fmpr_t c)
{
long m = fmpz_bits(fmpr_manref(c)) + fmpz_get_si(fmpr_expref(c));
long m = fmpr_abs_bound_lt_2exp_si(c);
return FLINT_MAX(m, 0);
}

20
gamma/stirling_choose_param.c
View file

@ -29,22 +29,6 @@
/* tuning factor */
#define GAMMA_STIRLING_BETA 0.27
static __inline__ long
_fmpr_mag(const fmpr_t x)
{
if (fmpr_is_special(x))
{
if (fmpr_is_zero(x))
return -FMPR_PREC_EXACT;
else
return FMPR_PREC_EXACT;
}
else
{
return fmpz_bits(fmpr_manref(x)) + *fmpr_expref(x) - 1;
}
}
#define PI 3.1415926535897932385
static long
@ -133,8 +117,8 @@ choose_large(int * reflect, long * r, long * n,
long ab, bb;
double log2z, argz;
ab = _fmpr_mag(a);
bb = _fmpr_mag(b);
ab = fmpr_abs_bound_lt_2exp_si(a);
bb = fmpr_abs_bound_lt_2exp_si(b);
log2z = FLINT_MAX(ab, bb);