From 05a6be124dcf57df28b39db36b4aaa1200d5f068 Mon Sep 17 00:00:00 2001 From: WizardHowl Date: Wed, 10 May 2023 07:39:18 +0000 Subject: [PATCH] libmpc: backport some bugfix patches 1.Fix bug in mpc_pow_fr 2.Fix bug in bug fix of mpc_pow_fr 3.Correct signs of 0 parts in exact Karatsuba multiplication Signed-off-by: WizardHowl --- backport-0002-Fix-bug-in-mpc_pow_fr.patch | 231 ++++++++++++++++++ ...003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch | 43 ++++ ...-0-parts-in-exact-Karatsuba-multipli.patch | 70 ++++++ libmpc.spec | 14 +- 4 files changed, 357 insertions(+), 1 deletion(-) create mode 100644 backport-0002-Fix-bug-in-mpc_pow_fr.patch create mode 100644 backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch create mode 100644 backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch diff --git a/backport-0002-Fix-bug-in-mpc_pow_fr.patch b/backport-0002-Fix-bug-in-mpc_pow_fr.patch new file mode 100644 index 0000000..d27e0b2 --- /dev/null +++ b/backport-0002-Fix-bug-in-mpc_pow_fr.patch @@ -0,0 +1,231 @@ +From 99595c0f98241bc27ddc63756e4e8ec932c1b9da Mon Sep 17 00:00:00 2001 +From: Andreas Enge +Date: Mon, 28 Nov 2022 12:48:56 +0100 +Subject: [PATCH 1/3] Fix bug in mpc_pow_fr. + +Reference:https://gitlab.inria.fr/mpc/mpc/-/commit/f04c4ccf618dfb01fe6d878f602006a643f13071 +Conflict:Context differences in tests/pow_fr.dat, NEWS and +doc/algorithms.tex(doesn't exist) + +* 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. +--- + src/pow.c | 146 ++++++++++++++++++++++++++--------------------- + tests/pow_fr.dat | 8 ++- + 2 files changed, 89 insertions(+), 65 deletions(-) + +diff --git a/src/pow.c b/src/pow.c +index 4fc90ae..185ef60 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 +@@ -765,78 +767,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 0816c13..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. + # +@@ -72,3 +72,9 @@ + # (+0 0.75)^7 = (-0 -0.13348388671875) is rounded to (-0 -0.125) + 0 + 2 -0 2 -0x1p-3 2 +0 2 0x3p-2 3 7 N N + 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 -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.36.1 + diff --git a/backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch b/backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch new file mode 100644 index 0000000..676ad63 --- /dev/null +++ b/backport-0003-Fix-bug-in-bug-fix-of-mpc_pow_fr.patch @@ -0,0 +1,43 @@ +From 957da63f50effb6e892948bb1e850d6544c06458 Mon Sep 17 00:00:00 2001 +From: Andreas Enge +Date: Wed, 30 Nov 2022 14:54:38 +0100 +Subject: [PATCH 2/3] Fix bug in bug fix of mpc_pow_fr. + +Reference:https://gitlab.inria.fr/mpc/mpc/-/commit/04ac91353e83962bfa99c4097b736a0efee4fea0 +Conflict:NA + +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 185ef60..31b8c77 100644 +--- a/src/pow.c ++++ b/src/pow.c +@@ -831,7 +831,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.36.1 + diff --git a/backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch b/backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch new file mode 100644 index 0000000..7e2671f --- /dev/null +++ b/backport-0004-Correct-signs-of-0-parts-in-exact-Karatsuba-multipli.patch @@ -0,0 +1,70 @@ +From a7928271cd02820faae17608bfb99170a174406d Mon Sep 17 00:00:00 2001 +From: Andreas Enge +Date: Wed, 30 Nov 2022 18:01:50 +0100 +Subject: [PATCH 3/3] Correct signs of 0 parts in exact Karatsuba + multiplication. + +Reference:https://gitlab.inria.fr/mpc/mpc/-/commit/153108f96fad542bf772b02abb95c748743b27da +Conflict:Context differences in src/mul.c for Copyright record + +* 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 61092e5..6a2c79f 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 INRIA ++Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2016, 2022 INRIA + + This file is part of GNU MPC. + +@@ -416,6 +416,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; + +@@ -621,6 +622,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 3f36e7c..51345fd 100644 +--- a/tests/mul.dat ++++ b/tests/mul.dat +@@ -182,3 +182,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.36.1 + diff --git a/libmpc.spec b/libmpc.spec index eeaa99e..b54ad8b 100644 --- a/libmpc.spec +++ b/libmpc.spec @@ -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 @@ -30,6 +33,9 @@ Header files and shared object symlinks for MPC is a C library. %prep %setup -q -n mpc-%{version} %patch6000 -p1 +%patch6001 -p1 +%patch6002 -p1 +%patch6003 -p1 %build %configure --disable-static @@ -72,6 +78,12 @@ fi %{_libdir}/libmpc.so %changelog +* Wed May 10 2023 Wenyu Liu - 1.2.0-6 +- backport some bugfix patches: + 1.Fix bug in mpc_pow_fr + 2.Fix bug in bug fix of mpc_pow_fr + 3.Correct signs of 0 parts in exact Karatsuba multiplication + * Fri Jan 6 2023 Wenyu Liu - 1.2.0-5 - Specify license version -- Gitee