mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
square roots and some tweaks
This commit is contained in:
Fredrik Johansson
parent
0dff48c4e7
commit
9b992a306c
7 changed files with 218 additions and 20 deletions
5
arb.h
5
arb.h
|
|
@ -556,6 +556,11 @@ arb_set_fmpq(arb_t y, const fmpq_t x, long prec)
|
|||
arb_fmpz_div_fmpz(y, fmpq_numref(x), fmpq_denref(x), prec);
|
||||
}
|
||||
|
||||
void arb_sqrt(arb_t z, const arb_t x, long prec);
|
||||
void arb_sqrt_arf(arb_t z, const arf_t x, long prec);
|
||||
void arb_sqrt_fmpz(arb_t z, const fmpz_t x, long prec);
|
||||
void arb_sqrt_ui(arb_t z, ulong x, long prec);
|
||||
|
||||
/* vector functions */
|
||||
|
||||
static __inline__ void
|
||||
|
|
|
|||
12
arb/addmul.c
12
arb/addmul.c
|
|
@ -36,20 +36,16 @@ arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, long prec)
|
|||
inexact = arf_addmul(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
{
|
||||
if (ARB_IS_LAGOM(z)) /* todo: arf_mag_fast_add_ulp -- or make arf_mag_add_ulp fast? */
|
||||
mag_fast_add_2exp_si(arb_radref(z), arb_radref(z), ARF_EXP(arb_midref(z)) - prec);
|
||||
else
|
||||
arf_mag_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
arf_mag_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
else if (ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
|
||||
{
|
||||
mag_fast_init_set_arf(ym, y);
|
||||
mag_fast_addmul(arb_radref(z), ym, arb_radref(x));
|
||||
inexact = arf_addmul(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(arb_radref(z), arb_radref(z), ARF_EXP(arb_midref(z)) - prec);
|
||||
arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -92,7 +88,7 @@ arb_addmul(arb_t z, const arb_t x, const arb_t y, long prec)
|
|||
prec, ARF_RND_DOWN);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(zr, zr, ARF_EXP(ARB_MIDREF(z)) - prec);
|
||||
arf_mag_fast_add_ulp(zr, zr, arb_midref(z), prec);
|
||||
|
||||
*ARB_RADREF(z) = *zr;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -50,7 +50,7 @@ arb_mul_arf(arb_t z, const arb_t x, const arf_t y, long prec)
|
|||
inexact = arf_mul(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(zr, zr, ARF_EXP(arb_midref(z)) - prec);
|
||||
arf_mag_fast_add_ulp(zr, zr, arb_midref(z), prec);
|
||||
|
||||
*arb_radref(z) = *zr;
|
||||
}
|
||||
|
|
@ -101,7 +101,7 @@ arb_mul(arb_t z, const arb_t x, const arb_t y, long prec)
|
|||
inexact = arf_mul(arb_midref(z), arb_midref(x), arb_midref(y), prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(zr, zr, ARF_EXP(arb_midref(z)) - prec);
|
||||
arf_mag_fast_add_ulp(zr, zr, arb_midref(z), prec);
|
||||
|
||||
*arb_radref(z) = *zr;
|
||||
}
|
||||
|
|
|
|||
113
arb/sqrt.c
Normal file
113
arb/sqrt.c
Normal file
|
|
@ -0,0 +1,113 @@
|
|||
/*=============================================================================
|
||||
|
||||
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_sqrt_ui(arb_t z, ulong x, long prec)
|
||||
{
|
||||
arf_t t;
|
||||
arf_init_set_ui(t, x); /* no need to free */
|
||||
arb_sqrt_arf(z, t, prec);
|
||||
}
|
||||
|
||||
void
|
||||
arb_sqrt_fmpz(arb_t z, const fmpz_t x, long prec)
|
||||
{
|
||||
arf_t t;
|
||||
arf_init(t);
|
||||
arf_set_fmpz(t, x);
|
||||
arb_sqrt_arf(z, t, prec);
|
||||
arf_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
arb_sqrt_arf(arb_t z, const arf_t x, long prec)
|
||||
{
|
||||
if (arf_sgn(x) < 0)
|
||||
{
|
||||
arb_indeterminate(z);
|
||||
}
|
||||
else
|
||||
{
|
||||
int inexact;
|
||||
|
||||
inexact = arf_sqrt(arb_midref(z), x, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
arf_mag_set_ulp(arb_radref(z), arb_midref(z), prec);
|
||||
else
|
||||
mag_zero(arb_radref(z));
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
arb_sqrt(arb_t z, const arb_t x, long prec)
|
||||
{
|
||||
mag_t rx, zr;
|
||||
int inexact;
|
||||
|
||||
if (arb_is_exact(x))
|
||||
{
|
||||
arb_sqrt_arf(z, arb_midref(x), prec);
|
||||
}
|
||||
else if (arb_contains_negative(x))
|
||||
{
|
||||
arb_indeterminate(z);
|
||||
}
|
||||
else
|
||||
{
|
||||
mag_init(zr);
|
||||
mag_init(rx);
|
||||
|
||||
/* rx = upper bound for x / r */
|
||||
arf_get_mag_lower(rx, arb_midref(x));
|
||||
mag_div(rx, rx, arb_radref(x));
|
||||
|
||||
inexact = arf_sqrt(arb_midref(z), arb_midref(x), prec, ARB_RND);
|
||||
|
||||
/* zr = upper bound for sqrt(x) */
|
||||
arf_get_mag(zr, arb_midref(z));
|
||||
if (inexact)
|
||||
arf_mag_add_ulp(zr, zr, arb_midref(z), prec);
|
||||
|
||||
/* propagated error: sqrt(x) - sqrt(x-r)
|
||||
= sqrt(x) * [1 - sqrt(1 - r/x)]
|
||||
<= sqrt(x) * 0.5 * (t + t^2), t = r/x */
|
||||
mag_addmul(rx, rx, rx);
|
||||
mag_mul(zr, zr, rx);
|
||||
mag_mul_2exp_si(zr, zr, -1);
|
||||
|
||||
/* add the rounding error */
|
||||
if (inexact)
|
||||
arf_mag_add_ulp(arb_radref(z), zr, arb_midref(z), prec);
|
||||
else
|
||||
mag_swap(arb_radref(z), zr);
|
||||
|
||||
mag_clear(zr);
|
||||
mag_clear(rx);
|
||||
}
|
||||
}
|
||||
|
||||
12
arb/submul.c
12
arb/submul.c
|
|
@ -36,20 +36,16 @@ arb_submul_arf(arb_t z, const arb_t x, const arf_t y, long prec)
|
|||
inexact = arf_submul(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
{
|
||||
if (ARB_IS_LAGOM(z)) /* todo: arf_mag_fast_add_ulp -- or make arf_mag_add_ulp fast? */
|
||||
mag_fast_add_2exp_si(arb_radref(z), arb_radref(z), ARF_EXP(arb_midref(z)) - prec);
|
||||
else
|
||||
arf_mag_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
arf_mag_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
else if (ARB_IS_LAGOM(x) && ARF_IS_LAGOM(y) && ARB_IS_LAGOM(z))
|
||||
{
|
||||
mag_fast_init_set_arf(ym, y);
|
||||
mag_fast_addmul(arb_radref(z), ym, arb_radref(x));
|
||||
inexact = arf_submul(arb_midref(z), arb_midref(x), y, prec, ARB_RND);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(arb_radref(z), arb_radref(z), ARF_EXP(arb_midref(z)) - prec);
|
||||
arf_mag_fast_add_ulp(arb_radref(z), arb_radref(z), arb_midref(z), prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -92,7 +88,7 @@ arb_submul(arb_t z, const arb_t x, const arb_t y, long prec)
|
|||
prec, ARF_RND_DOWN);
|
||||
|
||||
if (inexact)
|
||||
mag_fast_add_2exp_si(zr, zr, ARF_EXP(ARB_MIDREF(z)) - prec);
|
||||
arf_mag_fast_add_ulp(zr, zr, arb_midref(z), prec);
|
||||
|
||||
*ARB_RADREF(z) = *zr;
|
||||
}
|
||||
|
|
|
|||
78
arb/test/t-sqrt.c
Normal file
78
arb/test/t-sqrt.c
Normal file
|
|
@ -0,0 +1,78 @@
|
|||
/*=============================================================================
|
||||
|
||||
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 "arb.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
long iter;
|
||||
flint_rand_t state;
|
||||
|
||||
printf("sqrt....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 100000; iter++)
|
||||
{
|
||||
arb_t a, b, c;
|
||||
long prec = 2 + n_randint(state, 200);
|
||||
|
||||
arb_init(a);
|
||||
arb_init(b);
|
||||
arb_init(c);
|
||||
|
||||
arb_randtest_special(a, state, 1 + n_randint(state, 200), 100);
|
||||
|
||||
arb_sqrt(b, a, prec);
|
||||
arb_mul(c, b, b, prec);
|
||||
|
||||
if (!arb_contains(c, a))
|
||||
{
|
||||
printf("FAIL: containment\n\n");
|
||||
printf("a = "); arb_print(a); printf("\n\n");
|
||||
printf("b = "); arb_print(b); printf("\n\n");
|
||||
printf("c = "); arb_print(c); printf("\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_sqrt(a, a, prec);
|
||||
|
||||
if (!arb_equal(a, b))
|
||||
{
|
||||
printf("FAIL: aliasing\n\n");
|
||||
abort();
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
arb_clear(b);
|
||||
arb_clear(c);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
14
arf.h
14
arf.h
|
|
@ -1167,6 +1167,13 @@ mag_fast_init_set_arf(mag_t y, const arf_t x)
|
|||
}
|
||||
}
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ void
|
||||
arf_mag_fast_add_ulp(mag_t z, const mag_t x, const arf_t y, long prec)
|
||||
{
|
||||
mag_fast_add_2exp_si(z, x, ARF_EXP(y) - prec);
|
||||
}
|
||||
|
||||
/* TODO: document */
|
||||
static __inline__ void
|
||||
arf_mag_add_ulp(mag_t z, const mag_t x, const arf_t y, long prec)
|
||||
|
|
@ -1176,12 +1183,15 @@ arf_mag_add_ulp(mag_t z, const mag_t x, const arf_t y, long prec)
|
|||
printf("error: ulp error not defined for special value!\n");
|
||||
abort();
|
||||
}
|
||||
else if (MAG_IS_LAGOM(z) && MAG_IS_LAGOM(x) && ARF_IS_LAGOM(y))
|
||||
{
|
||||
arf_mag_fast_add_ulp(z, x, y, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* todo: speed up case r is special here */
|
||||
fmpz_t e;
|
||||
fmpz_init(e);
|
||||
_fmpz_add_fast(e, ARF_EXPREF(y), -prec);
|
||||
fmpz_sub_ui(e, ARF_EXPREF(y), prec);
|
||||
mag_add_2exp_fmpz(z, x, e);
|
||||
fmpz_clear(e);
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue