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

exp and log

Browse source
This commit is contained in:
Fredrik Johansson 2012-09-07 12:19:27 +02:00
parent 7ed98fc3fb
commit 0847d21f66
9 changed files with 524 additions and 34 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

24
fmpr.h
View file

@ -410,6 +410,10 @@ long fmpr_sqrt_ui(fmpr_t z, ulong x, long prec, fmpr_rnd_t rnd);
long fmpr_sqrt_fmpz(fmpr_t z, const fmpz_t x, long prec, fmpr_rnd_t rnd);
long fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
long fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
long fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
long fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd);
static __inline__ void
fmpr_neg(fmpr_t y, const fmpr_t x)
@ -529,5 +533,25 @@ void fmpr_get_fmpq(fmpq_t y, const fmpr_t x);
long fmpr_set_fmpq(fmpr_t x, const fmpq_t y, long prec, fmpr_rnd_t rnd);
#define CALL_MPFR_FUNC(r, func, y, x, prec, rnd) \
do { \
mpfr_t t, u; \
long r; \
mpfr_rnd_t rrnd; \
if (rnd == FMPR_RND_DOWN) rrnd = MPFR_RNDZ; \
else if (rnd == FMPR_RND_UP) rrnd = MPFR_RNDA; \
else if (rnd == FMPR_RND_FLOOR) rrnd = MPFR_RNDD; \
else if (rnd == FMPR_RND_CEIL) rrnd = MPFR_RNDU; \
else rrnd = MPFR_RNDN; \
mpfr_init2(t, 2 + fmpz_bits(fmpr_manref(x))); \
mpfr_init2(u, FLINT_MAX(2, prec)); \
fmpr_get_mpfr(t, x, MPFR_RNDD); \
func(u, t, rrnd); \
fmpr_set_mpfr(y, u); \
r = prec - fmpz_bits(fmpr_manref(y)); \
mpfr_clear(t); \
mpfr_clear(u); \
} while (0);
#endif

74
fmpr/exp.c Normal file
View file

@ -0,0 +1,74 @@
/*=============================================================================
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 "fmpr.h"
long
fmpr_exp(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
{
if (fmpr_is_special(x))
{
if (fmpr_is_zero(x))
fmpr_one(y);
else if (fmpr_is_pos_inf(x))
fmpr_pos_inf(y);
else if (fmpr_is_neg_inf(x))
fmpr_zero(y);
else
fmpr_nan(y);
return FMPR_RESULT_EXACT;
}
else
{
long r;
CALL_MPFR_FUNC(r, mpfr_exp, y, x, prec, rnd);
return r;
}
}
long
fmpr_expm1(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
{
if (fmpr_is_special(x))
{
if (fmpr_is_zero(x))
fmpr_zero(y);
else if (fmpr_is_pos_inf(x))
fmpr_pos_inf(y);
else if (fmpr_is_neg_inf(x))
fmpr_set_si(y, -1L);
else
fmpr_nan(y);
return FMPR_RESULT_EXACT;
}
else
{
long r;
CALL_MPFR_FUNC(r, mpfr_expm1, y, x, prec, rnd);
return r;
}
}

57
fmpr/log.c
View file

@ -25,16 +25,6 @@
#include "fmpr.h"
static mpfr_rnd_t fmpr_rnd_to_mpfr(fmpr_rnd_t rnd)
{
if (rnd == FMPR_RND_NEAR) return MPFR_RNDN;
if (rnd == FMPR_RND_FLOOR) return MPFR_RNDD;
if (rnd == FMPR_RND_CEIL) return MPFR_RNDU;
if (rnd == FMPR_RND_DOWN) return MPFR_RNDZ;
if (rnd == FMPR_RND_UP) return MPFR_RNDA;
abort();
}
long
fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
{
@ -49,37 +39,42 @@ fmpr_log(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
return FMPR_RESULT_EXACT;
}
if (fmpr_sgn(x) < 0)
else if (fmpr_sgn(x) < 0)
{
fmpr_nan(y);
return FMPR_RESULT_EXACT;
}
if (fmpr_is_one(x))
else if (fmpr_is_one(x))
{
fmpr_zero(y);
return FMPR_RESULT_EXACT;
}
else
{
mpfr_t t, u;
long r;
mpfr_init2(t, 1 + fmpz_bits(fmpr_manref(x)));
mpfr_init2(u, FLINT_MAX(2, prec));
fmpr_get_mpfr(t, x, MPFR_RNDD);
mpfr_log(u, t, fmpr_rnd_to_mpfr(rnd));
fmpr_set_mpfr(y, u);
r = prec - fmpz_bits(fmpr_manref(y));
mpfr_clear(t);
mpfr_clear(u);
CALL_MPFR_FUNC(r, mpfr_log, y, x, prec, rnd);
return r;
}
}
long
fmpr_log1p(fmpr_t y, const fmpr_t x, long prec, fmpr_rnd_t rnd)
{
if (fmpr_is_special(x))
{
if (fmpr_is_zero(x))
fmpr_zero(y);
else if (fmpr_is_pos_inf(x))
fmpr_pos_inf(y);
else
fmpr_nan(y);
return FMPR_RESULT_EXACT;
}
else
{
long r;
CALL_MPFR_FUNC(r, mpfr_log1p, y, x, prec, rnd);
return r;
}
}

103
fmpr/test/t-exp.c Normal file
View file

@ -0,0 +1,103 @@
/*=============================================================================
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 "fmpr.h"
int main()
{
long iter;
flint_rand_t state;
printf("exp....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
long bits;
fmpr_t x, z, w;
mpfr_t X, Z;
bits = 2 + n_randint(state, 200);
fmpr_init(x);
fmpr_init(z);
fmpr_init(w);
mpfr_init2(X, bits + 100);
mpfr_init2(Z, bits);
fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10);
fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10);
fmpr_get_mpfr(X, x, MPFR_RNDN);
switch (n_randint(state, 4))
{
case 0:
mpfr_exp(Z, X, MPFR_RNDZ);
fmpr_exp(z, x, bits, FMPR_RND_DOWN);
break;
case 1:
mpfr_exp(Z, X, MPFR_RNDA);
fmpr_exp(z, x, bits, FMPR_RND_UP);
break;
case 2:
mpfr_exp(Z, X, MPFR_RNDD);
fmpr_exp(z, x, bits, FMPR_RND_FLOOR);
break;
case 3:
mpfr_exp(Z, X, MPFR_RNDU);
fmpr_exp(z, x, bits, FMPR_RND_CEIL);
break;
}
fmpr_set_mpfr(w, Z);
if (!fmpr_equal(z, w))
{
printf("FAIL\n\n");
printf("bits = %ld\n", bits);
printf("x = "); fmpr_print(x); printf("\n\n");
printf("z = "); fmpr_print(z); printf("\n\n");
printf("w = "); fmpr_print(w); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(z);
fmpr_clear(w);
mpfr_clear(X);
mpfr_clear(Z);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

103
fmpr/test/t-expm1.c Normal file
View file

@ -0,0 +1,103 @@
/*=============================================================================
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 "fmpr.h"
int main()
{
long iter;
flint_rand_t state;
printf("expm1....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
long bits;
fmpr_t x, z, w;
mpfr_t X, Z;
bits = 2 + n_randint(state, 200);
fmpr_init(x);
fmpr_init(z);
fmpr_init(w);
mpfr_init2(X, bits + 100);
mpfr_init2(Z, bits);
fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10);
fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10);
fmpr_get_mpfr(X, x, MPFR_RNDN);
switch (n_randint(state, 4))
{
case 0:
mpfr_expm1(Z, X, MPFR_RNDZ);
fmpr_expm1(z, x, bits, FMPR_RND_DOWN);
break;
case 1:
mpfr_expm1(Z, X, MPFR_RNDA);
fmpr_expm1(z, x, bits, FMPR_RND_UP);
break;
case 2:
mpfr_expm1(Z, X, MPFR_RNDD);
fmpr_expm1(z, x, bits, FMPR_RND_FLOOR);
break;
case 3:
mpfr_expm1(Z, X, MPFR_RNDU);
fmpr_expm1(z, x, bits, FMPR_RND_CEIL);
break;
}
fmpr_set_mpfr(w, Z);
if (!fmpr_equal(z, w))
{
printf("FAIL\n\n");
printf("bits = %ld\n", bits);
printf("x = "); fmpr_print(x); printf("\n\n");
printf("z = "); fmpr_print(z); printf("\n\n");
printf("w = "); fmpr_print(w); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(z);
fmpr_clear(w);
mpfr_clear(X);
mpfr_clear(Z);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

103
fmpr/test/t-log1p.c Normal file
View file

@ -0,0 +1,103 @@
/*=============================================================================
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 "fmpr.h"
int main()
{
long iter;
flint_rand_t state;
printf("log1p....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 100000; iter++)
{
long bits;
fmpr_t x, z, w;
mpfr_t X, Z;
bits = 2 + n_randint(state, 200);
fmpr_init(x);
fmpr_init(z);
fmpr_init(w);
mpfr_init2(X, bits + 100);
mpfr_init2(Z, bits);
fmpr_randtest_special(x, state, bits + n_randint(state, 100), 10);
fmpr_randtest_special(z, state, bits + n_randint(state, 100), 10);
fmpr_get_mpfr(X, x, MPFR_RNDN);
switch (n_randint(state, 4))
{
case 0:
mpfr_log1p(Z, X, MPFR_RNDZ);
fmpr_log1p(z, x, bits, FMPR_RND_DOWN);
break;
case 1:
mpfr_log1p(Z, X, MPFR_RNDA);
fmpr_log1p(z, x, bits, FMPR_RND_UP);
break;
case 2:
mpfr_log1p(Z, X, MPFR_RNDD);
fmpr_log1p(z, x, bits, FMPR_RND_FLOOR);
break;
case 3:
mpfr_log1p(Z, X, MPFR_RNDU);
fmpr_log1p(z, x, bits, FMPR_RND_CEIL);
break;
}
fmpr_set_mpfr(w, Z);
if (!fmpr_equal(z, w))
{
printf("FAIL\n\n");
printf("bits = %ld\n", bits);
printf("x = "); fmpr_print(x); printf("\n\n");
printf("z = "); fmpr_print(z); printf("\n\n");
printf("w = "); fmpr_print(w); printf("\n\n");
abort();
}
fmpr_clear(x);
fmpr_clear(z);
fmpr_clear(w);
mpfr_clear(X);
mpfr_clear(Z);
}
flint_randclear(state);
_fmpz_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

2
fmprb.h
View file

@ -173,6 +173,8 @@ void fmprb_log(fmprb_t z, const fmprb_t x, long prec);
void fmprb_log_ui(fmprb_t z, ulong x, long prec);
void fmprb_log_fmpz(fmprb_t z, const fmpz_t x, long prec);
void fmprb_exp(fmprb_t z, const fmprb_t x, long prec);
void fmprb_fac_ui(fmprb_t x, ulong n, long prec);
void fmprb_bin_ui(fmprb_t x, const fmprb_t n, ulong k, long prec);

61
fmprb/exp.c Normal file
View file

@ -0,0 +1,61 @@
/*=============================================================================
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 "fmprb.h"
void
fmprb_exp(fmprb_t z, const fmprb_t x, long prec)
{
long r;
if (fmprb_is_exact(x))
{
r = fmpr_exp(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
}
else
{
/*
exp(a+b) - exp(a) = exp(a) * (exp(b)-1)
*/
fmpr_t t, u;
fmpr_init(t);
fmpr_init(u);
fmpr_exp(t, fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_expm1(u, fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_mul(t, t, u, FMPRB_RAD_PREC, FMPR_RND_UP);
r = fmpr_exp(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_add_error_result(fmprb_radref(z), t, fmprb_midref(z), r,
FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_clear(t);
fmpr_clear(u);
}
fmprb_adjust(z);
}

31
fmprb/log.c
View file

@ -30,12 +30,37 @@ fmprb_log(fmprb_t z, const fmprb_t x, long prec)
{
long r;
if (!fmprb_is_exact(x))
if (fmprb_contains_zero(x) || fmpr_sgn(fmprb_midref(x)) < 0)
{
printf("log: interval contains zero or negative numbers\n");
abort();
}
if (fmprb_is_exact(x))
{
r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
}
else
{
/*
Let the input be [a-b, a+b]. We require a > b >= 0 (otherwise the
interval contains zero or a negative number and the logarithm is not
defined). The error is largest at a-b, and we have
log(a) - log(a-b) = log(1 + b/(a-b)).
*/
fmpr_t err;
fmpr_init(err);
fmpr_sub(err, fmprb_midref(x), fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_DOWN);
fmpr_div(err, fmprb_radref(x), err, FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_log1p(err, err, FMPRB_RAD_PREC, FMPR_RND_UP);
r = fmpr_log(fmprb_midref(z), fmprb_midref(x), prec, FMPR_RND_DOWN);
fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r);
fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r,
FMPRB_RAD_PREC, FMPR_RND_UP);
fmpr_clear(err);
}
fmprb_adjust(z);
}