1/* mpc_cmp -- Compare two complex numbers.
2
3Copyright (C) 2016 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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)) */
24int
25mpc_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