diff options
| author | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2003-01-16 09:54:54 +0000 |
|---|---|---|
| committer | zimmerma <zimmerma@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2003-01-16 09:54:54 +0000 |
| commit | d8857148e690936e54f1fbc8d0193a8ee0d4dcc1 (patch) | |
| tree | 691cbfe393b09240bcedcb8df55ae19df1528bfd | |
| parent | a9117d598d31cf572d78f9d8a060df82647522fc (diff) | |
| download | mpc-d8857148e690936e54f1fbc8d0193a8ee0d4dcc1.tar.gz | |
Initial revision
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@2 211d60ee-9f03-0410-a15a-8952a2c7a4e4
| -rw-r--r-- | COPYING.LIB | 504 | ||||
| -rw-r--r-- | INSTALL | 56 | ||||
| -rw-r--r-- | abs.c | 62 | ||||
| -rw-r--r-- | add.c | 36 | ||||
| -rw-r--r-- | add_fr.c | 36 | ||||
| -rw-r--r-- | add_ui.c | 36 | ||||
| -rw-r--r-- | clear.c | 31 | ||||
| -rw-r--r-- | cmp.c | 36 | ||||
| -rw-r--r-- | conj.c | 35 | ||||
| -rw-r--r-- | div.c | 153 | ||||
| -rw-r--r-- | div_2exp.c | 35 | ||||
| -rw-r--r-- | div_fr.c | 35 | ||||
| -rw-r--r-- | div_ui.c | 35 | ||||
| -rw-r--r-- | exp.c | 94 | ||||
| -rw-r--r-- | init.c | 30 | ||||
| -rw-r--r-- | init2.c | 31 | ||||
| -rw-r--r-- | init3.c | 31 | ||||
| -rw-r--r-- | inp_str.c | 65 | ||||
| -rw-r--r-- | makefile | 92 | ||||
| -rw-r--r-- | mpc-impl.h | 53 | ||||
| -rw-r--r-- | mpc.h | 180 | ||||
| -rw-r--r-- | mpc.texi | 829 | ||||
| -rw-r--r-- | mul.c | 382 | ||||
| -rw-r--r-- | mul_2exp.c | 35 | ||||
| -rw-r--r-- | mul_fr.c | 35 | ||||
| -rw-r--r-- | mul_ui.c | 35 | ||||
| -rw-r--r-- | neg.c | 35 | ||||
| -rw-r--r-- | norm.c | 61 | ||||
| -rw-r--r-- | out_str.c | 38 | ||||
| -rw-r--r-- | random.c | 41 | ||||
| -rw-r--r-- | random2.c | 31 | ||||
| -rw-r--r-- | set.c | 35 | ||||
| -rw-r--r-- | set_d_d.c | 35 | ||||
| -rw-r--r-- | set_dfl_prec.c | 39 | ||||
| -rw-r--r-- | set_prec.c | 31 | ||||
| -rw-r--r-- | set_si_si.c | 35 | ||||
| -rw-r--r-- | set_ui_fr.c | 36 | ||||
| -rw-r--r-- | set_ui_ui.c | 36 | ||||
| -rw-r--r-- | sqr.c | 161 | ||||
| -rw-r--r-- | sqrt.c | 147 | ||||
| -rw-r--r-- | sub.c | 35 | ||||
| -rw-r--r-- | sub_ui.c | 36 | ||||
| -rw-r--r-- | tdiv.c | 279 | ||||
| -rw-r--r-- | test.c | 222 | ||||
| -rw-r--r-- | texp.c | 104 | ||||
| -rw-r--r-- | tmul.c | 334 | ||||
| -rw-r--r-- | tsqr.c | 220 | ||||
| -rw-r--r-- | ui_div.c | 50 |
48 files changed, 5023 insertions, 0 deletions
diff --git a/COPYING.LIB b/COPYING.LIB new file mode 100644 index 0000000..b1e3f5a --- /dev/null +++ b/COPYING.LIB @@ -0,0 +1,504 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 2.1, February 1999 + + Copyright (C) 1991, 1999 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + +[This is the first released version of the Lesser GPL. It also counts + as the successor of the GNU Library Public License, version 2, hence + the version number 2.1.] + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +Licenses are intended to guarantee your freedom to share and change +free software--to make sure the software is free for all its users. + + This license, the Lesser General Public License, applies to some +specially designated software packages--typically libraries--of the +Free Software Foundation and other authors who decide to use it. You +can use it too, but we suggest you first think carefully about whether +this license or the ordinary General Public License is the better +strategy to use in any particular case, based on the explanations below. + + When we speak of free software, we are referring to freedom of use, +not price. Our General Public Licenses are designed to make sure that +you have the freedom to distribute copies of free software (and charge +for this service if you wish); that you receive source code or can get +it if you want it; that you can change the software and use pieces of +it in new free programs; and that you are informed that you can do +these things. + + To protect your rights, we need to make restrictions that forbid +distributors to deny you these rights or to ask you to surrender these +rights. These restrictions translate to certain responsibilities for +you if you distribute copies of the library or if you modify it. + + For example, if you distribute copies of the library, whether gratis +or for a fee, you must give the recipients all the rights that we gave +you. You must make sure that they, too, receive or can get the source +code. If you link other code with the library, you must provide +complete object files to the recipients, so that they can relink them +with the library after making changes to the library and recompiling +it. And you must show them these terms so they know their rights. + + We protect your rights with a two-step method: (1) we copyright the +library, and (2) we offer you this license, which gives you legal +permission to copy, distribute and/or modify the library. + + To protect each distributor, we want to make it very clear that +there is no warranty for the free library. Also, if the library is +modified by someone else and passed on, the recipients should know +that what they have is not the original version, so that the original +author's reputation will not be affected by problems that might be +introduced by others. + + Finally, software patents pose a constant threat to the existence of +any free program. We wish to make sure that a company cannot +effectively restrict the users of a free program by obtaining a +restrictive license from a patent holder. Therefore, we insist that +any patent license obtained for a version of the library must be +consistent with the full freedom of use specified in this license. + + Most GNU software, including some libraries, is covered by the +ordinary GNU General Public License. This license, the GNU Lesser +General Public License, applies to certain designated libraries, and +is quite different from the ordinary General Public License. We use +this license for certain libraries in order to permit linking those +libraries into non-free programs. + + When a program is linked with a library, whether statically or using +a shared library, the combination of the two is legally speaking a +combined work, a derivative of the original library. The ordinary +General Public License therefore permits such linking only if the +entire combination fits its criteria of freedom. The Lesser General +Public License permits more lax criteria for linking other code with +the library. + + We call this license the "Lesser" General Public License because it +does Less to protect the user's freedom than the ordinary General +Public License. It also provides other free software developers Less +of an advantage over competing non-free programs. These disadvantages +are the reason we use the ordinary General Public License for many +libraries. However, the Lesser license provides advantages in certain +special circumstances. + + For example, on rare occasions, there may be a special need to +encourage the widest possible use of a certain library, so that it becomes +a de-facto standard. To achieve this, non-free programs must be +allowed to use the library. A more frequent case is that a free +library does the same job as widely used non-free libraries. In this +case, there is little to gain by limiting the free library to free +software only, so we use the Lesser General Public License. + + In other cases, permission to use a particular library in non-free +programs enables a greater number of people to use a large body of +free software. For example, permission to use the GNU C Library in +non-free programs enables many more people to use the whole GNU +operating system, as well as its variant, the GNU/Linux operating +system. + + Although the Lesser General Public License is Less protective of the +users' freedom, it does ensure that the user of a program that is +linked with the Library has the freedom and the wherewithal to run +that program using a modified version of the Library. + + The precise terms and conditions for copying, distribution and +modification follow. Pay close attention to the difference between a +"work based on the library" and a "work that uses the library". The +former contains code derived from the library, whereas the latter must +be combined with the library in order to run. + + GNU LESSER GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License Agreement applies to any software library or other +program which contains a notice placed by the copyright holder or +other authorized party saying it may be distributed under the terms of +this Lesser General Public License (also called "this License"). +Each licensee is addressed as "you". + + A "library" means a collection of software functions and/or data +prepared so as to be conveniently linked with application programs +(which use some of those functions and data) to form executables. + + The "Library", below, refers to any such software library or work +which has been distributed under these terms. A "work based on the +Library" means either the Library or any derivative work under +copyright law: that is to say, a work containing the Library or a +portion of it, either verbatim or with modifications and/or translated +straightforwardly into another language. (Hereinafter, translation is +included without limitation in the term "modification".) + + "Source code" for a work means the preferred form of the work for +making modifications to it. For a library, complete source code means +all the source code for all modules it contains, plus any associated +interface definition files, plus the scripts used to control compilation +and installation of the library. + + Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running a program using the Library is not restricted, and output from +such a program is covered only if its contents constitute a work based +on the Library (independent of the use of the Library in a tool for +writing it). Whether that is true depends on what the Library does +and what the program that uses the Library does. + + 1. You may copy and distribute verbatim copies of the Library's +complete source code as you receive it, in any medium, provided that +you conspicuously and appropriately publish on each copy an +appropriate copyright notice and disclaimer of warranty; keep intact +all the notices that refer to this License and to the absence of any +warranty; and distribute a copy of this License along with the +Library. + + You may charge a fee for the physical act of transferring a copy, +and you may at your option offer warranty protection in exchange for a +fee. + + 2. You may modify your copy or copies of the Library or any portion +of it, thus forming a work based on the Library, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) The modified work must itself be a software library. + + b) You must cause the files modified to carry prominent notices + stating that you changed the files and the date of any change. + + c) You must cause the whole of the work to be licensed at no + charge to all third parties under the terms of this License. + + d) If a facility in the modified Library refers to a function or a + table of data to be supplied by an application program that uses + the facility, other than as an argument passed when the facility + is invoked, then you must make a good faith effort to ensure that, + in the event an application does not supply such function or + table, the facility still operates, and performs whatever part of + its purpose remains meaningful. + + (For example, a function in a library to compute square roots has + a purpose that is entirely well-defined independent of the + application. Therefore, Subsection 2d requires that any + application-supplied function or table used by this function must + be optional: if the application does not supply it, the square + root function must still compute square roots.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Library, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Library, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote +it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Library. + +In addition, mere aggregation of another work not based on the Library +with the Library (or with a work based on the Library) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may opt to apply the terms of the ordinary GNU General Public +License instead of this License to a given copy of the Library. To do +this, you must alter all the notices that refer to this License, so +that they refer to the ordinary GNU General Public License, version 2, +instead of to this License. (If a newer version than version 2 of the +ordinary GNU General Public License has appeared, then you can specify +that version instead if you wish.) Do not make any other change in +these notices. + + Once this change is made in a given copy, it is irreversible for +that copy, so the ordinary GNU General Public License applies to all +subsequent copies and derivative works made from that copy. + + This option is useful when you wish to copy part of the code of +the Library into a program that is not a library. + + 4. You may copy and distribute the Library (or a portion or +derivative of it, under Section 2) in object code or executable form +under the terms of Sections 1 and 2 above provided that you accompany +it with the complete corresponding machine-readable source code, which +must be distributed under the terms of Sections 1 and 2 above on a +medium customarily used for software interchange. + + If distribution of object code is made by offering access to copy +from a designated place, then offering equivalent access to copy the +source code from the same place satisfies the requirement to +distribute the source code, even though third parties are not +compelled to copy the source along with the object code. + + 5. A program that contains no derivative of any portion of the +Library, but is designed to work with the Library by being compiled or +linked with it, is called a "work that uses the Library". Such a +work, in isolation, is not a derivative work of the Library, and +therefore falls outside the scope of this License. + + However, linking a "work that uses the Library" with the Library +creates an executable that is a derivative of the Library (because it +contains portions of the Library), rather than a "work that uses the +library". The executable is therefore covered by this License. +Section 6 states terms for distribution of such executables. + + When a "work that uses the Library" uses material from a header file +that is part of the Library, the object code for the work may be a +derivative work of the Library even though the source code is not. +Whether this is true is especially significant if the work can be +linked without the Library, or if the work is itself a library. The +threshold for this to be true is not precisely defined by law. + + If such an object file uses only numerical parameters, data +structure layouts and accessors, and small macros and small inline +functions (ten lines or less in length), then the use of the object +file is unrestricted, regardless of whether it is legally a derivative +work. (Executables containing this object code plus portions of the +Library will still fall under Section 6.) + + Otherwise, if the work is a derivative of the Library, you may +distribute the object code for the work under the terms of Section 6. +Any executables containing that work also fall under Section 6, +whether or not they are linked directly with the Library itself. + + 6. As an exception to the Sections above, you may also combine or +link a "work that uses the Library" with the Library to produce a +work containing portions of the Library, and distribute that work +under terms of your choice, provided that the terms permit +modification of the work for the customer's own use and reverse +engineering for debugging such modifications. + + You must give prominent notice with each copy of the work that the +Library is used in it and that the Library and its use are covered by +this License. You must supply a copy of this License. If the work +during execution displays copyright notices, you must include the +copyright notice for the Library among them, as well as a reference +directing the user to the copy of this License. Also, you must do one +of these things: + + a) Accompany the work with the complete corresponding + machine-readable source code for the Library including whatever + changes were used in the work (which must be distributed under + Sections 1 and 2 above); and, if the work is an executable linked + with the Library, with the complete machine-readable "work that + uses the Library", as object code and/or source code, so that the + user can modify the Library and then relink to produce a modified + executable containing the modified Library. (It is understood + that the user who changes the contents of definitions files in the + Library will not necessarily be able to recompile the application + to use the modified definitions.) + + b) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (1) uses at run time a + copy of the library already present on the user's computer system, + rather than copying library functions into the executable, and (2) + will operate properly with a modified version of the library, if + the user installs one, as long as the modified version is + interface-compatible with the version that the work was made with. + + c) Accompany the work with a written offer, valid for at + least three years, to give the same user the materials + specified in Subsection 6a, above, for a charge no more + than the cost of performing this distribution. + + d) If distribution of the work is made by offering access to copy + from a designated place, offer equivalent access to copy the above + specified materials from the same place. + + e) Verify that the user has already received a copy of these + materials or that you have already sent this user a copy. + + For an executable, the required form of the "work that uses the +Library" must include any data and utility programs needed for +reproducing the executable from it. However, as a special exception, +the materials to be distributed need not include anything that is +normally distributed (in either source or binary form) with the major +components (compiler, kernel, and so on) of the operating system on +which the executable runs, unless that component itself accompanies +the executable. + + It may happen that this requirement contradicts the license +restrictions of other proprietary libraries that do not normally +accompany the operating system. Such a contradiction means you cannot +use both them and the Library together in an executable that you +distribute. + + 7. 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 not covered by this License, and distribute such a combined +library, provided that the separate distribution of the work based on +the Library and of the other library facilities is otherwise +permitted, and provided that you do these two things: + + a) Accompany the combined library with a copy of the same work + based on the Library, uncombined with any other library + facilities. This must be distributed under the terms of the + Sections above. + + b) Give prominent notice with the combined library of the fact + that part of it is a work based on the Library, and explaining + where to find the accompanying uncombined form of the same work. + + 8. You may not copy, modify, sublicense, link with, or distribute +the Library except as expressly provided under this License. Any +attempt otherwise to copy, modify, sublicense, link with, or +distribute the Library is void, and will automatically terminate your +rights under this License. However, parties who have received copies, +or rights, from you under this License will not have their licenses +terminated so long as such parties remain in full compliance. + + 9. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Library or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Library (or any work based on the +Library), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Library or works based on it. + + 10. Each time you redistribute the Library (or any work based on the +Library), the recipient automatically receives a license from the +original licensor to copy, distribute, link with or modify the Library +subject to these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties with +this License. + + 11. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Library at all. For example, if a patent +license would not permit royalty-free redistribution of the Library by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Library. + +If any portion of this section is held invalid or unenforceable under any +particular circumstance, the balance of the section is intended to apply, +and the section as a whole is intended to apply in other circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 12. If the distribution and/or use of the Library is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Library under this License may add +an explicit geographical distribution limitation excluding those countries, +so that distribution is permitted only in or among countries not thus +excluded. In such case, this License incorporates the limitation as if +written in the body of this License. + + 13. The Free Software Foundation may publish revised and/or new +versions of the 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 +specifies a version number of this License which applies to it and +"any later version", you have the option of following the terms and +conditions either of that version or of any later version published by +the Free Software Foundation. If the Library does not specify a +license version number, you may choose any version ever published by +the Free Software Foundation. + + 14. If you wish to incorporate parts of the Library into other free +programs whose distribution conditions are incompatible with these, +write to the author to ask for permission. For software which is +copyrighted by the Free Software Foundation, write to the Free +Software Foundation; we sometimes make exceptions for this. Our +decision will be guided by the two goals of preserving the free status +of all derivatives of our free software and of promoting the sharing +and reuse of software generally. + + NO WARRANTY + + 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO +WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR +OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY +KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE +LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME +THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN +WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY +AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU +FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR +CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE +LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING +RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A +FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF +SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Libraries + + If you develop a new library, and you want it to be of the greatest +possible use to the public, we recommend making it free software that +everyone can redistribute and change. You can do so by permitting +redistribution under these terms (or, alternatively, under the terms of the +ordinary General Public License). + + To apply these terms, attach the following notices to the library. It is +safest to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least the +"copyright" line and a pointer to where the full notice is found. + + <one line to give the library's name and a brief idea of what it does.> + Copyright (C) <year> <name of author> + + This library 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 2.1 of the License, or (at your option) any later version. + + This library 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 library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Also add information on how to contact you by electronic and paper mail. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the library, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the + library `Frob' (a library for tweaking knobs) written by James Random Hacker. + + <signature of Ty Coon>, 1 April 1990 + Ty Coon, President of Vice + +That's all there is to it! + + @@ -0,0 +1,56 @@ + Installing MPC + ============== + +0. You first need to install GMP. See <http://www.swox.com/gmp/>. + MPC requires GMP version 4.1 or later. + + You need to configure GMP with --enable-mpfr, to install the MPFR + library, too. + +1. In the MPC build directory, type + + make + + This assumes that GMP is installed into /usr/local. + + Otherwise, type + + make GMP=<gmp_install_dir> + + where <gmp_install_dir> is the directory where you installed GMP. The GMP + header files like gmp.h are then expected to be in <gmp_install_dir>/include, + and the GMP library files like libgmp.a in <gmp_install_dir>/lib. + + It is also possible to specify a separate directory for MPFR by + + make GMP=<gmp_install_dir> MPFR=<mpfr_install_dir> + + which allows to use a newer version of MPFR. + +2. Type + + make check + + resp. + + make check GMP=<gmp_install_dir> + + resp. + + make check GMP=<gmp_install_dir> MPFR=<mpfr_install_dir> + + to produce and run the test files (please check that <gmp_install_dir>/lib + is first in your LD_LIBRARY_PATH environment variable if you use the dynamic + GMP libraries, or compile the test files with the -static option). + +3. Type + + make mpc.dvi + + or + make mpc.ps + + to produce the documentation in DVI or Postscript format. + +In case of difficulties, send a description of the problem to +<enge@lix.polytechnique.fr, zimmerma@loria.fr>. @@ -0,0 +1,62 @@ +/* mpc_abs -- Absolute value of a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +int +mpc_abs (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + mpfr_t u, v; + mp_prec_t prec; + int inexact; + + prec = MPFR_PREC(a); + + mpfr_init (u); + mpfr_init (v); + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 3; + + mpfr_set_prec (u, prec); + mpfr_set_prec (v, prec); + + /* first compute norm(b)^2 */ + inexact = mpfr_mul (u, MPC_RE(b), MPC_RE(b), GMP_RNDD); /* err<=1ulp */ + inexact |= mpfr_mul (v, MPC_IM(b), MPC_IM(b), GMP_RNDD); /* err<=1ulp*/ + + inexact |= mpfr_add (u, u, v, GMP_RNDD); /* err <= 3 ulps */ + inexact |= mpfr_sqrt (u, u, GMP_RNDD); /* err <= 4 ulps */ + } + while (inexact != 0 && + mpfr_can_round (u, prec - 2, GMP_RNDD, rnd, MPFR_PREC(a)) == 0); + + inexact |= mpfr_set (a, u, rnd); + + mpfr_clear (u); + mpfr_clear (v); + + return inexact; +} @@ -0,0 +1,36 @@ +/* mpc_add -- Add two complex numbers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +/* return 0 iff both the real and imaginary parts are exact */ +int +mpc_add (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_add (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); + inex_im = mpfr_add (MPC_IM(a), MPC_IM(b), MPC_IM(c), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/add_fr.c b/add_fr.c new file mode 100644 index 0000000..7575a08 --- /dev/null +++ b/add_fr.c @@ -0,0 +1,36 @@ +/* mpc_add_fr -- Add a complex number and a floating-point number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.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, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_add (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/add_ui.c b/add_ui.c new file mode 100644 index 0000000..fd7f3ac --- /dev/null +++ b/add_ui.c @@ -0,0 +1,36 @@ +/* mpc_add_ui -- Add a complex number and an unsigned long int. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.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, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_add_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,31 @@ +/* mpc_clear -- Clear a complex variable. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_clear (mpc_t x) +{ + mpfr_clear (MPC_RE(x)); + mpfr_clear (MPC_IM(x)); +} @@ -0,0 +1,36 @@ +/* mpc_cmp -- Compare two complex numbers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.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_RE(a), MPC_RE(b)); + cmp_im = mpfr_cmp (MPC_IM(a), MPC_IM(b)); + + return MPC_INEX(cmp_re, cmp_im); +} @@ -0,0 +1,35 @@ +/* mpc_conj -- Conjugate of a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_conj (mpc_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); + inex_im = mpfr_neg (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,153 @@ +/* mpc_div -- Divide two complex numbers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +int +mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int ok_re=0, ok_im=0; + mpc_t res, c_conj; + mpfr_t q; + mp_prec_t prec; + int inexact_prod, inexact_norm, inexact_re, inexact_im; + + prec = MPC_MAX_PREC(a); + + mpc_init (res); + mpfr_init (q); + + /* create the conjugate of c in c_conj without allocating new memory */ + MPC_RE (c_conj)[0] = MPC_RE (c)[0]; + MPC_IM (c_conj)[0] = MPC_IM (c)[0]; + MPFR_CHANGE_SIGN (MPC_IM (c_conj)); + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 5; + + mpc_set_prec (res, prec); + mpfr_set_prec (q, prec); + + /* first compute norm(c)^2 */ + inexact_norm = mpc_norm (q, c, GMP_RNDD); + + /* now compute b*conjugate(c) */ + /* We need rounding away from zero for both the real and the imagin- */ + /* ary part; then the final result is also rounded away from zero. */ + /* The error is less than 1 ulp. Since this is not implemented in */ + /* mpfr, we round towards zero and add 1 ulp to the absolute values */ + /* if they are not exact. */ + inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ); + inexact_re = MPC_INEX_RE (inexact_prod); + inexact_im = MPC_INEX_IM (inexact_prod); + if (inexact_re != 0) + mpfr_add_one_ulp (MPC_RE (res), GMP_RNDN); + if (inexact_im != 0) + mpfr_add_one_ulp (MPC_IM (res), GMP_RNDN); + + /* divide the product by the norm*/ + if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) + { + /* The divison 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. */ + if (MPFR_SIGN (MPC_RE (res)) > 0) + { + inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDU); + ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU, + MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a))); + } + else + { + inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDD); + ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD, + MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a))); + } + + if (ok_re || !inexact_re) /* compute imaginary part */ + { + if (MPFR_SIGN (MPC_IM (res)) > 0) + { + inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDU); + ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU, + MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a))); + } + else + { + inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDD); + ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD, + MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a))); + } + } + } + else + { + /* The division is inexact, so for efficiency reasons we invert q */ + /* only once and multiply by the inverse. */ + /* We do not decide about the sign of the difference. */ + inexact_re = 1; + inexact_im = 1; + mpfr_ui_div (q, 1, q, GMP_RNDU); + if (MPFR_SIGN (MPC_RE (res)) > 0) + { + mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDU); + ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU, + MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a))); + } + else + { + mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDD); + ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD, + MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a))); + } + + if (ok_re) /* compute imaginary part */ + { + if (MPFR_SIGN (MPC_IM (res)) > 0) + { + mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDU); + ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU, + MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a))); + } + else + { + mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDD); + ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD, + MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a))); + } + } + } + } + while ((!ok_re && inexact_re) || (!ok_im && inexact_im)); + + mpc_set (a, res, rnd); + + mpc_clear (res); + mpfr_clear (q); + + return (MPC_INEX (inexact_re, inexact_im)); + /* Only exactness vs. inexactness is tested, not the sign. */ +} diff --git a/div_2exp.c b/div_2exp.c new file mode 100644 index 0000000..0aa4604 --- /dev/null +++ b/div_2exp.c @@ -0,0 +1,35 @@ +/* mpc_div_2exp -- Divide a complex number by 2^e. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_div_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_div_2exp (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_2exp (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/div_fr.c b/div_fr.c new file mode 100644 index 0000000..02a673d --- /dev/null +++ b/div_fr.c @@ -0,0 +1,35 @@ +/* mpc_div_fr -- Divide a complex number by a floating-point number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_div_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_div (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/div_ui.c b/div_ui.c new file mode 100644 index 0000000..2bd560d --- /dev/null +++ b/div_ui.c @@ -0,0 +1,35 @@ +/* mpc_div_ui -- Divide a complex number by a nonnegative integer. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_div_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_div_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_ui (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,94 @@ +/* mpc_exp -- exponential of a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +void +mpc_exp (mpc_ptr rop, mpc_srcptr op, mp_rnd_t rnd) +{ + mpfr_t x, y, z; + mp_prec_t prec; + int ok = 0; + + /* let op = a + i*b, then exp(op) = exp(a)*[cos(b) + i*sin(b)] + = exp(a)*cos(b) + i*exp(a)*sin(b). + + We use the following algorithm (same for the imaginary part): + + (1) x = o(exp(a)) rounded towards +infinity: + (2) y = o(cos(b)) rounded to nearest + (3) r = o(x*y) + then the error on r for the real part is at most 4 ulps: + |r - exp(a)*cos(b)| <= ulp(r) + |x*y - exp(a)*cos(b)| + <= ulp(r) + |x*y - exp(a)*y| + exp(a) * |y - cos(b)| + <= ulp(r) + |y| ulp(x) + 1/2 * x * ulp(y) + <= ulp(r) + 2 * ulp(x*y) + ulp(x*y) [Rule 4] + <= 4 * ulp(r) [Rule 8] + */ + + /* special case when the input is real */ + if (mpfr_cmp_ui (MPC_IM(op), 0) == 0) + { + mpfr_exp (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd)); + mpfr_set_ui (MPC_IM(rop), 0, MPC_RND_IM(rnd)); + return; + } + + prec = MPC_MAX_PREC(rop); + + mpfr_init2 (x, 2); + mpfr_init2 (y, 2); + mpfr_init2 (z, 2); + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 5; + + mpfr_set_prec (x, prec); + mpfr_set_prec (y, prec); + mpfr_set_prec (z, prec); + + mpfr_exp (x, MPC_RE(op), GMP_RNDU); + mpfr_sin_cos (z, y, MPC_IM(op), GMP_RNDN); + mpfr_mul (y, y, x, GMP_RNDN); + ok = mpfr_can_round (y, prec - 2, GMP_RNDN, MPC_RND_RE(rnd), + MPFR_PREC(MPC_RE(rop))); + if (ok) /* compute imaginary part */ + { + mpfr_mul (z, z, x, GMP_RNDN); + ok = mpfr_can_round (z, prec - 2, GMP_RNDN, MPC_RND_IM(rnd), + MPFR_PREC(MPC_IM(rop))); + } + } + while (ok == 0); + + mpfr_set (MPC_RE(rop), y, MPC_RND_RE(rnd)); + mpfr_set (MPC_IM(rop), z, MPC_RND_IM(rnd)); + + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); + + +} @@ -0,0 +1,30 @@ +/* mpc_init -- Initialize a complex variable. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_init (mpc_t x) +{ + mpc_init2 (x, __mpc_default_fp_bit_precision); +} @@ -0,0 +1,31 @@ +/* mpc_init2 -- Initialize a complex variable with a given precision. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_init2 (mpc_t x, mp_prec_t prec) +{ + mpfr_init2 (MPC_RE(x), prec); + mpfr_init2 (MPC_IM(x), prec); +} @@ -0,0 +1,31 @@ +/* mpc_init3 -- Initialize a complex variable with given precisions. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_init3 (mpc_t x, mp_prec_t prec_re, mp_prec_t prec_im) +{ + mpfr_init2 (MPC_RE(x), prec_re); + mpfr_init2 (MPC_IM(x), prec_im); +} diff --git a/inp_str.c b/inp_str.c new file mode 100644 index 0000000..8d0c8e8 --- /dev/null +++ b/inp_str.c @@ -0,0 +1,65 @@ +/* mpc_inp_str -- Input a complex number from a given stream. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <ctype.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +size_t +mpc_inp_str (mpc_ptr rop, FILE *stream, int base, mp_rnd_t rnd_mode) +{ + size_t size, size_im; + int c; + + /* input real part */ + size = mpfr_inp_str (MPC_RE(rop), stream, base, MPC_RND_RE(rnd_mode)); + if (size == 0) /* error while reading the real part */ + return 0; + + /* get 'i' or 'I' followed by '*' */ + do + { + c = getc (stream); + size ++; + } + while (isspace (c)); + + if (c != '+') + return 0; /* error */ + + c = getc (stream); + size ++; + if (c != 'I') + return 0; /* error */ + + c = getc (stream); + size ++; + if (c != '*') + return 0; /* error */ + + size_im = mpfr_inp_str (MPC_IM(rop), stream, base, MPC_RND_IM(rnd_mode)); + if (size_im == 0) /* error while reading the imaginary part */ + return 0; + + return size + size_im; +} diff --git a/makefile b/makefile new file mode 100644 index 0000000..fa2b263 --- /dev/null +++ b/makefile @@ -0,0 +1,92 @@ +# Makefile for the MPC library. +# +# Copyright (C) 2002 Andreas Enge, Paul Zimmermann +# +# This file is part of the MPC Library. +# +# The MPC Library 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 2.1 of the License, or (at your +# option) any later version. +# +# The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +# MA 02111-1307, USA. + +# directory where GMP is installed (you need GMP 4.1 or later, configured +# with MPFR installed, i.e. ./configure --enable-mpfr): +# headers are expected in $(GMP)/include, lib in $(GMP)/lib +GMP=/usr/local +MPFR=$(GMP) + +######################## do not edit below this line ########################## + +VERSION=0.4.1 + +.SUFFIXES: .c .o + +OBJECTS= abs.o add.o add_fr.o add_ui.o clear.o cmp.o conj.o div.o div_2exp.o div_fr.o div_ui.o exp.o init.o init2.o init3.o inp_str.o mul.o mul_2exp.o mul_fr.o mul_ui.o neg.o norm.o out_str.o random.o random2.o set.o set_d_d.o set_dfl_prec.o set_prec.o set_ui_fr.o set_si_si.o set_ui_ui.o sqr.o sqrt.o sub.o sub_ui.o ui_div.o +SOURCES= abs.c add.c add_fr.c add_ui.c clear.c cmp.c conj.c div.c div_2exp.c div_fr.c div_ui.c exp.c init.c init2.c init3.c inp_str.c mul.c mul_2exp.c mul_fr.c mul_ui.c neg.c norm.c out_str.c random.c random2.c set.c set_d_d.c set_dfl_prec.c set_prec.c set_ui_fr.c set_si_si.c set_ui_ui.c sqr.c sqrt.c sub.c sub_ui.c ui_div.c +TESTS= test.c tmul.c tsqr.c tdiv.c texp.c +DIST= $(SOURCES) $(TESTS) makefile mpc.h mpc-impl.h COPYING.LIB mpc.texi INSTALL + +CFLAGS= -g -O2 -Wall -Wmissing-prototypes -ansi -pedantic + +.c.o: + gcc -I$(MPFR)/include -I$(GMP)/include -I. $(CFLAGS) -c $< + +libmpc.a: $(OBJECTS) + ar cru libmpc.a $(OBJECTS) + ranlib libmpc.a + +check: test tmul tsqr tdiv texp + @echo Testing all functions + rm -f mpc_test + ./test + rm -f mpc_test + @echo Testing mpc_mul + ./tmul + @echo Testing mpc_sqr + ./tsqr + @echo Testing mpc_div + ./tdiv + @echo Testing mpc_exp + ./texp + +test: test.c libmpc.a + gcc $(CFLAGS) -I$(MPFR)/include -I$(GMP)/include -I. -L. test.c -o test -lmpc $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a + +tmul: tmul.c libmpc.a + gcc $(CFLAGS) -I$(MPFR)/include -I$(GMP)/include -I. -L. tmul.c -o tmul -lmpc $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a + +tsqr: tsqr.c libmpc.a + gcc $(CFLAGS) -I$(MPFR)/include -I$(GMP)/include -I. -L. tsqr.c -o tsqr -lmpc $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a + +tdiv: tdiv.c libmpc.a + gcc $(CFLAGS) -I$(MPFR)/include -I$(GMP)/include -I. -L. tdiv.c -o tdiv -lmpc $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a + +texp: texp.c libmpc.a + gcc $(CFLAGS) -I$(MPFR)/include -I$(GMP)/include -I. -L. texp.c -o texp -lmpc $(MPFR)/lib/libmpfr.a $(GMP)/lib/libgmp.a + +clean: + rm -f *.o *~ libmpc.a test tmul tsqr tdiv texp mpc-$(VERSION).tar.gz mpc.aux mpc.cp mpc.cps mpc.dvi mpc.fn mpc.fns mpc.ky mpc.log mpc.pg mpc.ps mpc.toc mpc.tp mpc.vr mpc.vrs + +dist: $(DIST) + rm -fr mpc-$(VERSION) + mkdir mpc-$(VERSION) + cp $(DIST) mpc-$(VERSION) + tar cf mpc-$(VERSION).tar mpc-$(VERSION) + gzip --best mpc-$(VERSION).tar + rm -fr mpc-$(VERSION) + +mpc.dvi: mpc.texi + texi2dvi mpc.texi + +mpc.ps: mpc.dvi + dvips mpc -o diff --git a/mpc-impl.h b/mpc-impl.h new file mode 100644 index 0000000..8c6e774 --- /dev/null +++ b/mpc-impl.h @@ -0,0 +1,53 @@ +/* mpc-impl.h -- Internal include file for mpc. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +/* stolen from gmp-impl.h */ +#define MAX(h,i) ((h) > (i) ? (h) : (i)) + +#ifndef MUL_KARATSUBA_THRESHOLD +#define MUL_KARATSUBA_THRESHOLD 23 +#endif + +#ifndef BITS_PER_MP_LIMB +#define BITS_PER_MP_LIMB mp_bits_per_limb +#endif + + +/* avoids to include mpfr-impl.h */ +/* This unsigned type must correspond to the signed one defined in gmp.h */ +#if defined (_CRAY) && ! defined (_CRAYMPP) +typedef unsigned int mp_size_unsigned_t; +#else +typedef unsigned long int mp_size_unsigned_t; +#endif +#define MPFR_PREC(x) ((x)->_mpfr_prec) +#define MPFR_SIZE(x) ((x)->_mpfr_size) +#define MPFR_MANT(x) ((x)->_mpfr_d) +#define MPFR_EXP(x) ((x)->_mpfr_exp) +#define MPFR_CHANGE_SIGN(x) (MPFR_SIZE(x) ^= (((mp_size_unsigned_t) 1) << 31)) +#define MPFR_IS_ZERO(x) \ + (MPFR_MANT(x)[(MPFR_PREC(x)-1)/BITS_PER_MP_LIMB] == (mp_limb_t) 0) +extern long _mpfr_ceil_log2 _PROTO((double)); + +#define MPC_MAX_PREC(x) MAX(MPFR_PREC(MPC_RE(x)), MPFR_PREC(MPC_IM(x))) + +/* forgotten in mpfr.h from GMP 4.1 */ +int mpfr_cmp_abs _PROTO ((mpfr_srcptr, mpfr_srcptr)); @@ -0,0 +1,180 @@ +/* mpc.h -- Include file for mpc. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#ifndef __MPC_H +#define __MPC_H + +/* check if stdio.h is included */ +#if defined (FILE) || defined (H_STDIO) || defined (_H_STDIO) \ + || defined (_STDIO_H) || defined (_STDIO_H_) || defined (__STDIO_H__) \ + || defined (_STDIO_INCLUDED) || defined (__dj_include_stdio_h_) +#define _MPC_H_HAVE_FILE 1 +#endif + +/* Definition of rounding modes */ + +/* a complex rounding mode is just a pair of two real rounding modes + we reserve three bits for a real rounding mode, to have space to + add the "away from zero" rounding mode. + */ +#define RNDC(r1,r2) ((r1) + (r2 << 3)) +#define MPC_RND_RE(x) ((x) & 0x07) +#define MPC_RND_IM(x) ((x) >> 3) + +#define GMP_RNDA 4 /* round away from zero */ + +#define MPC_RNDNN RNDC(GMP_RNDN,GMP_RNDN) +#define MPC_RNDNZ RNDC(GMP_RNDN,GMP_RNDZ) +#define MPC_RNDNU RNDC(GMP_RNDN,GMP_RNDU) +#define MPC_RNDND RNDC(GMP_RNDN,GMP_RNDD) +#define MPC_RNDNA RNDC(GMP_RNDN,GMP_RNDA) + +#define MPC_RNDZN RNDC(GMP_RNDZ,GMP_RNDN) +#define MPC_RNDZZ RNDC(GMP_RNDZ,GMP_RNDZ) +#define MPC_RNDZU RNDC(GMP_RNDZ,GMP_RNDU) +#define MPC_RNDZD RNDC(GMP_RNDZ,GMP_RNDD) +#define MPC_RNDZA RNDC(GMP_RNDZ,GMP_RNDA) + +#define MPC_RNDUN RNDC(GMP_RNDU,GMP_RNDN) +#define MPC_RNDUZ RNDC(GMP_RNDU,GMP_RNDZ) +#define MPC_RNDUU RNDC(GMP_RNDU,GMP_RNDU) +#define MPC_RNDUD RNDC(GMP_RNDU,GMP_RNDD) +#define MPC_RNDUA RNDC(GMP_RNDU,GMP_RNDA) + +#define MPC_RNDDN RNDC(GMP_RNDD,GMP_RNDN) +#define MPC_RNDDZ RNDC(GMP_RNDD,GMP_RNDZ) +#define MPC_RNDDU RNDC(GMP_RNDD,GMP_RNDU) +#define MPC_RNDDD RNDC(GMP_RNDD,GMP_RNDD) +#define MPC_RNDDA RNDC(GMP_RNDD,GMP_RNDA) + +#define MPC_RNDAN RNDC(GMP_RNDA,GMP_RNDN) +#define MPC_RNDAZ RNDC(GMP_RNDA,GMP_RNDZ) +#define MPC_RNDAU RNDC(GMP_RNDA,GMP_RNDU) +#define MPC_RNDAD RNDC(GMP_RNDA,GMP_RNDD) +#define MPC_RNDAA RNDC(GMP_RNDA,GMP_RNDA) + +/* return values */ +/* transform negative to 2, positive to 1, leave 0 unchanged */ +#define MPC_INEX_POS(inex) (((inex) < 0) ? 2 : ((inex) == 0) ? 0 : 1) +/* transform 2 to negative, 1 to positive, leave 0 unchanged */ +#define MPC_INEX_NEG(inex) (((inex) == 2) ? -1 : ((inex) == 0) ? 0 : 1) + +#define MPC_INEX(inex_re, inex_im) \ + (MPC_INEX_POS(inex_re) | (MPC_INEX_POS(inex_im) << 2)) +#define MPC_INEX_RE(inex) MPC_INEX_NEG((inex) & 3) +#define MPC_INEX_IM(inex) MPC_INEX_NEG((inex) >> 2) + +/* Definitions of types and their semantics */ + +/* typedef int mp_rnd_t; */ /* imported from mpfr.h, + only 2*3 = 6 bits are used */ + +typedef struct { + mpfr_t re; + mpfr_t im; +} +__mpc_struct; + +#define MPC_RE(x) ((x)->re) +#define MPC_IM(x) ((x)->im) + +typedef __mpc_struct mpc_t[1]; +typedef __mpc_struct *mpc_ptr; +typedef __gmp_const __mpc_struct *mpc_srcptr; + +/* Prototypes */ + +#ifndef _PROTO +#if defined (__STDC__) || defined (__cplusplus) +#define _PROTO(x) x +#else +#define _PROTO(x) () +#endif +#endif + +#if defined (__cplusplus) +extern "C" { +#endif + +extern mp_prec_t __mpc_default_fp_bit_precision; + +int mpc_add _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_add_fr _PROTO ((mpc_ptr, mpc_srcptr, mpfr_srcptr, mp_rnd_t)); +int mpc_add_ui _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_sub _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_sub_ui _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_mul _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_mul_naive _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_mul_karatsuba _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_mul_fr _PROTO ((mpc_ptr, mpc_srcptr, mpfr_srcptr, mp_rnd_t)); +int mpc_mul_ui _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_sqr _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_div _PROTO ((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); +int mpc_div_fr _PROTO ((mpc_ptr, mpc_srcptr, mpfr_srcptr, mp_rnd_t)); +int mpc_div_ui _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_ui_div _PROTO ((mpc_ptr, unsigned long int, mpc_srcptr, mp_rnd_t)); +int mpc_div_2exp _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_mul_2exp _PROTO ((mpc_ptr, mpc_srcptr, unsigned long int, mp_rnd_t)); +int mpc_conj _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_neg _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_norm _PROTO ((mpfr_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_abs _PROTO ((mpfr_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_sqrt _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_set _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +int mpc_set_d_d _PROTO ((mpc_ptr, double, double, mp_rnd_t)); +int mpc_set_ui_fr _PROTO ((mpc_ptr, unsigned long int, mpfr_srcptr, mp_rnd_t)); +int mpc_set_ui_ui _PROTO ((mpc_ptr, unsigned long int, unsigned long int, mp_rnd_t)); +int mpc_set_si_si _PROTO ((mpc_ptr, long int, long int, mp_rnd_t)); +int mpc_cmp _PROTO ((mpc_srcptr, mpc_srcptr)); +void mpc_exp _PROTO ((mpc_ptr, mpc_srcptr, mp_rnd_t)); +void mpc_clear _PROTO ((mpc_ptr)); +void mpc_init _PROTO ((mpc_ptr)); +void mpc_random _PROTO ((mpc_ptr)); +void mpc_random2 _PROTO ((mpc_ptr, mp_size_t, mp_exp_t)); +void mpc_init2 _PROTO ((mpc_ptr, mp_prec_t)); +void mpc_init3 _PROTO ((mpc_ptr, mp_prec_t, mp_prec_t)); +void mpc_set_prec _PROTO ((mpc_ptr, mp_prec_t)); +void mpc_set_default_prec _PROTO ((mp_prec_t)); +mp_prec_t mpc_get_default_prec _PROTO ((void)); +#ifdef _MPC_H_HAVE_FILE +size_t mpc_inp_str _PROTO ((mpc_ptr, FILE *, int, mp_rnd_t)); +size_t mpc_out_str _PROTO ((FILE *, int, size_t, mpc_srcptr, mp_rnd_t)); +#endif + +#if defined (__cplusplus) +} +#endif + +#define mpc_set_d(x, y, rnd) mpc_set_d_d(x, y, 0.0, rnd) +#define mpc_set_ui(x, y, rnd) mpc_set_ui_ui(x, y, 0, rnd) +#define mpc_set_si(x, y, rnd) mpc_set_si_si(x, y, 0, rnd) +#define mpc_init_set(x, y, rnd) \ + ( mpc_init(x), mpc_set((x), (y), (rnd)) ) +#define mpc_init_set_ui_fr(x, y, z, rnd) \ + ( mpc_init(x), mpc_set_ui_fr((x), (y), (z), (rnd)) ) +#define mpc_init_set_ui_ui(x, y, z, rnd) \ + ( mpc_init(x), mpc_set_ui_ui((x), (y), (z), (rnd)) ) +#define mpc_init_set_si_si(x, y, z, rnd) \ + ( mpc_init(x), mpc_set_si_si((x), (y), (z), (rnd)) ) +#define mpc_add_si(x, y, z, rnd) \ + ( (z) >= 0 ? mpc_add_ui ((x), (y), (z), (rnd)) : mpc_sub_ui ((x), (y), -(z), (rnd)) ) + +#endif /* ifndef __MPC_H */ diff --git a/mpc.texi b/mpc.texi new file mode 100644 index 0000000..a1392d6 --- /dev/null +++ b/mpc.texi @@ -0,0 +1,829 @@ +\input texinfo @c -*-texinfo-*- +@c %**start of header +@setfilename mpc.info +@settitle MPC @value{VERSION} +@synindex tp fn +@iftex +@afourpaper +@end iftex +@comment %**end of header + +@set VERSION 0.4.1 +@set DATE {October 2002} + +@ifinfo +@format +START-INFO-DIR-ENTRY +* mpc: (mpc.info). Multiple Precision Complex Library. +END-INFO-DIR-ENTRY +@end format +@end ifinfo + +@c smallbook + +@iftex +@finalout +@end iftex + + +@ifinfo +This file documents MPC, a library for multiple precision complex arithmetic + +Copyright (C) 2002, Andreas Enge, Paul Zimmermann + +Permission is granted to make and distribute verbatim copies of +this manual provided the copyright notice and this permission notice +are preserved on all copies. + +@ignore +Permission is granted to process this file through TeX and print the +results, provided the printed document carries copying permission +notice identical to this one except for the removal of this paragraph +(this paragraph not being relevant to the printed manual). + +@end ignore +Permission is granted to copy and distribute modified versions of this +manual under the conditions for verbatim copying, provided that the entire +resulting derived work is distributed under the terms of a permission +notice identical to this one. + +Permission is granted to copy and distribute translations of this manual +into another language, under the above conditions for modified versions, +except that this permission notice may be stated in a translation approved +by the Foundation. +@end ifinfo + +@setchapternewpage on +@titlepage +@c use the new format for titles + +@title MPC +@subtitle The Multiple Precision Complex Library +@subtitle Edition @value{VERSION} +@subtitle @value{DATE} + +@author Andreas Enge, Paul Zimmermann + +@c Include the Distribution inside the titlepage so +@c that headings are turned off. + +@tex +\global\parindent=0pt +\global\parskip=8pt +\global\baselineskip=13pt +@end tex + +@page +@vskip 0pt plus 1filll +Copyright @copyright{} 2002 Andreas Enge, Paul Zimmermann + +@sp 2 + +@c Published by the Free Software Foundation @* +@c 59 Temple Place - Suite 330 @* +@c Boston, MA 02111-1307, USA @* + +Permission is granted to make and distribute verbatim copies of +this manual provided the copyright notice and this permission notice +are preserved on all copies. + +Permission is granted to copy and distribute modified versions of this +manual under the conditions for verbatim copying, provided that the entire +resulting derived work is distributed under the terms of a permission +notice identical to this one. + +Permission is granted to copy and distribute translations of this manual +into another language, under the above conditions for modified versions, +except that this permission notice may be stated in a translation approved +by the Foundation. +@end titlepage +@headings double + +@ifinfo +@node Top, Copying, (dir), (dir) + +@top MPC + +This manual documents how to install and use the Multiple Precision +Complex Library, version @value{VERSION} +@end ifinfo + +@menu +* Copying:: GMP Copying Conditions (LGPL). +* Introduction to MPC :: Brief introduction to MPC. +* Installing MPC:: How to configure and compile the MPC library. +* MPC Basics:: What every MPC user should now. +* Reporting Bugs:: How to usefully report bugs. +* Complex Functions:: Functions for arithmetic on complex numbers. + +* Contributors:: +* References:: +* Concept Index:: +* Function Index:: +@end menu + +@node Copying, Introduction to MPC, Top, Top +@comment node-name, next, previous, up +@unnumbered MPC Copying Conditions +@cindex Copying conditions +@cindex Conditions for copying MPC + +This library is @dfn{free}; this means that everyone is free to use it and +free to redistribute it on a free basis. The library is not in the public +domain; it is copyrighted and there are restrictions on its distribution, but +these restrictions are designed to permit everything that a good cooperating +citizen would want to do. What is not allowed is to try to prevent others +from further sharing any version of this library that they might get from +you.@refill + +Specifically, we want to make sure that you have the right to give away copies +of the library, that you receive source code or else can get it if you want +it, that you can change this library or use pieces of it in new free programs, +and that you know you can do these things.@refill + +To make sure that everyone has such rights, we have to forbid you to deprive +anyone else of these rights. For example, if you distribute copies of the +MPC library, you must give the recipients all the rights that you have. You +must make sure that they, too, receive or can get the source code. And you +must tell them their rights.@refill + +Also, for our own protection, we must make certain that everyone finds out +that there is no warranty for the MPC library. If it is modified by +someone else and passed on, we want their recipients to know that what they +have is not what we distributed, so that any problems introduced by others +will not reflect on our reputation.@refill + +The precise conditions of the license for the MPC library are found in the +Lesser General Public License that accompanies the source code. +See the file COPYING.LIB.@refill + +@node Introduction to MPC, Installing MPC, Copying, Top +@comment node-name, next, previous, up +@chapter Introduction to MPC + + +MPC is a portable library written in C for arbitrary precision arithmetic +on complex numbers. It is based on the GNU MP library. + +This version of MPC is released under the GNU Lesser General Public +License. +It is permitted to link MPC to non-free programs, as long as when +distributing them the MPC source code and a means to re-link with a +modified MPC is provided. + +@section How to use this Manual + +Everyone should read @ref{MPC Basics}. If you need to install the library +yourself, you need to read @ref{Installing MPC}, too. + +The rest of the manual can be used for later reference, although it is +probably a good idea to glance through it. + +@node Installing MPC, Reporting Bugs, Introduction to MPC, Top +@comment node-name, next, previous, up +@chapter Installing MPC +@cindex Installation + +To build MPC, you first have to install GNU MP +(version 4.1 or higher) on your computer. +You need to configure GNU MP with @code{--enable-mpfr}, to install +also the MPFR sub-library. +You need a C compiler, preferably GCC, but any reasonable compiler should +work. And you need a standard Unix @samp{make} program, plus some other +standard Unix utility programs. + +Here are the steps needed to install the library on Unix systems: + +@enumerate +@item +@samp{make} + +if GMP is installed into the standard directory @samp{/usr/local}. +This will compile MPC, and create a library archive file @file{libmpc.a} +in the working directory. + +@samp{make GMP=<gmp_install_dir>}. + +is used to indicate a different location of GMP. + +@samp{make GMP=<gmp_install_dir> MPFR=<mpfr_install_dir>}. + +allows to use a different version of MPFR than that provided with GMP. + + +@item +@samp{make check} resp. + +@samp {make check GMP=<gmp_install_dir>} resp. + +@samp {make check GMP=<gmp_install_dir> MPFR=<mpfr_install_dir>} + +This will make sure MPC was built correctly. +If you get error messages, please +report them to @samp{enge@@lix.polytechnique.fr,zimmerma@@loria.fr}. +(@xref{Reporting Bugs}, for +information on what to include in useful bug reports.) + +@end enumerate + +There are some other useful make targets: + +@itemize @bullet +@item +@samp{make mpc.dvi} + +Create a DVI version of the manual, in @file{mpc.dvi}. + +@item +@samp{make mpc.ps} + +Create a Postscript version of the manual, in @file{mpc.ps}. + +@item +@samp{make clean} + +Delete all object and archive files. + +@end itemize + + +@section Known Build Problems + +MPC suffers from all bugs from the GNU MP and MPFR libraries, +plus many many more. + +Please report other problems to +@samp{enge@@lix.polytechnique.fr,zimmerma@@loria.fr}. +@xref{Reporting Bugs}. + +@node Reporting Bugs, MPC Basics, Installing MPC, Top +@comment node-name, next, previous, up +@chapter Reporting Bugs +@cindex Reporting bugs + +If you think you have found a bug in the MPC library, +please investigate +and report it. We have made this library available to you, and it is not to ask +too much from you, to ask you to report the bugs that you find. + +There are a few things you should think about when you put your bug report +together. + +You have to send us a test case that makes it possible for us to reproduce the +bug. Include instructions on how to run the test case. + +You also have to explain what is wrong; if you get a crash, or if the results +printed are incorrect and in that case, in what way. + +Please include compiler version information +in your bug report. This can be extracted using @samp{cc -V} on some +machines, or, +if you are using gcc, @samp{gcc -v}. Also, include the output from @samp{uname +-a}. + +If your bug report is good, we will do our best to help you to get a corrected +version of the library; if the bug report is poor, we will not do anything about +it (aside of chiding you to send better bug reports). + +Send your bug report to: @samp{enge@@lix.polytechnique.fr,zimmerma@@loria.fr}. + +If you think something in this manual is unclear, or downright incorrect, or if +the language needs to be improved, please send a note to the same address. + +@node MPC Basics, Complex Functions, Reporting Bugs, Top +@comment node-name, next, previous, up +@chapter MPC Basics + + +@cindex @file{mpc.h} +All declarations needed to use MPC are collected in the include file +@file{mpc.h}. It is designed to work with both C and C++ compilers. +You should include that file in any program using the MPC library. +As MPC is based on the MPFR library, itself based on GNU MP, you +should also include the corresponding header files in the following order: +@example + #include "gmp.h" + #include "mpfr.h" + #include "mpc.h" +@end example + +@section Nomenclature and Types + +@cindex Complex number +@tindex @code{mpc_t} +@noindent +@dfn{Complex number} or @dfn{Complex} for short, is a pair of two +arbitrary precision floating-point numbers (for the real and imaginary parts). +The C data type for such objects is @code{mpc_t}. + +@cindex Precision +@tindex @code{mp_prec_t} +@noindent +The @dfn{Precision} is the number of bits used to represent the mantissa +of the real and imaginary parts; +the corresponding C data type is @code{mp_prec_t}. +See the MPFR documentation for more details on the allowed precision range. + +@cindex Rounding Mode +@tindex @code{mp_rnd_t} +@noindent +The @dfn{rounding mode} specifies the way to round the result of a +complex operation, in case the exact result can not be represented +exactly in the destination mantissa; +the corresponding C data type is @code{mp_rnd_t}. +A complex rounding mode is a pair of two rounding modes: one for the real +part, one for the imaginary part. + +@section Function Classes + +There is only one class of functions in the MPC library, namely functions for +complex arithmetic. The function names begin with @code{mpc_}. The +associated type is @code{mpc_t}. + + +@section MPC Variable Conventions + +As a general rule, all MPC functions expect output arguments before input +arguments. This notation is based on an analogy with the assignment operator. + +MPC allows you to use the same variable for both input and output in the same +expression. For example, the main function for floating-point multiplication, +@code{mpc_mul}, can be used like this: @code{mpc_mul (x, x, x, rnd_mode)}. +This +computes the square of @var{x} with rounding mode @code{rnd_mode} +and puts the result back in @var{x}. + +Before you can assign to an MPC variable, you need to initialize it by calling +one of the special initialization functions. When you are done with a +variable, you need to clear it out, using one of the functions for that +purpose. + +A variable should only be initialized once, or at least cleared out between +each initialization. After a variable has been initialized, it may be +assigned to any number of times. + +For efficiency reasons, avoid to initialize and clear out a variable in loops. +Instead, initialize it before entering the loop, and clear it out after the +loop has exited. + +You do not need to be concerned about allocating additional space for MPC +variables, since each of its real and imaginary part +has a mantissa of fixed size. +Hence unless you change its precision, or clear and reinitialize it, +a complex variable will have the same allocated space during all its +life. + +@node Complex Functions, Contributors, MPC Basics, Top +@comment node-name, next, previous, up +@chapter Complex Functions +@cindex Complex functions + +The complex functions expect arguments of type @code{mpc_t}. + +The MPC floating-point functions have an interface that is similar to the +GNU MP +integer functions. The function prefix for operations on complex numbers is +@code{mpc_}. + +@cindex User-defined precision +The precision of a computation is defined as follows: Compute the requested +operation exactly (with ``infinite precision''), and round the result to +the destination variable precision with the given rounding mode. + +The MPC complex functions are intended to be a smooth extension +of the IEEE P754 arithmetic. The results obtained on one +computer should not differ from the results obtained on a computer with a +different word size. + +@menu +* Rounding Modes:: +* Initializing Complex Numbers:: +* Assigning Complex Numbers:: +* Simultaneous Init & Assign:: +* Complex Arithmetic:: +* Complex Comparison:: +* Special Complex Functions:: +* I/O of Complex Numbers:: +* Miscellaneous Complex Functions:: +* Internals:: +@end menu + +@node Rounding Modes, Initializing Complex Numbers, , Complex Functions +@comment node-name, next, previous, up +@cindex Rounding modes +@section Rounding Modes + +A complex rounding mode is of the form @code{MPC_RNDxy} where +@code{x} and @code{y} are one of @code{N} (to nearest), @code{Z} (towards +zero), @code{U} (towards plus infinity), @code{D} (towards minus infinity). +The first letter refers to the rounding mode for the real part, +and the second one for the imaginary part. +For example @code{MPC_RNDZU} indicates to round the real part towards zero, +and the imaginary part towards plus infinity. + +The @samp{round to nearest} mode works as in the IEEE P754 standard: in case +the number to be rounded lies exactly in the middle of two representable +numbers, it is rounded to the one with the least significant bit set to zero. +For example, the number 5, which is represented by (101) in binary, is rounded +to (100)=4 with a precision of two bits, and not to (110)=6. + +Most MPC functions have a return value of type @code{int}, which is used +to indicate the position of the rounded real or imaginary parts with respect +to the exact (infinite precision) values. +If this integer is @code{i}, the macros @code{MPC_INEX_RE(i)} and +@code{MPC_INEX_IM(i)} give 0 if the corresponding rounded value is exact, +a negative value if the rounded value is less than the exact one, +and a positive value if it is greater than the exact one. +However, some functions do not completely fulfill this: +in some cases the sign is not guaranteed, and in some cases a non-zero value +is returned although the result is exact; in these cases the +function documentation explains the exact meaning of the return value. +However, the return value never wrongly indicates an exact computation. + +@node Initializing Complex Numbers, Assigning Complex Numbers, Rounding Modes, Complex Functions +@comment node-name, next, previous, up +@section Initialization Functions + +@deftypefun void mpc_set_default_prec (mp_prec_t @var{prec}) +Set the default precision to be @strong{exactly} @var{prec} bits. +All subsequent calls to @code{mpc_init} will use this precision, but previously +initialized variables are unaffected. +This default precision is set to 53 bits initially. +It is valid for the real and the imaginary parts alike. +@end deftypefun + +@deftypefun mp_prec_t mpc_get_default_prec () +Returns the default MPC precision in bits. +@end deftypefun + +An @code{mpc_t} object must be initialized before storing the first value in +it. The functions @code{mpc_init}, @code{mpc_init2} and @code{mpc_init3} +are used for that purpose. + +@deftypefun void mpc_init (mpc_t @var{z}) +Initialize @var{z}, and set its real and imaginary parts to NaN. + +Normally, a variable should be initialized once only +or at least be cleared, using @code{mpc_clear}, between initializations. The +precision of @var{x} is the default precision, which can be changed +by a call to @code{mpc_set_default_prec}. +@end deftypefun + +@deftypefun void mpc_init2 (mpc_t @var{z}, mp_prec_t @var{prec}) +Initialize @var{z}, set its precision to be +@var{prec} bits, and set its real and imaginary parts to NaN. +@end deftypefun + +@deftypefun void mpc_init3 (mpc_t @var{z}, mp_prec_t @var{prec_r}, mp_prec_t @var{prec_i}) +Initialize @var{z}, set the precision of its real part to @var{prec_r} bits, +the precision of its imaginary part to @var{prec_i} bits, +and set its real and imaginary parts to NaN. +@end deftypefun + +@deftypefun void mpc_clear (mpc_t @var{z}) +Free the space occupied by @var{z}. Make sure to call this function for all +@code{mpc_t} variables when you are done with them. +@end deftypefun + +@need 2000 +Here is an example on how to initialize complex variables: +@example +@{ + mpc_t x, y, z; + mpc_init (x); /* use default precision */ + mpc_init2 (y, 256); /* precision @emph{exactly} 256 bits */ + mpc_init3 (z, 100, 50); /* 100/50 bits for the real/imaginary part */ + @dots{} + /* Unless the program is about to exit, do ... */ + mpc_clear (x); + mpc_clear (y); + mpc_clear (z); +@} +@end example + +The following function is useful for changing the precision during a +calculation. A typical use would be for adjusting the precision gradually in +iterative algorithms like Newton-Raphson, making the computation precision +closely match the actual accurate part of the numbers. + +@deftypefun void mpc_set_prec (mpc_t @var{x}, mp_prec_t @var{prec}) +Reset the precision of @var{x} to be @strong{exactly} @var{prec} bits, +and set its real/imaginary parts to NaN. +The previous value stored in @var{x} is lost. It is equivalent to +a call to @code{mpc_clear(x)} followed by a call to +@code{mpc_init2(x, prec)}, but more efficient as no allocation is done in +case the current allocated space for the mantissa of @var{x} is sufficient. +@end deftypefun + +@node Assigning Complex Numbers, Simultaneous Init & Assign, Initializing Complex Numbers, Complex Functions +@comment node-name, next, previous, up +@section Assignment Functions +@cindex Complex assignment functions + +These functions assign new values to already initialized complex numbers +(@pxref{Initializing Complex Numbers}). + +@deftypefn Function int mpc_set (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_set_ui (mpc_t @var{rop}, unsigned long int @var{op}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_set_si (mpc_t @var{rop}, long int @var{op}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_set_d (mpc_t @var{rop}, double @var{op}, mp_rnd_t @var{rnd}) +Set the value of @var{rop} from @var{op}, rounded to the precision of @var{rop} +with the given rounding mode @var{rnd}. +Except for @code{mpc_set}, the argument @var{op} is interpreted as real, +so the imaginary part of @var{rop} is set to zero. +Please note that even a @code{long int} may have to be rounded, +if the destination precision is less than the machine word width. +For @code{mpc_set_d}, be careful that the input number @var{op} +may not be exactly representable as a double-precision number (this happens for +0.1 for instance), in which case it is first +rounded by the C compiler to a double-precision number, and then only +to a complex number. +@end deftypefn + +@deftypefun int mpc_set_d_d (mpc_t @var{rop}, double @var{op1}, double @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_set_ui_ui (mpc_t @var{rop}, unsigned long int @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_set_si_si (mpc_t @var{rop}, long int @var{op1}, long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_set_ui_fr (mpc_t @var{rop}, unsigned long int @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) +Set the real part of @var{rop} from @var{op1}, and its imaginary part from +@var{op2}, according to the rounding mode @var{rnd}. +@end deftypefun + + +@node Simultaneous Init & Assign, Complex Arithmetic, Assigning Complex Numbers, Complex Functions +@comment node-name, next, previous, up +@section Combined Initialization and Assignment Functions +@cindex Initialization and assignment functions + +@deftypefn Macro int mpc_init_set (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_init_set_ui (mpc_t @var{rop}, unsigned long int @var{op}, mp_rnd_t @var{rnd}) +Initialize @var{rop} and set its value from @var{op}, rounded with the +rounding mode @var{rnd}. +The precision of @var{rop} will be taken from the active default precision, +as set by @code{mpc_set_default_prec}. +@end deftypefn + +@deftypefn Macro int mpc_init_set_ui_ui (mpc_t @var{rop}, unsigned long int @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_init_set_ui_fr (mpc_t @var{rop}, unsigned long int @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefnx Macro int mpc_init_set_si_si (mpc_t @var{rop}, long int @var{op1}, long int @var{op2}, mp_rnd_t @var{rnd}) +Initialize @var{rop}, set its real part from @var{op1}, and its imaginary +part from @var{op2}, rounded with the rounding mode @var{rnd}. +The precision of @var{rop} will be taken from the active default precision, +as set by @code{mpc_set_default_prec}. +@end deftypefn + +@node Complex Arithmetic, Complex Comparison, Simultaneous Init & Assign, Complex Functions +@comment node-name, next, previous, up +@section Basic Arithmetic Functions +@cindex Complex arithmetic functions +@cindex Arithmetic functions + +All the following functions are designed in such a way that, when working +with real numbers instead of complex numbers, their complexity should +essentially be the same as with the MPFR library, with only a marginal +overhead due to the MPC layer. + +@deftypefun int mpc_add (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_add_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_add_fr (mpc_t @var{rop}, mpc_t @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1} @math{+} @var{op2} rounded according to @var{rnd}. +@end deftypefun + +@deftypefun int mpc_sub (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_sub_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1} @minus{} @var{op2} rounded according to @var{rnd}. +@end deftypefun + +@deftypefun int mpc_mul (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_mul_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_mul_fr (mpc_t @var{rop}, mpc_t @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1} times @var{op2} rounded according to @var{rnd}. +@end deftypefun + +@deftypefun int mpc_sqr (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set @var{rop} to the square of @var{op} rounded according to @var{rnd}. +@end deftypefun + +@deftypefun int mpc_div (mpc_t @var{rop}, mpc_t @var{op1}, mpc_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_div_ui (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_ui_div (mpc_t @var{rop}, unsigned long int @var{op1}, mpc_t @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpc_div_fr (mpc_t @var{rop}, mpc_t @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1}/@var{op2} rounded according to @var{rnd}. +For @code{mpc_div} and @code{mpc_ui_div}, the return value may fail to +recognize some exact results, and its sign is not significant. +@end deftypefun + +@deftypefun int mpc_sqrt (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set @var{rop} to the square root of @var{op} rounded according to @var{rnd}. +Here, when the return value is 0, it means the result is exact, but the +contrary may be false. +@end deftypefun + +@deftypefun int mpc_neg (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set @var{rop} to @minus{}@var{op} rounded according to @var{rnd}. +Just changes the sign if @var{rop} and @var{op} are the same variable. +@end deftypefun + +@deftypefun int mpc_conj (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set @var{rop} to the conjugate of @var{op} rounded according to @var{rnd}. +Just changes the sign of the imaginary part +if @var{rop} and @var{op} are the same variable. +@end deftypefun + +@deftypefun int mpc_abs (mpfr_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set the floating-point number @var{rop} to the absolute value of @var{op}, +rounded in the direction @var{rnd}. +The returned value is zero iff the result is exact. +Note the destination is of type @code{mpfr_t}, not @code{mpc_t}. +@end deftypefun + +@deftypefun int mpc_norm (mpfr_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set the floating-point number @var{rop} to the norm of @var{op} +(i.e. the square of its absolute value), +rounded in the direction @var{rnd}. +The returned value is zero iff the result is exact. +Note that the destination is of type @code{mpfr_t}, not @code{mpc_t}. +@end deftypefun + +@deftypefun int mpc_mul_2exp (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1} times 2 raised to @var{op2} +rounded according to @var{rnd}. Just increases the exponents +of the real and imaginary parts by @var{op2} +when @var{rop} and @var{op1} are identical. +@end deftypefun + +@deftypefun int mpc_div_2exp (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) +Set @var{rop} to @var{op1} divided by 2 raised to @var{op2} +rounded according to @var{rnd}. Just decreases the exponents +of the real and imaginary parts by @var{op2} +when @var{rop} and @var{op1} are identical. +@end deftypefun + +@node Complex Comparison, Special Complex Functions, Complex Arithmetic, Complex Functions +@comment node-name, next, previous, up +@section Comparison Functions +@cindex Complex comparisons functions +@cindex Comparison functions + +@deftypefun int mpc_cmp (mpc_t @var{op1}, mpc_t @var{op2}) +Compare @var{op1} and @var{op2}. +The return value @var{c} can be decomposed into @code{x = MPC_INEX_RE(c)} +and @code{y = MPC_INEX_IM(c)}, such that @var{x} is +positive if the real part of @var{op1} is greater than that of @var{op2}, +zero if both real parts are equal, and negative if the real part of @var{op1} +is less than that of @var{op2}, and likewise for @var{y}. +Both @var{op1} and @var{op2} are considered to their full own precision, +which may differ. +It is not allowed that one of the operands has a NaN (Not-a-Number) part. + +The storage of the return value is such that equality can be simply checked +with @code{mpc_cmp (op1, op2) == 0}. +@end deftypefun + +@node Special Complex Functions, I/O of Complex Numbers, Complex Comparison, Complex Functions +@comment node-name, next, previous, up +@section Special Functions +@cindex Special functions +@cindex Special complex functions + +@deftypefun void mpc_exp (mpc_t @var{rop}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Set @var{rop} to the exponential of @var{op}, +rounded according to @var{rnd} with the precision of @var{rop}. +@end deftypefun + +@node I/O of Complex Numbers, Miscellaneous Complex Functions, Special Complex Functions, Complex Functions +@comment node-name, next, previous, up +@section Input and Output Functions +@cindex Complex input and output functions +@cindex Input functions +@cindex Output functions +@cindex I/O functions + +Functions that perform input from a standard input/output +stream, and functions that output to +a standard input/output stream. +Passing a null pointer for a @var{stream} argument to any of +these functions will make them read from @code{stdin} and write to +@code{stdout}, respectively. + +When using any of these functions, it is a good idea to include @file{stdio.h} +before @file{mpc.h}, since that will allow @file{mpc.h} to define prototypes +for these functions. + +@deftypefun size_t mpc_out_str (FILE *@var{stream}, int @var{base}, size_t @var{n_digits}, mpc_t @var{op}, mp_rnd_t @var{rnd}) +Output @var{op} on stdio stream @var{stream}, in +base @var{base}, rounded according to @var{rnd}. +First the real part is printed, then @code{+I*} followed by the imaginary part. +The base may vary from 2 to 36. Print at most +@var{n_digits} significant digits for each part, +or if @var{n_digits} is 0, the maximum +number of digits accurately representable by @var{op}. + +In addition to the significant digits, a decimal point at the right of the +first digit and a +trailing exponent, in the form @samp{eNNN}, are printed. If @var{base} +is greater than 10, @samp{@@} will be used instead of @samp{e} as +exponent delimiter. + +Return the number of characters written. +@end deftypefun + +@deftypefun size_t mpc_inp_str (mpc_t @var{rop}, FILE *@var{stream}, int @var{base}, mp_rnd_t @var{rnd}) +Input a string in base @var{base} from stdio stream @var{stream}, +rounded according to @var{rnd}, and put the +read complex in @var{rop}. +Each of the real and imaginary part should be of the form @samp{M@@N} or, +if the base is 10 or less, alternatively @samp{MeN} or @samp{MEN}. +@samp{M} is the mantissa and +@samp{N} is the exponent. The mantissa is always in the specified base. The +exponent is always read in decimal. +This function first reads the real part, then @code{+} followed by +@code{I*} and the imaginary part. + +The argument @var{base} may be in the range 2 to 36. + +Return the number of bytes read, or if an error occurred, return 0. +@end deftypefun + +@node Miscellaneous Complex Functions, Internals, I/O of Complex Numbers, Complex Functions +@comment node-name, next, previous, up +@section Miscellaneous Functions +@cindex Miscellaneous complex functions + +@deftypefun void mpc_random (mpc_t @var{rop}) +Generate a random complex, with real and imaginary parts +uniformly distributed in the interval -1 < X < 1. +@end deftypefun + +@deftypefun void mpc_random2 (mpc_t @var{rop}, mp_size_t @var{max_size}, mp_exp_t @var{max_exp}) +Generate a random complex, with real and imaginary part +of at most @var{max_size} limbs, with long strings of +zeros and ones in the binary representation. The exponent of the +real (resp. imaginary) part is in +the interval @minus{}@var{exp} to @var{exp}. +This function is useful for +testing functions and algorithms, since this kind of random numbers have +proven to be more likely to trigger corner-case bugs. +Negative parts are generated when @var{max_size} is negative. +@end deftypefun + +@node Internals, , Miscellaneous Complex Functions, Complex Functions +@comment node-name, next, previous, up +@section Internals + +These types and +functions were mainly designed for the implementation of MPC, +but may be useful for users too. +However no upward compatibility is guaranteed. +You need to include @code{mpc-impl.h} to use them. + +The @code{mpc_t} type consists of two fields of type @code{mpfr_t}, +one for the real part, one for the imaginary part. +These fields can be accessed through @code{MPC_RE(z)} and @code{MPC_IM(z)}. + +Normally the real and imaginary part have the same precision, but +the function @code{mpc_init3} enables one to have different precisions, +and the user may also use @code{mpfr_set_prec} to change their precision. +The macro @code{MPC_MAX_PREC(z)} gives the maximum of the precisions +of the real and imaginary parts. + +@node Contributors, References, Complex Functions, Top +@comment node-name, next, previous, up +@unnumbered Contributors + +The developers of the MPC library are Andreas Enge and +Paul Zimmermann. + +@node References, Concept Index, Contributors, Top +@comment node-name, next, previous, up +@unnumbered References + +@itemize @bullet + +@item +Torbjorn Granlund, "GNU MP: The GNU Multiple Precision Arithmetic Library", + version 4.1, 2002. + +@item +IEEE standard for binary floating-point arithmetic, Technical Report +ANSI-IEEE Standard 754-1985, New York, 1985. +Approved March 21, 1985: IEEE Standards Board; approved July 26, + 1985: American National Standards Institute, 18 pages. + +@item +Donald E. Knuth, "The Art of Computer Programming", vol 2, +"Seminumerical Algorithms", 2nd edition, Addison-Wesley, 1981. + +@end itemize + +@node Concept Index, Function Index, References, Top +@comment node-name, next, previous, up +@unnumbered Concept Index +@printindex cp + +@node Function Index, , Concept Index, Top +@comment node-name, next, previous, up +@unnumbered Function and Type Index +@printindex fn + +@contents +@bye @@ -0,0 +1,382 @@ +/* mpc_mul -- Multiply two complex numbers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" +#include <stdio.h> + +#ifdef MPFR202 +#define MPFR_CMP_ABS mpfr_cmp_abs +#else /* mpfr-2.0.1 doesn't allow 0 in mpfr_cmp_abs */ +#define MPFR_CMP_ABS(a,b) \ + ((MPFR_IS_ZERO (a)) ? 0 : ((MPFR_IS_ZERO (b)) ? 1 : mpfr_cmp_abs(a, b))) +#endif +#define INV_RND(r) \ + (((r) == GMP_RNDU) ? GMP_RNDD : (((r) == GMP_RNDD) ? GMP_RNDU : (r))) +#define SWAP(a,b) { mpfr_srcptr tmp; tmp = a; a = b; b = tmp; } + +/* return inex such that MPC_INEX_RE(inex) = -1, 0, 1 + MPC_INEX_IM(inex) = -1, 0, 1 + depending on the signs of the real/imaginary parts of the result +*/ +int +mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + int overlap; + mpc_t result; + + /* first check for real multiplication */ + if (MPFR_IS_ZERO(MPC_IM(b))) /* b * (x+i*y) = b*x + i *(b*y) */ + { + /* first start with imaginary part in case a=b */ + inex_im = mpfr_mul (MPC_IM(a), MPC_RE(b), MPC_IM(c), MPC_RND_IM(rnd)); + inex_re = mpfr_mul (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); + return MPC_INEX(inex_re, inex_im); + } + if (MPFR_IS_ZERO(MPC_IM(c))) + { + inex_im = mpfr_mul (MPC_IM(a), MPC_RE(c), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_mul (MPC_RE(a), MPC_RE(c), MPC_RE(b), MPC_RND_RE(rnd)); + return MPC_INEX(inex_re, inex_im); + } + /* check for purely complex multiplication */ + if (MPFR_IS_ZERO(MPC_RE(b))) /* i*b * (x+i*y) = -b*y + i*(b*x) */ + { + overlap = (a == b || a == c); + if (overlap) + mpc_init3 (result, MPFR_PREC (MPC_RE (a)), + MPFR_PREC (MPC_IM (a))); + else + result [0] = a [0]; + inex_re = -mpfr_mul (MPC_RE(result), MPC_IM(b), MPC_IM(c), + INV_RND(MPC_RND_RE(rnd))); + mpfr_neg (MPC_RE (result), MPC_RE (result), GMP_RNDN); + inex_im = mpfr_mul (MPC_IM(result), MPC_IM(b), MPC_RE(c), + MPC_RND_IM(rnd)); + mpc_set (a, result, MPC_RNDNN); + if (overlap) + mpc_clear (result); + return MPC_INEX(inex_re, inex_im); + } + if (MPFR_IS_ZERO(MPC_RE(c))) + { + overlap = (a == b || a == c); + if (overlap) + mpc_init3 (result, MPFR_PREC (MPC_RE (a)), + MPFR_PREC (MPC_IM (a))); + else + result [0] = a [0]; + inex_re = -mpfr_mul (MPC_RE(result), MPC_IM(c), MPC_IM(b), + INV_RND(MPC_RND_RE(rnd))); + mpfr_neg (MPC_RE (result), MPC_RE (result), GMP_RNDN); + inex_im = mpfr_mul (MPC_IM(result), MPC_IM(c), MPC_RE(b), + MPC_RND_IM(rnd)); + mpc_set (a, result, MPC_RNDNN); + if (overlap) + mpc_clear (result); + return MPC_INEX(inex_re, inex_im); + } + + return ((MPC_MAX_PREC(a) <= MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB) + ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd); +} + +int +mpc_mul_naive (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int overlap, inex_re, inex_im; + mpfr_t u, v, t; + mp_prec_t prec; + + overlap = (a == b) || (a == c); + + prec = MPC_MAX_PREC(b) + MPC_MAX_PREC(c); + + mpfr_init2 (u, prec); + mpfr_init2 (v, prec); + + /* Re(a) = Re(b)*Re(c) - Im(b)*Im(c) */ + mpfr_mul (u, MPC_RE(b), MPC_RE(c), MPC_RNDNN); /* exact */ + mpfr_mul (v, MPC_IM(b), MPC_IM(c), MPC_RNDNN); /* exact */ + + if (overlap) + { + mpfr_init2 (t, MPFR_PREC(MPC_RE(a))); + inex_re = mpfr_sub (t, u, v, MPC_RND_RE(rnd)); + } + else + inex_re = mpfr_sub (MPC_RE(a), u, v, MPC_RND_RE(rnd)); + + /* Im(a) = Re(b)*Im(c) + Im(b)*Re(c) */ + mpfr_mul (u, MPC_RE(b), MPC_IM(c), MPC_RNDNN); /* exact */ + if (b == c) /* square case */ + inex_im = mpfr_mul_2exp (MPC_IM(a), u, 1, MPC_RND_IM(rnd)); + else + { + mpfr_mul (v, MPC_IM(b), MPC_RE(c), MPC_RNDNN); /* exact */ + inex_im = mpfr_add (MPC_IM(a), u, v, MPC_RND_IM(rnd)); + } + + mpfr_clear (u); + mpfr_clear (v); + + if (overlap) + { + mpfr_set (MPC_RE(a), t, GMP_RNDN); /* exact */ + mpfr_clear (t); + } + + return MPC_INEX(inex_re, inex_im); +} + + +/* Karatsuba multiplication, with 3 real multiplies */ +int +mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mp_rnd_t rnd) +{ + mpfr_srcptr a, b, c, d; + int mul_i, ok, inexact, mul_a, mul_c, inex_re, inex_im, sign_x, sign_u; + mpfr_t u, v, w, x; + mp_prec_t prec, prec_re; + mp_rnd_t rnd_re, rnd_u, rnd_x; + int overlap; + /* true if rop == op1 or rop == op2 */ + mpc_t result; + /* overlap is quite difficult to handle, because we have to tentatively + round the variable u in the end to either the real or the imaginary + part of rop (it is not possible to tell now whether the real or + 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. */ + + overlap = (rop == op1) || (rop == op2); + if (overlap) + mpc_init3 (result, MPFR_PREC (MPC_RE (rop)), + MPFR_PREC (MPC_IM (rop))); + else + result [0] = rop [0]; + + a = MPC_RE(op1); + b = MPC_IM(op1); + c = MPC_RE(op2); + d = MPC_IM(op2); + + /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */ + + mul_i = 0; /* number of multiplications by i */ + mul_a = 1; /* implicit factor for a */ + mul_c = 1; /* implicit factor for c */ + + if (MPFR_CMP_ABS (a, b) < 0) + { + SWAP(a, b); + mul_i ++; + mul_a = -1; /* consider i * (a+i*b) = -b + i*a */ + } + + if (MPFR_CMP_ABS (c, d) < 0) + { + SWAP(c, d); + mul_i ++; + mul_c = -1; /* consider -d + i*c instead of c + i*d */ + } + + /* find the precision and rounding mode for the new real part. + */ + if (mul_i % 2) + { + prec_re = MPFR_PREC(MPC_IM(rop)); + rnd_re = MPC_RND_IM(rnd); + } + else /* mul_i = 0 or 2 */ + { + prec_re = MPFR_PREC(MPC_RE(rop)); + rnd_re = MPC_RND_RE(rnd); + } + + if (mul_i) + rnd_re = INV_RND(rnd_re); + + /* now |a| >= |b| and |c| >= |d| */ + prec = MPC_MAX_PREC(rop); + + mpfr_init2 (u, 2); + mpfr_init2 (v, MPFR_PREC(a) + MPFR_PREC(d)); + mpfr_init2 (w, MPFR_PREC(b) + MPFR_PREC(c)); + mpfr_init2 (x, 2); + + mpfr_mul (v, a, d, GMP_RNDN); /* exact */ + if (mul_a == -1) + mpfr_neg (v, v, GMP_RNDN); + + mpfr_mul (w, b, c, GMP_RNDN); /* exact */ + if (mul_c == -1) + mpfr_neg (w, w, GMP_RNDN); + + /* compute sign(v-w) */ + sign_x = MPFR_CMP_ABS (v, w); + if (sign_x > 0) + sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w); + else if (sign_x == 0) + sign_x = mpfr_sgn (v) - mpfr_sgn (w); + else + sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w); + + sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c); + + if (sign_x * sign_u < 0) + { /* swap inputs */ + SWAP (a, c); + SWAP (b, d); + mpfr_swap (v, w); + { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; } + sign_x = - sign_x; + } + + /* now sign_x * sign_u >= 0 */ + do + { + do + { + /* the following should give failures with prob. <= 1/prec */ + prec += _mpfr_ceil_log2 ((double) prec) + 3; + + mpfr_set_prec (u, prec); + mpfr_set_prec (x, prec); + + /* first compute away(b +/- a) and store it in u */ + rnd_u = (mpfr_sgn (a) > 0) ? GMP_RNDU : GMP_RNDD; + if (mul_a == -1) + rnd_u = INV_RND(rnd_u); + inexact = ((mul_a == -1) ? mpfr_sub : mpfr_add) (u, b, a, rnd_u); + + /* then compute away(+/-c - d) and store it in x */ + rnd_x = (mpfr_sgn (c) > 0) ? GMP_RNDU : GMP_RNDD; + inexact |= ((mul_c == -1) ? mpfr_add : mpfr_sub) (x, c, d, rnd_x); + if (mul_c == -1) + mpfr_neg (x, x, GMP_RNDN); + + /* compute away(u*x) and store it in u */ + rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD; + inexact |= mpfr_mul (u, u, x, rnd_u); /* (a+b)*(c-d) */ + + inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */ + + /* in case u=0, ensure that rnd_u rounds x away from zero */ + if (mpfr_sgn (u) == 0) + rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD; + inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */ + + ok = inexact == 0 || + mpfr_can_round (u, prec - 3, rnd_u, rnd_re, prec_re); + } + while (ok == 0); + + if (mul_i == 0) + { + inex_re = mpfr_set (MPC_RE(result), u, MPC_RND_RE(rnd)); + if (inex_re == 0) + { + inex_re = inexact; /* u is rounded away from 0 */ + inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd)); + } + else if (MPC_RND_RE (rnd) == GMP_RNDN && inexact != 0 + && MPC_INEX_POS (inex_re) == MPC_INEX_POS (-MPFR_SIGN (u)) + && !(MPFR_SIGN (u) > 0 + ? mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDU, prec_re) + : mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDD, prec_re))) + /* rounding did work, but we do not know whether we are larger + or smaller than the correct result */ + { + inex_im = 0; /* just to avoid the gcc warning message */ + ok = 0; + } + else + inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd)); + } + else if (mul_i == 1) /* (x+i*y)/i = y - i*x */ + { + inex_im = mpfr_neg (MPC_IM(result), u, MPC_RND_IM(rnd)); + if (inex_im == 0) + { + inex_im = -inexact; /* u is rounded away from 0 */ + inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd)); + } + else if (MPC_RND_IM (rnd) == GMP_RNDN && inexact != 0 + && MPC_INEX_POS (inex_im) == MPC_INEX_POS (MPFR_SIGN (u)) + && !(MPFR_SIGN (u) > 0 + ? mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDU, prec_re) + : mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDD, prec_re))) + /* rounding did work, but we do not know whether we are larger + or smaller than the correct result */ + { + inex_re = 0; /* just to avoid the gcc warning message */ + ok = 0; + } + else + inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd)); + + } + else /* mul_i = 2, z/i^2 = -z */ + { + inex_re = mpfr_neg (MPC_RE(result), u, MPC_RND_RE(rnd)); + if (inex_re == 0) + { + inex_re = -inexact; /* u is rounded away from 0 */ + inex_im = -mpfr_add (MPC_IM(result), v, w, + INV_RND(MPC_RND_IM(rnd))); + mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd)); + } + else if (MPC_RND_RE (rnd) == GMP_RNDN && inexact != 0 + && MPC_INEX_POS (inex_re) == MPC_INEX_POS (MPFR_SIGN (u)) + && !(MPFR_SIGN (u) > 0 + ? mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDU, prec_re) + : mpfr_can_round (u, prec - 3, rnd_u, GMP_RNDD, prec_re))) + /* rounding did work, but we do not know whether we are larger + or smaller than the correct result */ + { + inex_im = 0; /* just to avoid the gcc warning message */ + ok = 0; + } + else + { + inex_im = -mpfr_add (MPC_IM(result), v, w, + INV_RND(MPC_RND_IM(rnd))); + mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd)); + } + } + } + while (ok == 0); + + mpc_set (rop, result, MPC_RNDNN); + + mpfr_clear (u); + mpfr_clear (v); + mpfr_clear (w); + mpfr_clear (x); + if (overlap) + mpc_clear (result); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/mul_2exp.c b/mul_2exp.c new file mode 100644 index 0000000..c70dc5c --- /dev/null +++ b/mul_2exp.c @@ -0,0 +1,35 @@ +/* mpc_mul_2exp -- Multiply a complex number by 2^e. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_mul_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_mul_2exp (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_2exp (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/mul_fr.c b/mul_fr.c new file mode 100644 index 0000000..7fc5f78 --- /dev/null +++ b/mul_fr.c @@ -0,0 +1,35 @@ +/* mpc_mul_fr -- Multiply a complex number by a floating-point number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_mul_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_mul (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/mul_ui.c b/mul_ui.c new file mode 100644 index 0000000..7ebfa6d --- /dev/null +++ b/mul_ui.c @@ -0,0 +1,35 @@ +/* mpc_mul_ui -- Multiply a complex number by a nonnegative integer. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_mul_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_mul_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_ui (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,35 @@ +/* mpc_neg -- Negate a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_neg (mpc_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_neg (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); + inex_im = mpfr_neg (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,61 @@ +/* mpc_norm -- Square of the norm of a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +int +mpc_norm (mpfr_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + mpfr_t u, v; + mp_prec_t prec; + int inexact; + + prec = MPFR_PREC(a); + + mpfr_init (u); + mpfr_init (v); + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 3; + + mpfr_set_prec (u, prec); + mpfr_set_prec (v, prec); + + /* first compute norm(b)^2 */ + inexact = mpfr_mul (u, MPC_RE(b), MPC_RE(b), GMP_RNDN); /* err<=1/2ulp */ + inexact |= mpfr_mul (v, MPC_IM(b), MPC_IM(b), GMP_RNDN); /* err<=1/2ulp*/ + + inexact |= mpfr_add (u, u, v, GMP_RNDN); /* err <= 3/2 ulps */ + } + while (inexact != 0 && + mpfr_can_round (u, prec - 2, GMP_RNDN, rnd, MPFR_PREC(a)) == 0); + + inexact |= mpfr_set (a, u, rnd); + + mpfr_clear (u); + mpfr_clear (v); + + return inexact; +} diff --git a/out_str.c b/out_str.c new file mode 100644 index 0000000..783049a --- /dev/null +++ b/out_str.c @@ -0,0 +1,38 @@ +/* mpc_out_str -- Output a complex number on a given stream. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <ctype.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +size_t +mpc_out_str (FILE *stream, int base, size_t n, mpc_srcptr op, mp_rnd_t rnd) +{ + size_t size; + + size = mpfr_out_str (stream, base, n, MPC_RE(op), MPC_RND_RE(rnd)); + size += fprintf (stream, "+I*"); + size += mpfr_out_str (stream, base, n, MPC_IM(op), MPC_RND_IM(rnd)); + + return size; +} diff --git a/random.c b/random.c new file mode 100644 index 0000000..feb17d7 --- /dev/null +++ b/random.c @@ -0,0 +1,41 @@ +/* mpc_random -- Generate a random complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +void +mpc_random (mpc_ptr a) +{ + mpfr_set_prec (MPC_RE (a), MPFR_PREC (MPC_RE (a)) + 1); + mpfr_random (MPC_RE(a)); + mpfr_mul_2ui (MPC_RE (a), MPC_RE (a), 1, GMP_RNDN); + mpfr_sub_ui (MPC_RE (a), MPC_RE (a), 1, GMP_RNDN); + mpfr_round_prec (MPC_RE (a), GMP_RNDZ, MPFR_PREC (MPC_RE (a)) - 1); + + mpfr_set_prec (MPC_IM (a), MPFR_PREC (MPC_IM (a)) + 1); + mpfr_random (MPC_IM(a)); + mpfr_mul_2ui (MPC_IM (a), MPC_IM (a), 1, GMP_RNDN); + mpfr_sub_ui (MPC_IM (a), MPC_IM (a), 1, GMP_RNDN); + mpfr_round_prec (MPC_IM (a), GMP_RNDZ, MPFR_PREC (MPC_IM (a)) - 1); +} diff --git a/random2.c b/random2.c new file mode 100644 index 0000000..bb08f04 --- /dev/null +++ b/random2.c @@ -0,0 +1,31 @@ +/* mpc_random2 -- Generate a random complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_random2 (mpc_ptr a, mp_size_t size, mp_exp_t exp) +{ + mpfr_random2 (MPC_RE(a), size, exp); + mpfr_random2 (MPC_IM(a), size, exp); +} @@ -0,0 +1,35 @@ +/* mpc_set -- Set a complex number from another complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_set (mpc_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); + inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/set_d_d.c b/set_d_d.c new file mode 100644 index 0000000..d8f14a1 --- /dev/null +++ b/set_d_d.c @@ -0,0 +1,35 @@ +/* mpc_set_dd -- Set a complex number from two doubles. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_set_d_d (mpc_ptr a, double b, double c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set_d (MPC_RE(a), b, MPC_RND_RE(rnd)); + inex_im = mpfr_set_d (MPC_IM(a), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/set_dfl_prec.c b/set_dfl_prec.c new file mode 100644 index 0000000..920e8a1 --- /dev/null +++ b/set_dfl_prec.c @@ -0,0 +1,39 @@ +/* mpc_set_default_prec, mpc_get_default_prec -- set/get default precision + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +/* default is 53 bits */ +mp_prec_t __mpc_default_fp_bit_precision = 53; + +void +mpc_set_default_prec (mp_prec_t prec) +{ + __mpc_default_fp_bit_precision = prec; +} + +mp_prec_t +mpc_get_default_prec (void) +{ + return __mpc_default_fp_bit_precision; +} diff --git a/set_prec.c b/set_prec.c new file mode 100644 index 0000000..1d254be --- /dev/null +++ b/set_prec.c @@ -0,0 +1,31 @@ +/* mpc_set_prec -- reset the precision of a complex variable. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +void +mpc_set_prec (mpc_t x, mp_prec_t prec) +{ + mpfr_set_prec (MPC_RE(x), prec); + mpfr_set_prec (MPC_IM(x), prec); +} diff --git a/set_si_si.c b/set_si_si.c new file mode 100644 index 0000000..fa22c39 --- /dev/null +++ b/set_si_si.c @@ -0,0 +1,35 @@ +/* mpc_set_si_si -- Set a complex number from two signed long integers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_set_si_si (mpc_ptr a, long int b, long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set_si (MPC_RE(a), b, MPC_RND_RE(rnd)); + inex_im = mpfr_set_si (MPC_IM(a), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/set_ui_fr.c b/set_ui_fr.c new file mode 100644 index 0000000..7cb795c --- /dev/null +++ b/set_ui_fr.c @@ -0,0 +1,36 @@ +/* mpc_set_ui_fr -- Set a complex number from an unsigned long integer + and a floating-point number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_set_ui_fr (mpc_ptr a, unsigned long int b, mpfr_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set_ui (MPC_RE(a), b, MPC_RND_RE(rnd)); + inex_im = mpfr_set (MPC_IM(a), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/set_ui_ui.c b/set_ui_ui.c new file mode 100644 index 0000000..ab70b98 --- /dev/null +++ b/set_ui_ui.c @@ -0,0 +1,36 @@ +/* mpc_set_ui_ui -- Set a complex number from two unsigned long integers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_set_ui_ui (mpc_ptr a, unsigned long int b, unsigned long int c, + mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_set_ui (MPC_RE(a), b, MPC_RND_RE(rnd)); + inex_im = mpfr_set_ui (MPC_IM(a), c, MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,161 @@ +/* mpc_sqr -- Square a complex number. + +Copyright (C) 2002 Free Software Foundation, Inc. + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +#define INV_RND(r) \ + (((r) == GMP_RNDU) ? GMP_RNDD : (((r) == GMP_RNDD) ? GMP_RNDU : (r))) + +int +mpc_sqr (mpc_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + int ok; + mpfr_t u, v; + mpfr_t x; + /* a temporary variable to hold the real part of b, + needed in the case a==b */ + mp_prec_t prec; + int inex_re, inex_im, inexact; + + prec = MPC_MAX_PREC(a); + + /* first check for real resp. purely imaginary number */ + if (MPFR_IS_ZERO (MPC_IM(b))) + { + inex_re = mpfr_mul (MPC_RE(a), MPC_RE(b), MPC_RE(b), MPC_RND_RE(rnd)); + inex_im = mpfr_set_ui (MPC_IM(a), 0, GMP_RNDN); + return MPC_INEX(inex_re, inex_im); + } + if (MPFR_IS_ZERO (MPC_RE(b))) + { + inex_re = -mpfr_mul (MPC_RE(a), MPC_IM(b), MPC_IM(b), + INV_RND (MPC_RND_RE(rnd))); + mpfr_neg (MPC_RE(a), MPC_RE(a), GMP_RNDN); + inex_im = mpfr_set_ui (MPC_IM(a), 0, GMP_RNDN); + return MPC_INEX(inex_re, inex_im); + } + + mpfr_init (u); + mpfr_init (v); + + if (a == b) + { + mpfr_init2 (x, MPFR_PREC (b->re)); + mpfr_set (x, b->re, GMP_RNDN); + } + else + x [0] = b->re [0]; + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 5; + + mpfr_set_prec (u, prec); + mpfr_set_prec (v, prec); + + /* Let b = x + iy. We need u = x+y and v = x-y, rounded away. */ + /* As this is not implemented in mpfr, we round towards zero and */ + /* add one ulp if the result is not exact. The error is bounded */ + /* above by 1 ulp. */ + /* We first let inexact be 1 if the real part is not computed */ + /* exactly and determine the sign later. */ + inexact = 0; + if (mpfr_add (u, x, MPC_IM (b), GMP_RNDZ)) + /* we have to use x here instead of MPC_RE (b), as MPC_RE (a) + may be modified later in the attempt to assign u to it */ + { + mpfr_add_one_ulp (u, GMP_RNDN); + inexact = 1; + } + if (mpfr_sub (v, x, MPC_IM (b), GMP_RNDZ)) + { + mpfr_add_one_ulp (v, GMP_RNDN); + inexact = 1; + } + + /* compute the real part as u*v, rounded away */ + /* determine also the sign of inex_re */ + if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) + { + /* as we have rounded away, the result is exact */ + mpfr_set_ui (u, 0, GMP_RNDN); + mpfr_set_ui (MPC_RE (a), 0, GMP_RNDN); + inex_re = 0; + ok = 1; + } + else if (mpfr_sgn (u) * mpfr_sgn (v) > 0) + { + inexact |= mpfr_mul (u, u, v, GMP_RNDU); /* error 5 */ + ok = !inexact | mpfr_can_round (u, prec - 3, GMP_RNDU, + MPC_RND_RE (rnd), MPFR_PREC (MPC_RE (a))); + if (ok) + { + inex_re = mpfr_set (MPC_RE (a), u, MPC_RND_RE (rnd)); + if (inex_re == 0) + /* remember that u was already rounded */ + inex_re = inexact; + /* even if rounding did work, we might not know whether the + result is too large, too small or exact */ + else if (inexact && MPC_RND_RE (rnd) == GMP_RNDN + && inex_re < 0 + && !mpfr_can_round (u, prec - 3, GMP_RNDU, GMP_RNDU, + MPFR_PREC (MPC_RE (a)))) + ok = 0; + } + } + else + { + inexact |= mpfr_mul (u, u, v, GMP_RNDD); /* error 5 */ + ok = !inexact | mpfr_can_round (u, prec - 3, GMP_RNDD, + MPC_RND_RE (rnd), MPFR_PREC (MPC_RE (a))); + if (ok) + { + inex_re = mpfr_set (MPC_RE (a), u, MPC_RND_RE (rnd)); + if (inex_re == 0) + inex_re = inexact; + /* even if rounding did work, we might not know whether the + result is too large, too small or exact */ + else if (inexact && MPC_RND_RE (rnd) == GMP_RNDN + && inex_re > 0 + && !mpfr_can_round (u, prec - 3, GMP_RNDD, GMP_RNDD, + MPFR_PREC (MPC_RE (a)))) + ok = 0; + } + } + } + while (!ok); + + /* compute the imaginary part as 2*x*y, which is always possible */ + inex_im = mpfr_mul (MPC_IM (a), x, MPC_IM (b), + MPC_RND_IM (rnd)); + MPFR_EXP (MPC_IM (a))++; + + mpfr_clear (u); + mpfr_clear (v); + + if (a == b) + mpfr_clear (x); + + return MPC_INEX (inex_re, inex_im); +} @@ -0,0 +1,147 @@ +/* mpc_sqrt -- Take the square root of a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +int +mpc_sqrt (mpc_ptr a, mpc_srcptr b, mp_rnd_t rnd) +{ + int ok=0; + mpfr_t w, t; + mp_rnd_t rnd_w, rnd_t; + mp_prec_t prec_w, prec_t; + /* the rounding mode and the precision required for w and t, which can */ + /* be either the real or the imaginary part of a */ + int re_cmp, im_cmp; + /* comparison of the real/imaginary part of b with 0 */ + mp_prec_t prec; + int inexact; + + im_cmp = mpfr_cmp_ui (MPC_IM (b), 0); + re_cmp = mpfr_cmp_ui (MPC_RE (b), 0); + if (im_cmp == 0) + { + if (re_cmp == 0) + { + mpc_set_ui_ui (a, 0, 0, MPC_RNDNN); + return 0; + } + else if (re_cmp > 0) + { + inexact = mpfr_sqrt (MPC_RE (a), MPC_RE (b), MPC_RND_RE (rnd)); + mpfr_set_ui (MPC_IM (a), 0, GMP_RNDN); + return inexact; + } + else + { + mpfr_init2 (w, MPFR_PREC (MPC_RE (b))); + inexact = mpfr_neg (w, MPC_RE (b), GMP_RNDN); + inexact = mpfr_sqrt (MPC_IM (a), w, MPC_RND_IM (rnd)); + mpfr_set_ui (MPC_RE (a), 0, GMP_RNDN); + mpfr_clear (w); + return inexact; + } + } + + prec = MPC_MAX_PREC(a); + + mpfr_init (w); + mpfr_init (t); + + if (re_cmp >= 0) + { + rnd_w = MPC_RND_RE (rnd); + prec_w = MPFR_PREC (MPC_RE (a)); + rnd_t = MPC_RND_IM(rnd); + prec_t = MPFR_PREC (MPC_IM (a)); + } + else + { + rnd_w = MPC_RND_IM(rnd); + prec_w = MPFR_PREC (MPC_IM (a)); + rnd_t = MPC_RND_RE(rnd); + prec_t = MPFR_PREC (MPC_RE (a)); + if (im_cmp < 0) + { + if (rnd_w == GMP_RNDD) + rnd_w = GMP_RNDU; + else if (rnd_w == GMP_RNDU) + rnd_w = GMP_RNDD; + if (rnd_t == GMP_RNDD) + rnd_t = GMP_RNDU; + if (rnd_t == GMP_RNDU) + rnd_t = GMP_RNDD; + } + } + + do + { + prec += _mpfr_ceil_log2 ((double) prec) + 4; + mpfr_set_prec (w, prec); + mpfr_set_prec (t, prec); + /* let b = x + iy */ + /* w = sqrt ((|x| + sqrt (x^2 + y^2)) / 2), rounded down */ + /* total error bounded by 3 ulps */ + inexact = mpc_abs (w, b, GMP_RNDD); + if (re_cmp < 0) + inexact |= mpfr_sub (w, w, MPC_RE (b), GMP_RNDD); + else + inexact |= mpfr_add (w, w, MPC_RE (b), GMP_RNDD); + inexact |= mpfr_div_2ui (w, w, 1, GMP_RNDD); + inexact |= mpfr_sqrt (w, w, GMP_RNDD); + + ok = mpfr_can_round (w, prec - 2, GMP_RNDD, rnd_w, prec_w); + if (ok) + { + /* t = y / 2w, rounded up */ + /* total error bounded by 7 ulps */ + inexact |= mpfr_div (t, MPC_IM (b), w, GMP_RNDU); + inexact |= mpfr_div_2ui (t, t, 1, GMP_RNDU); + ok = mpfr_can_round (t, prec - 3, GMP_RNDU, rnd_t, prec_t); + } + } + while (inexact != 0 && ok == 0); + + if (re_cmp > 0) + { + inexact |= mpfr_set (MPC_RE (a), w, rnd_w); + inexact |= mpfr_set (MPC_IM (a), t, rnd_t); + } + else + { + inexact |= mpfr_set (MPC_RE(a), t, MPC_RND_RE(rnd)); + inexact |= mpfr_set (MPC_IM(a), w, MPC_RND_IM(rnd)); + if (im_cmp < 0) + { + inexact |= mpfr_neg (MPC_RE (a), MPC_RE (a), MPC_RND_RE (rnd)); + inexact |= mpfr_neg (MPC_IM (a), MPC_IM (a), MPC_RND_IM (rnd)); + } + } + + mpfr_clear (w); + mpfr_clear (t); + + return inexact; +} @@ -0,0 +1,35 @@ +/* mpc_sub -- Subtract two complex numbers. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +int +mpc_sub (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_sub (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); + inex_im = mpfr_sub (MPC_IM(a), MPC_IM(b), MPC_IM(c), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} diff --git a/sub_ui.c b/sub_ui.c new file mode 100644 index 0000000..6ecc4db --- /dev/null +++ b/sub_ui.c @@ -0,0 +1,36 @@ +/* mpc_sub_ui -- Add a complex number and an unsigned long int. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +/* return 0 iff both the real and imaginary parts are exact */ +int +mpc_sub_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mp_rnd_t rnd) +{ + int inex_re, inex_im; + + inex_re = mpfr_sub_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + + return MPC_INEX(inex_re, inex_im); +} @@ -0,0 +1,279 @@ +/* tdiv -- test file for mpc_div. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +#define OUT(x) { mpc_out_str (stdout, 2, 0, x, MPC_RNDNN); putchar ('\n'); } +#define ERR(x) { mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); \ + fprintf (stderr, "\n"); } + +#define INV_RND(r) \ + (((r) == GMP_RNDU) ? GMP_RNDD : (((r) == GMP_RNDD) ? GMP_RNDU : (r))) + +int mpc_div_ref _PROTO((mpc_ptr, mpc_srcptr, mpc_srcptr, mp_rnd_t)); + +int +mpc_div_ref (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mp_rnd_t rnd) +{ + int ok = 0, inexact_q, inex_re, inex_im = 0, cancel = 0, err, sgn; + mpfr_t u, v, q, t; + mp_prec_t prec; + mp_rnd_t rnd_re, rnd_im; + + if (mpfr_cmp_ui (MPC_IM(c), 0) == 0) /* c is real */ + { + /* compute imaginary part first in case a=b or a=c */ + inex_im = mpfr_div (MPC_IM(a), MPC_IM(b), MPC_RE(c), MPC_RND_IM(rnd)); + inex_re = mpfr_div (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); + return MPC_INEX(inex_re, inex_im); + } + + prec = MPC_MAX_PREC(a); + + /* try to guess the signs of the real and imaginary parts */ + sgn = mpfr_sgn (MPC_RE(b)) * mpfr_sgn (MPC_RE(c)) + + mpfr_sgn (MPC_IM(b)) * mpfr_sgn (MPC_IM(c)); + rnd_re = (sgn > 0) ? GMP_RNDU : GMP_RNDD; + sgn = mpfr_sgn (MPC_RE(b)) * mpfr_sgn (MPC_IM(c)) + - mpfr_sgn (MPC_IM(b)) * mpfr_sgn (MPC_RE(c)); + rnd_im = (sgn > 0) ? GMP_RNDU : GMP_RNDD; + + /* first try with only one real division */ + prec += _mpfr_ceil_log2 ((double) prec) + 3; + + mpfr_init2 (u, prec); + mpfr_init2 (v, prec); + mpfr_init2 (q, prec); + mpfr_init2 (t, prec); + + /* first compute 1/norm(c)^2 rounded away */ + inexact_q = mpfr_mul (u, MPC_RE(c), MPC_RE(c), GMP_RNDZ); + inexact_q |= mpfr_mul (v, MPC_IM(c), MPC_IM(c), GMP_RNDZ); + inexact_q |= mpfr_add (q, u, v, GMP_RNDZ); /* 3 ulps */ + inexact_q |= mpfr_ui_div (q, 1, q, GMP_RNDU); /* 7 ulps */ + + /* now compute b*conjugate(c) */ + + /* real part */ + inex_re = mpfr_mul (u, MPC_RE(b), MPC_RE(c), rnd_re); /* 1 */ + inex_re |= mpfr_mul (v, MPC_IM(b), MPC_IM(c), rnd_re); /* 1 */ + inex_re |= mpfr_add (t, u, v, rnd_re); + cancel = MAX(MPFR_EXP(u), MPFR_EXP(v)) - MPFR_EXP(t); + /* err(t) <= [1 + 2*2^cancel] ulp(t) */ + inex_re |= mpfr_mul (t, t, q, rnd_re) || inexact_q; + /* err(t) <= [1 + 3*(1 + 2*2^cancel) + 14] ulp(t) + */ + err = (cancel <= 1) ? 5 : ((cancel == 2) ? 6 : + ((cancel <= 4) ? 7 : cancel + 3)); + ok = (inex_re == 0) || ((err < prec) && mpfr_can_round (t, prec - err, + rnd_re, MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)))); + if ((rnd_re == GMP_RNDU && mpfr_sgn (t) < 0) + || (rnd_re == GMP_RNDD && mpfr_sgn (t) > 0)) + { + ok = 0; + rnd_re = INV_RND(rnd_re); + } + + if (ok) /* compute imaginary part */ + { + inex_im = mpfr_mul (u, MPC_RE(b), MPC_IM(c), INV_RND(rnd_im)); + inex_im |= mpfr_mul (v, MPC_IM(b), MPC_RE(c), rnd_im); + cancel = MAX(MPFR_EXP(u), MPFR_EXP(v)); + inex_im |= mpfr_sub (u, v, u, rnd_im); + cancel = cancel - MPFR_EXP(u); + inex_im |= mpfr_mul (u, u, q, rnd_im) || inexact_q; + err = (cancel <= 1) ? 5 : ((cancel == 2) ? 6 : + ((cancel <= 4) ? 7 : cancel + 3)); + ok = (inex_im == 0) || ((err < prec) && mpfr_can_round (u, + prec - err, GMP_RNDN, MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)))); + if ((rnd_im == GMP_RNDU && mpfr_sgn (u) < 0) + || (rnd_im == GMP_RNDD && mpfr_sgn (u) > 0)) + { + ok = 0; + rnd_im = INV_RND(rnd_im); + } + } + + /* now loop with two real divisions, to detect exact quotients */ + while (ok == 0) + { + prec += _mpfr_ceil_log2 ((double) prec) + 3; + + mpfr_set_prec (u, prec); + mpfr_set_prec (v, prec); + mpfr_set_prec (q, prec); + mpfr_set_prec (t, prec); + + /* first compute 1/norm(c)^2 rounded away */ + inexact_q = mpfr_mul (u, MPC_RE(c), MPC_RE(c), GMP_RNDZ); + inexact_q |= mpfr_mul (v, MPC_IM(c), MPC_IM(c), GMP_RNDZ); + inexact_q |= mpfr_add (q, u, v, GMP_RNDZ); /* 3 ulps */ + + /* now compute b*conjugate(c) */ + + /* real part */ + inex_re = mpfr_mul (u, MPC_RE(b), MPC_RE(c), rnd_re); /* 1 */ + inex_re |= mpfr_mul (v, MPC_IM(b), MPC_IM(c), rnd_re); /* 1 */ + inex_re |= mpfr_add (t, u, v, rnd_re); + cancel = MAX(MPFR_EXP(u), MPFR_EXP(v)) - MPFR_EXP(t); + /* err(t) <= [1 + 2*2^cancel] ulp(t) */ + inex_re |= mpfr_div (t, t, q, rnd_re) || inexact_q; + /* err(t) <= [1 + 2*(1 + 2*2^cancel) + 6] ulp(t) + */ + err = (cancel <= 0) ? 4 : cancel + 3; + ok = (inex_re == 0) || ((err < prec) && mpfr_can_round (t, prec - err, + rnd_re, MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)))); + if ((rnd_re == GMP_RNDU && mpfr_sgn (t) < 0) + || (rnd_re == GMP_RNDD && mpfr_sgn (t) > 0)) + { + ok = 0; + rnd_re = INV_RND(rnd_re); + } + + if (ok) /* compute imaginary part */ + { + inex_im = mpfr_mul (u, MPC_RE(b), MPC_IM(c), INV_RND(rnd_im)); + inex_im |= mpfr_mul (v, MPC_IM(b), MPC_RE(c), rnd_im); + cancel = MAX(MPFR_EXP(u), MPFR_EXP(v)); + inex_im |= mpfr_sub (u, v, u, rnd_im); + cancel = cancel - MPFR_EXP(u); + inex_im |= mpfr_div (u, u, q, rnd_im) || inexact_q; + err = (cancel <= 0) ? 4 : cancel + 3; + ok = (inex_im == 0) || ((err < prec) && mpfr_can_round (u, + prec - err, GMP_RNDN, MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)))); + if ((rnd_im == GMP_RNDU && mpfr_sgn (u) < 0) + || (rnd_im == GMP_RNDD && mpfr_sgn (u) > 0)) + { + ok = 0; + rnd_im = INV_RND(rnd_im); + } + } + } + + inex_re = mpfr_set (MPC_RE(a), t, MPC_RND_RE(rnd)) || inex_re; + inex_im = mpfr_set (MPC_IM(a), u, MPC_RND_IM(rnd)) || inex_im; + + mpfr_clear (u); + mpfr_clear (v); + mpfr_clear (q); + mpfr_clear (t); + + return MPC_INEX(inex_re, inex_im); +} + +int +main() +{ + mpc_t b, c, q, q_ref; + int inex, i; + mp_prec_t prec; + mp_rnd_t rnd_re, rnd_im, rnd; + + mpc_init (b); + mpc_init (c); + mpc_init (q); + mpc_init (q_ref); + + mpc_set_prec (b, 10); + mpc_set_prec (c, 10); + mpc_set_prec (q, 10); + + mpc_set_ui_ui (b, 973, 964, MPC_RNDNN); + mpc_set_ui_ui (c, 725, 745, MPC_RNDNN); + + inex = mpc_div (q, b, c, MPC_RNDZZ); + mpc_set_si_si (b, 43136, -787, MPC_RNDNN); + mpc_div_2exp (b, b, 15, MPC_RNDNN); + if (mpc_cmp (q, b)) + { + fprintf (stderr, "mpc_div failed for (973+I*964)/(725+I*745)\n"); + exit (1); + } + + mpc_set_si_si (b, -837, 637, MPC_RNDNN); + mpc_set_si_si (c, 63, -5, MPC_RNDNN); + inex = mpc_div (q, b, c, MPC_RNDZN); + + mpc_set_prec (b, 2); + mpc_set_prec (c, 2); + mpc_set_prec (q, 2); + + mpc_set_ui_ui (b, 4, 3, MPC_RNDNN); + mpc_set_ui_ui (c, 1, 2, MPC_RNDNN); + + inex = mpc_div (q, b, c, MPC_RNDNN); + + for (prec = 2; prec < 1000; prec++) + { + mpc_set_prec (b, prec); + mpc_set_prec (c, prec); + mpc_set_prec (q, prec); + mpc_set_prec (q_ref, prec); + + for (i = 0; i < 1000/prec; i++) + { + mpc_random (b); + /* generate a non-zero divisor */ + do + { + mpc_random (c); + } + while (mpfr_sgn (MPC_RE(c)) == 0 && mpfr_sgn (MPC_IM(c)) == 0); + + for (rnd_re = 0; rnd_re < 4; rnd_re++) + for (rnd_im = 0; rnd_im < 4; rnd_im++) + { + rnd = RNDC(rnd_re, rnd_im); + mpc_div (q, b, c, rnd); + mpc_div_ref (q_ref, b, c, rnd); + if (mpc_cmp (q, q_ref)) + { + fprintf (stderr, "mpc_div and mpc_div_ref differ" + "for prec=%lu rnd=(%s,%s)\n", prec, + mpfr_print_rnd_mode (rnd_re), + mpfr_print_rnd_mode (rnd_im)); + fprintf (stderr, "b="); + ERR(b); + fprintf (stderr, "c="); + ERR(c); + fprintf (stderr, "q="); + ERR(q); + fprintf (stderr, "q_ref="); + ERR(q_ref); + exit (1); + } + } + } + + } + + mpc_clear (b); + mpc_clear (c); + mpc_clear (q); + mpc_clear (q_ref); + + return 0; +} @@ -0,0 +1,222 @@ +/* test file for mpc. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + +#define PRINT(x) {} +#define OUT(x) { \ + mpfr_print_binary (MPC_RE(x)); \ + printf ("+I*"); \ + mpfr_print_binary (MPC_IM(x)); \ + printf ("\n"); \ +} + +int +main() +{ + mpc_t x, y, z; + mp_prec_t prec; + mpfr_t f, g; + FILE *file; + int nread; + const char *filename = "mpc_test"; + + mpc_init (x); + mpc_init2 (y, 2); + mpc_init3 (z, 2, 2); + mpfr_init (f); + mpfr_init (g); + + for (prec = 2; prec <= 1000; prec++) + { + mpc_set_prec (x, prec); + mpc_set_prec (y, prec); + mpc_set_prec (z, prec); + mpfr_set_prec (f, prec); + mpfr_set_prec (g, prec); + + mpc_set_ui (x, 1, MPC_RNDNN); + mpc_mul_2exp (x, x, prec, MPC_RNDNN); + mpc_set_ui (y, 1, MPC_RNDNN); + if (mpc_add (z, x, y, MPC_RNDNN) == 0) + { + fprintf (stderr, "Error in mpc_add: 2^(-prec)+1 cannot be exact\n"); + exit (1); + } + + mpc_random (x); + + mpc_random2 (y, 1 + (prec - 1) / mp_bits_per_limb, 0); + + PRINT ("Testing mpc_add\n"); + mpc_add (z, x, y, MPC_RNDNN); + if (mpc_add (z, z, z, MPC_RNDNN)) + { + fprintf (stderr, "Error in mpc_add: z+z should be exact\n"); + fprintf (stderr, "Precision %li\n", prec); + exit (1); + } + + PRINT ("Testing mpc_add_fr\n"); + mpfr_random (f); + mpc_add_fr (z, x, f, MPC_RNDNZ); + mpc_set_ui (z, 1, MPC_RNDNN); + mpfr_set_ui (f, 1, GMP_RNDN); + if (mpc_add_fr (z, z, f, MPC_RNDNZ)) + { + fprintf (stderr, "Error in mpc_add_fr: 1+1 should be exact\n"); + exit (1); + } + mpc_set_ui (z, 1, MPC_RNDNN); + mpc_mul_2exp (z, z, prec, MPC_RNDNN); + if (mpc_add_fr (z, z, f, MPC_RNDNN) == 0) + { + fprintf (stderr, "Error in mpc_add_fr: 2^prec+1 cannot be exact\n"); + exit (1); + } + + PRINT ("Testing mpc_add_ui\n"); + mpc_add_ui (z, x, 17, MPC_RNDNU); + mpc_set_ui (z, 1, MPC_RNDNN); + if (mpc_add_ui (z, z, 1, MPC_RNDNZ)) + { + fprintf (stderr, "Error in mpc_add_ui: 1+1 should be exact\n"); + exit (1); + } + mpc_set_ui (z, 1, MPC_RNDNN); + mpc_mul_2exp (z, z, prec, MPC_RNDNN); + if (mpc_add_ui (z, z, 1, MPC_RNDNN) == 0) + { + fprintf (stderr, "Error in mpc_add_ui: 2^prec+1 cannot be exact\n"); + exit (1); + } + + PRINT ("Testing mpc_conj\n"); + mpc_conj (z, z, MPC_RNDND); + + PRINT ("Testing mpc_div\n"); + mpc_div (z, x, y, MPC_RNDZN); + + PRINT ("Testing mpc_div_2exp\n"); + mpc_div_2exp (z, x, 1, MPC_RNDZZ); + + PRINT ("Testing mpc_div_fr\n"); + mpc_div_fr (z, x, f, MPC_RNDZU); + + PRINT ("Testing mpc_div_ui\n"); + mpc_div_ui (z, x, 17, MPC_RNDZD); + + PRINT ("Testing mpc_inp_str\n"); + if (!(file = fopen (filename, "w"))) + { + fprintf (stderr, "Could not open file %s\n", filename); + exit (1); + }; + fprintf (file, "1+I*1\n"); + fclose (file); + if (!(file = fopen (filename, "r"))) + { + fprintf (stderr, "Could not open file %s\n", filename); + exit (1); + }; + nread = mpc_inp_str (z, file, 10, MPC_RNDUZ); + fclose (file); + mpc_set_si_si (x, 1, 1, MPC_RNDNN); + if (nread && mpc_cmp (z, x)) + { + fprintf (stderr, "inp_str o out_str <> Id\n"); + mpc_out_str (stderr, 10, 0, z, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + + PRINT ("Testing mpc_mul\n"); + mpc_mul (z, x, y, MPC_RNDUU); + + PRINT ("Testing mpc_mul_2exp\n"); + mpc_mul_2exp (z, x, 17, MPC_RNDUD); + + PRINT ("Testing mpc_mul_fr\n"); + mpc_mul_fr (z, x, f, MPC_RNDDN); + + PRINT ("Testing mpc_mul_ui\n"); + mpc_mul_ui (z, x, 17, MPC_RNDDZ); + + PRINT ("Testing mpc_neg\n"); + mpc_neg (z, x, MPC_RNDDU); + + PRINT ("Testing mpc_norm\n"); + mpc_norm (f, x, MPC_RNDDD); + + PRINT ("Testing mpc_out_str\n"); + if (!(file = fopen (filename, "w"))) + { + fprintf (stderr, "Could not open file %s\n", filename); + exit (1); + }; + mpc_out_str (file, 10, 0, z, MPC_RNDNN); + fclose (file); + + PRINT ("Testing mpc_set\n"); + mpc_set (z, x, MPC_RNDNN); + + PRINT ("Testing mpc_set_d\n"); + mpc_set_d (z, 1.23456789, MPC_RNDNN); + + PRINT ("Testing mpc_set_d_d\n"); + mpc_set_d_d (z, 1.23456789, 1.23456789, MPC_RNDNN); + + PRINT ("Testing mpc_set_default_prec\n"); + mpc_set_default_prec (prec); + + PRINT ("Testing mpc_get_default_prec\n"); + if (mpc_get_default_prec () != prec) + abort (); + + PRINT ("Testing mpc_set_ui\n"); + mpc_set_ui (z, 17, MPC_RNDNN); + + PRINT ("Testing mpc_set_ui_fr\n"); + mpc_set_ui_fr (z, 17, f, MPC_RNDNN); + + PRINT ("Testing mpc_set_ui_ui\n"); + mpc_set_ui_ui (z, 17, 17, MPC_RNDNN); + + PRINT ("Testing mpc_sub\n"); + mpc_sub (z, x, y, MPC_RNDNN); + + PRINT ("Testing mpc_ui_div\n"); + mpc_ui_div (z, 17, x, MPC_RNDNN); + } + + + mpc_clear (x); + mpc_clear (y); + mpc_clear (z); + mpfr_clear (f); + mpfr_clear (g); + + return 0; +} @@ -0,0 +1,104 @@ +/* test file for mpc_exp. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" + + +int +main() +{ + mpc_t x, z, t, u; + mpfr_t f, g; + mp_prec_t prec; + + mpc_init (x); + mpc_init (z); + mpc_init (t); + mpc_init (u); + mpfr_init (f); + mpfr_init (g); + + for (prec = 2; prec <= 1000; prec++) + { + mpc_set_prec (x, prec); + mpc_set_prec (z, prec); + mpc_set_prec (t, prec); + mpc_set_prec (u, 4*prec); + mpfr_set_prec (f, prec); + mpfr_set_prec (g, prec); + + /* check that exp(I*b) = cos(b) + I*sin(b) */ + mpfr_set_ui (MPC_RE (x), 0, GMP_RNDN); + mpfr_random (MPC_IM (x)); + + mpc_exp (z, x, MPC_RNDNN); + mpfr_sin_cos (f, g, MPC_IM(x), GMP_RNDN); + if (mpfr_cmp (g, MPC_RE(z)) || mpfr_cmp (f, MPC_IM(z))) + { + fprintf (stderr, "Error in mpc_exp: exp(I*x) <> cos(x)+I*sin(x)\n" + "got "); + mpc_out_str (stderr, 10, 0, z, MPC_RNDNN); + fprintf (stderr, "\nexpected "); + mpfr_set (MPC_RE(z), g, GMP_RNDN); + mpfr_set (MPC_IM(z), f, GMP_RNDN); + mpc_out_str (stderr, 10, 0, z, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + } + + + /* We also compute the result with four times the precision and check */ + /* whether the rounding is correct. Error reports in this part of the */ + /* algorithm might still be wrong, though, since there are two */ + /* consecutive roundings. */ + mpc_random (x); + mpc_exp (z, x, MPC_RNDNN); + mpc_exp (t, x, MPC_RNDNN); + mpc_set (u, t, MPC_RNDNN); + + if (mpc_cmp (z, t)) + { + fprintf (stderr, "rounding in exp might be incorrect for\nx="); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nmpc_exp gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_exp quadruple precision gives "); + mpc_out_str (stderr, 2, 0, u, MPC_RNDNN); + fprintf (stderr, "\nand is rounded to "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_clear (x); + mpc_clear (z); + mpc_clear (t); + mpc_clear (u); + mpfr_clear (f); + mpfr_clear (g); + + return 0; +} @@ -0,0 +1,334 @@ +/* tmul -- test file for mpc_mul. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> +#include <sys/times.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +void cmpmul _PROTO((mpc_srcptr, mpc_srcptr, mp_rnd_t)); +void testmul _PROTO((long, long, long, long, mp_prec_t, mp_rnd_t)); +void special _PROTO((void)); +void timemul _PROTO((void)); + + +void cmpmul (mpc_srcptr x, mpc_srcptr y, mp_rnd_t rnd) + /* computes the product of x and y with the naive and Karatsuba methods */ + /* using the rounding mode rnd and compares the results and return */ + /* values. */ + /* In our current test suite, the real and imaginary parts of x and y */ + /* all have the same precision, and we use this precision also for the */ + /* result. */ + /* Furthermore, we check whether the multiplication with one of the */ + /* input arguments being also the output variable yields the same */ + /* result. */ + /* We also compute the result with four times the precision and check */ + /* whether the rounding is correct. Error reports in this part of the */ + /* algorithm might still be wrong, though, since there are two */ + /* consecutive roundings. */ +{ + mpc_t z, t, u; + int inexact_z, inexact_t; + + mpc_init2 (z, MPC_MAX_PREC (x)); + mpc_init2 (t, MPC_MAX_PREC (x)); + mpc_init2 (u, 4 * MPC_MAX_PREC (x)); + + inexact_z = mpc_mul_naive (z, x, y, rnd); + inexact_t = mpc_mul_karatsuba (t, x, y, rnd); + + if (mpc_cmp (z, t)) + { + fprintf (stderr, "mul and mul2 differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul_naive gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul_karatsuba gives "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + if (inexact_z != inexact_t) + { + fprintf (stderr, "The return values of mul and mul2 differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nand x*y="); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul_naive gives %i", inexact_z); + fprintf (stderr, "\nmpc_mul_karatsuba gives %i", inexact_t); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_set (t, x, MPC_RNDNN); + inexact_t = mpc_mul (t, t, y, rnd); + if (mpc_cmp (z, t)) + { + fprintf (stderr, "mul and mul with the first variable in place differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul in place gives "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + if (inexact_z != inexact_t) + { + fprintf (stderr, "The return values of mul and mul with the first variable in place differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nand x*y="); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives %i", inexact_z); + fprintf (stderr, "\nmpc_mul in place gives %i", inexact_t); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_set (t, y, MPC_RNDNN); + inexact_t = mpc_mul (t, x, t, rnd); + if (mpc_cmp (z, t)) + { + fprintf (stderr, "mul and mul with the second variable in place differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul in place gives "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + if (inexact_z != inexact_t) + { + fprintf (stderr, "The return values of mul and mul with the second variable in place differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nand x*y="); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives %i", inexact_z); + fprintf (stderr, "\nmpc_mul in place gives %i", inexact_t); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_mul (u, x, y, rnd); + mpc_set (t, u, rnd); + if (mpc_cmp (z, t)) + { + fprintf (stderr, "rounding in mul might be incorrect for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nand y="); + mpc_out_str (stderr, 2, 0, y, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul quadruple precision gives "); + mpc_out_str (stderr, 2, 0, u, MPC_RNDNN); + fprintf (stderr, "\nand is rounded to "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_clear (z); + mpc_clear (t); + mpc_clear (u); +} + + +void +testmul (long a, long b, long c, long d, mp_prec_t prec, mp_rnd_t rnd) +{ + mpc_t x, y; + + mpc_init2 (x, prec); + mpc_init2 (y, prec); + + mpc_set_si_si (x, a, b, rnd); + mpc_set_si_si (y, c, d, rnd); + + cmpmul (x, y, rnd); + + mpc_clear (x); + mpc_clear (y); +} + + +void +special () +{ + mpc_t x, y, z, t; + int inexact; + + mpc_init (x); + mpc_init (y); + mpc_init (z); + mpc_init (t); + + mpc_set_prec (x, 8); + mpc_set_prec (y, 8); + mpc_set_prec (z, 8); + mpc_set_si_si (x, 4, 3, MPC_RNDNN); + mpc_set_si_si (y, 1, -2, MPC_RNDNN); + inexact = mpc_mul (z, x, y, MPC_RNDNN); + if (MPC_INEX_RE(inexact) || MPC_INEX_IM(inexact)) + { + fprintf (stderr, "Error: (4+3*I)*(1-2*I) should be exact with prec=8\n"); + exit (1); + } + + mpc_set_prec (x, 27); + mpc_set_prec (y, 27); + mpfr_set_str_raw (MPC_RE(x), "1.11111011011000010101000000e-2"); + mpfr_set_str_raw (MPC_IM(x), "1.11010001010110111001110001e-3"); + mpfr_set_str_raw (MPC_RE(y), "1.10100101110110011011100100e-1"); + mpfr_set_str_raw (MPC_IM(y), "1.10111100011000001100110011e-1"); + cmpmul (x, y, 0); + + mpc_clear (x); + mpc_clear (y); +} + + +void +timemul () +{ + /* measures the time needed with different precisions for naive and */ + /* Karatsuba multiplication */ + + mpc_t x, y, z; + unsigned long int i, j; + const unsigned long int tests = 10000; + struct tms time_old, time_new; + double passed1, passed2; + + mpc_init (x); + mpc_init (y); + mpc_init_set_ui_ui (z, 1, 0, MPC_RNDNN); + + for (i = 1; i < 50; i++) + { + mpc_set_prec (x, i * BITS_PER_MP_LIMB); + mpc_set_prec (y, i * BITS_PER_MP_LIMB); + mpc_set_prec (z, i * BITS_PER_MP_LIMB); + mpc_random (x); + mpc_random (y); + + times (&time_old); + for (j = 0; j < tests; j++) + mpc_mul_naive (z, x, y, MPC_RNDNN); + times (&time_new); + passed1 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100; + + times (&time_old); + for (j = 0; j < tests; j++) + mpc_mul_karatsuba (z, x, y, MPC_RNDNN); + times (&time_new); + passed2 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100; + + printf ("Time for %3li limbs naive/Karatsuba: %5.2f %5.2f\n", i, + passed1, passed2); + } + + mpc_clear (x); + mpc_clear (y); + mpc_clear (z); +} + + +int +main() +{ + mpc_t x, y; + mp_rnd_t rnd_re, rnd_im; + mp_prec_t prec; + int i; + +/* + timemul (); +*/ + + special (); + + testmul (247, -65, -223, 416, 8, 24); + testmul (5, -896, 5, -32, 3, 2); + testmul (-3, -512, -1, -1, 2, 16); + testmul (266013312, 121990769, 110585572, 116491059, 27, 0); + testmul (170, 9, 450, 251, 8, 0); + testmul (768, 85, 169, 440, 8, 16); + testmul (145, 1816, 848, 169, 8, 24); + testmul (0, 1816, 848, 169, 8, 24); + testmul (145, 0, 848, 169, 8, 24); + testmul (145, 1816, 0, 169, 8, 24); + testmul (145, 1816, 848, 0, 8, 24); + + mpc_init (x); + mpc_init (y); + + for (prec = 2; prec < 1000; prec++) { + + mpc_set_prec (x, prec); + mpc_set_prec (y, prec); + + for (i = 0; i < 1000/prec; i++) + { + + mpc_random (x); + mpc_random (y); + + for (rnd_re = 0; rnd_re < 4; rnd_re ++) + for (rnd_im = 0; rnd_im < 4; rnd_im ++) + cmpmul (x, y, RNDC(rnd_re, rnd_im)); + } + } + + mpc_clear (x); + mpc_clear (y); + + return 0; +} @@ -0,0 +1,220 @@ +/* tsqr -- test file for mpc_sqr. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <stdlib.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +void cmpsqr _PROTO((mpc_srcptr, mp_rnd_t)); +void testsqr _PROTO((long, long, mp_prec_t, mp_rnd_t)); +void special _PROTO((void)); + + +void cmpsqr (mpc_srcptr x, mp_rnd_t rnd) + /* computes the square of x with the specific function or by simple */ + /* multiplication using the rounding mode rnd and compares the results */ + /* and return values. */ + /* In our current test suite, the real and imaginary parts of x have */ + /* the same precision, and we use this precision also for the result. */ + /* Furthermore, we check whether computing the square in the same */ + /* place yields the same result. */ + /* We also compute the result with four times the precision and check */ + /* whether the rounding is correct. Error reports in this part of the */ + /* algorithm might still be wrong, though, since there are two */ + /* consecutive roundings. */ +{ + mpc_t z, t, u; + int inexact_z, inexact_t; + + mpc_init2 (z, MPC_MAX_PREC (x)); + mpc_init2 (t, MPC_MAX_PREC (x)); + mpc_init2 (u, 4 * MPC_MAX_PREC (x)); + + inexact_z = mpc_sqr (z, x, rnd); + inexact_t = mpc_mul (t, x, x, rnd); + + if (mpc_cmp (z, t)) + { + fprintf (stderr, "sqr and mul differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_mul gives "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + if (inexact_z != inexact_t) + { + fprintf (stderr, "The return values of sqr and mul differ for rnd=(%s,%s) \nx= ", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nx^2="); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr gives %i", inexact_z); + fprintf (stderr, "\nmpc_mul gives %i", inexact_t); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_set (t, x, MPC_RNDNN); + inexact_t = mpc_sqr (t, t, rnd); + if (mpc_cmp (z, t)) + { + fprintf (stderr, "sqr and sqr in place differ for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr in place gives "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + if (inexact_z != inexact_t) + { + fprintf (stderr, "The return values of sqr and sqr in place differ for rnd=(%s,%s) \nx= ", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nx^2="); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr gives %i", inexact_z); + fprintf (stderr, "\nmpc_sqr in place gives %i", inexact_t); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_sqr (u, x, rnd); + mpc_set (t, u, rnd); + if (mpc_cmp (z, t)) + { + fprintf (stderr, "rounding in sqr might be incorrect for rnd=(%s,%s) \nx=", + mpfr_print_rnd_mode(MPC_RND_RE(rnd)), + mpfr_print_rnd_mode(MPC_RND_IM(rnd))); + mpc_out_str (stderr, 2, 0, x, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr gives "); + mpc_out_str (stderr, 2, 0, z, MPC_RNDNN); + fprintf (stderr, "\nmpc_sqr quadruple precision gives "); + mpc_out_str (stderr, 2, 0, u, MPC_RNDNN); + fprintf (stderr, "\nand is rounded to "); + mpc_out_str (stderr, 2, 0, t, MPC_RNDNN); + fprintf (stderr, "\n"); + exit (1); + } + + mpc_clear (z); + mpc_clear (t); +} + + +void +testsqr (long a, long b, mp_prec_t prec, mp_rnd_t rnd) +{ + mpc_t x; + + mpc_init2 (x, prec); + + mpc_set_si_si (x, a, b, rnd); + + cmpsqr (x, rnd); + + mpc_clear (x); +} + + +void +special () +{ + mpc_t x, z; + int inexact; + + mpc_init (x); + mpc_init (z); + + mpc_set_prec (x, 8); + mpc_set_prec (z, 8); + mpc_set_si_si (x, 4, 3, MPC_RNDNN); + inexact = mpc_sqr (z, x, MPC_RNDNN); + if (MPC_INEX_RE(inexact) || MPC_INEX_IM(inexact)) + { + fprintf (stderr, "Error: (4+3*I)^2 should be exact with prec=8\n"); + exit (1); + } + + mpc_set_prec (x, 27); + mpfr_set_str_raw (MPC_RE(x), "1.11111011011000010101000000e-2"); + mpfr_set_str_raw (MPC_IM(x), "1.11010001010110111001110001e-3"); + + cmpsqr (x, 0); + + mpc_clear (x); + mpc_clear (z); +} + + +int +main() +{ + mpc_t x; + mp_rnd_t rnd_re, rnd_im; + mp_prec_t prec; + int i; + + special (); + + testsqr (247, -65, 8, 24); + testsqr (5, -896, 3, 2); + testsqr (-3, -512, 2, 16); + testsqr (266013312, 121990769, 27, 0); + testsqr (170, 9, 8, 0); + testsqr (768, 85, 8, 16); + testsqr (145, 1816, 8, 24); + testsqr (0, 1816, 8, 24); + testsqr (145, 0, 8, 24); + + mpc_init (x); + + for (prec = 2; prec < 1000; prec++) + { + mpc_set_prec (x, prec); + + for (i = 0; i < 1000/prec; i++) + { + mpc_random (x); + + for (rnd_re = 0; rnd_re < 4; rnd_re ++) + for (rnd_im = 0; rnd_im < 4; rnd_im ++) + cmpsqr (x, RNDC(rnd_re, rnd_im)); + } + } + + mpc_clear (x); + + return 0; +} diff --git a/ui_div.c b/ui_div.c new file mode 100644 index 0000000..7949314 --- /dev/null +++ b/ui_div.c @@ -0,0 +1,50 @@ +/* mpc_ui_div -- Divide an unsigned long int by a complex number. + +Copyright (C) 2002 Andreas Enge, Paul Zimmermann + +This file is part of the MPC Library. + +The MPC Library 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 2.1 of the License, or (at your +option) any later version. + +The MPC Library 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 the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <limits.h> +#include "gmp.h" +#include "mpfr.h" +#include "mpc.h" +#include "mpc-impl.h" + +int +mpc_ui_div (mpc_ptr a, unsigned long int b, mpc_srcptr c, mp_rnd_t rnd) +{ + int inex; + mpc_t bb; + + if (mpfr_cmp_ui (MPC_IM(c), 0) == 0) /* c is real */ + { + inex = mpfr_ui_div (MPC_RE(a), b, MPC_RE(c), MPC_RND_RE(rnd)); + mpfr_set_ui (MPC_IM(a), 0, MPC_RND_IM(rnd)); + return MPC_INEX(inex, 0); + } + + mpc_init2 (bb, sizeof(unsigned long int) * CHAR_BIT); + + mpc_set_ui (bb, b, rnd); /* exact */ + + inex = mpc_div (a, bb, c, rnd); + + mpc_clear (bb); + + return inex; +} |