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

add mag_log1p

Browse source
This commit is contained in:
Fredrik Johansson 2014-06-18 16:50:15 +02:00
parent 498197bae9
commit 270a65569b
3 changed files with 197 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
mag.h
View file

@ -648,6 +648,9 @@ static __inline__ void mag_set_d(mag_t z, double x)
double mag_d_log_upper_bound(double x);
double mag_d_log_lower_bound(double x);
/* TODO: document */
void mag_log1p(mag_t z, const mag_t x);
/* TODO: document/test */
void mag_pow_ui(mag_t z, const mag_t x, ulong e);
void mag_pow_ui_lower(mag_t z, const mag_t x, ulong e);

96
mag/log1p.c Normal file
View file

@ -0,0 +1,96 @@
/*=============================================================================
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 "mag.h"
#include "double_extras.h"
void
mag_log1p(mag_t z, const mag_t x)
{
if (mag_is_special(x))
{
if (mag_is_zero(x))
mag_zero(z);
else
mag_inf(z);
}
else
{
fmpz exp = MAG_EXP(x);
if (!COEFF_IS_MPZ(exp))
{
/* Quick bound by x */
if (exp < -10)
{
mag_set(z, x);
return;
}
else if (exp < 1000)
{
double t;
t = ldexp(MAG_MAN(x), exp - MAG_BITS);
t = (1.0 + t) * (1 + 1e-14);
t = mag_d_log_upper_bound(t);
mag_set_d(z, t);
return;
}
}
else if (fmpz_sgn(MAG_EXPREF(x)) < 0)
{
/* Quick bound by x */
mag_set(z, x);
return;
}
/* Now we must have x >= 2^1000 */
/* Use log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v) */
{
double t;
fmpz_t b;
mag_t u;
mag_init(u);
fmpz_init(b);
/* incrementing the mantissa gives an upper bound for x+1 */
t = ldexp(MAG_MAN(x) + 1, 1 - MAG_BITS);
t = mag_d_log_upper_bound(t);
mag_set_d(u, t);
/* log(2) < 744261118/2^30 */
_fmpz_add_fast(b, MAG_EXPREF(x), -1);
fmpz_mul_ui(b, b, 744261118);
mag_set_fmpz(z, b);
_fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(z), -30);
mag_add(z, z, u);
mag_clear(u);
fmpz_clear(b);
}
}
}

98
mag/test/t-log1p.c Normal file
View file

@ -0,0 +1,98 @@
/*=============================================================================
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 "mag.h"
int main()
{
long iter;
flint_rand_t state;
printf("log1p....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
fmpr_t x, y, z, z2;
mag_t xb, yb;
fmpr_init(x);
fmpr_init(y);
fmpr_init(z);
fmpr_init(z2);
mag_init(xb);
mag_init(yb);
mag_randtest_special(xb, state, 25);
mag_randtest_special(yb, state, 25);
mag_log1p(yb, xb);
mag_get_fmpr(x, xb);
mag_get_fmpr(y, yb);
fmpr_log1p(z, x, MAG_BITS, FMPR_RND_UP);
fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP);
fmpr_mul_2exp_si(z2, z2, -10);
MAG_CHECK_BITS(xb)
MAG_CHECK_BITS(yb)
if (!(fmpr_cmpabs(z, y) <= 0 && fmpr_cmpabs(y, z2) <= 0))
{
printf("FAIL\n\n");
printf("x = "); fmpr_print(x); printf("\n\n");
printf("y = "); fmpr_print(y); printf("\n\n");
printf("z = "); fmpr_print(z); printf("\n\n");
printf("z2 = "); fmpr_print(z2); printf("\n\n");
abort();
}
mag_log1p(xb, xb);
if (!mag_equal(xb, yb))
{
printf("FAIL (aliasing)\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(y);
fmpr_clear(z);
fmpr_clear(z2);
mag_clear(xb);
mag_clear(yb);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}