/* mpc_cmp -- Compare two complex numbers.

Copyright (C) 2016 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://www.gnu.org/licenses/ .
*/

#include "mpc-impl.h"

/* return mpfr_cmp (mpc_abs (a), mpc_abs (b)) */
int
mpc_cmp_abs (mpc_srcptr a, mpc_srcptr b)
{
   mpfr_t nan;
   mpc_t z1, z2;
   mpfr_t n1, n2;
   mpfr_prec_t prec;
   int inex1, inex2, cmp;

   /* Handle numbers containing one NaN as mpfr_cmp. */
   mpfr_set_nan (nan); /* undocumented feature: can be used uninitialised */
   if (   mpfr_nan_p (mpc_realref (a)) || mpfr_nan_p (mpc_imagref (a))
       || mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)))
      return (mpfr_cmp (nan, nan));

   /* Handle infinities. */
   if (mpc_inf_p (a))
      if (mpc_inf_p (b))
         return (0);
      else
         return (1);
   else if (mpc_inf_p (b))
      return (-1);

   /* Replace all parts of a and b by their absolute values, then order
      them by size. */
   z1 [0] = a [0];
   z2 [0] = b [0];
   if (mpfr_signbit (mpc_realref (a)))
      MPFR_CHANGE_SIGN (mpc_realref (z1));
   if (mpfr_signbit (mpc_imagref (a)))
      MPFR_CHANGE_SIGN (mpc_imagref (z1));
   if (mpfr_signbit (mpc_realref (b)))
      MPFR_CHANGE_SIGN (mpc_realref (z2));
   if (mpfr_signbit (mpc_imagref (b)))
      MPFR_CHANGE_SIGN (mpc_imagref (z2));
   if (mpfr_cmp (mpc_realref (z1), mpc_imagref (z1)) > 0)
      mpfr_swap (mpc_realref (z1), mpc_imagref (z1));
   if (mpfr_cmp (mpc_realref (z2), mpc_imagref (z2)) > 0)
      mpfr_swap (mpc_realref (z2), mpc_imagref (z2));

   /* Handle cases in which only one part differs. */
   if (mpfr_cmp (mpc_realref (z1), mpc_realref (z2)) == 0)
      return (mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2)));
   if (mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2)) == 0)
      return (mpfr_cmp (mpc_realref (z1), mpc_realref (z2)));

   /* Implement the algorithm in algorithms.tex. */
   mpfr_init (n1);
   mpfr_init (n2);
   prec = MPC_MAX (50, MPC_MAX (MPC_MAX_PREC (z1), MPC_MAX_PREC (z2)) / 100) / 2;
   do {
      prec *= 2;
      mpfr_set_prec (n1, prec);
      mpfr_set_prec (n2, prec);
      inex1 = mpc_norm (n1, z1, MPFR_RNDD);
      inex2 = mpc_norm (n2, z2, MPFR_RNDD);
      cmp = mpfr_cmp (n1, n2);
      if (cmp != 0)
         return (cmp);
      else
         if (!inex1)
            if (inex2)
               return (-1);
            else
               return (0);
         else
            if (!inex1)
               return (1);
   } while (1);
   mpfr_clear (n1);
   mpfr_clear (n2);
}