Skip to content
Browse files

Merge branch 'vendor/MPC'

  • Loading branch information...
2 parents 2ec8bbd + d30dc8c commit bff3a498f59c557e506d7f7840e0f0bd0fb3787f @jrmarino jrmarino committed Sep 29, 2012
Showing with 9,007 additions and 0 deletions.
  1. +6 −0 contrib/mpc/AUTHORS
  2. +165 −0 contrib/mpc/COPYING.LESSER
  3. +11 −0 contrib/mpc/README
  4. +28 −0 contrib/mpc/src/abs.c
  5. +228 −0 contrib/mpc/src/acos.c
  6. +76 −0 contrib/mpc/src/acosh.c
  7. +33 −0 contrib/mpc/src/add.c
  8. +33 −0 contrib/mpc/src/add_fr.c
  9. +32 −0 contrib/mpc/src/add_si.c
  10. +33 −0 contrib/mpc/src/add_ui.c
  11. +27 −0 contrib/mpc/src/arg.c
  12. +226 −0 contrib/mpc/src/asin.c
  13. +55 −0 contrib/mpc/src/asinh.c
  14. +367 −0 contrib/mpc/src/atan.c
  15. +52 −0 contrib/mpc/src/atanh.c
  16. +28 −0 contrib/mpc/src/clear.c
  17. +33 −0 contrib/mpc/src/cmp.c
  18. +34 −0 contrib/mpc/src/cmp_si_si.c
  19. +32 −0 contrib/mpc/src/conj.c
  20. +27 −0 contrib/mpc/src/cos.c
  21. +35 −0 contrib/mpc/src/cosh.c
  22. +449 −0 contrib/mpc/src/div.c
  23. +32 −0 contrib/mpc/src/div_2si.c
  24. +32 −0 contrib/mpc/src/div_2ui.c
  25. +39 −0 contrib/mpc/src/div_fr.c
  26. +32 −0 contrib/mpc/src/div_ui.c
  27. +202 −0 contrib/mpc/src/exp.c
  28. +191 −0 contrib/mpc/src/fma.c
  29. +39 −0 contrib/mpc/src/fr_div.c
  30. +33 −0 contrib/mpc/src/fr_sub.c
  31. +28 −0 contrib/mpc/src/get_prec.c
  32. +29 −0 contrib/mpc/src/get_prec2.c
  33. +48 −0 contrib/mpc/src/get_version.c
  34. +236 −0 contrib/mpc/src/get_x.c
  35. +27 −0 contrib/mpc/src/imag.c
  36. +28 −0 contrib/mpc/src/init2.c
  37. +28 −0 contrib/mpc/src/init3.c
  38. +239 −0 contrib/mpc/src/inp_str.c
  39. +217 −0 contrib/mpc/src/log.c
  40. +297 −0 contrib/mpc/src/log10.c
  41. +147 −0 contrib/mpc/src/logging.c
  42. +46 −0 contrib/mpc/src/mem.c
  43. +194 −0 contrib/mpc/src/mpc-impl.h
  44. +51 −0 contrib/mpc/src/mpc-log.h
  45. +269 −0 contrib/mpc/src/mpc.h
  46. +639 −0 contrib/mpc/src/mul.c
  47. +32 −0 contrib/mpc/src/mul_2si.c
  48. +32 −0 contrib/mpc/src/mul_2ui.c
  49. +43 −0 contrib/mpc/src/mul_fr.c
  50. +80 −0 contrib/mpc/src/mul_i.c
  51. +32 −0 contrib/mpc/src/mul_si.c
  52. +32 −0 contrib/mpc/src/mul_ui.c
  53. +32 −0 contrib/mpc/src/neg.c
  54. +182 −0 contrib/mpc/src/norm.c
  55. +39 −0 contrib/mpc/src/out_str.c
  56. +819 −0 contrib/mpc/src/pow.c
  57. +38 −0 contrib/mpc/src/pow_d.c
  58. +37 −0 contrib/mpc/src/pow_fr.c
  59. +38 −0 contrib/mpc/src/pow_ld.c
  60. +30 −0 contrib/mpc/src/pow_si.c
  61. +169 −0 contrib/mpc/src/pow_ui.c
  62. +47 −0 contrib/mpc/src/pow_z.c
  63. +34 −0 contrib/mpc/src/proj.c
  64. +27 −0 contrib/mpc/src/real.c
  65. +32 −0 contrib/mpc/src/set.c
  66. +28 −0 contrib/mpc/src/set_prec.c
  67. +42 −0 contrib/mpc/src/set_str.c
  68. +104 −0 contrib/mpc/src/set_x.c
  69. +78 −0 contrib/mpc/src/set_x_x.c
  70. +27 −0 contrib/mpc/src/sin.c
  71. +402 −0 contrib/mpc/src/sin_cos.c
  72. +47 −0 contrib/mpc/src/sinh.c
  73. +324 −0 contrib/mpc/src/sqr.c
  74. +364 −0 contrib/mpc/src/sqrt.c
  75. +89 −0 contrib/mpc/src/strtoc.c
  76. +32 −0 contrib/mpc/src/sub.c
  77. +34 −0 contrib/mpc/src/sub_fr.c
  78. +33 −0 contrib/mpc/src/sub_ui.c
  79. +29 −0 contrib/mpc/src/swap.c
  80. +284 −0 contrib/mpc/src/tan.c
  81. +47 −0 contrib/mpc/src/tanh.c
  82. +33 −0 contrib/mpc/src/uceil_log2.c
  83. +36 −0 contrib/mpc/src/ui_div.c
  84. +34 −0 contrib/mpc/src/ui_ui_sub.c
  85. +32 −0 contrib/mpc/src/urandom.c
View
6 contrib/mpc/AUTHORS
@@ -0,0 +1,6 @@
+Main authors:
+ Andreas Enge
+ Philippe Théveny
+ Paul Zimmermann
+
+Mickaël Gastineau has contributed the file Makefile.vc.
View
165 contrib/mpc/COPYING.LESSER
@@ -0,0 +1,165 @@
+ GNU LESSER GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+ This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+ 0. Additional Definitions.
+
+ As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+ "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+ An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+ A "Combined Work" is a work produced by combining or linking an
+Application with the Library. The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+ The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+ The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+ 1. Exception to Section 3 of the GNU GPL.
+
+ You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+ 2. Conveying Modified Versions.
+
+ If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+ a) under this License, provided that you make a good faith effort to
+ ensure that, in the event an Application does not supply the
+ function or data, the facility still operates, and performs
+ whatever part of its purpose remains meaningful, or
+
+ b) under the GNU GPL, with none of the additional permissions of
+ this License applicable to that copy.
+
+ 3. Object Code Incorporating Material from Library Header Files.
+
+ The object code form of an Application may incorporate material from
+a header file that is part of the Library. You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+ a) Give prominent notice with each copy of the object code that the
+ Library is used in it and that the Library and its use are
+ covered by this License.
+
+ b) Accompany the object code with a copy of the GNU GPL and this license
+ document.
+
+ 4. Combined Works.
+
+ You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+ a) Give prominent notice with each copy of the Combined Work that
+ the Library is used in it and that the Library and its use are
+ covered by this License.
+
+ b) Accompany the Combined Work with a copy of the GNU GPL and this license
+ document.
+
+ c) For a Combined Work that displays copyright notices during
+ execution, include the copyright notice for the Library among
+ these notices, as well as a reference directing the user to the
+ copies of the GNU GPL and this license document.
+
+ d) Do one of the following:
+
+ 0) Convey the Minimal Corresponding Source under the terms of this
+ License, and the Corresponding Application Code in a form
+ suitable for, and under terms that permit, the user to
+ recombine or relink the Application with a modified version of
+ the Linked Version to produce a modified Combined Work, in the
+ manner specified by section 6 of the GNU GPL for conveying
+ Corresponding Source.
+
+ 1) Use a suitable shared library mechanism for linking with the
+ Library. A suitable mechanism is one that (a) uses at run time
+ a copy of the Library already present on the user's computer
+ system, and (b) will operate properly with a modified version
+ of the Library that is interface-compatible with the Linked
+ Version.
+
+ e) Provide Installation Information, but only if you would otherwise
+ be required to provide such information under section 6 of the
+ GNU GPL, and only to the extent that such information is
+ necessary to install and execute a modified version of the
+ Combined Work produced by recombining or relinking the
+ Application with a modified version of the Linked Version. (If
+ you use option 4d0, the Installation Information must accompany
+ the Minimal Corresponding Source and Corresponding Application
+ Code. If you use option 4d1, you must provide the Installation
+ Information in the manner specified by section 6 of the GNU GPL
+ for conveying Corresponding Source.)
+
+ 5. Combined Libraries.
+
+ You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+ a) Accompany the combined library with a copy of the same work based
+ on the Library, uncombined with any other library facilities,
+ conveyed under the terms of this License.
+
+ b) Give prominent notice with the combined library that part of it
+ is a work based on the Library, and explaining where to find the
+ accompanying uncombined form of the same work.
+
+ 6. Revised Versions of the GNU Lesser General Public License.
+
+ The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+ If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.
View
11 contrib/mpc/README
@@ -0,0 +1,11 @@
+Copyright (C) INRIA 2003, 2005, 2008, 2009, 2011
+
+Copying and distribution of this file, with or without modification,
+are permitted in any medium without royalty provided the copyright
+notice and this notice are preserved. This file is offered as-is,
+without any warranty.
+
+
+GNU MPC is a complex floating-point library with exact rounding.
+It is based on the GNU MPFR floating-point library (http://www.mpfr.org/),
+which is itself based on the GNU MP library (http://gmplib.org/).
View
28 contrib/mpc/src/abs.c
@@ -0,0 +1,28 @@
+/* mpc_abs -- Absolute value of a complex number.
+
+Copyright (C) 2008, 2009, 2011 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"
+
+/* the rounding mode is mpfr_rnd_t here since we return an mpfr number */
+int
+mpc_abs (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
+{
+ return mpfr_hypot (a, mpc_realref(b), mpc_imagref(b), rnd);
+}
View
228 contrib/mpc/src/acos.c
@@ -0,0 +1,228 @@
+/* mpc_acos -- arccosine of a complex number.
+
+Copyright (C) 2009, 2010, 2011, 2012 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 <stdio.h> /* for MPC_ASSERT */
+#include "mpc-impl.h"
+
+int
+mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im, inex;
+ mpfr_prec_t p_re, p_im, p;
+ mpc_t z1;
+ mpfr_t pi_over_2;
+ mpfr_exp_t e1, e2;
+ mpfr_rnd_t rnd_im;
+ mpc_rnd_t rnd1;
+
+ inex_re = 0;
+ inex_im = 0;
+
+ /* special values */
+ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
+ {
+ if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1);
+ mpfr_set_nan (mpc_realref (rop));
+ }
+ else if (mpfr_zero_p (mpc_realref (op)))
+ {
+ inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+ else
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+
+ return MPC_INEX (inex_re, 0);
+ }
+
+ if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ if (mpfr_inf_p (mpc_realref (op)))
+ {
+ if (mpfr_inf_p (mpc_imagref (op)))
+ {
+ if (mpfr_sgn (mpc_realref (op)) > 0)
+ {
+ inex_re =
+ set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
+ mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN);
+ }
+ else
+ {
+
+ /* the real part of the result is 3*pi/4
+ a = o(pi) error(a) < 1 ulp(a)
+ b = o(3*a) error(b) < 2 ulp(b)
+ c = b/4 exact
+ thus 1 bit is lost */
+ mpfr_t x;
+ mpfr_prec_t prec;
+ int ok;
+ mpfr_init (x);
+ prec = mpfr_get_prec (mpc_realref (rop));
+ p = prec;
+
+ do
+ {
+ p += mpc_ceil_log2 (p);
+ mpfr_set_prec (x, p);
+ mpfr_const_pi (x, GMP_RNDD);
+ mpfr_mul_ui (x, x, 3, GMP_RNDD);
+ ok =
+ mpfr_can_round (x, p - 1, GMP_RNDD, MPC_RND_RE (rnd),
+ prec+(MPC_RND_RE (rnd) == GMP_RNDN));
+
+ } while (ok == 0);
+ inex_re =
+ mpfr_div_2ui (mpc_realref (rop), x, 2, MPC_RND_RE (rnd));
+ mpfr_clear (x);
+ }
+ }
+ else
+ {
+ if (mpfr_sgn (mpc_realref (op)) > 0)
+ mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
+ else
+ inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd));
+ }
+ }
+ else
+ inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
+
+ mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1);
+
+ return MPC_INEX (inex_re, 0);
+ }
+
+ /* pure real argument */
+ if (mpfr_zero_p (mpc_imagref (op)))
+ {
+ int s_im;
+ s_im = mpfr_signbit (mpc_imagref (op));
+
+ if (mpfr_cmp_ui (mpc_realref (op), 1) > 0)
+ {
+ if (s_im)
+ inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
+ MPC_RND_IM (rnd));
+ else
+ inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
+ INV_RND (MPC_RND_IM (rnd)));
+
+ mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
+ }
+ else if (mpfr_cmp_si (mpc_realref (op), -1) < 0)
+ {
+ mpfr_t minus_op_re;
+ minus_op_re[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (minus_op_re);
+
+ if (s_im)
+ inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re,
+ MPC_RND_IM (rnd));
+ else
+ inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re,
+ INV_RND (MPC_RND_IM (rnd)));
+ inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd));
+ }
+ else
+ {
+ inex_re = mpfr_acos (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
+ mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd));
+ }
+
+ if (!s_im)
+ mpc_conj (rop, rop, MPC_RNDNN);
+
+ return MPC_INEX (inex_re, inex_im);
+ }
+
+ /* pure imaginary argument */
+ if (mpfr_zero_p (mpc_realref (op)))
+ {
+ inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd));
+ inex_im = -mpfr_asinh (mpc_imagref (rop), mpc_imagref (op),
+ INV_RND (MPC_RND_IM (rnd)));
+ mpc_conj (rop,rop, MPC_RNDNN);
+
+ return MPC_INEX (inex_re, inex_im);
+ }
+
+ /* regular complex argument: acos(z) = Pi/2 - asin(z) */
+ p_re = mpfr_get_prec (mpc_realref(rop));
+ p_im = mpfr_get_prec (mpc_imagref(rop));
+ p = p_re;
+ mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im,
+ with rounding mode opposite to rnd_im */
+ rnd_im = MPC_RND_IM(rnd);
+ /* the imaginary part of asin(z) has the same sign as Im(z), thus if
+ Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf
+ so that -Im(asin(z)) is rounded to zero */
+ if (rnd_im == GMP_RNDZ)
+ rnd_im = mpfr_sgn (mpc_imagref(op)) > 0 ? GMP_RNDD : GMP_RNDU;
+ else
+ rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD
+ : rnd_im == GMP_RNDD ? GMP_RNDU
+ : rnd_im; /* both RNDZ and RNDA map to themselves for -asin(z) */
+ rnd1 = MPC_RND (GMP_RNDN, rnd_im);
+ mpfr_init2 (pi_over_2, p);
+ for (;;)
+ {
+ p += mpc_ceil_log2 (p) + 3;
+
+ mpfr_set_prec (mpc_realref(z1), p);
+ mpfr_set_prec (pi_over_2, p);
+
+ set_pi_over_2 (pi_over_2, +1, GMP_RNDN);
+ e1 = 1; /* Exp(pi_over_2) */
+ inex = mpc_asin (z1, op, rnd1); /* asin(z) */
+ MPC_ASSERT (mpfr_sgn (mpc_imagref(z1)) * mpfr_sgn (mpc_imagref(op)) > 0);
+ inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */
+ e2 = mpfr_get_exp (mpc_realref(z1));
+ mpfr_sub (mpc_realref(z1), pi_over_2, mpc_realref(z1), GMP_RNDN);
+ if (!mpfr_zero_p (mpc_realref(z1)))
+ {
+ /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) +
+ 2^(e2-p-1) */
+ e1 = e1 >= e2 ? e1 + 1 : e2 + 1;
+ /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */
+ e1 -= mpfr_get_exp (mpc_realref(z1));
+ /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */
+ e1 = e1 <= 0 ? 0 : e1;
+ /* the error on x is bounded by 2^e1 * ulp(x) */
+ mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); /* exact */
+ inex_im = -inex_im;
+ if (mpfr_can_round (mpc_realref(z1), p - e1, GMP_RNDN, GMP_RNDZ,
+ p_re + (MPC_RND_RE(rnd) == GMP_RNDN)))
+ break;
+ }
+ }
+ inex = mpc_set (rop, z1, rnd);
+ inex_re = MPC_INEX_RE(inex);
+ mpc_clear (z1);
+ mpfr_clear (pi_over_2);
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
76 contrib/mpc/src/acosh.c
@@ -0,0 +1,76 @@
+/* mpc_acosh -- inverse hyperbolic cosine of a complex number.
+
+Copyright (C) 2009, 2011, 2012 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"
+
+int
+mpc_acosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ /* acosh(z) =
+ NaN + i*NaN, if z=0+i*NaN
+ -i*acos(z), if sign(Im(z)) = -
+ i*acos(z), if sign(Im(z)) = +
+ http://functions.wolfram.com/ElementaryFunctions/ArcCosh/27/02/03/01/01/
+ */
+ mpc_t a;
+ mpfr_t tmp;
+ int inex;
+
+ if (mpfr_zero_p (mpc_realref (op)) && mpfr_nan_p (mpc_imagref (op)))
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ return 0;
+ }
+
+ /* Note reversal of precisions due to later multiplication by i or -i */
+ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop));
+
+ if (mpfr_signbit (mpc_imagref (op)))
+ {
+ inex = mpc_acos (a, op,
+ MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd)));
+
+ /* change a to -i*a, i.e., -y+i*x to x+i*y */
+ tmp[0] = mpc_realref (a)[0];
+ mpc_realref (a)[0] = mpc_imagref (a)[0];
+ mpc_imagref (a)[0] = tmp[0];
+ MPFR_CHANGE_SIGN (mpc_imagref (a));
+ inex = MPC_INEX (MPC_INEX_IM (inex), -MPC_INEX_RE (inex));
+ }
+ else
+ {
+ inex = mpc_acos (a, op,
+ MPC_RND (MPC_RND_IM (rnd), INV_RND(MPC_RND_RE (rnd))));
+
+ /* change a to i*a, i.e., y-i*x to x+i*y */
+ tmp[0] = mpc_realref (a)[0];
+ mpc_realref (a)[0] = mpc_imagref (a)[0];
+ mpc_imagref (a)[0] = tmp[0];
+ MPFR_CHANGE_SIGN (mpc_realref (a));
+ inex = MPC_INEX (-MPC_INEX_IM (inex), MPC_INEX_RE (inex));
+ }
+
+ mpc_set (rop, a, rnd);
+
+ mpc_clear (a);
+
+ return inex;
+}
View
33 contrib/mpc/src/add.c
@@ -0,0 +1,33 @@
+/* mpc_add -- Add two complex numbers.
+
+Copyright (C) 2002, 2009, 2011 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 0 iff both the real and imaginary parts are exact */
+int
+mpc_add (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_add (mpc_realref(a), mpc_realref(b), mpc_realref(c), MPC_RND_RE(rnd));
+ inex_im = mpfr_add (mpc_imagref(a), mpc_imagref(b), mpc_imagref(c), MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
33 contrib/mpc/src/add_fr.c
@@ -0,0 +1,33 @@
+/* mpc_add_fr -- Add a complex number and a floating-point number.
+
+Copyright (C) 2002, 2009, 2011 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 0 iff both the real and imaginary parts are exact */
+int
+mpc_add_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_add (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
32 contrib/mpc/src/add_si.c
@@ -0,0 +1,32 @@
+/* mpc_add_si -- Add a complex number and a signed long int.
+
+Copyright (C) 2011 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"
+
+int
+mpc_add_si (mpc_ptr rop, mpc_srcptr op1, long int op2, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_add_si (mpc_realref (rop), mpc_realref (op1), op2, MPC_RND_RE (rnd));
+ inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op1), MPC_RND_IM (rnd));
+
+ return MPC_INEX (inex_re, inex_im);
+}
View
33 contrib/mpc/src/add_ui.c
@@ -0,0 +1,33 @@
+/* mpc_add_ui -- Add a complex number and an unsigned long int.
+
+Copyright (C) 2002, 2009, 2011 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 0 iff both the real and imaginary parts are exact */
+int
+mpc_add_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_add_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
27 contrib/mpc/src/arg.c
@@ -0,0 +1,27 @@
+/* mpc_arg -- Get the argument of a complex number.
+
+Copyright (C) 2008, 2009, 2011 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"
+
+int
+mpc_arg (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
+{
+ return mpfr_atan2 (a, mpc_imagref (b), mpc_realref (b), rnd);
+}
View
226 contrib/mpc/src/asin.c
@@ -0,0 +1,226 @@
+/* mpc_asin -- arcsine of a complex number.
+
+Copyright (C) 2009, 2010, 2011 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"
+
+int
+mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ mpfr_prec_t p, p_re, p_im, incr_p = 0;
+ mpfr_rnd_t rnd_re, rnd_im;
+ mpc_t z1;
+ int inex;
+
+ /* special values */
+ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
+ {
+ if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? -1 : +1);
+ }
+ else if (mpfr_zero_p (mpc_realref (op)))
+ {
+ mpfr_set (mpc_realref (rop), mpc_realref (op), GMP_RNDN);
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+ else
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+
+ return 0;
+ }
+
+ if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ int inex_re;
+ if (mpfr_inf_p (mpc_realref (op)))
+ {
+ int inf_im = mpfr_inf_p (mpc_imagref (op));
+
+ inex_re = set_pi_over_2 (mpc_realref (rop),
+ (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
+ mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1));
+
+ if (inf_im)
+ mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN);
+ }
+ else
+ {
+ mpfr_set_zero (mpc_realref (rop), (mpfr_signbit (mpc_realref (op)) ? -1 : 1));
+ inex_re = 0;
+ mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1));
+ }
+
+ return MPC_INEX (inex_re, 0);
+ }
+
+ /* pure real argument */
+ if (mpfr_zero_p (mpc_imagref (op)))
+ {
+ int inex_re;
+ int inex_im;
+ int s_im;
+ s_im = mpfr_signbit (mpc_imagref (op));
+
+ if (mpfr_cmp_ui (mpc_realref (op), 1) > 0)
+ {
+ if (s_im)
+ inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
+ INV_RND (MPC_RND_IM (rnd)));
+ else
+ inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
+ MPC_RND_IM (rnd));
+ inex_re = set_pi_over_2 (mpc_realref (rop),
+ (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
+ if (s_im)
+ mpc_conj (rop, rop, MPC_RNDNN);
+ }
+ else if (mpfr_cmp_si (mpc_realref (op), -1) < 0)
+ {
+ mpfr_t minus_op_re;
+ minus_op_re[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (minus_op_re);
+
+ if (s_im)
+ inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re,
+ INV_RND (MPC_RND_IM (rnd)));
+ else
+ inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re,
+ MPC_RND_IM (rnd));
+ inex_re = set_pi_over_2 (mpc_realref (rop),
+ (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
+ if (s_im)
+ mpc_conj (rop, rop, MPC_RNDNN);
+ }
+ else
+ {
+ inex_im = mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd));
+ if (s_im)
+ mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN);
+ inex_re = mpfr_asin (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
+ }
+
+ return MPC_INEX (inex_re, inex_im);
+ }
+
+ /* pure imaginary argument */
+ if (mpfr_zero_p (mpc_realref (op)))
+ {
+ int inex_im;
+ int s;
+ s = mpfr_signbit (mpc_realref (op));
+ mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
+ if (s)
+ mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN);
+ inex_im = mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
+
+ return MPC_INEX (0, inex_im);
+ }
+
+ /* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */
+ p_re = mpfr_get_prec (mpc_realref(rop));
+ p_im = mpfr_get_prec (mpc_imagref(rop));
+ rnd_re = MPC_RND_RE(rnd);
+ rnd_im = MPC_RND_IM(rnd);
+ p = p_re >= p_im ? p_re : p_im;
+ mpc_init2 (z1, p);
+ while (1)
+ {
+ mpfr_exp_t ex, ey, err;
+
+ p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */
+ incr_p = p / 2;
+ mpfr_set_prec (mpc_realref(z1), p);
+ mpfr_set_prec (mpc_imagref(z1), p);
+
+ /* z1 <- z^2 */
+ mpc_sqr (z1, op, MPC_RNDNN);
+ /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */
+ /* z1 <- 1-z1 */
+ ex = mpfr_get_exp (mpc_realref(z1));
+ mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), GMP_RNDN);
+ mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN);
+ ex = ex - mpfr_get_exp (mpc_realref(z1));
+ ex = (ex <= 0) ? 0 : ex;
+ /* err(x) <= 2^ex * ulp(x) */
+ ex = ex + mpfr_get_exp (mpc_realref(z1)) - p;
+ /* err(x) <= 2^ex */
+ ey = mpfr_get_exp (mpc_imagref(z1)) - p - 1;
+ /* err(y) <= 2^ey */
+ ex = (ex >= ey) ? ex : ey; /* err(x), err(y) <= 2^ex, i.e., the norm
+ of the error is bounded by |h|<=2^(ex+1/2) */
+ /* z1 <- sqrt(z1): if z1 = z + h, then sqrt(z1) = sqrt(z) + h/2/sqrt(t) */
+ ey = mpfr_get_exp (mpc_realref(z1)) >= mpfr_get_exp (mpc_imagref(z1))
+ ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
+ /* we have |z1| >= 2^(ey-1) thus 1/|z1| <= 2^(1-ey) */
+ mpc_sqrt (z1, z1, MPC_RNDNN);
+ ex = (2 * ex + 1) - 2 - (ey - 1); /* |h^2/4/|t| <= 2^ex */
+ ex = (ex + 1) / 2; /* ceil(ex/2) */
+ /* express ex in terms of ulp(z1) */
+ ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1))
+ ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
+ ex = ex - ey + p;
+ /* take into account the rounding error in the mpc_sqrt call */
+ err = (ex <= 0) ? 1 : ex + 1;
+ /* err(x) <= 2^err * ulp(x), err(y) <= 2^err * ulp(y) */
+ /* z1 <- i*z + z1 */
+ ex = mpfr_get_exp (mpc_realref(z1));
+ ey = mpfr_get_exp (mpc_imagref(z1));
+ mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), GMP_RNDN);
+ mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), GMP_RNDN);
+ if (mpfr_cmp_ui (mpc_realref(z1), 0) == 0 || mpfr_cmp_ui (mpc_imagref(z1), 0) == 0)
+ continue;
+ ex -= mpfr_get_exp (mpc_realref(z1)); /* cancellation in x */
+ ey -= mpfr_get_exp (mpc_imagref(z1)); /* cancellation in y */
+ ex = (ex >= ey) ? ex : ey; /* maximum cancellation */
+ err += ex;
+ err = (err <= 0) ? 1 : err + 1; /* rounding error in sub/add */
+ /* z1 <- log(z1): if z1 = z + h, then log(z1) = log(z) + h/t with
+ |t| >= min(|z1|,|z|) */
+ ex = mpfr_get_exp (mpc_realref(z1));
+ ey = mpfr_get_exp (mpc_imagref(z1));
+ ex = (ex >= ey) ? ex : ey;
+ err += ex - p; /* revert to absolute error <= 2^err */
+ mpc_log (z1, z1, GMP_RNDN);
+ err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */
+ /* express err in terms of ulp(z1) */
+ ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1))
+ ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
+ err = err - ey + p;
+ /* take into account the rounding error in the mpc_log call */
+ err = (err <= 0) ? 1 : err + 1;
+ /* z1 <- -i*z1 */
+ mpfr_swap (mpc_realref(z1), mpc_imagref(z1));
+ mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN);
+ if (mpfr_can_round (mpc_realref(z1), p - err, GMP_RNDN, GMP_RNDZ,
+ p_re + (rnd_re == GMP_RNDN)) &&
+ mpfr_can_round (mpc_imagref(z1), p - err, GMP_RNDN, GMP_RNDZ,
+ p_im + (rnd_im == GMP_RNDN)))
+ break;
+ }
+
+ inex = mpc_set (rop, z1, rnd);
+ mpc_clear (z1);
+
+ return inex;
+}
View
55 contrib/mpc/src/asinh.c
@@ -0,0 +1,55 @@
+/* mpc_asinh -- inverse hyperbolic sine of a complex number.
+
+Copyright (C) 2009, 2011, 2012 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"
+
+int
+mpc_asinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ /* asinh(op) = -i*asin(i*op) */
+ int inex;
+ mpc_t z, a;
+ mpfr_t tmp;
+
+ /* z = i*op */
+ mpc_realref (z)[0] = mpc_imagref (op)[0];
+ mpc_imagref (z)[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (mpc_realref (z));
+
+ /* Note reversal of precisions due to later multiplication by -i */
+ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop));
+
+ inex = mpc_asin (a, z,
+ MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd)));
+
+ /* if a = asin(i*op) = x+i*y, and we want y-i*x */
+
+ /* change a to -i*a */
+ tmp[0] = mpc_realref (a)[0];
+ mpc_realref (a)[0] = mpc_imagref (a)[0];
+ mpc_imagref (a)[0] = tmp[0];
+ MPFR_CHANGE_SIGN (mpc_imagref (a));
+
+ mpc_set (rop, a, MPC_RNDNN); /* exact */
+
+ mpc_clear (a);
+
+ return MPC_INEX (MPC_INEX_IM (inex), -MPC_INEX_RE (inex));
+}
View
367 contrib/mpc/src/atan.c
@@ -0,0 +1,367 @@
+/* mpc_atan -- arctangent of a complex number.
+
+Copyright (C) 2009, 2010, 2011, 2012 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 <stdio.h>
+#include "mpc-impl.h"
+
+/* set rop to
+ -pi/2 if s < 0
+ +pi/2 else
+ rounded in the direction rnd
+*/
+int
+set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd)
+{
+ int inex;
+
+ inex = mpfr_const_pi (rop, s < 0 ? INV_RND (rnd) : rnd);
+ mpfr_div_2ui (rop, rop, 1, GMP_RNDN);
+ if (s < 0)
+ {
+ inex = -inex;
+ mpfr_neg (rop, rop, GMP_RNDN);
+ }
+
+ return inex;
+}
+
+int
+mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ int s_re;
+ int s_im;
+ int inex_re;
+ int inex_im;
+ int inex;
+
+ inex_re = 0;
+ inex_im = 0;
+ s_re = mpfr_signbit (mpc_realref (op));
+ s_im = mpfr_signbit (mpc_imagref (op));
+
+ /* special values */
+ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
+ {
+ if (mpfr_nan_p (mpc_realref (op)))
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ if (mpfr_zero_p (mpc_imagref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
+ if (s_im)
+ mpc_conj (rop, rop, MPC_RNDNN);
+ }
+ else
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+ else
+ {
+ if (mpfr_inf_p (mpc_realref (op)))
+ {
+ inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
+ mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
+ }
+ else
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ }
+ }
+ return MPC_INEX (inex_re, 0);
+ }
+
+ if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
+ {
+ inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
+
+ mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
+ if (s_im)
+ mpc_conj (rop, rop, GMP_RNDN);
+
+ return MPC_INEX (inex_re, 0);
+ }
+
+ /* pure real argument */
+ if (mpfr_zero_p (mpc_imagref (op)))
+ {
+ inex_re = mpfr_atan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
+
+ mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN);
+ if (s_im)
+ mpc_conj (rop, rop, GMP_RNDN);
+
+ return MPC_INEX (inex_re, 0);
+ }
+
+ /* pure imaginary argument */
+ if (mpfr_zero_p (mpc_realref (op)))
+ {
+ int cmp_1;
+
+ if (s_im)
+ cmp_1 = -mpfr_cmp_si (mpc_imagref (op), -1);
+ else
+ cmp_1 = mpfr_cmp_ui (mpc_imagref (op), +1);
+
+ if (cmp_1 < 0)
+ {
+ /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1
+ atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */
+
+ mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
+ if (s_re)
+ mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN);
+
+ inex_im = mpfr_atanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
+ }
+ else if (cmp_1 == 0)
+ {
+ /* atan(+/-0+i) = NaN +i*inf
+ atan(+/-0-i) = NaN -i*inf */
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_inf (mpc_imagref (rop), s_im ? -1 : +1);
+ }
+ else
+ {
+ /* atan(+0+iy) = +pi/2 +i*atanh(1/y), if |y| > 1
+ atan(-0+iy) = -pi/2 +i*atanh(1/y), if |y| > 1 */
+ mpfr_rnd_t rnd_im, rnd_away;
+ mpfr_t y;
+ mpfr_prec_t p, p_im;
+ int ok;
+
+ rnd_im = MPC_RND_IM (rnd);
+ mpfr_init (y);
+ p_im = mpfr_get_prec (mpc_imagref (rop));
+ p = p_im;
+
+ /* a = o(1/y) with error(a) < 1 ulp(a)
+ b = o(atanh(a)) with error(b) < (1+2^{1+Exp(a)-Exp(b)}) ulp(b)
+
+ As |atanh (1/y)| > |1/y| we have Exp(a)-Exp(b) <=0 so, at most,
+ 2 bits of precision are lost.
+
+ We round atanh(1/y) away from 0.
+ */
+ do
+ {
+ p += mpc_ceil_log2 (p) + 2;
+ mpfr_set_prec (y, p);
+ rnd_away = s_im == 0 ? GMP_RNDU : GMP_RNDD;
+ inex_im = mpfr_ui_div (y, 1, mpc_imagref (op), rnd_away);
+ /* FIXME: should we consider the case with unreasonably huge
+ precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be
+ representable while 1/Im(op) underflows ?
+ This corresponds to |y| = 0.5*2^emin, in which case the
+ result may be wrong. */
+
+ /* atanh cannot underflow: |atanh(x)| > |x| for |x| < 1 */
+ inex_im |= mpfr_atanh (y, y, rnd_away);
+
+ ok = inex_im == 0
+ || mpfr_can_round (y, p - 2, rnd_away, GMP_RNDZ,
+ p_im + (rnd_im == GMP_RNDN));
+ } while (ok == 0);
+
+ inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd));
+ inex_im = mpfr_set (mpc_imagref (rop), y, rnd_im);
+ mpfr_clear (y);
+ }
+ return MPC_INEX (inex_re, inex_im);
+ }
+
+ /* regular number argument */
+ {
+ mpfr_t a, b, x, y;
+ mpfr_prec_t prec, p;
+ mpfr_exp_t err, expo;
+ int ok = 0;
+ mpfr_t minus_op_re;
+ mpfr_exp_t op_re_exp, op_im_exp;
+ mpfr_rnd_t rnd1, rnd2;
+
+ mpfr_inits2 (MPFR_PREC_MIN, a, b, x, y, (mpfr_ptr) 0);
+
+ /* real part: Re(arctan(x+i*y)) = [arctan2(x,1-y) - arctan2(-x,1+y)]/2 */
+ minus_op_re[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (minus_op_re);
+ op_re_exp = mpfr_get_exp (mpc_realref (op));
+ op_im_exp = mpfr_get_exp (mpc_imagref (op));
+
+ prec = mpfr_get_prec (mpc_realref (rop)); /* result precision */
+
+ /* a = o(1-y) error(a) < 1 ulp(a)
+ b = o(atan2(x,a)) error(b) < [1+2^{3+Exp(x)-Exp(a)-Exp(b)}] ulp(b)
+ = kb ulp(b)
+ c = o(1+y) error(c) < 1 ulp(c)
+ d = o(atan2(-x,c)) error(d) < [1+2^{3+Exp(x)-Exp(c)-Exp(d)}] ulp(d)
+ = kd ulp(d)
+ e = o(b - d) error(e) < [1 + kb*2^{Exp(b}-Exp(e)}
+ + kd*2^{Exp(d)-Exp(e)}] ulp(e)
+ error(e) < [1 + 2^{4+Exp(x)-Exp(a)-Exp(e)}
+ + 2^{4+Exp(x)-Exp(c)-Exp(e)}] ulp(e)
+ because |atan(u)| < |u|
+ < [1 + 2^{5+Exp(x)-min(Exp(a),Exp(c))
+ -Exp(e)}] ulp(e)
+ f = e/2 exact
+ */
+
+ /* p: working precision */
+ p = (op_im_exp > 0 || prec > SAFE_ABS (mpfr_prec_t, op_im_exp)) ? prec
+ : (prec - op_im_exp);
+ rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? GMP_RNDD : GMP_RNDU;
+ rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? GMP_RNDU : GMP_RNDD;
+
+ do
+ {
+ p += mpc_ceil_log2 (p) + 2;
+ mpfr_set_prec (a, p);
+ mpfr_set_prec (b, p);
+ mpfr_set_prec (x, p);
+
+ /* x = upper bound for atan (x/(1-y)). Since atan is increasing, we
+ need an upper bound on x/(1-y), i.e., a lower bound on 1-y for
+ x positive, and an upper bound on 1-y for x negative */
+ mpfr_ui_sub (a, 1, mpc_imagref (op), rnd1);
+ if (mpfr_sgn (a) == 0) /* y is near 1, thus 1+y is near 2, and
+ expo will be 1 or 2 below */
+ {
+ MPC_ASSERT (mpfr_cmp_ui (mpc_imagref(op), 1) == 0);
+ /* check for intermediate underflow */
+ err = 2; /* ensures err will be expo below */
+ }
+ else
+ err = mpfr_get_exp (a); /* err = Exp(a) with the notations above */
+ mpfr_atan2 (x, mpc_realref (op), a, GMP_RNDU);
+
+ /* b = lower bound for atan (-x/(1+y)): for x negative, we need a
+ lower bound on -x/(1+y), i.e., an upper bound on 1+y */
+ mpfr_add_ui (a, mpc_imagref(op), 1, rnd2);
+ /* if a is exactly zero, i.e., Im(op) = -1, then the error on a is 0,
+ and we can simply ignore the terms involving Exp(a) in the error */
+ if (mpfr_sgn (a) == 0)
+ {
+ MPC_ASSERT (mpfr_cmp_si (mpc_imagref(op), -1) == 0);
+ /* check for intermediate underflow */
+ expo = err; /* will leave err unchanged below */
+ }
+ else
+ expo = mpfr_get_exp (a); /* expo = Exp(c) with the notations above */
+ mpfr_atan2 (b, minus_op_re, a, GMP_RNDD);
+
+ err = err < expo ? err : expo; /* err = min(Exp(a),Exp(c)) */
+ mpfr_sub (x, x, b, GMP_RNDU);
+
+ err = 5 + op_re_exp - err - mpfr_get_exp (x);
+ /* error is bounded by [1 + 2^err] ulp(e) */
+ err = err < 0 ? 1 : err + 1;
+
+ mpfr_div_2ui (x, x, 1, GMP_RNDU);
+
+ /* Note: using RND2=RNDD guarantees that if x is exactly representable
+ on prec + ... bits, mpfr_can_round will return 0 */
+ ok = mpfr_can_round (x, p - err, GMP_RNDU, GMP_RNDD,
+ prec + (MPC_RND_RE (rnd) == GMP_RNDN));
+ } while (ok == 0);
+
+ /* Imaginary part
+ Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)] */
+ prec = mpfr_get_prec (mpc_imagref (rop)); /* result precision */
+
+ /* a = o(1+y) error(a) < 1 ulp(a)
+ b = o(a^2) error(b) < 5 ulp(b)
+ c = o(x^2) error(c) < 1 ulp(c)
+ d = o(b+c) error(d) < 7 ulp(d)
+ e = o(log(d)) error(e) < [1 + 7*2^{2-Exp(e)}] ulp(e) = ke ulp(e)
+ f = o(1-y) error(f) < 1 ulp(f)
+ g = o(f^2) error(g) < 5 ulp(g)
+ h = o(c+f) error(h) < 7 ulp(h)
+ i = o(log(h)) error(i) < [1 + 7*2^{2-Exp(i)}] ulp(i) = ki ulp(i)
+ j = o(e-i) error(j) < [1 + ke*2^{Exp(e)-Exp(j)}
+ + ki*2^{Exp(i)-Exp(j)}] ulp(j)
+ error(j) < [1 + 2^{Exp(e)-Exp(j)} + 2^{Exp(i)-Exp(j)}
+ + 7*2^{3-Exp(j)}] ulp(j)
+ < [1 + 2^{max(Exp(e),Exp(i))-Exp(j)+1}
+ + 7*2^{3-Exp(j)}] ulp(j)
+ k = j/4 exact
+ */
+ err = 2;
+ p = prec; /* working precision */
+
+ do
+ {
+ p += mpc_ceil_log2 (p) + err;
+ mpfr_set_prec (a, p);
+ mpfr_set_prec (b, p);
+ mpfr_set_prec (y, p);
+
+ /* a = upper bound for log(x^2 + (1+y)^2) */
+ ROUND_AWAY (mpfr_add_ui (a, mpc_imagref (op), 1, MPFR_RNDA), a);
+ mpfr_sqr (a, a, GMP_RNDU);
+ mpfr_sqr (y, mpc_realref (op), GMP_RNDU);
+ mpfr_add (a, a, y, GMP_RNDU);
+ mpfr_log (a, a, GMP_RNDU);
+
+ /* b = lower bound for log(x^2 + (1-y)^2) */
+ mpfr_ui_sub (b, 1, mpc_imagref (op), GMP_RNDZ); /* round to zero */
+ mpfr_sqr (b, b, GMP_RNDZ);
+ /* we could write mpfr_sqr (y, mpc_realref (op), GMP_RNDZ) but it is
+ more efficient to reuse the value of y (x^2) above and subtract
+ one ulp */
+ mpfr_nextbelow (y);
+ mpfr_add (b, b, y, GMP_RNDZ);
+ mpfr_log (b, b, GMP_RNDZ);
+
+ mpfr_sub (y, a, b, GMP_RNDU);
+
+ if (mpfr_zero_p (y))
+ /* FIXME: happens when x and y have very different magnitudes;
+ could be handled more efficiently */
+ ok = 0;
+ else
+ {
+ expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b));
+ expo = expo - mpfr_get_exp (y) + 1;
+ err = 3 - mpfr_get_exp (y);
+ /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */
+ if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */
+ err = (err < 0) ? 1 : err + 2;
+ else
+ err = (expo < 0) ? 1 : expo + 2;
+
+ mpfr_div_2ui (y, y, 2, GMP_RNDN);
+ MPC_ASSERT (!mpfr_zero_p (y));
+ /* FIXME: underflow. Since the main term of the Taylor series
+ in y=0 is 1/(x^2+1) * y, this means that y is very small
+ and/or x very large; but then the mpfr_zero_p (y) above
+ should be true. This needs a proof, or better yet,
+ special code. */
+
+ ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
+ prec + (MPC_RND_IM (rnd) == GMP_RNDN));
+ }
+ } while (ok == 0);
+
+ inex = mpc_set_fr_fr (rop, x, y, rnd);
+
+ mpfr_clears (a, b, x, y, (mpfr_ptr) 0);
+ return inex;
+ }
+}
View
52 contrib/mpc/src/atanh.c
@@ -0,0 +1,52 @@
+/* mpc_atanh -- inverse hyperbolic tangent of a complex number.
+
+Copyright (C) 2009, 2011, 2012 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"
+
+int
+mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ /* atanh(op) = -i*atan(i*op) */
+ int inex;
+ mpfr_t tmp;
+ mpc_t z, a;
+
+ mpc_realref (z)[0] = mpc_imagref (op)[0];
+ mpc_imagref (z)[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (mpc_realref (z));
+
+ /* Note reversal of precisions due to later multiplication by -i */
+ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop));
+
+ inex = mpc_atan (a, z,
+ MPC_RND (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd)));
+
+ /* change a to -i*a, i.e., x+i*y to y-i*x */
+ tmp[0] = mpc_realref (a)[0];
+ mpc_realref (a)[0] = mpc_imagref (a)[0];
+ mpc_imagref (a)[0] = tmp[0];
+ MPFR_CHANGE_SIGN (mpc_imagref (a));
+
+ mpc_set (rop, a, rnd);
+
+ mpc_clear (a);
+
+ return MPC_INEX (MPC_INEX_IM (inex), -MPC_INEX_RE (inex));
+}
View
28 contrib/mpc/src/clear.c
@@ -0,0 +1,28 @@
+/* mpc_clear -- Clear a complex variable.
+
+Copyright (C) 2002, 2009, 2011 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"
+
+void
+mpc_clear (mpc_t x)
+{
+ mpfr_clear (mpc_realref(x));
+ mpfr_clear (mpc_imagref(x));
+}
View
33 contrib/mpc/src/cmp.c
@@ -0,0 +1,33 @@
+/* mpc_cmp -- Compare two complex numbers.
+
+Copyright (C) 2002, 2009, 2010, 2011 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 0 iff a = b */
+int
+mpc_cmp (mpc_srcptr a, mpc_srcptr b)
+{
+ int cmp_re, cmp_im;
+
+ cmp_re = mpfr_cmp (mpc_realref(a), mpc_realref(b));
+ cmp_im = mpfr_cmp (mpc_imagref(a), mpc_imagref(b));
+
+ return MPC_INEX(cmp_re, cmp_im);
+}
View
34 contrib/mpc/src/cmp_si_si.c
@@ -0,0 +1,34 @@
+/* mpc_cmp_si_si -- Compare a complex number to a number of the form
+ b+c*i with b and c signed integers.
+
+Copyright (C) 2005, 2009, 2011 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 0 iff a = b */
+int
+mpc_cmp_si_si (mpc_srcptr a, long int b, long int c)
+{
+ int cmp_re, cmp_im;
+
+ cmp_re = mpfr_cmp_si (mpc_realref(a), b);
+ cmp_im = mpfr_cmp_si (mpc_imagref(a), c);
+
+ return MPC_INEX(cmp_re, cmp_im);
+}
View
32 contrib/mpc/src/conj.c
@@ -0,0 +1,32 @@
+/* mpc_conj -- Conjugate of a complex number.
+
+Copyright (C) 2002, 2009, 2011 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"
+
+int
+mpc_conj (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_set (mpc_realref(a), mpc_realref(b), MPC_RND_RE(rnd));
+ inex_im = mpfr_neg (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
27 contrib/mpc/src/cos.c
@@ -0,0 +1,27 @@
+/* mpc_cos -- cosine of a complex number.
+
+Copyright (C) 2010, 2011 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"
+
+int
+mpc_cos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ return MPC_INEX2 (mpc_sin_cos (NULL, rop, op, 0, rnd));
+}
View
35 contrib/mpc/src/cosh.c
@@ -0,0 +1,35 @@
+/* mpc_cosh -- hyperbolic cosine of a complex number.
+
+Copyright (C) 2008, 2009, 2011 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"
+
+int
+mpc_cosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ /* cosh(op) = cos(i*op) */
+ mpc_t z;
+
+ /* z = i*op without copying significand */
+ mpc_realref (z)[0] = mpc_imagref (op)[0];
+ mpc_imagref (z)[0] = mpc_realref (op)[0];
+ MPFR_CHANGE_SIGN (mpc_realref (z));
+
+ return mpc_cos (rop, z, rnd);
+}
View
449 contrib/mpc/src/div.c
@@ -0,0 +1,449 @@
+/* mpc_div -- Divide two complex numbers.
+
+Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012 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"
+
+/* this routine deals with the case where w is zero */
+static int
+mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
+/* Assumes w==0, implementation according to C99 G.5.1.8 */
+{
+ int sign = MPFR_SIGNBIT (mpc_realref (w));
+ mpfr_t infty;
+
+ mpfr_init2 (infty, MPFR_PREC_MIN);
+ mpfr_set_inf (infty, sign);
+ mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd));
+ mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd));
+ mpfr_clear (infty);
+ return MPC_INEX (0, 0); /* exact */
+}
+
+/* this routine deals with the case where z is infinite and w finite */
+static int
+mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
+/* Assumes w finite and non-zero and z infinite; implementation
+ according to C99 G.5.1.8 */
+{
+ int a, b, x, y;
+
+ a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0);
+ b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0);
+
+ /* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite
+ b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */
+
+ /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */
+ /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */
+ if (a == 0 || b == 0) {
+ /* only one of a or b can be zero, since z is infinite */
+ x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w));
+ y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w));
+ }
+ else {
+ /* Both parts of z are infinite; x could be determined by sign
+ considerations and comparisons. Since operations with non-finite
+ numbers are not considered time-critical, we let mpfr do the work. */
+ mpfr_t sign;
+
+ mpfr_init2 (sign, 2);
+ /* This is enough to determine the sign of sums and differences. */
+
+ if (a == 1)
+ if (b == 1) {
+ mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ x = MPC_MPFR_SIGN (sign);
+ mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ y = MPC_MPFR_SIGN (sign);
+ }
+ else { /* b == -1 */
+ mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ x = MPC_MPFR_SIGN (sign);
+ mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ y = -MPC_MPFR_SIGN (sign);
+ }
+ else /* a == -1 */
+ if (b == 1) {
+ mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN);
+ x = MPC_MPFR_SIGN (sign);
+ mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ y = MPC_MPFR_SIGN (sign);
+ }
+ else { /* b == -1 */
+ mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ x = -MPC_MPFR_SIGN (sign);
+ mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN);
+ y = MPC_MPFR_SIGN (sign);
+ }
+ mpfr_clear (sign);
+ }
+
+ if (x == 0)
+ mpfr_set_nan (mpc_realref (rop));
+ else
+ mpfr_set_inf (mpc_realref (rop), x);
+ if (y == 0)
+ mpfr_set_nan (mpc_imagref (rop));
+ else
+ mpfr_set_inf (mpc_imagref (rop), y);
+
+ return MPC_INEX (0, 0); /* exact */
+}
+
+
+/* this routine deals with the case where z if finite and w infinite */
+static int
+mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
+/* Assumes z finite and w infinite; implementation according to
+ C99 G.5.1.8 */
+{
+ mpfr_t c, d, a, b, x, y, zero;
+
+ mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */
+ mpfr_init2 (d, 2);
+ mpfr_init2 (x, 2);
+ mpfr_init2 (y, 2);
+ mpfr_init2 (zero, 2);
+ mpfr_set_ui (zero, 0ul, GMP_RNDN);
+ mpfr_init2 (a, mpfr_get_prec (mpc_realref (z)));
+ mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z)));
+
+ mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN);
+ MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN);
+ mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN);
+ MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN);
+
+ mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */
+ mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN);
+ mpfr_add (x, a, b, GMP_RNDN);
+
+ mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN);
+ mpfr_mul (a, mpc_realref (z), d, GMP_RNDN);
+ mpfr_sub (y, b, a, GMP_RNDN);
+
+ MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN);
+ MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN);
+
+ mpfr_clear (c);
+ mpfr_clear (d);
+ mpfr_clear (x);
+ mpfr_clear (y);
+ mpfr_clear (zero);
+ mpfr_clear (a);
+ mpfr_clear (b);
+
+ return MPC_INEX (0, 0); /* exact */
+}
+
+
+static int
+mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
+/* Assumes z finite and w finite and non-zero, with imaginary part
+ of w a signed zero. */
+{
+ int inex_re, inex_im;
+ /* save signs of operands in case there are overlaps */
+ int zrs = MPFR_SIGNBIT (mpc_realref (z));
+ int zis = MPFR_SIGNBIT (mpc_imagref (z));
+ int wrs = MPFR_SIGNBIT (mpc_realref (w));
+ int wis = MPFR_SIGNBIT (mpc_imagref (w));
+
+ /* warning: rop may overlap with z,w so treat the imaginary part first */
+ inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd));
+ inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd));
+
+ /* correct signs of zeroes if necessary, which does not affect the
+ inexact flags */
+ if (mpfr_zero_p (mpc_realref (rop)))
+ mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
+ GMP_RNDN); /* exact */
+ if (mpfr_zero_p (mpc_imagref (rop)))
+ mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
+ GMP_RNDN);
+
+ return MPC_INEX(inex_re, inex_im);
+}
+
+
+static int
+mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
+/* Assumes z finite and w finite and non-zero, with real part
+ of w a signed zero. */
+{
+ int inex_re, inex_im;
+ int overlap = (rop == z) || (rop == w);
+ int imag_z = mpfr_zero_p (mpc_realref (z));
+ mpfr_t wloc;
+ mpc_t tmprop;
+ mpc_ptr dest = (overlap) ? tmprop : rop;
+ /* save signs of operands in case there are overlaps */
+ int zrs = MPFR_SIGNBIT (mpc_realref (z));
+ int zis = MPFR_SIGNBIT (mpc_imagref (z));
+ int wrs = MPFR_SIGNBIT (mpc_realref (w));
+ int wis = MPFR_SIGNBIT (mpc_imagref (w));
+
+ if (overlap)
+ mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
+
+ wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */
+ inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd));
+ mpfr_neg (wloc, wloc, GMP_RNDN);
+ /* changes the sign only in wloc, not in w; no need to correct later */
+ inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd));
+
+ if (overlap) {
+ /* Note: we could use mpc_swap here, but this might cause problems
+ if rop and tmprop have been allocated using different methods, since
+ it will swap the significands of rop and tmprop. See
+ http://lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
+ mpc_set (rop, tmprop, MPC_RNDNN); /* exact */
+ mpc_clear (tmprop);
+ }
+
+ /* correct signs of zeroes if necessary, which does not affect the
+ inexact flags */
+ if (mpfr_zero_p (mpc_realref (rop)))
+ mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
+ GMP_RNDN); /* exact */
+ if (imag_z)
+ mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
+ GMP_RNDN);
+
+ return MPC_INEX(inex_re, inex_im);
+}
+
+
+int
+mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
+{
+ int ok_re = 0, ok_im = 0;
+ mpc_t res, c_conj;
+ mpfr_t q;
+ mpfr_prec_t prec;
+ int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
+ int underflow_norm, overflow_norm, underflow_prod, overflow_prod;
+ int underflow_re = 0, overflow_re = 0, underflow_im = 0, overflow_im = 0;
+ mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd);
+ int saved_underflow, saved_overflow;
+ int tmpsgn;
+
+ /* According to the C standard G.3, there are three types of numbers: */
+ /* finite (both parts are usual real numbers; contains 0), infinite */
+ /* (at least one part is a real infinity) and all others; the latter */
+ /* are numbers containing a nan, but no infinity, and could reasonably */
+ /* be called nan. */
+ /* By G.5.1.4, infinite/finite=infinite; finite/infinite=0; */
+ /* all other divisions that are not finite/finite return nan+i*nan. */
+ /* Division by 0 could be handled by the following case of division by */
+ /* a real; we handle it separately instead. */
+ if (mpc_zero_p (c))
+ return mpc_div_zero (a, b, c, rnd);
+ else if (mpc_inf_p (b) && mpc_fin_p (c))
+ return mpc_div_inf_fin (a, b, c);
+ else if (mpc_fin_p (b) && mpc_inf_p (c))
+ return mpc_div_fin_inf (a, b, c);
+ else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
+ mpc_set_nan (a);
+ return MPC_INEX (0, 0);
+ }
+ else if (mpfr_zero_p(mpc_imagref(c)))
+ return mpc_div_real (a, b, c, rnd);
+ else if (mpfr_zero_p(mpc_realref(c)))
+ return mpc_div_imag (a, b, c, rnd);
+
+ prec = MPC_MAX_PREC(a);
+
+ mpc_init2 (res, 2);
+ mpfr_init (q);
+
+ /* create the conjugate of c in c_conj without allocating new memory */
+ mpc_realref (c_conj)[0] = mpc_realref (c)[0];
+ mpc_imagref (c_conj)[0] = mpc_imagref (c)[0];
+ MPFR_CHANGE_SIGN (mpc_imagref (c_conj));
+
+ /* save the underflow or overflow flags from MPFR */
+ saved_underflow = mpfr_underflow_p ();
+ saved_overflow = mpfr_overflow_p ();
+
+ do {
+ loops ++;
+ prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2;
+
+ mpc_set_prec (res, prec);
+ mpfr_set_prec (q, prec);
+
+ /* first compute norm(c) */
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_norm = mpc_norm (q, c, GMP_RNDU);
+ underflow_norm = mpfr_underflow_p ();
+ overflow_norm = mpfr_overflow_p ();
+ if (underflow_norm)
+ mpfr_set_ui (q, 0ul, GMP_RNDN);
+ /* to obtain divisions by 0 later on */
+
+ /* now compute b*conjugate(c) */
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ);
+ inexact_re = MPC_INEX_RE (inexact_prod);
+ inexact_im = MPC_INEX_IM (inexact_prod);
+ underflow_prod = mpfr_underflow_p ();
+ overflow_prod = mpfr_overflow_p ();
+ /* unfortunately, does not distinguish between under-/overflow
+ in real or imaginary parts
+ hopefully, the side-effects of mpc_mul do indeed raise the
+ mpfr exceptions */
+ if (overflow_prod) {
+ int isinf = 0;
+ tmpsgn = mpfr_sgn (mpc_realref(res));
+ if (tmpsgn > 0)
+ {
+ mpfr_nextabove (mpc_realref(res));
+ isinf = mpfr_inf_p (mpc_realref(res));
+ mpfr_nextbelow (mpc_realref(res));
+ }
+ else if (tmpsgn < 0)
+ {
+ mpfr_nextbelow (mpc_realref(res));
+ isinf = mpfr_inf_p (mpc_realref(res));
+ mpfr_nextabove (mpc_realref(res));
+ }
+ if (isinf)
+ {
+ mpfr_set_inf (mpc_realref(res), tmpsgn);
+ overflow_re = 1;
+ }
+ tmpsgn = mpfr_sgn (mpc_imagref(res));
+ isinf = 0;
+ if (tmpsgn > 0)
+ {
+ mpfr_nextabove (mpc_imagref(res));
+ isinf = mpfr_inf_p (mpc_imagref(res));
+ mpfr_nextbelow (mpc_imagref(res));
+ }
+ else if (tmpsgn < 0)
+ {
+ mpfr_nextbelow (mpc_imagref(res));
+ isinf = mpfr_inf_p (mpc_imagref(res));
+ mpfr_nextabove (mpc_imagref(res));
+ }
+ if (isinf)
+ {
+ mpfr_set_inf (mpc_imagref(res), tmpsgn);
+ overflow_im = 1;
+ }
+ mpc_set (a, res, rnd);
+ goto end;
+ }
+
+ /* divide the product by the norm */
+ if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) {
+ /* The division has good chances to be exact in at least one part. */
+ /* Since this can cause problems when not rounding to the nearest, */
+ /* we use the division code of mpfr, which handles the situation. */
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ);
+ underflow_re = mpfr_underflow_p ();
+ overflow_re = mpfr_overflow_p ();
+ ok_re = !inexact_re || underflow_re || overflow_re
+ || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN,
+ GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
+
+ if (ok_re) /* compute imaginary part */ {
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ);
+ underflow_im = mpfr_underflow_p ();
+ overflow_im = mpfr_overflow_p ();
+ ok_im = !inexact_im || underflow_im || overflow_im
+ || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN,
+ GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
+ }
+ }
+ else {
+ /* The division is inexact, so for efficiency reasons we invert q */
+ /* only once and multiply by the inverse. */
+ if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) {
+ /* if 1/q is inexact, the approximations of the real and
+ imaginary part below will be inexact, unless RE(res)
+ or IM(res) is zero */
+ inexact_re |= ~mpfr_zero_p (mpc_realref (res));
+ inexact_im |= ~mpfr_zero_p (mpc_imagref (res));
+ }
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ);
+ underflow_re = mpfr_underflow_p ();
+ overflow_re = mpfr_overflow_p ();
+ ok_re = !inexact_re || underflow_re || overflow_re
+ || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN,
+ GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN));
+
+ if (ok_re) /* compute imaginary part */ {
+ mpfr_clear_underflow ();
+ mpfr_clear_overflow ();
+ inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ);
+ underflow_im = mpfr_underflow_p ();
+ overflow_im = mpfr_overflow_p ();
+ ok_im = !inexact_im || underflow_im || overflow_im
+ || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN,
+ GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN));
+ }
+ }
+ } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm
+ && !underflow_prod && !overflow_prod);
+
+ inex = mpc_set (a, res, rnd);
+ inexact_re = MPC_INEX_RE (inex);
+ inexact_im = MPC_INEX_IM (inex);
+
+ end:
+ /* fix values and inexact flags in case of overflow/underflow */
+ /* FIXME: heuristic, certainly does not cover all cases */
+ if (overflow_re || (underflow_norm && !underflow_prod)) {
+ mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res)));
+ inexact_re = mpfr_sgn (mpc_realref (res));
+ }
+ else if (underflow_re || (overflow_norm && !overflow_prod)) {
+ inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1;
+ mpfr_set_zero (mpc_realref (a), -inexact_re);
+ }
+ if (overflow_im || (underflow_norm && !underflow_prod)) {
+ mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res)));
+ inexact_im = mpfr_sgn (mpc_imagref (res));
+ }
+ else if (underflow_im || (overflow_norm && !overflow_prod)) {
+ inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1;
+ mpfr_set_zero (mpc_imagref (a), -inexact_im);
+ }
+
+ mpc_clear (res);
+ mpfr_clear (q);
+
+ /* restore underflow and overflow flags from MPFR */
+ if (saved_underflow)
+ mpfr_set_underflow ();
+ if (saved_overflow)
+ mpfr_set_overflow ();
+
+ return MPC_INEX (inexact_re, inexact_im);
+}
View
32 contrib/mpc/src/div_2si.c
@@ -0,0 +1,32 @@
+/* mpc_div_2si -- Divide a complex number by 2^e.
+
+Copyright (C) 2012 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"
+
+int
+mpc_div_2si (mpc_ptr a, mpc_srcptr b, long int c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_div_2si (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_div_2si (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
32 contrib/mpc/src/div_2ui.c
@@ -0,0 +1,32 @@
+/* mpc_div_2ui -- Divide a complex number by 2^e.
+
+Copyright (C) 2002, 2009, 2011, 2012 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"
+
+int
+mpc_div_2ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_div_2ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_div_2ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
39 contrib/mpc/src/div_fr.c
@@ -0,0 +1,39 @@
+/* mpc_div_fr -- Divide a complex number by a floating-point number.
+
+Copyright (C) 2002, 2008, 2009, 2010, 2011 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"
+
+int
+mpc_div_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+ mpfr_t real;
+
+ /* We have to use temporary variable in case c=mpc_realref (a). */
+ mpfr_init2 (real, MPC_PREC_RE (a));
+
+ inex_re = mpfr_div (real, mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_div (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd));
+ mpfr_set (mpc_realref (a), real, GMP_RNDN);
+
+ mpfr_clear (real);
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
32 contrib/mpc/src/div_ui.c
@@ -0,0 +1,32 @@
+/* mpc_div_ui -- Divide a complex number by a nonnegative integer.
+
+Copyright (C) 2002, 2009, 2011 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"
+
+int
+mpc_div_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd)
+{
+ int inex_re, inex_im;
+
+ inex_re = mpfr_div_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd));
+ inex_im = mpfr_div_ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd));
+
+ return MPC_INEX(inex_re, inex_im);
+}
View
202 contrib/mpc/src/exp.c
@@ -0,0 +1,202 @@
+/* mpc_exp -- exponential of a complex number.
+
+Copyright (C) 2002, 2009, 2010, 2011 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"
+
+int
+mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
+{
+ mpfr_t x, y, z;
+ mpfr_prec_t prec;
+ int ok = 0;
+ int inex_re, inex_im;
+ int saved_underflow, saved_overflow;
+
+ /* special values */
+ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
+ /* NaNs
+ exp(nan +i*y) = nan -i*0 if y = -0,
+ nan +i*0 if y = +0,
+ nan +i*nan otherwise
+ exp(x+i*nan) = +/-0 +/-i*0 if x=-inf,
+ +/-inf +i*nan if x=+inf,
+ nan +i*nan otherwise */
+ {
+ if (mpfr_zero_p (mpc_imagref (op)))
+ return mpc_set (rop, op, MPC_RNDNN);
+
+ if (mpfr_inf_p (mpc_realref (op)))
+ {
+ if (mpfr_signbit (mpc_realref (op)))
+ return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN);
+ else
+ {
+ mpfr_set_inf (mpc_realref (rop), +1);
+ mpfr_set_nan (mpc_imagref (rop));
+ return MPC_INEX(0, 0); /* Inf/NaN are exact */
+ }
+ }
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ return MPC_INEX(0, 0); /* NaN is exact */
+ }
+
+
+ if (mpfr_zero_p (mpc_imagref(op)))
+ /* special case when the input is real
+ exp(x-i*0) = exp(x) -i*0, even if x is NaN
+ exp(x+i*0) = exp(x) +i*0, even if x is NaN */
+ {
+ inex_re = mpfr_exp (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd));
+ inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref(op), MPC_RND_IM(rnd));
+ return MPC_INEX(inex_re, inex_im);
+ }
+
+ if (mpfr_zero_p (mpc_realref (op)))
+ /* special case when the input is imaginary */
+ {
+ inex_re = mpfr_cos (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE(rnd));
+ inex_im = mpfr_sin (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM(rnd));
+ return MPC_INEX(inex_re, inex_im);
+ }
+
+
+ if (mpfr_inf_p (mpc_realref (op)))
+ /* real part is an infinity,
+ exp(-inf +i*y) = 0*(cos y +i*sin y)
+ exp(+inf +i*y) = +/-inf +i*nan if y = +/-inf
+ +inf*(cos y +i*sin y) if 0 < |y| < inf */
+ {
+ mpfr_t n;
+
+ mpfr_init2 (n, 2);
+ if (mpfr_signbit (mpc_realref (op)))
+ mpfr_set_ui (n, 0, GMP_RNDN);
+ else
+ mpfr_set_inf (n, +1);
+
+ if (mpfr_inf_p (mpc_imagref (op)))
+ {
+ inex_re = mpfr_set (mpc_realref (rop), n, GMP_RNDN);
+ if (mpfr_signbit (mpc_realref (op)))
+ inex_im = mpfr_set (mpc_imagref (rop), n, GMP_RNDN);
+ else
+ {
+ mpfr_set_nan (mpc_imagref (rop));
+ inex_im = 0; /* NaN is exact */
+ }
+ }
+ else
+ {
+ mpfr_t c, s;
+ mpfr_init2 (c, 2);
+ mpfr_init2 (s, 2);
+
+ mpfr_sin_cos (s, c, mpc_imagref (op), GMP_RNDN);
+ inex_re = mpfr_copysign (mpc_realref (rop), n, c, GMP_RNDN);
+ inex_im = mpfr_copysign (mpc_imagref (rop), n, s, GMP_RNDN);
+
+ mpfr_clear (s);
+ mpfr_clear (c);
+ }
+
+ mpfr_clear (n);
+ return MPC_INEX(inex_re, inex_im);
+ }
+
+ if (mpfr_inf_p (mpc_imagref (op)))
+ /* real part is finite non-zero number, imaginary part is an infinity */
+ {
+ mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (mpc_imagref (rop));
+ return MPC_INEX(0, 0); /* NaN is exact */
+ }
+
+
+ /* from now on, both parts of op are regular numbers */