1*39f28e1eSmrg /* mpc_cmp -- Compare two complex numbers.
2*39f28e1eSmrg
3*39f28e1eSmrg Copyright (C) 2016 INRIA
4*39f28e1eSmrg
5*39f28e1eSmrg This file is part of GNU MPC.
6*39f28e1eSmrg
7*39f28e1eSmrg GNU MPC is free software; you can redistribute it and/or modify it under
8*39f28e1eSmrg the terms of the GNU Lesser General Public License as published by the
9*39f28e1eSmrg Free Software Foundation; either version 3 of the License, or (at your
10*39f28e1eSmrg option) any later version.
11*39f28e1eSmrg
12*39f28e1eSmrg GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13*39f28e1eSmrg WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14*39f28e1eSmrg FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15*39f28e1eSmrg more details.
16*39f28e1eSmrg
17*39f28e1eSmrg You should have received a copy of the GNU Lesser General Public License
18*39f28e1eSmrg along with this program. If not, see http://www.gnu.org/licenses/ .
19*39f28e1eSmrg */
20*39f28e1eSmrg
21*39f28e1eSmrg #include "mpc-impl.h"
22*39f28e1eSmrg
23*39f28e1eSmrg /* return mpfr_cmp (mpc_abs (a), mpc_abs (b)) */
24*39f28e1eSmrg int
mpc_cmp_abs(mpc_srcptr a,mpc_srcptr b)25*39f28e1eSmrg mpc_cmp_abs (mpc_srcptr a, mpc_srcptr b)
26*39f28e1eSmrg {
27*39f28e1eSmrg mpc_t z1, z2;
28*39f28e1eSmrg mpfr_t n1, n2;
29*39f28e1eSmrg mpfr_prec_t prec;
30*39f28e1eSmrg int inex1, inex2, ret;
31*39f28e1eSmrg
32*39f28e1eSmrg /* Handle numbers containing one NaN as mpfr_cmp. */
33*39f28e1eSmrg if ( mpfr_nan_p (mpc_realref (a)) || mpfr_nan_p (mpc_imagref (a))
34*39f28e1eSmrg || mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)))
35*39f28e1eSmrg {
36*39f28e1eSmrg mpfr_t nan;
37*39f28e1eSmrg mpfr_init (nan);
38*39f28e1eSmrg mpfr_set_nan (nan);
39*39f28e1eSmrg ret = mpfr_cmp (nan, nan);
40*39f28e1eSmrg mpfr_clear (nan);
41*39f28e1eSmrg return ret;
42*39f28e1eSmrg }
43*39f28e1eSmrg
44*39f28e1eSmrg /* Handle infinities. */
45*39f28e1eSmrg if (mpc_inf_p (a))
46*39f28e1eSmrg if (mpc_inf_p (b))
47*39f28e1eSmrg return 0;
48*39f28e1eSmrg else
49*39f28e1eSmrg return 1;
50*39f28e1eSmrg else if (mpc_inf_p (b))
51*39f28e1eSmrg return -1;
52*39f28e1eSmrg
53*39f28e1eSmrg /* Replace all parts of a and b by their absolute values, then order
54*39f28e1eSmrg them by size. */
55*39f28e1eSmrg z1 [0] = a [0];
56*39f28e1eSmrg z2 [0] = b [0];
57*39f28e1eSmrg if (mpfr_signbit (mpc_realref (a)))
58*39f28e1eSmrg MPFR_CHANGE_SIGN (mpc_realref (z1));
59*39f28e1eSmrg if (mpfr_signbit (mpc_imagref (a)))
60*39f28e1eSmrg MPFR_CHANGE_SIGN (mpc_imagref (z1));
61*39f28e1eSmrg if (mpfr_signbit (mpc_realref (b)))
62*39f28e1eSmrg MPFR_CHANGE_SIGN (mpc_realref (z2));
63*39f28e1eSmrg if (mpfr_signbit (mpc_imagref (b)))
64*39f28e1eSmrg MPFR_CHANGE_SIGN (mpc_imagref (z2));
65*39f28e1eSmrg if (mpfr_cmp (mpc_realref (z1), mpc_imagref (z1)) > 0)
66*39f28e1eSmrg mpfr_swap (mpc_realref (z1), mpc_imagref (z1));
67*39f28e1eSmrg if (mpfr_cmp (mpc_realref (z2), mpc_imagref (z2)) > 0)
68*39f28e1eSmrg mpfr_swap (mpc_realref (z2), mpc_imagref (z2));
69*39f28e1eSmrg
70*39f28e1eSmrg /* Handle cases in which only one part differs. */
71*39f28e1eSmrg if (mpfr_cmp (mpc_realref (z1), mpc_realref (z2)) == 0)
72*39f28e1eSmrg return mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2));
73*39f28e1eSmrg if (mpfr_cmp (mpc_imagref (z1), mpc_imagref (z2)) == 0)
74*39f28e1eSmrg return mpfr_cmp (mpc_realref (z1), mpc_realref (z2));
75*39f28e1eSmrg
76*39f28e1eSmrg /* Implement the algorithm in algorithms.tex. */
77*39f28e1eSmrg mpfr_init (n1);
78*39f28e1eSmrg mpfr_init (n2);
79*39f28e1eSmrg prec = MPC_MAX (50, MPC_MAX (MPC_MAX_PREC (z1), MPC_MAX_PREC (z2)) / 100);
80*39f28e1eSmrg do {
81*39f28e1eSmrg mpfr_set_prec (n1, prec);
82*39f28e1eSmrg mpfr_set_prec (n2, prec);
83*39f28e1eSmrg inex1 = mpc_norm (n1, z1, MPFR_RNDD);
84*39f28e1eSmrg inex2 = mpc_norm (n2, z2, MPFR_RNDD);
85*39f28e1eSmrg ret = mpfr_cmp (n1, n2);
86*39f28e1eSmrg if (ret != 0)
87*39f28e1eSmrg goto end;
88*39f28e1eSmrg else
89*39f28e1eSmrg if (inex1 == 0) /* n1 = norm(z1) */
90*39f28e1eSmrg if (inex2) /* n2 < norm(z2) */
91*39f28e1eSmrg {
92*39f28e1eSmrg ret = -1;
93*39f28e1eSmrg goto end;
94*39f28e1eSmrg }
95*39f28e1eSmrg else /* n2 = norm(z2) */
96*39f28e1eSmrg {
97*39f28e1eSmrg ret = 0;
98*39f28e1eSmrg goto end;
99*39f28e1eSmrg }
100*39f28e1eSmrg else /* n1 < norm(z1) */
101*39f28e1eSmrg if (inex2 == 0)
102*39f28e1eSmrg {
103*39f28e1eSmrg ret = 1;
104*39f28e1eSmrg goto end;
105*39f28e1eSmrg }
106*39f28e1eSmrg prec *= 2;
107*39f28e1eSmrg } while (1);
108*39f28e1eSmrg end:
109*39f28e1eSmrg mpfr_clear (n1);
110*39f28e1eSmrg mpfr_clear (n2);
111*39f28e1eSmrg return ret;
112*39f28e1eSmrg }
113*39f28e1eSmrg
114