1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
/* 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)
{
mpc_t z1, z2;
mpfr_t n1, n2;
mpfr_prec_t prec;
int inex1, inex2, ret;
/* Handle numbers containing one NaN as mpfr_cmp. */
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)))
{
mpfr_t nan;
mpfr_init (nan);
mpfr_set_nan (nan);
ret = mpfr_cmp (nan, nan);
mpfr_clear (nan);
return ret;
}
/* 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);
do {
mpfr_set_prec (n1, prec);
mpfr_set_prec (n2, prec);
inex1 = mpc_norm (n1, z1, MPFR_RNDD);
inex2 = mpc_norm (n2, z2, MPFR_RNDD);
ret = mpfr_cmp (n1, n2);
if (ret != 0)
goto end;
else
if (inex1 == 0) /* n1 = norm(z1) */
if (inex2) /* n2 < norm(z2) */
{
ret = -1;
goto end;
}
else /* n2 = norm(z2) */
{
ret = 0;
goto end;
}
else /* n1 < norm(z1) */
if (inex2 == 0)
{
ret = 1;
goto end;
}
prec *= 2;
} while (1);
end:
mpfr_clear (n1);
mpfr_clear (n2);
return ret;
}
|
