packages/libmpc
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

backport some bugfix patches

Browse Source
This commit is contained in:
volcanodragon 2023-05-05 10:12:47 +08:00
parent 7f2b1326e0
commit 5ef54d5a53
4 changed files with 375 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

262
backport-0002-Fix-bug-in-mpc_pow_fr.patch Normal file
View File

@ -0,0 +1,262 @@
From f04c4ccf618dfb01fe6d878f602006a643f13071 Mon Sep 17 00:00:00 2001
From: Andreas Enge <andreas.enge@inria.fr>
Date: Mon, 28 Nov 2022 12:48:56 +0100
Subject: [PATCH] Fix bug in mpc_pow_fr.
* tests/pow_fr.dat: Correct test and add more tests.
* src/pow.c (mpc_pow): Correct sign of zero part in result for c*(1+I)
or c*(1-I) raised to an even positive integer.
* doc/algorithms.tex: Add comment concerning the sign.
* NEWS: Add entry.
---
NEWS | 3 +
doc/algorithms.tex | 9 +++
src/pow.c | 146 +++++++++++++++++++++++++--------------------
tests/pow_fr.dat | 7 ++-
4 files changed, 99 insertions(+), 66 deletions(-)
diff --git a/NEWS b/NEWS
index ba37c9d..6926338 100644
--- a/NEWS
+++ b/NEWS
@@ -6,6 +6,9 @@ Changes in version 1.3.0:
- New experimental function: mpc_eta_fund
- Bug fixes:
- mpc_asin for asin(z) with small |Re(z)| and tiny |Im(z)|
+ - mpc_pow_fr: sign of zero part of result when the base has up to sign
+ the same real and imaginary part, and the exponent is an even
+ positive integer
- mpc_fma
Changes in version 1.2.1, released in October 2020:
diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index aaff643..4cdbb09 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -2126,6 +2126,15 @@ the determined cases:
where $\sigma_2$ (resp $\rho_1$, $\rho_2$) is the sign of $x_2$ (resp. $y_1$,
$y_2$) and with the convention $0^0=+1$.
+FIXME: This misses the cases $(x_1 \pm x_1 i)^{y_1 + 0 i}$ for $y_1 \geq 2$
+even and $x_1 \neq 0$; the case (d)(ii) above applies when $x_2 > 0$, but
+misses $x_2 < 0$. In any case, the sign of the imaginary part
+(for $y_1$ divisible by~$4$)
+or the real part
+(for $y_1$ divisible by~$2$, but not by~$4$)
+cannot be determined.
+
+
\subsection {\texttt {mpc\_pow\_ui}}
\label {ssec:mpcpowui}
diff --git a/src/pow.c b/src/pow.c
index 0fc6932..2bab8b8 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -1,6 +1,6 @@
/* mpc_pow -- Raise a complex number to the power of another complex number.
-Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020 INRIA
+Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020, 2022 INRIA
This file is part of GNU MPC.
@@ -481,7 +481,8 @@ is_odd (mpfr_srcptr y, mpfr_exp_t k)
int
mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
{
- int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
+ int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0,
+ ramified = 0;
mpc_t t, u;
mpfr_prec_t p, pr, pi, maxprec;
int saved_underflow, saved_overflow;
@@ -640,6 +641,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
{
+ ramified = 1;
if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
z_imag = 1;
else
@@ -764,78 +766,94 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
if (z_real)
{
- /* When the result is real (see algorithm.tex for details),
+ /* When the result is real (see algorithm.tex for details) and
+ x=x1 * (1 \pm i), y a positive integer divisible by 4, then
+ Im(x^y) = 0i with a sign that cannot be determined (and is thus
+ chosen as _1). Otherwise,
Im(x^y) =
+ sign(imag(y))*0i, if |x| > 1
+ sign(imag(x))*sign(real(y))*0i, if |x| = 1
- sign(imag(y))*0i, if |x| < 1
*/
- mpfr_t n;
- int inex, cx1;
- int sign_zi, sign_rex, sign_imx;
- /* cx1 < 0 if |x| < 1
- cx1 = 0 if |x| = 1
- cx1 > 0 if |x| > 1
- */
-
- sign_rex = mpfr_signbit (mpc_realref (x));
- sign_imx = mpfr_signbit (mpc_imagref (x));
- mpfr_init (n);
- inex = mpc_norm (n, x, MPFR_RNDN);
- cx1 = mpfr_cmp_ui (n, 1);
- if (cx1 == 0 && inex != 0)
- cx1 = -inex;
-
- sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
- || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
- || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
-
- /* copy RE(y) to n since if z==y we will destroy Re(y) below */
- mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
- mpfr_set (n, mpc_realref (y), MPFR_RNDN);
- ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
- if (y_real && (x_real || x_imag))
- {
- /* FIXME: with y_real we assume Im(y) is really 0, which is the case
- for example when y comes from pow_fr, but in case Im(y) is +0 or
- -0, we might get different results */
- mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
- fix_sign (z, sign_rex, sign_imx, n);
- ret = MPC_INEX(ret, 0); /* imaginary part is exact */
- }
- else
- {
- inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
- ret = MPC_INEX (ret, inex_im);
- /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
- if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
- mpc_conj (z, z, MPC_RNDNN);
- }
+ if (ramified)
+ ret = MPC_INEX (
+ mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)),
+ mpfr_set_ui (mpc_imagref (z), 0, MPFR_RNDN));
+ else {
+ mpfr_t n;
+ int inex, cx1;
+ int sign_zi, sign_rex, sign_imx;
+ /* cx1 < 0 if |x| < 1
+ cx1 = 0 if |x| = 1
+ cx1 > 0 if |x| > 1
+ */
+
+ sign_rex = mpfr_signbit (mpc_realref (x));
+ sign_imx = mpfr_signbit (mpc_imagref (x));
+ mpfr_init (n);
+ inex = mpc_norm (n, x, MPFR_RNDN);
+ cx1 = mpfr_cmp_ui (n, 1);
+ if (cx1 == 0 && inex != 0)
+ cx1 = -inex;
+
+ sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
+ || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
+ || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
+
+ /* copy RE(y) to n since if z==y we will destroy Re(y) below */
+ mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
+ mpfr_set (n, mpc_realref (y), MPFR_RNDN);
+ ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
+ if (y_real && (x_real || x_imag))
+ {
+ /* FIXME: with y_real we assume Im(y) is really 0, which is the case
+ for example when y comes from pow_fr, but in case Im(y) is +0 or
+ -0, we might get different results */
+ mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
+ fix_sign (z, sign_rex, sign_imx, n);
+ ret = MPC_INEX(ret, 0); /* imaginary part is exact */
+ }
+ else
+ {
+ inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
+ ret = MPC_INEX (ret, inex_im);
+ /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
+ if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
+ mpc_conj (z, z, MPC_RNDNN);
+ }
- mpfr_clear (n);
+ mpfr_clear (n);
+ }
}
else if (z_imag)
{
- ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
- /* if z is imaginary and y real, then x cannot be real */
- if (y_real && x_imag)
- {
- int sign_rex = mpfr_signbit (mpc_realref (x));
-
- /* If z overlaps with y we set Re(z) before checking Re(y) below,
- but in that case y=0, which was dealt with above. */
- mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
- /* Note: fix_sign only does something when y is an integer,
- then necessarily y = 1 or 3 (mod 4), and in that case the
- sign of Im(x) is irrelevant. */
- fix_sign (z, sign_rex, 0, mpc_realref (y));
- ret = MPC_INEX(0, ret);
- }
+ if (ramified)
+ ret = MPC_INEX (
+ mpfr_set_ui (mpc_realref (z), 0, MPFR_RNDN),
+ mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_RE(rnd)));
else
- {
- inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd));
- ret = MPC_INEX(inex_re, ret);
- }
+ {
+ ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
+ /* if z is imaginary and y real, then x cannot be real */
+ if (y_real && x_imag)
+ {
+ int sign_rex = mpfr_signbit (mpc_realref (x));
+
+ /* If z overlaps with y we set Re(z) before checking Re(y) below,
+ but in that case y=0, which was dealt with above. */
+ mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
+ /* Note: fix_sign only does something when y is an integer,
+ then necessarily y = 1 or 3 (mod 4), and in that case the
+ sign of Im(x) is irrelevant. */
+ fix_sign (z, sign_rex, 0, mpc_realref (y));
+ ret = MPC_INEX(0, ret);
+ }
+ else
+ {
+ inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd));
+ ret = MPC_INEX(inex_re, ret);
+ }
+ }
}
else
ret = mpc_set (z, u, rnd);
diff --git a/tests/pow_fr.dat b/tests/pow_fr.dat
index 66162c9..632c39a 100644
--- a/tests/pow_fr.dat
+++ b/tests/pow_fr.dat
@@ -1,6 +1,6 @@
# Data file for mpc_pow_fr.
#
-# Copyright (C) 2011 INRIA
+# Copyright (C) 2011, 2022 INRIA
#
# This file is part of GNU MPC.
#
@@ -74,4 +74,7 @@
0 - 8 -0 8 -0x89p-10 2 +0 2 0x3p-2 3 7 N N
# issue revealed by random tests (with GMP_CHECK_RANDOMIZE=1669437260)
-- 0 2 -0x1p-28 2 +0 2 0x1.8p-8 2 0x1.8p-8 2 4 N N
+- 0 2 -0x3p-29 2 +0 2 0x1.8p-8 2 0x1.8p-8 2 4 N N
+- 0 2 -0x3p-29 2 +0 2 0x1.8p-8 2 -0x1.8p-8 2 4 N N
+0 - 2 +0 2 0x1p-14 2 0x1.8p-8 2 0x1.8p-8 2 2 N N
+0 + 2 +0 2 -0x1p-14 2 0x1.8p-8 2 -0x1.8p-8 2 2 N N
--
2.33.0

40
backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch Normal file
View File

@ -0,0 +1,40 @@
From 04ac91353e83962bfa99c4097b736a0efee4fea0 Mon Sep 17 00:00:00 2001
From: Andreas Enge <andreas.enge@inria.fr>
Date: Wed, 30 Nov 2022 14:54:38 +0100
Subject: [PATCH] Fix bug in bug fix of mpc_pow_fr.
This is a follow-up to commit f04c4ccf618dfb01fe6d878f602006a643f13071
* src/pow.c (mpc_pow): Correct a typo.
* tests/pow_fr.dat: Add a test.
---
src/pow.c | 2 +-
tests/pow_fr.dat | 2 ++
2 files changed, 3 insertions(+), 1 deletion(-)
diff --git a/src/pow.c b/src/pow.c
index 2bab8b8..f436942 100644
--- a/src/pow.c
+++ b/src/pow.c
@@ -830,7 +830,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
if (ramified)
ret = MPC_INEX (
mpfr_set_ui (mpc_realref (z), 0, MPFR_RNDN),
- mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_RE(rnd)));
+ mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)));
else
{
ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
diff --git a/tests/pow_fr.dat b/tests/pow_fr.dat
index 632c39a..14cd82c 100644
--- a/tests/pow_fr.dat
+++ b/tests/pow_fr.dat
@@ -78,3 +78,5 @@
- 0 2 -0x3p-29 2 +0 2 0x1.8p-8 2 -0x1.8p-8 2 4 N N
0 - 2 +0 2 0x1p-14 2 0x1.8p-8 2 0x1.8p-8 2 2 N N
0 + 2 +0 2 -0x1p-14 2 0x1.8p-8 2 -0x1.8p-8 2 2 N N
+# issue revealed by random tests (with GMP_CHECK_RANDOMIZE=1670627686)
+0 + 2 +0 2 -0x3p-18 2 -0x3p6 2 -0x3p6 2 -2 N Z
--
2.33.0

66
backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch Normal file
View File

@ -0,0 +1,66 @@
From 153108f96fad542bf772b02abb95c748743b27da Mon Sep 17 00:00:00 2001
From: Andreas Enge <andreas.enge@inria.fr>
Date: Wed, 30 Nov 2022 18:01:50 +0100
Subject: [PATCH] Correct signs of 0 parts in exact Karatsuba multiplication.
* src/mul.c (mpc_mul_karatsuba): Clear sign if a part of a result is an
exact 0.
* tests/mul.dat: Add test.
---
src/mul.c | 15 ++++++++++++++-
tests/mul.dat | 3 +++
2 files changed, 17 insertions(+), 1 deletion(-)
diff --git a/src/mul.c b/src/mul.c
index 9ca95e1..e1574a1 100644
--- a/src/mul.c
+++ b/src/mul.c
@@ -1,6 +1,6 @@
/* mpc_mul -- Multiply two complex numbers
-Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2016, 2020 INRIA
+Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2016, 2020, 2022 INRIA
This file is part of GNU MPC.
@@ -219,6 +219,7 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
imaginary part is used). If this fails, we have to start again and
need the correct values of op1 and op2.
So we just create a new variable for the result in this case. */
+ mpfr_ptr ref;
int loop;
const int MAX_MUL_LOOP = 1;
@@ -424,6 +425,18 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd)
}
}
+ /* Clear potential signs of 0. */
+ if (!inex_re) {
+ ref = mpc_realref (result);
+ if (mpfr_zero_p (ref) && mpfr_signbit (ref))
+ MPFR_CHANGE_SIGN (ref);
+ }
+ if (!inex_im) {
+ ref = mpc_imagref (result);
+ if (mpfr_zero_p (ref) && mpfr_signbit (ref))
+ MPFR_CHANGE_SIGN (ref);
+ }
+
mpc_set (rop, result, MPC_RNDNN);
}
diff --git a/tests/mul.dat b/tests/mul.dat
index ad2515a..b560098 100644
--- a/tests/mul.dat
+++ b/tests/mul.dat
@@ -184,3 +184,6 @@
+ 0 2 1 2 0x2p-536870913 2 1 2 0x1p-536870913 2 1 2 0x1p-536870913 N N
0 - 2 0 2 1 2 0x1p-536870913 2 1 2 1 2 0x1p-536870913 N N
+
+# error in sign of 0 with exact Karatsuba found on 2022-11-30
+0 0 1473 -50 1473 +0 20 1 20 2 20 -10 20 20 N N
--
2.33.0

8
libmpc.spec
View File

@ -1,12 +1,15 @@
Name: libmpc
Version: 1.2.0
Release: 5
Release: 6
Summary: C library for multiple precision complex arithmetic
License: LGPLv3+ and GFDL-1.3-only
URL: http://www.multiprecision.org/
Source0: https://ftp.gnu.org/gnu/mpc/mpc-%{version}.tar.gz
Patch6000: backport-0001-added-missing-include.patch
Patch6001: backport-0002-Fix-bug-in-mpc_pow_fr.patch
Patch6002: backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch
Patch6003: backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch
BuildRequires: gcc
BuildRequires: gmp-devel >= 5.0.0
@ -72,6 +75,9 @@ fi
%{_libdir}/libmpc.so
%changelog
* Thu Apr 20 2023 Hewenliang<hewenliang4@huawei.com> - 1.2.0-6
- backport some bugfix patches.
* Fri Jan 6 2023 Wenyu Liu<liuwenyu7@huawei.com> - 1.2.0-5
- Specify license version