add mag_get_d_log2_approx as a public method

This commit is contained in:
Fredrik Johansson 2017-02-22 11:21:58 +01:00
parent f0570449f9
commit f6d26a52e3
9 changed files with 59 additions and 48 deletions

View file

@ -11,8 +11,6 @@
#include "acb_hypgeom.h"
double mag_get_log2_d_approx(const mag_t x);
#define D_ABS(x) ((x) < 0.0 ? (-(x)) : (x))
int
@ -108,7 +106,7 @@ acb_hypgeom_pfq_choose_n_max(acb_srcptr a, slong p,
bim = bre + q;
acb_get_mag(zmag, z);
log2_z = mag_get_log2_d_approx(zmag);
log2_z = mag_get_d_log2_approx(zmag);
n_skip = 1;
n_min = 1;
@ -217,7 +215,7 @@ acb_hypgeom_pfq_series_choose_n(const acb_poly_struct * a, slong p,
bim = bre + q;
acb_get_mag(zmag, z->coeffs);
log2_z = mag_get_log2_d_approx(zmag);
log2_z = mag_get_d_log2_approx(zmag);
n_skip = 1;
n_min = 1;

View file

@ -11,8 +11,6 @@
#include "acb_modular.h"
double mag_get_log2_d_approx(const mag_t x);
static const int pentagonal_best_m[] = {
2, 5, 7, 11, 13, 17, 19, 23, 25, 35,
55, 65, 77, 91, 119, 133, 143, 175, 275, 325,
@ -227,7 +225,7 @@ acb_modular_eta_sum(acb_t eta, const acb_t q, slong prec)
q_is_real = arb_is_zero(acb_imagref(q));
acb_get_mag(qmag, q);
log2q_approx = mag_get_log2_d_approx(qmag);
log2q_approx = mag_get_d_log2_approx(qmag);
if (log2q_approx >= 0.0)
{

View file

@ -11,8 +11,6 @@
#include "acb_modular.h"
double mag_get_log2_d_approx(const mag_t x);
void
acb_modular_theta_const_sum(acb_t theta2, acb_t theta3, acb_t theta4,
const acb_t q, slong prec)
@ -26,7 +24,7 @@ acb_modular_theta_const_sum(acb_t theta2, acb_t theta3, acb_t theta4,
mag_init(err);
acb_get_mag(qmag, q);
log2q_approx = mag_get_log2_d_approx(qmag);
log2q_approx = mag_get_d_log2_approx(qmag);
is_real = arb_is_zero(acb_imagref(q));
is_real_or_imag = is_real || arb_is_zero(acb_realref(q));

View file

@ -11,8 +11,6 @@
#include "acb_modular.h"
double mag_get_log2_d_approx(const mag_t x);
void
acb_modular_theta_const_sum_basecase(acb_t theta2, acb_t theta3, acb_t theta4,
const acb_t q, slong N, slong prec)
@ -125,7 +123,7 @@ acb_modular_theta_const_sum_basecase(acb_t theta2, acb_t theta3, acb_t theta4,
acb_modular_fill_addseq(tab, N);
acb_get_mag(qmag, q);
log2q_approx = mag_get_log2_d_approx(qmag);
log2q_approx = mag_get_d_log2_approx(qmag);
for (k = 0; k < N; k++)
{

View file

@ -11,8 +11,6 @@
#include "acb_modular.h"
double mag_get_log2_d_approx(const mag_t x);
static const int square_best_m[] = {
2, 3, 4, 8, 12, 16, 32, 48, 80, 96,
112, 144, 240, 288, 336, 480, 560, 576, 720, 1008,
@ -95,7 +93,7 @@ acb_modular_theta_const_sum_rs(acb_t theta2, acb_t theta3, acb_t theta4,
mag_init(qmag);
acb_get_mag(qmag, q);
log2q_approx = mag_get_log2_d_approx(qmag);
log2q_approx = mag_get_d_log2_approx(qmag);
mag_clear(qmag);
acb_init(tmp1);

View file

@ -11,36 +11,6 @@
#include "acb_modular.h"
double
mag_get_log2_d_approx(const mag_t x)
{
if (mag_is_zero(x))
{
return COEFF_MIN;
}
else if (mag_is_inf(x))
{
return COEFF_MAX;
}
else if (COEFF_IS_MPZ(MAG_EXP(x)))
{
if (fmpz_sgn(MAG_EXPREF(x)) < 0)
return COEFF_MIN;
else
return COEFF_MAX;
}
else
{
slong e = MAG_EXP(x);
if (e < -20 || e > 20)
return e;
else
return e + 1.4426950408889634074 *
mag_d_log_upper_bound(MAG_MAN(x) * ldexp(1.0, -MAG_BITS));
}
}
void
acb_modular_theta_sum(acb_ptr theta1,
acb_ptr theta2,
@ -84,7 +54,7 @@ acb_modular_theta_sum(acb_ptr theta1,
acb_inv(v, w, prec);
acb_get_mag(qmag, q);
log2q_approx = mag_get_log2_d_approx(qmag);
log2q_approx = mag_get_d_log2_approx(qmag);
if (w_is_unit)
{
@ -97,7 +67,7 @@ acb_modular_theta_sum(acb_ptr theta1,
acb_get_mag(wmag, w);
acb_get_mag(vmag, v);
mag_max(wmag, wmag, vmag);
log2w_approx = mag_get_log2_d_approx(wmag);
log2w_approx = mag_get_d_log2_approx(wmag);
}
if (log2q_approx >= 0.0)

View file

@ -175,6 +175,12 @@ Conversions
Returns a *double* giving an upper bound for *x*.
.. function:: double mag_get_d_log2_approx(const mag_t x)
Returns a *double* approximating `\log_2(x)`, suitable for estimating
magnitudes (not for rigorous bounds).
The value is clamped between COEFF_MIN and COEFF_MAX.
.. function:: void mag_set_ui_lower(mag_t z, ulong x)
.. function:: void mag_set_fmpz_lower(mag_t z, const fmpz_t x)

2
mag.h
View file

@ -546,6 +546,8 @@ MAG_INLINE void mag_set_d(mag_t z, double x)
double mag_get_d(const mag_t z);
double mag_get_d_log2_approx(const mag_t x);
/* TODO: document */
double mag_d_log_upper_bound(double x);
double mag_d_log_lower_bound(double x);

43
mag/get_d_log2_approx.c Normal file
View file

@ -0,0 +1,43 @@
/*
Copyright (C) 2014 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "mag.h"
double
mag_get_d_log2_approx(const mag_t x)
{
if (mag_is_zero(x))
{
return COEFF_MIN;
}
else if (mag_is_inf(x))
{
return COEFF_MAX;
}
else if (COEFF_IS_MPZ(MAG_EXP(x)))
{
if (fmpz_sgn(MAG_EXPREF(x)) < 0)
return COEFF_MIN;
else
return COEFF_MAX;
}
else
{
slong e = MAG_EXP(x);
if (e < -20 || e > 20)
return e;
else
return e + 1.4426950408889634074 *
mag_d_log_upper_bound(MAG_MAN(x) * ldexp(1.0, -MAG_BITS));
}
}