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

port the div_2expm1_ui code completely to the arb type

Browse source
This commit is contained in:
Fredrik Johansson 2015-07-21 18:05:30 +02:00
parent b24d5d482d
commit 6c81cba454
6 changed files with 50 additions and 172 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

61
fmprb/div_2expm1.c → arb/div_2expm1_ui.c
View file

@ -19,61 +19,64 @@
=============================================================================*/ =============================================================================*/
/****************************************************************************** /******************************************************************************
Copyright (C) 2012 Fredrik Johansson Copyright (C) 2015 Fredrik Johansson
******************************************************************************/ ******************************************************************************/
#include "fmprb.h" #include "arb.h"
/* TODO: fix/optimize for n > prec */
void void
fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec) arb_div_2expm1_ui(arb_t y, const arb_t x, ulong n, long prec)
{ {
if (n < FLINT_BITS) if (n < FLINT_BITS)
{ {
fmprb_div_ui(y, x, (1UL << n) - 1, prec); arb_div_ui(y, x, (1UL << n) - 1, prec);
} }
else if (n < 1024 + prec / 32 || n > LONG_MAX / 2) else if (n < 1024 + prec / 32 || n > LONG_MAX / 4)
{ {
fmprb_t t; arb_t t;
fmprb_init(t); fmpz_t e;
fmprb_one(t);
fmpz_set_ui(fmpr_expref(fmprb_midref(t)), n); arb_init(t);
fmprb_sub_ui(t, t, 1, prec); fmpz_init_set_ui(e, n);
fmprb_div(y, x, t, prec);
fmprb_clear(t); arb_one(t);
arb_mul_2exp_fmpz(t, t, e);
arb_sub_ui(t, t, 1, prec);
arb_div(y, x, t, prec);
arb_clear(t);
fmpz_clear(e);
} }
else else
{ {
fmprb_t s, t; arb_t s, t;
long i, b; long i, b;
fmprb_init(s); arb_init(s);
fmprb_init(t); arb_init(t);
/* x / (2^n - 1) = sum_{k>=1} x * 2^(-k*n)*/ /* x / (2^n - 1) = sum_{k>=1} x * 2^(-k*n)*/
arb_mul_2exp_si(s, x, -n);
fmprb_mul_2exp_si(s, x, -n); arb_set(t, s);
fmprb_set(t, s);
b = 1; b = 1;
for (i = 2; i <= prec / n + 1; i++) for (i = 2; i <= prec / n + 1; i++)
{ {
fmprb_mul_2exp_si(t, t, -n); arb_mul_2exp_si(t, t, -n);
fmprb_add(s, s, t, prec); arb_add(s, s, t, prec);
b = i; b = i;
} }
/* error bound: sum_{k>b} x * 2^(-k*n) <= x * 2^(-b*n - (n-1)) */ /* error bound: sum_{k>b} x * 2^(-k*n) <= x * 2^(-b*n - (n-1)) */
fmprb_mul_2exp_si(t, x, -b*n - (n-1)); arb_mul_2exp_si(t, x, -b*n - (n-1));
fmprb_abs(t, t); arb_abs(t, t);
fmprb_add_error(s, t); arb_add_error(s, t);
fmprb_set(y, s); arb_set(y, s);
fmprb_clear(s); arb_clear(s);
fmprb_clear(t); arb_clear(t);
} }
} }

37
arb/fmprb_wrappers.c
View file

@ -1,37 +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) 2014 Fredrik Johansson
******************************************************************************/
#include "arb.h"
void arb_div_2expm1_ui(arb_t z, const arb_t x, ulong n, long prec)
{
fmprb_t t;
fmprb_init(t);
arb_get_fmprb(t, x);
fmprb_div_2expm1_ui(t, t, n, prec);
arb_set_fmprb(z, t);
fmprb_clear(t);
}

21
arb/test/t-div_2expm1_ui.c
View file

@ -49,12 +49,26 @@ int main()
prec = 2 + n_randint(state, 10000); prec = 2 + n_randint(state, 10000);
arb_randtest(a, state, 1 + n_randint(state, 10000), 10); arb_randtest(a, state, 1 + n_randint(state, 10000), 10);
n = 1 + (n_randtest(state) % (10 * prec));
if (n_randint(state, 2))
n = 1 + (n_randtest(state) % (10 * prec));
else
n = n_randtest(state);
arb_div_2expm1_ui(b, a, n, prec); arb_div_2expm1_ui(b, a, n, prec);
arb_one(c); arb_one(c);
arb_mul_2exp_si(c, c, n); if (n >= (1UL << (FLINT_BITS-1)))
{
arb_mul_2exp_si(c, c, (1UL << (FLINT_BITS-2)));
arb_mul_2exp_si(c, c, (1UL << (FLINT_BITS-2)));
arb_mul_2exp_si(c, c, n - (1UL << (FLINT_BITS-1)));
}
else
{
arb_mul_2exp_si(c, c, n);
}
arb_sub_ui(c, c, 1, prec); arb_sub_ui(c, c, 1, prec);
arb_div(c, a, c, prec); arb_div(c, a, c, prec);
@ -64,13 +78,14 @@ int main()
if (!arb_overlaps(b, c)) if (!arb_overlaps(b, c))
{ {
printf("FAIL: containment\n\n"); printf("FAIL: containment\n\n");
printf("n = %lu\n", n);
printf("a = "); arb_print(a); printf("\n\n"); printf("a = "); arb_print(a); printf("\n\n");
printf("b = "); arb_print(b); printf("\n\n"); printf("b = "); arb_print(b); printf("\n\n");
printf("c = "); arb_print(c); printf("\n\n"); printf("c = "); arb_print(c); printf("\n\n");
abort(); abort();
} }
if ((acc2 < FLINT_MIN(prec, acc1) - 10) && if (n > 0 && (acc2 < FLINT_MIN(prec, acc1) - 10) &&
!(acc1 == -ARF_PREC_EXACT && acc2 == -ARF_PREC_EXACT)) !(acc1 == -ARF_PREC_EXACT && acc2 == -ARF_PREC_EXACT))
{ {
printf("FAIL: poor accuracy\n\n"); printf("FAIL: poor accuracy\n\n");

4
doc/source/fmprb.rst
View file

@ -437,10 +437,6 @@ Arithmetic
where the triangle inequality has been applied to the numerator and where the triangle inequality has been applied to the numerator and
the reverse triangle inequality has been applied to the denominator. the reverse triangle inequality has been applied to the denominator.
.. function:: void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec)
Sets `y = x / (2^n - 1)`, rounded to *prec* bits.
.. function:: void fmprb_addmul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec) .. function:: void fmprb_addmul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)
.. function:: void fmprb_addmul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec) .. function:: void fmprb_addmul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec)

3
fmprb.h
View file

@ -292,9 +292,6 @@ fmprb_inv(fmprb_t y, const fmprb_t x, long prec)
fmprb_ui_div(y, 1, x, prec); fmprb_ui_div(y, 1, x, prec);
} }
void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec);
void fmprb_mul_fmpr_naive(fmprb_t z, const fmprb_t x, const fmpr_t y, long prec); void fmprb_mul_fmpr_naive(fmprb_t z, const fmprb_t x, const fmpr_t y, long prec);
void fmprb_mul_main_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec); void fmprb_mul_main_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec);
void fmprb_mul_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec); void fmprb_mul_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec);

96
fmprb/test/t-div_2expm1_ui.c
View file

@ -1,96 +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()
{
long iter;
flint_rand_t state;
printf("div_2expm1_ui....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmprb_t a, b, c;
ulong n;
long prec, acc1, acc2;
fmpz_t t;
fmprb_init(a);
fmprb_init(b);
fmprb_init(c);
fmpz_init(t);
prec = 2 + n_randint(state, 10000);
fmprb_randtest(a, state, 1 + n_randint(state, 10000), 10);
n = 1 + (n_randtest(state) % (10 * prec));
fmprb_div_2expm1_ui(b, a, n, prec);
fmprb_one(c);
fmpz_set_ui(fmpr_expref(fmprb_midref(c)), n);
fmprb_sub_ui(c, c, 1, prec);
fmprb_div(c, a, c, prec);
acc1 = fmprb_rel_accuracy_bits(a);
acc2 = fmprb_rel_accuracy_bits(b);
if (!fmprb_overlaps(b, c))
{
printf("FAIL: containment\n\n");
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
if ((acc2 < FLINT_MIN(prec, acc1) - 10) &&
!(acc1 == -FMPR_PREC_EXACT && acc2 == -FMPR_PREC_EXACT))
{
printf("FAIL: poor accuracy\n\n");
printf("prec=%ld, acc1=%ld, acc2=%ld\n\n", prec, acc1, acc2);
printf("n = %lu\n\n", n);
printf("a = "); fmprb_print(a); printf("\n\n");
printf("b = "); fmprb_print(b); printf("\n\n");
printf("c = "); fmprb_print(c); printf("\n\n");
abort();
}
fmprb_clear(a);
fmprb_clear(b);
fmprb_clear(c);
fmpz_clear(t);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}