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

remove some obsolete fmprb methods

Browse source
This commit is contained in:
Fredrik Johansson 2016-02-25 02:28:19 +01:00
parent 85808d0a46
commit 3bc3c9e895
7 changed files with 0 additions and 760 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

31
doc/source/fmprb.rst
View file

@ -225,18 +225,6 @@ Radius and interval operations
Sets *x* to a ball containing the interval `[a, b]`. We
require that `a \le b`.
.. function:: slong fmprb_rel_error_bits(const fmprb_t x)
Returns the effective relative error of *x* measured in bits, defined as
the difference between the position of the top bit in the radius
and the top bit in the midpoint, plus one.
The result is clamped between plus/minus *FMPR_PREC_EXACT*.
.. function:: slong fmprb_rel_accuracy_bits(const fmprb_t x)
Returns the effective relative accuracy of *x* measured in bits,
equal to the negative of the return value from *fmprb_rel_error_bits*.
.. function:: slong fmprb_bits(const fmprb_t x)
Returns the number of bits needed to represent the absolute value
@ -244,14 +232,6 @@ Radius and interval operations
sufficient to represent *x* exactly. Returns 0 if the midpoint
of *x* is a special value.
.. function:: void fmprb_trim(fmprb_t y, const fmprb_t x)
Sets *y* to a trimmed copy of *x*: rounds *x* to a number of bits
equal to the accuracy of *x* (as indicated by its radius),
plus a few guard bits. The resulting ball is guaranteed to
contain *x*, but is more economical if *x* has
less than full accuracy.
.. function:: int fmprb_get_unique_fmpz(fmpz_t z, const fmprb_t x)
If *x* contains a unique integer, sets *z* to that value and returns
@ -491,14 +471,3 @@ Powers and roots
Sets *z* to the reciprocal square root of *x*, rounded to *prec* bits.
At high precision, this is faster than computing a square root.
.. function:: void fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec)
Sets *z* to the *k*-th root of *x*, rounded to *prec* bits.
As currently implemented, this function is only fast for small
fixed *k*. For large *k* it is better to use :func:`fmprb_pow_fmpq`
or :func:`fmprb_pow`.
.. function:: void fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec)
Sets *z* to the arithmetic-geometric mean of *x* and *y*.

285
fmprb.h
View file

@ -75,27 +75,6 @@ fmprb_clear(fmprb_t x)
fmpr_clear(fmprb_radref(x));
}
static __inline__ fmprb_ptr
_fmprb_vec_init(slong n)
{
slong i;
fmprb_ptr v = (fmprb_ptr) flint_malloc(sizeof(fmprb_struct) * n);
for (i = 0; i < n; i++)
fmprb_init(v + i);
return v;
}
static __inline__ void
_fmprb_vec_clear(fmprb_ptr v, slong n)
{
slong i;
for (i = 0; i < n; i++)
fmprb_clear(v + i);
flint_free(v);
}
static __inline__ int
fmprb_is_exact(const fmprb_t x)
{
@ -165,8 +144,6 @@ fmprb_set(fmprb_t x, const fmprb_t y)
void fmprb_set_round(fmprb_t z, const fmprb_t x, slong prec);
void fmprb_trim(fmprb_t y, const fmprb_t x);
static __inline__ void
fmprb_swap(fmprb_t x, fmprb_t y)
{
@ -324,10 +301,6 @@ void fmprb_submul_ui(fmprb_t z, const fmprb_t x, ulong y, slong prec);
void fmprb_submul_si(fmprb_t z, const fmprb_t x, slong y, slong prec);
void fmprb_submul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, slong prec);
void fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec);
void fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec);
static __inline__ void
fmprb_print(const fmprb_t x)
{
@ -491,38 +464,6 @@ void fmprb_set_interval_fmpr(fmprb_t x, const fmpr_t a, const fmpr_t b, slong pr
void fmprb_union(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec);
static __inline__ slong
fmprb_rel_error_bits(const fmprb_t x)
{
fmpz_t midmag, radmag;
slong result;
if (fmpr_is_zero(fmprb_radref(x)))
return -FMPR_PREC_EXACT;
if (fmpr_is_special(fmprb_midref(x)) || fmpr_is_special(fmprb_radref(x)))
return FMPR_PREC_EXACT;
fmpz_init(midmag);
fmpz_init(radmag);
fmpr_abs_bound_lt_2exp_fmpz(midmag, fmprb_midref(x));
fmpr_abs_bound_lt_2exp_fmpz(radmag, fmprb_radref(x));
fmpz_add_ui(radmag, radmag, 1);
result = _fmpz_sub_small(radmag, midmag);
fmpz_clear(midmag);
fmpz_clear(radmag);
return result;
}
static __inline__ slong
fmprb_rel_accuracy_bits(const fmprb_t x)
{
return -fmprb_rel_error_bits(x);
}
static __inline__ slong
fmprb_bits(const fmprb_t x)
{
@ -542,232 +483,6 @@ void fmprb_randtest_special(fmprb_t x, flint_rand_t state, slong prec, slong mag
void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, slong bits);
#define DEF_CACHED_CONSTANT(name, comp_func) \
TLS_PREFIX slong name ## _cached_prec = 0; \
TLS_PREFIX fmprb_t name ## _cached_value; \
void name ## _cleanup(void) \
{ \
fmprb_clear(name ## _cached_value); \
name ## _cached_prec = 0; \
} \
void name(fmprb_t x, slong prec) \
{ \
if (name ## _cached_prec < prec) \
{ \
if (name ## _cached_prec == 0) \
{ \
fmprb_init(name ## _cached_value); \
flint_register_cleanup_function(name ## _cleanup); \
} \
comp_func(name ## _cached_value, prec); \
name ## _cached_prec = prec; \
} \
fmprb_set_round(x, name ## _cached_value, prec); \
}
/* vector functions */
static __inline__ void
_fmprb_vec_zero(fmprb_ptr A, slong n)
{
slong i;
for (i = 0; i < n; i++)
fmprb_zero(A + i);
}
static __inline__ int
_fmprb_vec_is_zero(fmprb_srcptr vec, slong len)
{
slong i;
for (i = 0; i < len; i++)
if (!fmprb_is_zero(vec + i))
return 0;
return 1;
}
static __inline__ void
_fmprb_vec_set(fmprb_ptr res, fmprb_srcptr vec, slong len)
{
slong i;
for (i = 0; i < len; i++)
fmprb_set(res + i, vec + i);
}
static __inline__ void
_fmprb_vec_set_round(fmprb_ptr res, fmprb_srcptr vec, slong len, slong prec)
{
slong i;
for (i = 0; i < len; i++)
fmprb_set_round(res + i, vec + i, prec);
}
static __inline__ void
_fmprb_vec_swap(fmprb_ptr res, fmprb_ptr vec, slong len)
{
slong i;
for (i = 0; i < len; i++)
fmprb_swap(res + i, vec + i);
}
static __inline__ void
_fmprb_vec_neg(fmprb_ptr B, fmprb_srcptr A, slong n)
{
slong i;
for (i = 0; i < n; i++)
fmprb_neg(B + i, A + i);
}
static __inline__ void
_fmprb_vec_sub(fmprb_ptr C, fmprb_srcptr A,
fmprb_srcptr B, slong n, slong prec)
{
slong i;
for (i = 0; i < n; i++)
fmprb_sub(C + i, A + i, B + i, prec);
}
static __inline__ void
_fmprb_vec_add(fmprb_ptr C, fmprb_srcptr A,
fmprb_srcptr B, slong n, slong prec)
{
slong i;
for (i = 0; i < n; i++)
fmprb_add(C + i, A + i, B + i, prec);
}
static __inline__ void
_fmprb_vec_scalar_mul(fmprb_ptr res, fmprb_srcptr vec,
slong len, const fmprb_t c, slong prec)
{
slong i;
for (i = 0; i < len; i++)
fmprb_mul(res + i, vec + i, c, prec);
}
static __inline__ void
_fmprb_vec_scalar_div(fmprb_ptr res, fmprb_srcptr vec,
slong len, const fmprb_t c, slong prec)
{
slong i;
for (i = 0; i < len; i++)
fmprb_div(res + i, vec + i, c, prec);
}
static __inline__ void
_fmprb_vec_scalar_mul_fmpz(fmprb_ptr res, fmprb_srcptr vec,
slong len, const fmpz_t c, slong prec)
{
slong i;
fmpr_t t;
fmpr_init(t);
fmpr_set_fmpz(t, c);
for (i = 0; i < len; i++)
fmprb_mul_fmpr(res + i, vec + i, t, prec);
fmpr_clear(t);
}
static __inline__ void
_fmprb_vec_scalar_mul_2exp_si(fmprb_ptr res, fmprb_srcptr src, slong len, slong c)
{
slong i;
for (i = 0; i < len; i++)
fmprb_mul_2exp_si(res + i, src + i, c);
}
static __inline__ void
_fmprb_vec_scalar_addmul(fmprb_ptr res, fmprb_srcptr vec,
slong len, const fmprb_t c, slong prec)
{
slong i;
for (i = 0; i < len; i++)
fmprb_addmul(res + i, vec + i, c, prec);
}
static __inline__ void
_fmprb_vec_get_abs_ubound_fmpr(fmpr_t bound, fmprb_srcptr vec,
slong len, slong prec)
{
fmpr_t t;
slong i;
if (len < 1)
{
fmpr_zero(bound);
}
else
{
fmprb_get_abs_ubound_fmpr(bound, vec, prec);
fmpr_init(t);
for (i = 1; i < len; i++)
{
fmprb_get_abs_ubound_fmpr(t, vec + i, prec);
if (fmpr_cmp(t, bound) > 0)
fmpr_set(bound, t);
}
fmpr_clear(t);
}
}
static __inline__ slong
_fmprb_vec_bits(fmprb_srcptr x, slong len)
{
slong i, b, c;
b = 0;
for (i = 0; i < len; i++)
{
c = fmprb_bits(x + i);
b = FLINT_MAX(b, c);
}
return b;
}
static __inline__ void
_fmprb_vec_set_powers(fmprb_ptr xs, const fmprb_t x, slong len, slong prec)
{
slong i;
for (i = 0; i < len; i++)
{
if (i == 0)
fmprb_one(xs + i);
else if (i == 1)
fmprb_set_round(xs + i, x, prec);
else if (i % 2 == 0)
fmprb_mul(xs + i, xs + i / 2, xs + i / 2, prec);
else
fmprb_mul(xs + i, xs + i - 1, x, prec);
}
}
static __inline__ void
_fmprb_vec_add_error_fmpr_vec(fmprb_ptr res, fmpr_srcptr err, slong len)
{
slong i;
for (i = 0; i < len; i++)
fmprb_add_error_fmpr(res + i, err + i);
}
static __inline__ void
_fmprb_vec_indeterminate(fmprb_ptr vec, slong len)
{
slong i;
for (i = 0; i < len; i++)
{
fmpr_nan(fmprb_midref(vec + i));
fmpr_pos_inf(fmprb_radref(vec + i));
}
}
static __inline__ void
_fmprb_vec_trim(fmprb_ptr res, fmprb_srcptr vec, slong len)
{
slong i;
for (i = 0; i < len; i++)
fmprb_trim(res + i, vec + i);
}
#ifdef __cplusplus
}
#endif

81
fmprb/agm.c
View file

@ -1,81 +0,0 @@
/*=============================================================================
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 "fmprb.h"
void
fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec)
{
fmprb_t t, u, v, w;
if (fmprb_contains_negative(x) || fmprb_contains_negative(y))
{
fmprb_indeterminate(z);
return;
}
if (fmprb_is_zero(x) || fmprb_is_zero(y))
{
fmprb_zero(z);
return;
}
fmprb_init(t);
fmprb_init(u);
fmprb_init(v);
fmprb_init(w);
fmprb_set(t, x);
fmprb_set(u, y);
while (!fmprb_overlaps(t, u) &&
!fmprb_contains_nonpositive(t) &&
!fmprb_contains_nonpositive(u))
{
fmprb_add(v, t, u, prec);
fmprb_mul_2exp_si(v, v, -1);
fmprb_mul(w, t, u, prec);
fmprb_sqrt(w, w, prec);
fmprb_swap(v, t);
fmprb_swap(w, u);
}
if (!fmprb_is_finite(t) || !fmprb_is_finite(u))
{
fmprb_indeterminate(z);
}
else
{
fmprb_union(z, t, u, prec);
}
fmprb_clear(t);
fmprb_clear(u);
fmprb_clear(v);
fmprb_clear(w);
}

102
fmprb/root.c
View file

@ -1,102 +0,0 @@
/*=============================================================================
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) 2012, 2013 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
void
fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec)
{
slong r;
if (k == 1)
{
fmprb_set_round(z, x, prec);
return;
}
if (k == 2)
{
fmprb_sqrt(z, x, prec);
return;
}
if (k == -2)
{
fmprb_rsqrt(z, x, prec);
return;
}
if (fmprb_contains_negative(x))
{
fmprb_indeterminate(z);
return;
}
if (fmprb_is_exact(x))
{
r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
}
else
{
fmpr_t t, u, err;
fmpr_init(t);
fmpr_init(u);
fmpr_init(err);
/* lower point */
fmpr_sub(err, fmprb_midref(x), fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_DOWN);
if (fmpr_is_zero(err))
{
fmpr_root(err, fmprb_midref(x), k, FMPRB_RAD_PREC, FMPR_RND_UP);
r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN);
fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r,
FMPRB_RAD_PREC, FMPR_RND_UP);
}
else
{
/* derivative x^(1/k) / (x k) at lower point */
fmpr_root(t, err, k, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul_ui(u, err, k, FMPRB_RAD_PREC, FMPR_RND_DOWN);
fmpr_divappr_abs_ubound(t, t, u, FMPRB_RAD_PREC);
/* multiplied by distance */
fmpr_mul(err, t, fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN);
fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r,
FMPRB_RAD_PREC, FMPR_RND_UP);
}
fmpr_clear(t);
fmpr_clear(u);
fmpr_clear(err);
}
fmprb_adjust(z);
}

102
fmprb/test/t-agm.c
View file

@ -1,102 +0,0 @@
/*=============================================================================
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) 2012 Fredrik Johansson
******************************************************************************/
#include "fmprb.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("agm....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmprb_t a, b, c;
fmpq_t q, r;
mpfr_t t, u;
slong prec = 2 + n_randint(state, 200);
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmpq_init(q);
fmpq_init(r);
mpfr_init2(t, prec + 100);
mpfr_init2(u, prec + 100);
fmprb_randtest(a, state, 1 + n_randint(state, 200), 3);
fmprb_randtest(b, state, 1 + n_randint(state, 200), 3);
fmprb_randtest(c, state, 1 + n_randint(state, 200), 3);
fmprb_agm(c, a, b, prec);
if (fmprb_equal(a, b))
{
if (!fmprb_contains(c, a))
{
flint_printf("FAIL: containment (identity)\n\n");
flint_printf("a = "); fmprb_print(a); flint_printf("\n\n");
flint_printf("b = "); fmprb_print(b); flint_printf("\n\n");
flint_printf("c = "); fmprb_print(c); flint_printf("\n\n");
abort();
}
}
else
{
fmprb_get_rand_fmpq(q, state, a, 1 + n_randint(state, 200));
fmprb_get_rand_fmpq(r, state, b, 1 + n_randint(state, 200));
fmpq_get_mpfr(t, q, MPFR_RNDN);
fmpq_get_mpfr(u, r, MPFR_RNDN);
mpfr_agm(t, t, u, MPFR_RNDN);
if (!fmprb_contains_mpfr(c, t))
{
flint_printf("FAIL: containment\n\n");
flint_printf("a = "); fmprb_print(a); flint_printf("\n\n");
flint_printf("b = "); fmprb_print(b); flint_printf("\n\n");
flint_printf("c = "); fmprb_print(c); flint_printf("\n\n");
abort();
}
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
fmpq_clear(q);
fmpq_clear(r);
mpfr_clear(t);
mpfr_clear(u);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

91
fmprb/test/t-trim.c
View file

@ -1,91 +0,0 @@
/*=============================================================================
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 "fmprb.h"
int main()
{
slong iter;
flint_rand_t state;
flint_printf("trim....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmprb_t x, y;
slong acc1, acc2;
int accuracy_ok;
fmprb_init(x);
fmprb_init(y);
fmprb_randtest_special(x, state, 1000, 100);
fmprb_trim(y, x);
if (!fmprb_contains(y, x))
{
flint_printf("FAIL (containment):\n\n");
flint_printf("x = "); fmprb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmprb_print(y); flint_printf("\n\n");
abort();
}
acc1 = fmprb_rel_accuracy_bits(x);
acc2 = fmprb_rel_accuracy_bits(y);
accuracy_ok = (acc1 < 0 && acc2 <= acc1) || (acc2 >= acc1 - 1);
if (!accuracy_ok)
{
flint_printf("FAIL (accuracy):\n\n");
flint_printf("x: %wd, y = %wd\n\n", acc1, acc2);
flint_printf("x = "); fmprb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmprb_print(y); flint_printf("\n\n");
abort();
}
fmprb_trim(x, x);
if (!fmprb_equal(y, x))
{
flint_printf("FAIL (aliasing):\n\n");
flint_printf("x = "); fmprb_print(x); flint_printf("\n\n");
flint_printf("y = "); fmprb_print(y); flint_printf("\n\n");
abort();
}
fmprb_clear(x);
fmprb_clear(y);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

68
fmprb/trim.c
View file

@ -1,68 +0,0 @@
/*=============================================================================
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 "fmprb.h"
#define TRIM_PADDING 16
void
fmprb_trim(fmprb_t y, const fmprb_t x)
{
if (fmpr_is_zero(fmprb_radref(x)) || fmpr_is_special(fmprb_midref(x)))
{
fmprb_set(y, x);
}
else if (fmpr_is_special(fmprb_radref(x)))
{
/* midpoint must be finite, so set to 0 +/- inf */
fmprb_zero_pm_inf(y);
}
else
{
slong bits, accuracy;
bits = fmprb_bits(x);
accuracy = fmprb_rel_accuracy_bits(x);
if (accuracy < -TRIM_PADDING)
{
/* set to 0 +/- rad */
if (fmpr_sgn(fmprb_midref(x)) < 0)
fmpr_sub(fmprb_radref(y), fmprb_radref(x), fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
else
fmpr_add(fmprb_radref(y), fmprb_radref(x), fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_zero(fmprb_midref(y));
}
else if (accuracy < bits - TRIM_PADDING)
{
fmprb_set_round(y, x, FLINT_MAX(0, accuracy) + TRIM_PADDING);
}
else
{
fmprb_set(y, x);
}
}
}