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

change arb/acb_mat_bound_inf_norm to output a mag_t; tweak arb/acb_mat_exp

Browse source
This commit is contained in:
Fredrik Johansson 2014-08-15 15:16:07 +02:00
parent 6880e8b223
commit 24b18ec104
12 changed files with 131 additions and 105 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
acb_mat.h
View file

@ -99,7 +99,7 @@ void acb_mat_one(acb_mat_t mat);
/* Norms */
void acb_mat_bound_inf_norm(arf_t b, const acb_mat_t A, long prec);
void acb_mat_bound_inf_norm(mag_t b, const acb_mat_t A);
/* Arithmetic */

22
acb_mat/bound_inf_norm.c
View file

@ -26,37 +26,37 @@
#include "acb_mat.h"
void
acb_mat_bound_inf_norm(arf_t b, const acb_mat_t A, long prec)
acb_mat_bound_inf_norm(mag_t b, const acb_mat_t A)
{
long i, j, r, c;
arf_t s, t;
mag_t s, t;
r = acb_mat_nrows(A);
c = acb_mat_ncols(A);
arf_zero(b);
mag_zero(b);
if (r == 0 || c == 0)
return;
arf_init(s);
arf_init(t);
mag_init(s);
mag_init(t);
for (i = 0; i < r; i++)
{
arf_zero(s);
mag_zero(s);
for (j = 0; j < c; j++)
{
acb_get_abs_ubound_arf(t, acb_mat_entry(A, i, j), prec);
arf_add(s, s, t, prec, ARF_RND_UP);
acb_get_mag(t, acb_mat_entry(A, i, j));
mag_add(s, s, t);
}
arf_max(b, b, s);
mag_max(b, b, s);
}
arf_clear(s);
arf_clear(t);
mag_clear(s);
mag_clear(t);
}

31
acb_mat/exp.c
View file

@ -26,9 +26,7 @@
#include "double_extras.h"
#include "acb_mat.h"
long _arb_mat_exp_choose_N(const arf_t norm, long prec);
void _arb_mat_exp_bound(arf_t err, const arf_t norm, long N);
long _arb_mat_exp_choose_N(const mag_t norm, long prec);
/* evaluates the truncated Taylor series (assumes no aliasing) */
void
@ -123,7 +121,7 @@ void
acb_mat_exp(acb_mat_t B, const acb_mat_t A, long prec)
{
long i, j, dim, wp, N, q, r;
arf_t norm, err;
mag_t norm, err;
acb_mat_t T;
dim = acb_mat_nrows(A);
@ -146,39 +144,38 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, long prec)
wp = prec + 3 * FLINT_BIT_COUNT(prec);
arf_init(norm);
arf_init(err);
mag_init(norm);
mag_init(err);
acb_mat_init(T, dim, dim);
acb_mat_bound_inf_norm(norm, A, MAG_BITS);
acb_mat_bound_inf_norm(norm, A);
if (arf_is_zero(norm))
if (mag_is_zero(norm))
{
acb_mat_one(B);
}
else
{
r = arf_abs_bound_lt_2exp_si(norm);
q = pow(wp, 0.25); /* wanted magnitude */
if (r > 2 * wp) /* too big */
if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */
r = 2 * wp;
else if (r < -q) /* tiny, no need to reduce */
else if (mag_cmp_2exp_si(norm, -q) < 0) /* tiny, no need to reduce */
r = 0;
else
r += q; /* reduce to magnitude 2^(-r) */
r = FLINT_MAX(0, q + MAG_EXP(norm)); /* reduce to magnitude 2^(-r) */
acb_mat_scalar_mul_2exp_si(T, A, -r);
arf_mul_2exp_si(norm, norm, -r);
mag_mul_2exp_si(norm, norm, -r);
N = _arb_mat_exp_choose_N(norm, wp);
_arb_mat_exp_bound(err, norm, N);
mag_exp_tail(err, norm, N);
_acb_mat_exp_taylor(B, T, N, wp);
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
acb_add_error_arf(acb_mat_entry(B, i, j), err);
acb_add_error_mag(acb_mat_entry(B, i, j), err);
for (i = 0; i < r; i++)
{
@ -192,8 +189,8 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, long prec)
acb_mat_entry(B, i, j), prec);
}
arf_clear(norm);
arf_clear(err);
mag_clear(norm);
mag_clear(err);
acb_mat_clear(T);
}

2
arb_mat.h
View file

@ -97,7 +97,7 @@ void arb_mat_one(arb_mat_t mat);
/* Norms */
void arb_mat_bound_inf_norm(arf_t b, const arb_mat_t A, long prec);
void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A);
/* Arithmetic */

22
arb_mat/bound_inf_norm.c
View file

@ -26,37 +26,37 @@
#include "arb_mat.h"
void
arb_mat_bound_inf_norm(arf_t b, const arb_mat_t A, long prec)
arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A)
{
long i, j, r, c;
arf_t s, t;
mag_t s, t;
r = arb_mat_nrows(A);
c = arb_mat_ncols(A);
arf_zero(b);
mag_zero(b);
if (r == 0 || c == 0)
return;
arf_init(s);
arf_init(t);
mag_init(s);
mag_init(t);
for (i = 0; i < r; i++)
{
arf_zero(s);
mag_zero(s);
for (j = 0; j < c; j++)
{
arb_get_abs_ubound_arf(t, arb_mat_entry(A, i, j), prec);
arf_add(s, s, t, prec, ARF_RND_UP);
arb_get_mag(t, arb_mat_entry(A, i, j));
mag_add(s, s, t);
}
arf_max(b, b, s);
mag_max(b, b, s);
}
arf_clear(s);
arf_clear(t);
mag_clear(s);
mag_clear(t);
}

69
arb_mat/exp.c
View file

@ -26,30 +26,30 @@
#include "double_extras.h"
#include "arb_mat.h"
long _fmprb_mat_exp_choose_N(const fmpr_t norm, long prec);
void _fmprb_mat_exp_bound(fmpr_t err, const fmpr_t norm, long N);
#define LOG2_OVER_E 0.25499459743395350926
long _arb_mat_exp_choose_N(const arf_t norm, long prec)
long
_arb_mat_exp_choose_N(const mag_t norm, long prec)
{
fmpr_t fnorm;
long r;
fmpr_init(fnorm);
arf_get_fmpr(fnorm, norm);
r = _fmprb_mat_exp_choose_N(fnorm, prec);
fmpr_clear(fnorm);
return r;
}
if (mag_is_special(norm) || mag_cmp_2exp_si(norm, 30) > 0 ||
mag_cmp_2exp_si(norm, -prec) < 0)
{
return 1;
}
else if (mag_cmp_2exp_si(norm, -300) < 0)
{
long N = -MAG_EXP(norm);
return (prec + N - 1) / N;
}
else
{
double c, t;
void _arb_mat_exp_bound(arf_t err, const arf_t norm, long N)
{
fmpr_t ferr, fnorm;
fmpr_init(ferr);
fmpr_init(fnorm);
arf_get_fmpr(fnorm, norm);
_fmprb_mat_exp_bound(ferr, fnorm, N);
arf_set_fmpr(err, ferr);
fmpr_clear(ferr);
fmpr_clear(fnorm);
c = mag_get_d(norm);
t = d_lambertw(prec * LOG2_OVER_E / c);
t = c * exp(t + 1.0);
return FLINT_MIN((long) (t + 1.0), 2 * prec);
}
}
/* evaluates the truncated Taylor series (assumes no aliasing) */
@ -145,7 +145,7 @@ void
arb_mat_exp(arb_mat_t B, const arb_mat_t A, long prec)
{
long i, j, dim, wp, N, q, r;
arf_t norm, err;
mag_t norm, err;
arb_mat_t T;
dim = arb_mat_nrows(A);
@ -168,39 +168,38 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, long prec)
wp = prec + 3 * FLINT_BIT_COUNT(prec);
arf_init(norm);
arf_init(err);
mag_init(norm);
mag_init(err);
arb_mat_init(T, dim, dim);
arb_mat_bound_inf_norm(norm, A, MAG_BITS);
arb_mat_bound_inf_norm(norm, A);
if (arf_is_zero(norm))
if (mag_is_zero(norm))
{
arb_mat_one(B);
}
else
{
r = arf_abs_bound_lt_2exp_si(norm);
q = pow(wp, 0.25); /* wanted magnitude */
if (r > 2 * wp) /* too big */
if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */
r = 2 * wp;
else if (r < -q) /* tiny, no need to reduce */
else if (mag_cmp_2exp_si(norm, -q) < 0) /* tiny, no need to reduce */
r = 0;
else
r += q; /* reduce to magnitude 2^(-r) */
r = FLINT_MAX(0, q + MAG_EXP(norm)); /* reduce to magnitude 2^(-r) */
arb_mat_scalar_mul_2exp_si(T, A, -r);
arf_mul_2exp_si(norm, norm, -r);
mag_mul_2exp_si(norm, norm, -r);
N = _arb_mat_exp_choose_N(norm, wp);
_arb_mat_exp_bound(err, norm, N);
mag_exp_tail(err, norm, N);
_arb_mat_exp_taylor(B, T, N, wp);
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
arb_add_error_arf(arb_mat_entry(B, i, j), err);
arb_add_error_mag(arb_mat_entry(B, i, j), err);
for (i = 0; i < r; i++)
{
@ -214,8 +213,8 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, long prec)
arb_mat_entry(B, i, j), prec);
}
arf_clear(norm);
arf_clear(err);
mag_clear(norm);
mag_clear(err);
arb_mat_clear(T);
}

5
doc/source/acb_mat.rst
View file

@ -111,11 +111,10 @@ Special matrices
Norms
-------------------------------------------------------------------------------
.. function:: void acb_mat_bound_inf_norm(arf_t b, const acb_mat_t A, long prec)
.. function:: void acb_mat_bound_inf_norm(mag_t b, const acb_mat_t A)
Sets *b* to an upper bound for the infinity norm (i.e. the largest
absolute value row sum) of *A*, computed using floating-point arithmetic
at *prec* bits with all operations rounded up.
absolute value row sum) of *A*.
Arithmetic

5
doc/source/arb_mat.rst
View file

@ -111,11 +111,10 @@ Special matrices
Norms
-------------------------------------------------------------------------------
.. function:: void arb_mat_bound_inf_norm(arf_t b, const arb_mat_t A, long prec)
.. function:: void arb_mat_bound_inf_norm(mag_t b, const arb_mat_t A)
Sets *b* to an upper bound for the infinity norm (i.e. the largest
absolute value row sum) of *A*, computed using floating-point arithmetic
at *prec* bits with all operations rounded up.
absolute value row sum) of *A*.
Arithmetic
-------------------------------------------------------------------------------

21
hypgeom/bound.c
View file

@ -243,27 +243,6 @@ mag_sub_lower(mag_t z, const mag_t x, const mag_t y)
}
}
double
mag_get_d(const mag_t z)
{
long exp;
if (mag_is_zero(z))
return 0.0;
else if (mag_is_inf(z))
return D_INF;
/* XXX: overflow/underflow properly */
exp = MAG_EXP(z);
if (exp < -1000)
return ldexp(1.0, -1000);
else if (exp > 1000)
return D_INF;
else
return ldexp(MAG_MAN(z), exp - MAG_BITS);
}
long
hypgeom_bound(mag_t error, int r,
long A, long B, long K, const mag_t TK, const mag_t z, long tol_2exp)

2
hypgeom/estimate_terms_d.c
View file

@ -39,8 +39,6 @@ static __inline__ double d_root(double x, int r)
return pow(x, 1. / r);
}
double mag_get_d(const mag_t z);
long
hypgeom_estimate_terms(const mag_t z, int r, long prec)
{

3
mag.h
View file

@ -523,6 +523,9 @@ static __inline__ void mag_set_d(mag_t z, double x)
fmpz_clear(e);
}
/* TODO: test/document */
double mag_get_d(const mag_t z);
/* TODO: document */
double mag_d_log_upper_bound(double x);
double mag_d_log_lower_bound(double x);

52
mag/get_d.c Normal file
View file

@ -0,0 +1,52 @@
/*=============================================================================
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 "double_extras.h"
#include "mag.h"
double
mag_get_d(const mag_t z)
{
if (mag_is_zero(z))
{
return 0.0;
}
else if (mag_is_inf(z))
{
return D_INF;
}
else if (MAG_EXP(z) < -1000 || MAG_EXP(z) > 1000)
{
if (fmpz_sgn(MAG_EXPREF(z)) < 0)
return ldexp(1.0, -1000);
else
return D_INF;
}
else
{
return ldexp(MAG_MAN(z), MAG_EXP(z) - MAG_BITS);
}
}