contrib/mpc/src/sin_cos.c

   1 /* mpc_sin_cos -- combined sine and cosine of a complex number.
   2
   3 Copyright (C) 2010, 2011 INRIA
   4
   5 This file is part of GNU MPC.
   6
   7 GNU MPC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
  13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
  15 more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with this program. If not, see http://www.gnu.org/licenses/ .
  19 */
  20
  21 #include "mpc-impl.h"
  22
  23 static int
  24 mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
  25    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
  26    /* assumes that op (that is, its real or imaginary part) is not finite */
  27 {
  28    int overlap;
  29    mpc_t op_loc;
  30
  31    overlap = (rop_sin == op || rop_cos == op);
  32    if (overlap) {
  33       mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
  34       mpc_set (op_loc, op, MPC_RNDNN);
  35    }
  36    else
  37       op_loc [0] = op [0];
  38
  39    if (rop_sin != NULL) {
  40       if (mpfr_nan_p (mpc_realref (op_loc)) || mpfr_nan_p (mpc_imagref (op_loc))) {
  41          mpc_set (rop_sin, op_loc, rnd_sin);
  42          if (mpfr_nan_p (mpc_imagref (op_loc))) {
  43             /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
  44             /* sin(-0 +i*NaN) = -0 +i*NaN */
  45             /* sin(+0 +i*NaN) = +0 +i*NaN */
  46             if (!mpfr_zero_p (mpc_realref (op_loc)))
  47                mpfr_set_nan (mpc_realref (rop_sin));
  48          }
  49          else /* op = NaN + i*y */
  50             if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
  51             /* sin(NaN -i*Inf) = NaN -i*Inf */
  52             /* sin(NaN -i*0) = NaN -i*0 */
  53             /* sin(NaN +i*0) = NaN +i*0 */
  54             /* sin(NaN +i*Inf) = NaN +i*Inf */
  55             /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */
  56             mpfr_set_nan (mpc_imagref (rop_sin));
  57       }
  58       else if (mpfr_inf_p (mpc_realref (op_loc))) {
  59          mpfr_set_nan (mpc_realref (rop_sin));
  60
  61          if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
  62             /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
  63             mpfr_set_nan (mpc_imagref (rop_sin));
  64          else
  65             /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
  66             /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
  67             /* sin(+/-Inf -i*0) = NaN -i*0 */
  68             /* sin(+/-Inf +i*0) = NaN +i*0 */
  69             mpfr_set (mpc_imagref (rop_sin), mpc_imagref (op_loc), MPC_RND_IM (rnd_sin));
  70       }
  71       else if (mpfr_zero_p (mpc_realref (op_loc))) {
  72          /* sin(-0 -i*Inf) = -0 -i*Inf */
  73          /* sin(+0 -i*Inf) = +0 -i*Inf */
  74          /* sin(-0 +i*Inf) = -0 +i*Inf */
  75          /* sin(+0 +i*Inf) = +0 +i*Inf */
  76          mpc_set (rop_sin, op_loc, rnd_sin);
  77       }
  78       else {
  79          /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
  80          /* sin(x +i*Inf) = +Inf*(sin(x) +i*cos(x)) */
  81          mpfr_t s, c;
  82          mpfr_init2 (s, 2);
  83          mpfr_init2 (c, 2);
  84          mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDZ);
  85          mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s));
  86          mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc)));
  87          mpfr_clear (s);
  88          mpfr_clear (c);
  89       }
  90    }
  91
  92    if (rop_cos != NULL) {
  93       if (mpfr_nan_p (mpc_realref (op_loc))) {
  94          /* cos(NaN + i * NaN) = NaN + i * NaN */
  95          /* cos(NaN - i * Inf) = +Inf + i * NaN */
  96          /* cos(NaN + i * Inf) = +Inf + i * NaN */
  97          /* cos(NaN - i * 0) = NaN - i * 0 */
  98          /* cos(NaN + i * 0) = NaN + i * 0 */
  99          /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */
 100          if (mpfr_inf_p (mpc_imagref (op_loc)))
 101             mpfr_set_inf (mpc_realref (rop_cos), +1);
 102          else
 103             mpfr_set_nan (mpc_realref (rop_cos));
 104
 105          if (mpfr_zero_p (mpc_imagref (op_loc)))
 106             mpfr_set (mpc_imagref (rop_cos), mpc_imagref (op_loc), MPC_RND_IM (rnd_cos));
 107          else
 108             mpfr_set_nan (mpc_imagref (rop_cos));
 109       }
 110       else if (mpfr_nan_p (mpc_imagref (op_loc))) {
 111           /* cos(-Inf + i * NaN) = NaN + i * NaN */
 112           /* cos(+Inf + i * NaN) = NaN + i * NaN */
 113           /* cos(-0 + i * NaN) = NaN - i * 0 */
 114           /* cos(+0 + i * NaN) = NaN + i * 0 */
 115           /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */
 116          if (mpfr_zero_p (mpc_realref (op_loc)))
 117             mpfr_set (mpc_imagref (rop_cos), mpc_realref (op_loc), MPC_RND_IM (rnd_cos));
 118          else
 119             mpfr_set_nan (mpc_imagref (rop_cos));
 120
 121          mpfr_set_nan (mpc_realref (rop_cos));
 122       }
 123       else if (mpfr_inf_p (mpc_realref (op_loc))) {
 124          /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */
 125          /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */
 126          /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */
 127          /* cos(-Inf +i*0) = cos(+Inf -i*0) = NaN +i*0 */
 128          /* cos(-Inf +i*y) = cos(+Inf +i*y) = NaN +i*NaN, when y != 0 */
 129
 130          const int same_sign =
 131             mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
 132
 133          if (mpfr_inf_p (mpc_imagref (op_loc)))
 134             mpfr_set_inf (mpc_realref (rop_cos), (same_sign ? -1 : +1));
 135          else
 136             mpfr_set_nan (mpc_realref (rop_cos));
 137
 138          if (mpfr_zero_p (mpc_imagref (op_loc)))
 139             mpfr_setsign (mpc_imagref (rop_cos), mpc_imagref (op_loc), same_sign,
 140                           MPC_RND_IM(rnd_cos));
 141          else
 142             mpfr_set_nan (mpc_imagref (rop_cos));
 143       }
 144       else if (mpfr_zero_p (mpc_realref (op_loc))) {
 145          /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */
 146          /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */
 147          const int same_sign =
 148             mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
 149
 150          mpfr_setsign (mpc_imagref (rop_cos), mpc_realref (op_loc), same_sign,
 151                        MPC_RND_IM (rnd_cos));
 152          mpfr_set_inf (mpc_realref (rop_cos), +1);
 153       }
 154       else {
 155          /* cos(x -i*Inf) = +Inf*cos(x) +i*Inf*sin(x), when x != 0 */
 156          /* cos(x +i*Inf) = +Inf*cos(x) -i*Inf*sin(x), when x != 0 */
 157          mpfr_t s, c;
 158          mpfr_init2 (c, 2);
 159          mpfr_init2 (s, 2);
 160          mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDN);
 161          mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c));
 162          mpfr_set_inf (mpc_imagref (rop_cos),
 163             (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1));
 164          mpfr_clear (s);
 165          mpfr_clear (c);
 166       }
 167    }
 168
 169    if (overlap)
 170       mpc_clear (op_loc);
 171
 172    return MPC_INEX12 (MPC_INEX (0,0), MPC_INEX (0,0));
 173       /* everything is exact */
 174 }
 175
 176
 177 static int
 178 mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
 179    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
 180    /* assumes that op is real */
 181 {
 182    int inex_sin_re = 0, inex_cos_re = 0;
 183       /* Until further notice, assume computations exact; in particular,
 184          by definition, for not computed values.                         */
 185    mpfr_t s, c;
 186    int inex_s, inex_c;
 187    int sign_im = mpfr_signbit (mpc_imagref (op));
 188
 189    /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
 190    /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
 191    if (rop_sin != 0)
 192       mpfr_init2 (s, MPC_PREC_RE (rop_sin));
 193    else
 194       mpfr_init2 (s, 2); /* We need only the sign. */
 195    if (rop_cos != NULL)
 196       mpfr_init2 (c, MPC_PREC_RE (rop_cos));
 197    else
 198       mpfr_init2 (c, 2);
 199    inex_s = mpfr_sin (s, mpc_realref (op), MPC_RND_RE (rnd_sin));
 200    inex_c = mpfr_cos (c, mpc_realref (op), MPC_RND_RE (rnd_cos));
 201       /* We cannot use mpfr_sin_cos since we may need two distinct rounding
 202          modes and the exact return values. If we need only the sign, an
 203          arbitrary rounding mode will work.                                 */
 204
 205    if (rop_sin != NULL) {
 206       mpfr_set (mpc_realref (rop_sin), s, GMP_RNDN); /* exact */
 207       inex_sin_re = inex_s;
 208       mpfr_set_zero (mpc_imagref (rop_sin),
 209          (     ( sign_im && !mpfr_signbit(c))
 210             || (!sign_im &&  mpfr_signbit(c)) ? -1 : 1));
 211    }
 212
 213    if (rop_cos != NULL) {
 214       mpfr_set (mpc_realref (rop_cos), c, GMP_RNDN); /* exact */
 215       inex_cos_re = inex_c;
 216       mpfr_set_zero (mpc_imagref (rop_cos),
 217          (     ( sign_im &&  mpfr_signbit(s))
 218             || (!sign_im && !mpfr_signbit(s)) ? -1 : 1));
 219    }
 220
 221    mpfr_clear (s);
 222    mpfr_clear (c);
 223
 224    return MPC_INEX12 (MPC_INEX (inex_sin_re, 0), MPC_INEX (inex_cos_re, 0));
 225 }
 226
 227
 228 static int
 229 mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
 230    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
 231    /* assumes that op is purely imaginary, but not zero */
 232 {
 233    int inex_sin_im = 0, inex_cos_re = 0;
 234       /* assume exact if not computed */
 235    int overlap;
 236    mpc_t op_loc;
 237
 238    overlap = (rop_sin == op || rop_cos == op);
 239    if (overlap) {
 240       mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
 241       mpc_set (op_loc, op, MPC_RNDNN);
 242    }
 243    else
 244       op_loc [0] = op [0];
 245
 246    if (rop_sin != NULL) {
 247       /* sin(+-O +i*y) = +-0 +i*sinh(y) */
 248       mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), GMP_RNDN);
 249       inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin));
 250    }
 251
 252    if (rop_cos != NULL) {
 253       /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0,
 254          cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0,
 255          where y > 0 */
 256       inex_cos_re = mpfr_cosh (mpc_realref (rop_cos), mpc_imagref (op_loc), MPC_RND_RE (rnd_cos));
 257
 258       mpfr_set_ui (mpc_imagref (rop_cos), 0ul, MPC_RND_IM (rnd_cos));
 259       if (mpfr_signbit (mpc_realref (op_loc)) ==  mpfr_signbit (mpc_imagref (op_loc)))
 260          MPFR_CHANGE_SIGN (mpc_imagref (rop_cos));
 261    }
 262
 263    if (overlap)
 264       mpc_clear (op_loc);
 265
 266    return MPC_INEX12 (MPC_INEX (0, inex_sin_im), MPC_INEX (inex_cos_re, 0));
 267 }
 268
 269
 270 int
 271 mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
 272    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
 273    /* Feature not documented in the texinfo file: One of rop_sin or
 274       rop_cos may be NULL, in which case it is not computed, and the
 275       corresponding ternary inexact value is set to 0 (exact).       */
 276 {
 277    if (!mpc_fin_p (op))
 278       return mpc_sin_cos_nonfinite (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
 279    else if (mpfr_zero_p (mpc_imagref (op)))
 280       return mpc_sin_cos_real (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
 281    else if (mpfr_zero_p (mpc_realref (op)))
 282       return mpc_sin_cos_imag (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
 283    else {
 284       /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b)
 285                            and  cos(op) = cos(a)*cosh(b) - i*sin(a)*sinh(b).
 286
 287          For Re(sin(op)) (and analogously, the other parts), we use the
 288          following algorithm, with rounding to nearest for all operations
 289          and working precision w:
 290
 291          (1) x = o(sin(a))
 292          (2) y = o(cosh(b))
 293          (3) r = o(x*y)
 294          then the error on r is at most 4 ulps, since we can write
 295          r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
 296          thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
 297          thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
 298       */
 299       mpfr_t s, c, sh, ch, sch, csh;
 300       mpfr_prec_t prec;
 301       int ok;
 302       int inex_re, inex_im, inex_sin, inex_cos;
 303
 304       prec = 2;
 305       if (rop_sin != NULL)
 306          prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
 307       if (rop_cos != NULL)
 308          prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos));
 309
 310       mpfr_init2 (s, 2);
 311       mpfr_init2 (c, 2);
 312       mpfr_init2 (sh, 2);
 313       mpfr_init2 (ch, 2);
 314       mpfr_init2 (sch, 2);
 315       mpfr_init2 (csh, 2);
 316
 317       do {
 318          ok = 1;
 319          prec += mpc_ceil_log2 (prec) + 5;
 320
 321          mpfr_set_prec (s, prec);
 322          mpfr_set_prec (c, prec);
 323          mpfr_set_prec (sh, prec);
 324          mpfr_set_prec (ch, prec);
 325          mpfr_set_prec (sch, prec);
 326          mpfr_set_prec (csh, prec);
 327
 328          mpfr_sin_cos (s, c, mpc_realref(op), GMP_RNDN);
 329          mpfr_sinh_cosh (sh, ch, mpc_imagref(op), GMP_RNDN);
 330
 331          if (rop_sin != NULL) {
 332             /* real part of sine */
 333             mpfr_mul (sch, s, ch, GMP_RNDN);
 334             ok = (!mpfr_number_p (sch))
 335                   || mpfr_can_round (sch, prec - 2, GMP_RNDN, GMP_RNDZ,
 336                         MPC_PREC_RE (rop_sin)
 337                         + (MPC_RND_RE (rnd_sin) == GMP_RNDN));
 338
 339             if (ok) {
 340                /* imaginary part of sine */
 341                mpfr_mul (csh, c, sh, GMP_RNDN);
 342                ok = (!mpfr_number_p (csh))
 343                      || mpfr_can_round (csh, prec - 2, GMP_RNDN, GMP_RNDZ,
 344                            MPC_PREC_IM (rop_sin)
 345                            + (MPC_RND_IM (rnd_sin) == GMP_RNDN));
 346             }
 347          }
 348
 349          if (rop_cos != NULL && ok) {
 350             /* real part of cosine */
 351             mpfr_mul (c, c, ch, GMP_RNDN);
 352             ok = (!mpfr_number_p (c))
 353                   || mpfr_can_round (c, prec - 2, GMP_RNDN, GMP_RNDZ,
 354                         MPC_PREC_RE (rop_cos)
 355                         + (MPC_RND_RE (rnd_cos) == GMP_RNDN));
 356
 357             if (ok) {
 358                /* imaginary part of cosine */
 359                mpfr_mul (s, s, sh, GMP_RNDN);
 360                mpfr_neg (s, s, GMP_RNDN);
 361                ok = (!mpfr_number_p (s))
 362                      || mpfr_can_round (s, prec - 2, GMP_RNDN, GMP_RNDZ,
 363                            MPC_PREC_IM (rop_cos)
 364                            + (MPC_RND_IM (rnd_cos) == GMP_RNDN));
 365             }
 366          }
 367       } while (ok == 0);
 368
 369       if (rop_sin != NULL) {
 370          inex_re = mpfr_set (mpc_realref (rop_sin), sch, MPC_RND_RE (rnd_sin));
 371          if (mpfr_inf_p (sch))
 372             inex_re = mpfr_sgn (sch);
 373          inex_im = mpfr_set (mpc_imagref (rop_sin), csh, MPC_RND_IM (rnd_sin));
 374          if (mpfr_inf_p (csh))
 375             inex_im = mpfr_sgn (csh);
 376          inex_sin = MPC_INEX (inex_re, inex_im);
 377       }
 378       else
 379          inex_sin = MPC_INEX (0,0); /* return exact if not computed */
 380
 381       if (rop_cos != NULL) {
 382          inex_re = mpfr_set (mpc_realref (rop_cos), c, MPC_RND_RE (rnd_cos));
 383          if (mpfr_inf_p (c))
 384             inex_re = mpfr_sgn (c);
 385          inex_im = mpfr_set (mpc_imagref (rop_cos), s, MPC_RND_IM (rnd_cos));
 386          if (mpfr_inf_p (s))
 387             inex_im = mpfr_sgn (s);
 388          inex_cos = MPC_INEX (inex_re, inex_im);
 389       }
 390       else
 391          inex_cos = MPC_INEX (0,0); /* return exact if not computed */
 392
 393       mpfr_clear (s);
 394       mpfr_clear (c);
 395       mpfr_clear (sh);
 396       mpfr_clear (ch);
 397       mpfr_clear (sch);
 398       mpfr_clear (csh);
 399
 400       return (MPC_INEX12 (inex_sin, inex_cos));
 401    }
 402 }