Home | History | Annotate | Line # | Download | only in src 
      1 /* mpc_cmp -- Compare two complex numbers.
      2 
      3 Copyright (C) 2016 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 /* return mpfr_cmp (mpc_abs (a), mpc_abs (b)) */
     24 int
     25 mpc_cmp_abs (mpc_srcptr a, mpc_srcptr b)
     26 {
     27    mpc_t z1, z2;
     28    mpfr_t n1, n2;
     29    mpfr_prec_t prec;
     30    int inex1, inex2, ret;
     31 
     32    /* Handle numbers containing one NaN as mpfr_cmp. */
     33    if (   mpfr_nan_p (mpc_realref (a)) || mpfr_nan_p (mpc_imagref (a))
     34        || mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)))
     35      {
     36        mpfr_t nan;
     37        mpfr_init (nan);
     38        mpfr_set_nan (nan);
     39        ret = mpfr_cmp (nan, nan);
     40        mpfr_clear (nan);
     41        return ret;
     42      }
     43 
     44    /* Handle infinities. */
     45    if (mpc_inf_p (a))
     46       if (mpc_inf_p (b))
     47          return 0;
     48       else
     49          return 1;
     50    else if (mpc_inf_p (b))
     51       return -1;
     52 
     53    /* Replace all parts of a and b by their absolute values, then order
     54       them by size. */
     55    z1 [0] = a [0];
     56    z2 [0] = b [0];
     57    if (mpfr_signbit (mpc_realref (a)))
     58       MPFR_CHANGE_SIGN (mpc_realref (z1));
     59    if (mpfr_signbit (mpc_imagref (a)))
     60       MPFR_CHANGE_SIGN (mpc_imagref (z1));
     61    if (mpfr_signbit (mpc_realref (b)))
     62       MPFR_CHANGE_SIGN (mpc_realref (z2));
     63    if (mpfr_signbit (mpc_imagref (b)))
     64       MPFR_CHANGE_SIGN (mpc_imagref (z2));
     65    if (mpfr_cmp (mpc_realref (z1), mpc_imagref (z1)) > 0)
     66       mpfr_swap (mpc_realref (z1), mpc_imagref (z1));
     67    if (mpfr_cmp (mpc_realref (z2), mpc_imagref (z2)) > 0)
     68       mpfr_swap (mpc_realref (z2), mpc_imagref (z2));
     69 
     70    /* Handle cases in which only one part differs. */
     71    if (mpfr_cmp (mpc_realref (z1), mpc_realref (z2)) == 0)
     72       return mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2));
     73    if (mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2)) == 0)
     74       return mpfr_cmp (mpc_realref (z1), mpc_realref (z2));
     75 
     76    /* Implement the algorithm in algorithms.tex. */
     77    mpfr_init (n1);
     78    mpfr_init (n2);
     79    prec = MPC_MAX (50, MPC_MAX (MPC_MAX_PREC (z1), MPC_MAX_PREC (z2)) / 100);
     80    do {
     81       mpfr_set_prec (n1, prec);
     82       mpfr_set_prec (n2, prec);
     83       inex1 = mpc_norm (n1, z1, MPFR_RNDD);
     84       inex2 = mpc_norm (n2, z2, MPFR_RNDD);
     85       ret = mpfr_cmp (n1, n2);
     86       if (ret != 0)
     87         goto end;
     88       else
     89          if (inex1 == 0) /* n1 = norm(z1) */
     90             if (inex2)   /* n2 < norm(z2) */
     91               {
     92                 ret = -1;
     93                 goto end;
     94               }
     95             else /* n2 = norm(z2) */
     96               {
     97                 ret = 0;
     98                 goto end;
     99               }
    100          else /* n1 < norm(z1) */
    101             if (inex2 == 0)
    102               {
    103                 ret = 1;
    104                 goto end;
    105               }
    106       prec *= 2;
    107    } while (1);
    108  end:
    109    mpfr_clear (n1);
    110    mpfr_clear (n2);
    111    return ret;
    112 }
    113 
    114