src/log10.c - master - niceSOFT/mpc

src/log10.c
159 строк6 KB
Paul Zimmermann
review increase of precision in mpc_log
10 сен 2024, 12:02
10 сен 2024, 12:024c62e36
/* mpc_log10 -- Take the base-10 logarithm of a complex number.

Copyright (C) 2012, 2020, 2024 INRIA

This file is part of GNU MPC.

GNU MPC is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

GNU MPC 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 Lesser General Public License for
more details.

You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see http://logw.gnu.org/licenses/ .
*/

#include <limits.h> /* for CHAR_BIT */
#include "mpc-impl.h"

static void
mpfr_const_log10 (mpfr_ptr log10)
{
   mpfr_set_ui (log10, 10, MPFR_RNDN); /* exact since prec >= 4 */
   mpfr_log (log10, log10, MPFR_RNDN); /* error <= 1/2 ulp */
}


int
mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
{
   int ok = 0, loop = 0, check_exact = 0, special_re, special_im,
       inex, inex_re, inex_im;
   mpfr_prec_t prec;
   mpfr_t log10;
   mpc_t log;
   mpfr_exp_t saved_emin, saved_emax;

   saved_emin = mpfr_get_emin ();
   saved_emax = mpfr_get_emax ();
   mpfr_set_emin (mpfr_get_emin_min ());
   mpfr_set_emax (mpfr_get_emax_max ());

   mpfr_init2 (log10, 2);
   mpc_init2 (log, 2);
   prec = MPC_MAX_PREC (rop);
   /* compute log(op)/log(10) */
   while (ok == 0) {
      MPC_LOOP_NEXT(loop, op, rop);
      prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
      mpfr_set_prec (log10, prec);
      mpc_set_prec (log, prec);

      inex = mpc_log (log, op, rnd); /* error <= 1 ulp */

      if (!mpfr_number_p (mpc_imagref (log))
         || mpfr_zero_p (mpc_imagref (log))) {
         /* no need to divide by log(10) */
         special_im = 1;
         ok = 1;
      }
      else {
         special_im = 0;
         mpfr_const_log10 (log10);
         mpfr_div (mpc_imagref (log), mpc_imagref (log), log10, MPFR_RNDN);

         ok = mpfr_can_round (mpc_imagref (log), prec - 2,
                  MPFR_RNDN, MPFR_RNDZ,
                  MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN));
      }

      if (ok) {
         if (!mpfr_number_p (mpc_realref (log))
            || mpfr_zero_p (mpc_realref (log)))
            special_re = 1;
         else {
            special_re = 0;
            if (special_im)
               /* log10 not yet computed */
               mpfr_const_log10 (log10);
            mpfr_div (mpc_realref (log), mpc_realref (log), log10, MPFR_RNDN);
               /* error <= 24/7 ulp < 4 ulp for prec >= 4, see algorithms.tex */

            ok = mpfr_can_round (mpc_realref (log), prec - 2,
                     MPFR_RNDN, MPFR_RNDZ,
                     MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == MPFR_RNDN));
         }

         /* Special code to deal with cases where the real part of log10(x+i*y)
            is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2
            this happens whenever x^2+y^2 is a nonnegative power of 10.
            Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0,
            since x^2+y^2 is rational, and 10^(a/2^b) is irrational.
            Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2
            is a rational with denominator a power of 2.
            Now let x^2+y^2 = 10^s. Without loss of generality we can assume
            x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e)
            thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily
            u = v = 0 mod 2^e, thus x and y are necessarily integers.
         */
         if (!ok && !check_exact && mpfr_integer_p (mpc_realref (op)) &&
            mpfr_integer_p (mpc_imagref (op))) {
            mpz_t x, y;
            unsigned long s, v;

            check_exact = 1;
            mpz_init (x);
            mpz_init (y);
            mpfr_get_z (x, mpc_realref (op), MPFR_RNDN); /* exact */
            mpfr_get_z (y, mpc_imagref (op), MPFR_RNDN); /* exact */
            mpz_mul (x, x, x);
            mpz_mul (y, y, y);
            mpz_add (x, x, y); /* x^2+y^2 */
            v = mpz_scan1 (x, 0);
            /* if x = 10^s then necessarily s = v */
            s = mpz_sizeinbase (x, 10);
            /* since s is either the number of digits of x or one more,
               then x = 10^(s-1) or 10^(s-2) */
            if (s == v + 1 || s == v + 2) {
               mpz_div_2exp (x, x, v);
               mpz_ui_pow_ui (y, 5, v);
               if (mpz_cmp (y, x) == 0) {
                  /* Re(log10(x+i*y)) is exactly v/2
                     we reset the precision of Re(log) so that v can be
                     represented exactly */
                  mpfr_set_prec (mpc_realref (log),
                                 sizeof(unsigned long)*CHAR_BIT);
                  mpfr_set_ui_2exp (mpc_realref (log), v, -1, MPFR_RNDN);
                     /* exact */
                  ok = 1;
               }
            }
            mpz_clear (x);
            mpz_clear (y);
         }
      }
   }

   inex_re = mpfr_set (mpc_realref(rop), mpc_realref (log), MPC_RND_RE (rnd));
   if (special_re)
      inex_re = MPC_INEX_RE (inex);
      /* recover flag from call to mpc_log above */
   inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (log), MPC_RND_IM (rnd));
   if (special_im)
      inex_im = MPC_INEX_IM (inex);
   mpfr_clear (log10);
   mpc_clear (log);

   /* restore the exponent range, and check the range of results */
   mpfr_set_emin (saved_emin);
   mpfr_set_emax (saved_emax);
   inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd));
   inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd));

   return MPC_INEX(inex_re, inex_im);
}