xref: /netbsd-src/external/lgpl3/mpc/dist/src/cmp_abs.c (revision 39f28e1e142c5bfb6be935a49cb55e2287fec7ea)
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
mpc_cmp_abs(mpc_srcptr a,mpc_srcptr b)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