diff options
| -rw-r--r-- | ChangeLog | 13 | ||||
| -rw-r--r-- | TODO | 3 | ||||
| -rw-r--r-- | doc/algorithms.tex | 15 | ||||
| -rw-r--r-- | doc/mpc.texi | 108 | ||||
| -rw-r--r-- | src/log.c | 2 | ||||
| -rw-r--r-- | tests/tadd.c | 2 | ||||
| -rw-r--r-- | tests/tdiv.c | 1 | ||||
| -rw-r--r-- | tests/tlog.c | 67 | ||||
| -rw-r--r-- | tests/tsub.c | 2 |
9 files changed, 157 insertions, 56 deletions
@@ -1,16 +1,19 @@ -2008-03-21 Philippe THEVENY <philippe DOT theveny AT loria DOT fr> +2008-04-11 Andreas Enge <enge@lix.polytechnique.fr> + * new function mpc_log + * bug fix in mpc_sqrt with directed rounding + * algorithms and error analysis document + +2008-03-27 Philippe THEVENY <philippe DOT theveny AT loria DOT fr> + * new function mpc_cos +2008-03-21 Philippe THEVENY <philippe DOT theveny AT loria DOT fr> * INSTALL: short description of the new compilation commands * mpc.texi: deeper explanations on the compilation options * configure.ac: add options selecting GMP/MPFR installation directory 2008-03-14 Philippe THEVENY <philippe DOT theveny AT loria DOT fr> - * NEWS: Added for autotools support * AUTHORS: Added for autotools support * Makefile.am: Added for autotools support * configure.ac: Added for autotools support - -Changes since version 0.4.6: - nothing worth noting for the moment @@ -29,8 +29,7 @@ New functions to implement: - from Andreas Enge and Philippe Théveny 9 April 2008 work plan for completing functionality 1) tan - 2) log and agm - 3) inverse trigonometric functions, require log + 2) inverse trigonometric functions, require log hyperbolic functions and cproj at any moment (relatively easy) - from Andreas Enge and Philippe Théveny 9 April 2008 correct handling of Nan and infinities: diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 0a2a53c..6310a77 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -6,11 +6,14 @@ \usepackage{amsmath,amssymb} \newcommand {\mpc}{\texttt {mpc}} +\newcommand {\mpfr}{\texttt {mpfr}} \newcommand {\ulp}[1]{#1~ulp} +\newcommand {\atantwo}{\operatorname {atan2}} + \title {MPC: Algorithms and Error Analysis} -\author {The {\mpc} team} -\date {April 10, 2008} +\author {Andreas Enge} +\date {April 11, 2008} \begin {document} @@ -47,4 +50,12 @@ t + i w & \text {if } x < 0, y > 0 \\ $w$ is rounded down. $\sqrt {x^2 + y^2}$ is computed with an error of \ulp{1}; $|x|$ is added with an error of \ulp{1}, since both terms are positive. The generic error of the real square root in the special case that the argument was rounded down is \ulp{1}, so that the total error in computing $w$ is \ulp{3}. $t$ is rounded up. The generic error of real division, applied to an error of \ulp{3} for $w$ and \ulp{0} for $y$ implies an error of \ulp{7}. + + +\subsection {\texttt {mpc\_log}} + +Let $z = x + i y$. Then $\log (z) = \frac {1}{2} \log (x^2 + y^2) + i \atantwo (y, x)$. The imaginary part is computed by a call to the corresponding {\mpfr} function. + +Let $w = \log (x^2 + y^2)$, rounded down. The error of the complex norm is \ulp{1}. The generic error of the real logarithm is then given by \ulp{$2^{2 - e_w} + 1$}, where $e_w$ is the exponent of $w$. For $e_w \geq 2$, this is bounded by \ulp{2} or 2~digits; otherwise, it is bounded by \ulp{$2^{3 - e_w}$} or $3 - e_w$ digits. + \end {document} diff --git a/doc/mpc.texi b/doc/mpc.texi index ddd9d41..158c81d 100644 --- a/doc/mpc.texi +++ b/doc/mpc.texi @@ -154,8 +154,8 @@ modified MPC is provided. Everyone should read @ref{MPC Basics}. If you need to install the library yourself, you need to read @ref{Installing MPC}, too. -The rest of the manual can be used for later reference, although it is -probably a good idea to glance through it. +The remainder of the manual can be used for later reference, although it is +probably a good idea to skim through it. @node Installing MPC, Reporting Bugs, Introduction to MPC, Top @comment node-name, next, previous, up @@ -424,14 +424,16 @@ of the IEEE P754 arithmetic. The results obtained on one computer should not differ from the results obtained on a computer with a different word size. + @menu * Rounding Modes:: * Initializing Complex Numbers:: * Assigning Complex Numbers:: * Simultaneous Init & Assign:: -* Complex Arithmetic:: * Complex Comparison:: -* Special Complex Functions:: +* Basic Arithmetic:: +* Power Functions and Logarithm:: +* Trigonometric Functions:: * I/O of Complex Numbers:: * Miscellaneous Complex Functions:: * Internals:: @@ -554,6 +556,7 @@ Returns the precision of the real part of @var{x} via @var{pr} and of its imagin via @var{pi}. @end deftypefun + @node Assigning Complex Numbers, Simultaneous Init & Assign, Initializing Complex Numbers, Complex Functions @comment node-name, next, previous, up @section Assignment Functions @@ -588,7 +591,7 @@ Set the real part of @var{rop} from @var{op1}, and its imaginary part from @end deftypefun -@node Simultaneous Init & Assign, Complex Arithmetic, Assigning Complex Numbers, Complex Functions +@node Simultaneous Init & Assign, Complex Comparison, Assigning Complex Numbers, Complex Functions @comment node-name, next, previous, up @section Combined Initialization and Assignment Functions @cindex Initialization and assignment functions @@ -610,7 +613,35 @@ The precision of @var{rop} will be taken from the active default precision, as set by @code{mpc_set_default_prec}. @end deftypefn -@node Complex Arithmetic, Complex Comparison, Simultaneous Init & Assign, Complex Functions + + +@node Complex Comparison, Basic Arithmetic, Simultaneous Init & Assign, Complex Functions +@comment node-name, next, previous, up +@section Comparison Functions +@cindex Complex comparisons functions +@cindex Comparison functions + +@deftypefn Function int mpc_cmp (mpc_t @var{op1}, mpc_t @var{op2}) +@deftypefnx Function int mpc_cmp_si_si (mpc_t @var{op1}, long int @var{op2r}, long int @var{op2i}) +@deftypefnx Macro int mpc_cmp_si (mpc_t @var{op1}, long int @var{op2}) + +Compare @var{op1} and @var{op2}, where in the case of @code{mpc_cmp_si_si}, +@var{op2} is taken to be @var{op2r} + i @var{op2i}. +The return value @var{c} can be decomposed into @code{x = MPC_INEX_RE(c)} +and @code{y = MPC_INEX_IM(c)}, such that @var{x} is +positive if the real part of @var{op1} is greater than that of @var{op2}, +zero if both real parts are equal, and negative if the real part of @var{op1} +is less than that of @var{op2}, and likewise for @var{y}. +Both @var{op1} and @var{op2} are considered to their full own precision, +which may differ. +It is not allowed that one of the operands has a NaN (Not-a-Number) part. + +The storage of the return value is such that equality can be simply checked +with @code{mpc_cmp (op1, op2) == 0}. +@end deftypefn + + +@node Basic Arithmetic, Power Functions and Logarithm, Complex Comparison, Complex Functions @comment node-name, next, previous, up @section Basic Arithmetic Functions @cindex Complex arithmetic functions @@ -661,12 +692,6 @@ For @code{mpc_div} and @code{mpc_ui_div}, the return value may fail to recognize some exact results, and its sign is not significant. @end deftypefun -@deftypefun int mpc_sqrt (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) -Set @var{rop} to the square root of @var{op} rounded according to @var{rnd}. -Here, when the return value is 0, it means the result is exact, but the -contrary may be false. -@end deftypefun - @deftypefun int mpc_neg (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) Set @var{rop} to @minus{}@var{op} rounded according to @var{rnd}. Just changes the sign if @var{rop} and @var{op} are the same variable. @@ -707,49 +732,50 @@ of the real and imaginary parts by @var{op2} when @var{rop} and @var{op1} are identical. @end deftypefun -@node Complex Comparison, Special Complex Functions, Complex Arithmetic, Complex Functions + +@node Power Functions and Logarithm, Trigonometric Functions, Basic Arithmetic, Complex Functions @comment node-name, next, previous, up -@section Comparison Functions -@cindex Complex comparisons functions -@cindex Comparison functions +@section Power Functions and Logarithm +@cindex Power functions +@cindex Logarithm -@deftypefn Function int mpc_cmp (mpc_t @var{op1}, mpc_t @var{op2}) -@deftypefnx Function int mpc_cmp_si_si (mpc_t @var{op1}, long int @var{op2r}, long int @var{op2i}) -@deftypefnx Macro int mpc_cmp_si (mpc_t @var{op1}, long int @var{op2}) +@deftypefun int mpc_sqrt (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the square root of @var{op} rounded according to @var{rnd}. +Here, when the return value is 0, it means the result is exact, but the +contrary may be false. +@end deftypefun -Compare @var{op1} and @var{op2}, where in the case of @code{mpc_cmp_si_si}, -@var{op2} is taken to be @var{op2r} + i @var{op2i}. -The return value @var{c} can be decomposed into @code{x = MPC_INEX_RE(c)} -and @code{y = MPC_INEX_IM(c)}, such that @var{x} is -positive if the real part of @var{op1} is greater than that of @var{op2}, -zero if both real parts are equal, and negative if the real part of @var{op1} -is less than that of @var{op2}, and likewise for @var{y}. -Both @var{op1} and @var{op2} are considered to their full own precision, -which may differ. -It is not allowed that one of the operands has a NaN (Not-a-Number) part. +@deftypefun void mpc_exp (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the exponential of @var{op}, +rounded according to @var{rnd} with the precision of @var{rop}. +@end deftypefun -The storage of the return value is such that equality can be simply checked -with @code{mpc_cmp (op1, op2) == 0}. -@end deftypefn +@deftypefun void mpc_log (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the logarithm of @var{op}, +rounded according to @var{rnd} with the precision of @var{rop}. +The principal branch is chosen, with the branch cut on the negative real axis, +so that the imaginary part of the result lies in +@math{]-\pi , \pi]}. +@end deftypefun -@node Special Complex Functions, I/O of Complex Numbers, Complex Comparison, Complex Functions +@node Trigonometric Functions, I/O of Complex Numbers, Power Functions and Logarithm, Complex Functions @comment node-name, next, previous, up -@section Special Functions -@cindex Special functions -@cindex Special complex functions +@section Trigonometric Functions +@cindex Trigonometric functions -@deftypefun void mpc_exp (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) -Set @var{rop} to the exponential of @var{op}, +@deftypefun void mpc_sin (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the sine of @var{op}, rounded according to @var{rnd} with the precision of @var{rop}. @end deftypefun -@deftypefun void mpc_sin (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) -Set @var{rop} to the sine of @var{op}, +@deftypefun void mpc_cos (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the cosine of @var{op}, rounded according to @var{rnd} with the precision of @var{rop}. @end deftypefun -@node I/O of Complex Numbers, Miscellaneous Complex Functions, Special Complex Functions, Complex Functions + +@node I/O of Complex Numbers, Miscellaneous Complex Functions, Trigonometric Functions, Complex Functions @comment node-name, next, previous, up @section Input and Output Functions @cindex Complex input and output functions @@ -102,7 +102,7 @@ mpc_log (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) mpc_norm (w, b, GMP_RNDD); /* error 1 ulp */ mpfr_log (w, w, GMP_RNDD); - /* generic error of log: (2^(2 - exp(w)) + 0.5) ulp */ + /* generic error of log: (2^(2 - exp(w)) + 1) ulp */ if (MPFR_EXP (w) >= 2) ok = mpfr_can_round (w, prec - 2, GMP_RNDD, MPC_RND_RE(rnd), MPC_PREC_RE(a)); diff --git a/tests/tadd.c b/tests/tadd.c index e9e2bab..709215c 100644 --- a/tests/tadd.c +++ b/tests/tadd.c @@ -1,6 +1,6 @@ /* test file for mpc_add. -Copyright (C) 2008 INRIA. +Copyright (C) 2008 Philippe Th\'eveny This file is part of the MPC Library. diff --git a/tests/tdiv.c b/tests/tdiv.c index 7568f3e..ef3ac3d 100644 --- a/tests/tdiv.c +++ b/tests/tdiv.c @@ -26,7 +26,6 @@ MA 02111-1307, USA. */ #include "mpc.h" #include "mpc-impl.h" -#define OUT(x) { mpc_out_str (stdout, 2, 0, x, MPC_RNDNN); putchar ('\n'); } #define ERR(x) { mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); \ fprintf (stderr, "\n"); } diff --git a/tests/tlog.c b/tests/tlog.c index 5c2f567..d60f312 100644 --- a/tests/tlog.c +++ b/tests/tlog.c @@ -31,7 +31,70 @@ MA 02111-1307, USA. */ int main() { - tgeneric (); + mpc_t z, z2, tmp; + mpfr_t twopi; + mp_prec_t prec; - return 0; + tgeneric (); + + mpc_init (z); + mpc_init (z2); + mpc_init (tmp); + mpfr_init2 (twopi, 1000); + mpfr_const_pi (twopi, GMP_RNDD); + mpfr_div_2ui (twopi, twopi, 1, GMP_RNDN); + + /* Test whether exp (log (z)) = z */ + for (prec = 2; prec <= 1000; prec++) + { + mpc_set_prec (z, prec); + mpc_set_prec (z2, prec); + mpc_set_prec (tmp, 4*prec); + + mpc_random (z); + mpc_log (tmp, z, MPC_RNDNN); + mpc_exp (z2, tmp, MPC_RNDNN); + + if (mpc_cmp (z, z2) != 0) + { + printf ("Possible error in log; difference between z and z2=exp(log(z)):\n"); + OUT (z); + OUT (tmp); + OUT (z2); + exit (1); + } + } + + + /* Test whether log (exp (z)) = z for purely imaginary z; then exp (x) */ + /* lies on the unit cercle, a critical case for the logarithm. */ + for (prec = 2; prec <= 1000; prec++) + { + mpc_set_prec (z, prec); + mpc_set_prec (z2, prec); + mpc_set_prec (tmp, 4*prec); + + mpfr_set_ui (MPC_RE (z), 0, GMP_RNDN); + mpfr_random (MPC_IM (z)); + mpfr_remainder (MPC_IM (z), MPC_IM (z), twopi, GMP_RNDZ); + mpc_exp (tmp, z, MPC_RNDNN); + mpc_log (z2, tmp, MPC_RNDNN); + + /* There is a tiny real part, do not care if it si sufficiently small. */ + if (mpfr_cmp (MPC_IM (z), MPC_IM (z2)) != 0 || MPFR_EXP (MPC_RE (z)) > -4*prec) + { + printf ("Possible error in log; difference between z and z2=log(exp(z)):\n"); + OUT (z); + OUT (tmp); + OUT (z2); + //exit (1); + } + } + + mpc_clear (z); + mpc_clear (z2); + mpc_clear (tmp); + mpfr_clear (twopi); + + return 0; } diff --git a/tests/tsub.c b/tests/tsub.c index b00f94f..99dc5b5 100644 --- a/tests/tsub.c +++ b/tests/tsub.c @@ -1,6 +1,6 @@ /* test file for mpc_sub. -Copyright (C) 2008 INRIA. +Copyright (C) 2008 Philippe Th\'eveny This file is part of the MPC Library. |