summaryrefslogtreecommitdiff
path: root/mpfr
diff options
context:
space:
mode:
authorKevin Ryde <user42@zip.com.au>2001-02-01 23:31:48 +0100
committerKevin Ryde <user42@zip.com.au>2001-02-01 23:31:48 +0100
commit4b133dd08ea286095823032802a0fe78b15e2715 (patch)
tree5a645599aa66bd9f76c9c2cec26cefb63c3c1f3f /mpfr
parent618be102b4c29f1e02e4a7d2244ae57b4bb544c9 (diff)
downloadgmp-4b133dd08ea286095823032802a0fe78b15e2715.tar.gz
* mpfr/*: Update to mpfr 2001.
Diffstat (limited to 'mpfr')
-rw-r--r--mpfr/AUTHORS17
-rw-r--r--mpfr/COPYING481
-rw-r--r--mpfr/NEWS36
-rw-r--r--mpfr/TODO18
-rw-r--r--mpfr/karadiv.c137
-rw-r--r--mpfr/karasqrt.c102
-rw-r--r--mpfr/tests/reuse.c386
-rw-r--r--mpfr/tests/tabs.c89
-rw-r--r--mpfr/tests/tadd.c532
-rw-r--r--mpfr/tests/tagm.c67
-rw-r--r--mpfr/tests/tcan_round.c12
-rw-r--r--mpfr/tests/tcmp.c108
-rw-r--r--mpfr/tests/tcmp2.c62
-rw-r--r--mpfr/tests/tcmp_ui.c36
-rw-r--r--mpfr/tests/tdiv.c204
-rw-r--r--mpfr/tests/tdiv_ui.c54
-rw-r--r--mpfr/tests/tdump.c (renamed from mpfr/set_dfl_rnd.c)27
-rw-r--r--mpfr/tests/teq.c76
-rw-r--r--mpfr/tests/texp.c171
-rw-r--r--mpfr/tests/tget_str.c86
-rw-r--r--mpfr/tests/tlog.c75
-rw-r--r--mpfr/tests/tlog2.c10
-rw-r--r--mpfr/tests/tmul.c176
-rw-r--r--mpfr/tests/tmul_2exp.c46
-rw-r--r--mpfr/tests/tmul_ui.c56
-rw-r--r--mpfr/tests/tout_str.c69
-rw-r--r--mpfr/tests/tpi.c10
-rw-r--r--mpfr/tests/tpow.c66
-rw-r--r--mpfr/tests/trandom.c182
-rw-r--r--mpfr/tests/tround.c11
-rw-r--r--mpfr/tests/tset_d.c23
-rw-r--r--mpfr/tests/tset_f.c35
-rw-r--r--mpfr/tests/tset_q.c78
-rw-r--r--mpfr/tests/tset_si.c49
-rw-r--r--mpfr/tests/tset_str.c80
-rw-r--r--mpfr/tests/tset_z.c14
-rw-r--r--mpfr/tests/tsin_cos.c101
-rw-r--r--mpfr/tests/tsqrt.c259
-rw-r--r--mpfr/tests/tsqrt_ui.c79
-rw-r--r--mpfr/tests/tswap.c42
-rw-r--r--mpfr/tests/ttrunc.c144
-rw-r--r--mpfr/tests/tui_div.c98
-rw-r--r--mpfr/tests/tui_sub.c99
-rw-r--r--mpfr/tests/tzeta.c70
-rw-r--r--mpfr/zeta.c116
45 files changed, 3560 insertions, 1129 deletions
diff --git a/mpfr/AUTHORS b/mpfr/AUTHORS
new file mode 100644
index 000000000..3c5c3016b
--- /dev/null
+++ b/mpfr/AUTHORS
@@ -0,0 +1,17 @@
+The MPFR library was developed by the PolKA project at INRIA Lorraine and
+LORIA. Authors and contributors are:
+
+Sylvie Boldo <sboldo@ens-lyon.fr> wrote the agm and log routines
+Guillaume Hanrot <hanrot@loria.fr> wrote part of the library
+Emmanuel Jeandel <E_Jeandel@mail.dotcom.fr> wrote the hypergeometric code
+Vincent Lefe`vre <lefevre@loria.fr> contributed many test cases for elementary
+ functions
+Thom Mulders <mulders@inf.ethz.ch> contributed the short multiplication code
+Fabrice Rouillier <rouillie@loria.fr> helped a lot for testing and porting
+Jean-Luc Re'my <remy@loria.fr> wrote a prototype of the Zeta function
+Paul Zimmermann <zimmerma@loria.fr> wrote part of the library
+
+All the authors are included in the MPFR mailing-list <mpfr@loria.fr>.
+This is the preferred way to contact us. For further information, please
+look at the MPFR web page <http://www.loria.fr/projets/mpfr/>.
+
diff --git a/mpfr/COPYING b/mpfr/COPYING
new file mode 100644
index 000000000..92b8903ff
--- /dev/null
+++ b/mpfr/COPYING
@@ -0,0 +1,481 @@
+ GNU LIBRARY GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1991 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 library GPL. It is
+ numbered 2 because it goes with version 2 of the ordinary GPL.]
+
+ 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 Library General Public License, applies to some
+specially designated Free Software Foundation software, and to any
+other libraries whose authors decide to use it. You can use it for
+your libraries, too.
+
+ When we speak of free software, we are referring to freedom, 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 or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the 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 a program 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.
+
+ Our method of protecting your rights has two steps: (1) copyright
+the library, and (2) offer you this license which gives you legal
+permission to copy, distribute and/or modify the library.
+
+ Also, for each distributor's protection, we want to make certain
+that everyone understands that there is no warranty for this free
+library. If the library is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original
+version, so that any problems introduced by others will not reflect on
+the original authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that companies distributing free
+software will individually obtain patent licenses, thus in effect
+transforming the program into proprietary software. To prevent this,
+we have made it clear that any patent must be licensed for everyone's
+free use or not licensed at all.
+
+ Most GNU software, including some libraries, is covered by the ordinary
+GNU General Public License, which was designed for utility programs. This
+license, the GNU Library General Public License, applies to certain
+designated libraries. This license is quite different from the ordinary
+one; be sure to read it in full, and don't assume that anything in it is
+the same as in the ordinary license.
+
+ The reason we have a separate public license for some libraries is that
+they blur the distinction we usually make between modifying or adding to a
+program and simply using it. Linking a program with a library, without
+changing the library, is in some sense simply using the library, and is
+analogous to running a utility program or application program. However, in
+a textual and legal sense, the linked executable is a combined work, a
+derivative of the original library, and the ordinary General Public License
+treats it as such.
+
+ Because of this blurred distinction, using the ordinary General
+Public License for libraries did not effectively promote software
+sharing, because most developers did not use the libraries. We
+concluded that weaker conditions might promote sharing better.
+
+ However, unrestricted linking of non-free programs would deprive the
+users of those programs of all benefit from the free status of the
+libraries themselves. This Library General Public License is intended to
+permit developers of non-free programs to use free libraries, while
+preserving your freedom as a user of such programs to change the free
+libraries that are incorporated in them. (We have not seen how to achieve
+this as regards changes in header files, but we have achieved it as regards
+changes in the actual functions of the Library.) The hope is that this
+will lead to faster development of free libraries.
+
+ 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, while the latter only
+works together with the library.
+
+ Note that it is possible for a library to be covered by the ordinary
+General Public License rather than by this special one.
+
+ GNU LIBRARY GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License Agreement applies to any software library which
+contains a notice placed by the copyright holder or other authorized
+party saying it may be distributed under the terms of this Library
+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 compile 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) 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.
+
+ c) 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.
+
+ d) 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 source code 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 to
+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 Library 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 Library General Public
+ License as published by the Free Software Foundation; either
+ version 2 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
+ Library General Public License for more details.
+
+ You should have received a copy of the GNU Library General Public
+ License along with this library; if not, write to the Free
+ 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!
diff --git a/mpfr/NEWS b/mpfr/NEWS
new file mode 100644
index 000000000..1792ad624
--- /dev/null
+++ b/mpfr/NEWS
@@ -0,0 +1,36 @@
+Changes from version 1.0 to version 2001:
+- the default installation does not provide any more access to machine
+ rounding mode, and as a consequence does not compare MPFR results with
+ precision=53 to machine results. Add option -DTEST if you want to have
+ access to machine rounding mode, and to check MPFR results against.
+- the MPFR files do not need <math.h> any more
+- the header file <mpfr.h> was split into <mpfr.h> for exported functions
+ and <mpfr-impl.h> for internal functions. The user should not use functions
+ or macros from <mpfr-impl.h>, since those may change in further releases.
+- <mpfr.h> was modified in order to make easy a C++ interface
+- MPFR now deals with infinities (+infinity and -infinity) and NaN
+- the missing function mpfr_swap is now available
+- mpfr_zeta was removed (was incomplete)
+- mpfr_init and mpfr_init2 now initialize the corresponding variable to 0
+ (like in other initialization functions from GNU MP)
+- in case memory allocation fails, an error message is output
+- several bugs of version 0.4 were fixed
+
+Changes from version 0.4 to version 1.0:
+
+- Version 1.0 now uses a standard configure/make installation.
+- Version 1.0 implements all functions that are available in the MPF class
+ from GMP 3.1 (except mpf_swap) and a header file mpf2mpfr.h is included in
+ the distribution for easy change from MPF to MPFR.
+- Version 1.0 implements new elementary functions: mpfr_sincos
+- Some functions and macros have been renamed: mpfr_log2 is now
+ mpfr_const_log2, mpfr_pi is now mpfr_const_pi, SIGN is now MPFR_SIGN.
+- Version 1.0 uses faster algorithms for mpfr_exp, mpfr_const_pi,
+ mpfr_const_log2. Compare the timings from version 1.0 and version 0.4.
+- Version 1.0 corrects some bugs of version 0.4.
+- The precision of MPFR variables is now named mpfr_prec, which makes it
+ easier to change it, to say unsigned long long. Same for the rounding mode
+ which is called mp_rnd_t.
+
+You'll find other news concerning the MPFR library on the web
+page <http://www.loria.fr/projets/mpfr/>.
diff --git a/mpfr/TODO b/mpfr/TODO
new file mode 100644
index 000000000..4c62eaf4c
--- /dev/null
+++ b/mpfr/TODO
@@ -0,0 +1,18 @@
+New functions to implement:
+
+- mpfr_mul_mpz, mpfr_div_mpz
+
+Miscellaneous:
+
+- implement a C++ wrapper [see with Torbjo"rn and/or Fabrice]
+
+- in all functions that do rounding, return an int indicating if the result
+ is exact (i.e. no bit was lost) or not (like in mpfr_div_ui)
+
+- detect overflow/underflow in exponent (from Ben Hinkle <bhinkle4@juno.com>)
+
+- specify exponent size (suggestion from Ben Hinkle <bhinkle4@juno.com>)
+
+- add mpfr_get_ld for 'long double' [asked by J-C Fauge`re] ?
+ (exists since K&R, but has not necessarily a greater precision than double)
+ cf http://anubis.dkuug.dk/jtc1/sc22/wg14/www/docs/n869/
diff --git a/mpfr/karadiv.c b/mpfr/karadiv.c
deleted file mode 100644
index 94db05353..000000000
--- a/mpfr/karadiv.c
+++ /dev/null
@@ -1,137 +0,0 @@
-/* mpn_divrem_n -- Karatsuba division in 2*K(n) limb operations
-
-Copyright (C) 1999-2000 PolKA project, Inria Lorraine and Loria
-
-This file is part of the MPFR Library.
-
-The MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Library General Public License as published by
-the Free Software Foundation; either version 2 of the License, or (at your
-option) any later version.
-
-The MPFR 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 Library General Public
-License for more details.
-
-You should have received a copy of the GNU Library General Public License
-along with the MPFR 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. */
-
-/*
-[1] Fast Recursive Division, by Christoph Burnikel and Joachim Ziegler,
- Technical report MPI-I-98-1-022, october 1998,
- cf http://www.mpi-sb.mpg.de/~ziegler/TechRep.ps.gz.
- Implemented by Paul Zimmermann, 1999.
-*/
-
-#include "gmp.h"
-#include "gmp-impl.h"
-#include "mpfr.h"
-
-extern void mpn_divrem_n2 (mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t, mp_limb_t *);
-extern void mpn_divrem_3by2 (mp_limb_t *, mp_limb_t *, mp_limb_t *, mp_size_t, mp_limb_t *);
-
-/* mpn_divrem_n(n) calls 2*mul(n/2)+2*div(n/2), thus to be faster then
- div(n)=4*div(n/2), we need mul(n/2) to be faster than the classic way,
- i.e. n/2 >= KARATSUBA_MUL_THRESHOLD */
-
-#define DIV_LIMIT (7*KARATSUBA_MUL_THRESHOLD)
-
-static void mpn_decr(mp_limb_t *Q)
-{
- while ((*Q++)-- == 0);
-}
-
-/* implements algorithm of page 8 in [1]: divides (A,2n) by (B,n) and puts the
- quotient in (Q,n), the remainder in (A,n).
- Returns most significant limb of the quotient, which is 0 or 1.
-*/
-mp_limb_t
-mpn_divrem_n(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B, mp_size_t n)
-{
- if (n<DIV_LIMIT) return mpn_divrem(Q, 0, A, 2*n, B, n);
- else {
- mp_limb_t cc=0;
- if (mpn_cmp(A+n, B, n)>=0) {
- cc=1;
- mpn_sub_n(A+n, A+n, B, n);
- }
- if (n%2) {
- /* divide (2n-2) most significant limbs from A by those (n-1) from B */
- mpn_divrem_n(Q+1, A+2, B+1, n-1);
- /* now (Q+1, n-1) contains the quotient of (A+2,2n-2) by (B+1,n-1)
- and (A+2, n-1) contains the remainder */
- if (mpn_sub_1(A+n, A+n, 1, mpn_submul_1(A+1, Q+1, n-1, B[0]))) {
- /* quotient two large */
- mpn_decr(Q+1);
- if (mpn_add_n(A+1, A+1, B, n)==0) {
- mpn_decr(Q+1); mpn_add_n(A+1, A+1, B, n);
- }
- }
- /* now divide (A,n+1) by (B,n) */
- mpn_divrem(Q, 0, A, n+1, B, n);
- }
- else {
- mp_limb_t *tmp; int n2=n/2;
- TMP_DECL (marker);
-
- TMP_MARK (marker);
- tmp = (mp_limb_t*) TMP_ALLOC(n*sizeof(mp_limb_t));
- mpn_divrem_3by2(Q+n2, A+n2, B, n2, tmp);
- mpn_divrem_3by2(Q, A, B, n2, tmp);
- TMP_FREE (marker);
- }
- return cc;
- }
-}
-
-/* inner procedure, with no memory allocation
- assumes mpn_cmp(A+n, B, n) < 0
-*/
-void mpn_divrem_n2(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B, mp_size_t n,
- mp_limb_t *tmp)
-{
- if (n%2) {
- /* divide (2n-2) most significant limbs from A by those (n-1) from B */
- mpn_divrem_n2(Q+1, A+2, B+1, n-1, tmp);
- /* now (Q+1, n-1) contains the quotient of (A+2,2n-2) by (B+1,n-1)
- and (A+2, n-1) contains the remainder */
- if (mpn_sub_1(A+n, A+n, 1, mpn_submul_1(A+1, Q+1, n-1, B[0]))) {
- /* quotient two large */
- mpn_decr(Q+1);
- if (mpn_add_n(A+1, A+1, B, n)==0) { /* still too large */
- mpn_decr(Q+1); mpn_add_n(A+1, A+1, B, n);
- }
- }
- /* now divide (A,n+1) by (B,n) */
- mpn_divrem(Q, 0, A, n+1, B, n);
- }
- else {
- int n2=n/2;
- mpn_divrem_3by2(Q+n2, A+n2, B, n2, tmp);
- mpn_divrem_3by2(Q, A, B, n2, tmp);
- }
-}
-
-/* divides (A,3n) by (B,2n) and puts the quotient in (Q,n),
- the remainder in (A,2n) */
-void
-mpn_divrem_3by2(mp_limb_t *Q, mp_limb_t *A, mp_limb_t *B,
- mp_size_t n, mp_limb_t *tmp)
-{
- int twon = n+n;
-
- if (n<DIV_LIMIT) mpn_divrem(Q, 0, A+n, twon, B+n, n);
- else mpn_divrem_n2(Q, A+n, B+n, n, tmp);
- /* q=(Q,n), c=(A+n,n) with the notations of [1] */
- mpn_mul_n(tmp, Q, B, n);
- if (mpn_sub_n(A, A, tmp, twon)) /* R=(A,2n) */ {
- mpn_decr(Q);
- if (mpn_add_n(A, A, B, twon)==0) { /* Q still too large */
- mpn_decr(Q); mpn_add_n(A, A, B, twon);
- }
- }
- return;
-}
diff --git a/mpfr/karasqrt.c b/mpfr/karasqrt.c
deleted file mode 100644
index 74a4bf805..000000000
--- a/mpfr/karasqrt.c
+++ /dev/null
@@ -1,102 +0,0 @@
-/* kara_sqrtrem -- Karatsuba square root
-
-Copyright (C) 1999-2000 PolKA project, Inria Lorraine and Loria
-
-This file is part of the MPFR Library.
-
-The MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Library General Public License as published by
-the Free Software Foundation; either version 2 of the License, or (at your
-option) any later version.
-
-The MPFR 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 Library General Public
-License for more details.
-
-You should have received a copy of the GNU Library General Public License
-along with the MPFR 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. */
-
-/* Reference: Karatsuba Square Root, Paul Zimmermann, Research Report 3805,
- INRIA, November 1999. */
-
-#include "gmp.h"
-#include "gmp-impl.h"
-#include "mpfr.h"
-
-#define SQRT_LIMIT KARATSUBA_MUL_THRESHOLD /* must be at least 3, should be
- near from optimal */
-
-/* n must be even */
-mp_size_t kara_sqrtrem(mp_limb_t *s, mp_limb_t *r, mp_limb_t *op, mp_size_t n)
-{
- if (n<SQRT_LIMIT) return mpn_sqrtrem(s, r, op, n);
- else {
- mp_size_t nn, rn, rrn, sn, qn; mp_limb_t *q, tmp;
- TMP_DECL (marker);
-
- TMP_MARK (marker);
- nn = n/4; /* block size 'b' corresponds to nn limbs */
- rn = kara_sqrtrem(s+nn, r+nn, op+2*nn, n-2*nn);
- /* rn <= ceil(n-2*nn, 2) + 1 <= ceil(2*nn+3, 2) + 1 <= nn+3 */
- /* to divide by 2*s', first divide by 2, to ensure the dividend is
- less than b^2 */
- sn=(n-2*nn+1)/2; /* sn >= nn */
- MPN_COPY(r, op+nn, nn); /* copy a_1 */
- tmp = mpn_rshift(r, r, nn+rn, 1);
- if (r[nn+rn-1]==0) rn--;
- q = (mp_limb_t*) TMP_ALLOC(2*(sn+1)*sizeof(mp_limb_t));
- if (nn+rn < 2*sn) MPN_ZERO(r+nn+rn, 2*sn-nn-rn);
- qn = sn; if (mpn_cmp(r+sn, s+nn, sn)>=0) {
- q[qn++]=1; mpn_sub_n(r+sn, r+sn, s+nn, sn);
- }
-#if 0
- mpn_divrem(q, 0, r, 2*sn, s+nn, sn);
-#else
- mpn_divrem_n(q, r, s+nn, sn);
-#endif
- while (qn>nn && q[qn-1]==0) qn--;
- MPN_COPY(s, q, nn);
- if (nn+rn > 2*sn) {
- tmp=mpn_add_n(s+sn, s+sn, q+sn, nn+rn-2*sn);
- if (tmp) mpn_add_1(s+nn+rn-sn, s+nn+rn-sn, (n+1)/2-nn-rn+sn, tmp);
- }
- /* multiply remainder by two and add low bit of a_1 */
- rrn = nn+sn; /* size of output remainder */
- rrn += mpn_lshift(r+nn, r, sn, 1);
- r[nn] |= (op[nn] & 1);
- sn += nn;
- if (qn>nn) {
- MPN_COPY(r, s+nn, qn-nn); /* save the qn-nn limbs from s */
- MPN_COPY(s+nn, q+nn, qn-nn); /* replace by those of q */
- }
- mpn_mul_n(q, s, s, qn);
- if (qn>nn) { /* restore the limbs from s, adding them to those of q */
- mp_limb_t cy;
-
- cy = mpn_add_n(s+nn, s+nn, r, qn-nn);
- if (qn<sn) cy = mpn_add_1(s+qn, s+qn, sn-qn, cy);
- if (cy) s[sn++]=1;
- }
- MPN_COPY(r, op, nn); /* copy a_0 */
- qn = 2*qn;
- if (qn<sn) MPN_ZERO(q+qn, sn-qn);
- if (rrn<sn) MPN_ZERO(r+rrn, sn-rrn);
- if (mpn_sub_n(r, r, q, sn) || (qn>sn)) {
- if (rrn>sn) rrn=sn;
- else {
- /* one shift and one add is faster than two add's */
- r[sn] = mpn_lshift(q, s, sn, 1) + mpn_add_n(r, r, q, sn)
- - mpn_sub_1(r, r, sn, 1) - 1;
- rrn = sn + r[sn];
- mpn_sub_1(s, s, sn, 1);
- }
- }
- else if (rrn>sn) r[sn]=1;
- TMP_FREE (marker);
- MPN_NORMALIZE(r, rrn);
- return rrn;
- }
-}
diff --git a/mpfr/tests/reuse.c b/mpfr/tests/reuse.c
new file mode 100644
index 000000000..d8214acb8
--- /dev/null
+++ b/mpfr/tests/reuse.c
@@ -0,0 +1,386 @@
+/* Test file for in-place operations.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 "mpfr-impl.h"
+
+void (*testfunc) () = NULL;
+void test3 (char *, mp_prec_t, mp_rnd_t);
+void test3a (char *, mp_prec_t, mp_rnd_t);
+void test2ui (char *, mp_prec_t, mp_rnd_t);
+void testui2 (char *, mp_prec_t, mp_rnd_t);
+void test2 (char *, mp_prec_t, mp_rnd_t);
+void test2a (char *, mp_prec_t);
+int mpfr_compare (mpfr_t, mpfr_t);
+
+/* same than mpfr_cmp, but returns 0 for both NaN's */
+int mpfr_compare (mpfr_t a, mpfr_t b)
+{
+ return (MPFR_IS_NAN(a)) ? !MPFR_IS_NAN(b) : mpfr_cmp(a, b);
+}
+
+void test3 (char *foo, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpfr_t ref1, ref2, ref3;
+ mpfr_t res1;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref2, prec);
+ mpfr_init2 (ref3, prec);
+ mpfr_init2 (res1, prec);
+
+ /* for each variable, consider each of the following 6 possibilities:
+ NaN, +Infinity, -Infinity, +0, -0 or a random number */
+ for (i=0; i<36; i++) {
+ if (i%6==0) MPFR_SET_NAN(ref2);
+ if (i%6==1) mpfr_set_d (ref2, 1.0/0.0, GMP_RNDN);
+ if (i%6==2) mpfr_set_d (ref2, -1.0/0.0, GMP_RNDN);
+ if (i%6==3) mpfr_set_d (ref2, 0.0, GMP_RNDN);
+ if (i%6==4) mpfr_set_d (ref2, -0.0, GMP_RNDN);
+ if (i%6==5) mpfr_random (ref2);
+
+ if (i/6==0) MPFR_SET_NAN(ref3);
+ if (i/6==1) mpfr_set_d (ref3, 1.0/0.0, GMP_RNDN);
+ if (i/6==2) mpfr_set_d (ref3, -1.0/0.0, GMP_RNDN);
+ if (i/6==3) mpfr_set_d (ref3, 0.0, GMP_RNDN);
+ if (i/6==4) mpfr_set_d (ref3, -0.0, GMP_RNDN);
+ if (i/6==5) mpfr_random (ref3);
+
+ /* reference call: foo(a, b, c) */
+ testfunc (ref1, ref2, ref3, rnd);
+
+ /* foo(a, a, c) */
+ mpfr_set (res1, ref2, rnd); /* exact operation */
+ testfunc (res1, res1, ref3, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, a, c) for a=%e, c=%e\n", foo,
+ mpfr_get_d (ref2), mpfr_get_d (ref3));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+
+ /* foo(a, b, a) */
+ mpfr_set (res1, ref3, rnd);
+ testfunc (res1, ref2, res1, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, b, a) for b=%e, a=%e\n", foo,
+ mpfr_get_d (ref2), mpfr_get_d (ref3));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+
+ /* foo(a, a, a) */
+ mpfr_set (ref3, ref2, rnd);
+ testfunc (ref1, ref2, ref3, rnd);
+ mpfr_set (res1, ref2, rnd);
+ testfunc (res1, res1, res1, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, a, a) for a=%e\n", foo,
+ mpfr_get_d (ref2));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref2);
+ mpfr_clear (ref3);
+ mpfr_clear (res1);
+}
+
+void test2ui (char *foo, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpfr_t ref1, ref2;
+ unsigned int ref3;
+ mp_limb_t c[1];
+ mpfr_t res1;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref2, prec);
+ mpfr_init2 (res1, prec);
+
+ /* ref2 can be NaN, +Inf, -Inf, +0, -0 or any number
+ ref3 can be 0 or any number */
+ for (i=0; i<12; i++) {
+ if (i%6==0) MPFR_SET_NAN(ref2);
+ if (i%6==1) mpfr_set_d (ref2, 1.0/0.0, GMP_RNDN);
+ if (i%6==2) mpfr_set_d (ref2, -1.0/0.0, GMP_RNDN);
+ if (i%6==3) mpfr_set_d (ref2, 0.0, GMP_RNDN);
+ if (i%6==4) mpfr_set_d (ref2, -0.0, GMP_RNDN);
+ if (i%6==5) mpfr_random (ref2);
+
+ if (i/6==0) ref3=0;
+ else {
+ mpn_random (c, 1);
+ ref3 = (unsigned int) c[0];
+ }
+
+ /* reference call: foo(a, b, c) */
+ testfunc (ref1, ref2, ref3, rnd);
+
+ /* foo(a, a, c) */
+ mpfr_set (res1, ref2, rnd); /* exact operation */
+ testfunc (res1, res1, ref3, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, a, c) for a=%e c=%u\n", foo,
+ mpfr_get_d (ref2), ref3);
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref2);
+ mpfr_clear (res1);
+}
+
+void testui2 (char *foo, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpfr_t ref1, ref3;
+ unsigned int ref2;
+ mp_limb_t c[1];
+ mpfr_t res1;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref3, prec);
+ mpfr_init2 (res1, prec);
+ mpfr_random (ref3);
+ mpn_random (c, 1);
+ ref2 = (unsigned int) c[0];
+
+ for (i=0; i<12; i++) {
+ if (i%6==0) MPFR_SET_NAN(ref3);
+ if (i%6==1) mpfr_set_d (ref3, 1.0/0.0, GMP_RNDN);
+ if (i%6==2) mpfr_set_d (ref3, -1.0/0.0, GMP_RNDN);
+ if (i%6==3) mpfr_set_d (ref3, 0.0, GMP_RNDN);
+ if (i%6==4) mpfr_set_d (ref3, -0.0, GMP_RNDN);
+ if (i%6==5) mpfr_random (ref3);
+
+ if (i/6==0) ref2=0;
+ else {
+ mpn_random (c, 1);
+ ref2 = (unsigned int) c[0];
+ }
+
+ /* reference call: foo(a, b, c) */
+ testfunc (ref1, ref2, ref3, rnd);
+
+ /* foo(a, b, a) */
+ mpfr_set (res1, ref3, rnd); /* exact operation */
+ testfunc (res1, ref2, res1, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, b, a) for b=%u a=%e\n", foo,
+ ref2, mpfr_get_d (ref3));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref3);
+ mpfr_clear (res1);
+}
+
+/* foo(mpfr_ptr, mpfr_srcptr, mp_rndt) */
+void test2 (char *foo, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpfr_t ref1, ref2;
+ mpfr_t res1;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref2, prec);
+ mpfr_init2 (res1, prec);
+ mpfr_random (ref2);
+
+ for (i=0; i<6; i++) {
+ if (i==0) MPFR_SET_NAN(ref2);
+ if (i==1) mpfr_set_d (ref2, 1.0/0.0, GMP_RNDN);
+ if (i==2) mpfr_set_d (ref2, -1.0/0.0, GMP_RNDN);
+ if (i==3) mpfr_set_d (ref2, 0.0, GMP_RNDN);
+ if (i==4) mpfr_set_d (ref2, -0.0, GMP_RNDN);
+ if (i==5) mpfr_random (ref2);
+
+ /* reference call: foo(a, b) */
+ testfunc (ref1, ref2, rnd);
+
+ /* foo(a, a) */
+ mpfr_set (res1, ref2, rnd); /* exact operation */
+ testfunc (res1, res1, rnd);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, a) for a=%e\n", foo, mpfr_get_d (ref2));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref2);
+ mpfr_clear (res1);
+}
+
+/* foo(mpfr_ptr, mpfr_srcptr) */
+void test2a (char *foo, mp_prec_t prec)
+{
+ mpfr_t ref1, ref2;
+ mpfr_t res1;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref2, prec);
+ mpfr_init2 (res1, prec);
+ mpfr_random (ref2);
+
+ for (i=0; i<6; i++) {
+ if (i==0) MPFR_SET_NAN(ref2);
+ if (i==1) mpfr_set_d (ref2, 1.0/0.0, GMP_RNDN);
+ if (i==2) mpfr_set_d (ref2, -1.0/0.0, GMP_RNDN);
+ if (i==3) mpfr_set_d (ref2, 0.0, GMP_RNDN);
+ if (i==4) mpfr_set_d (ref2, -0.0, GMP_RNDN);
+ if (i==5) mpfr_random (ref2);
+
+ /* reference call: foo(a, b) */
+ testfunc (ref1, ref2);
+
+ /* foo(a, a) */
+ mpfr_set (res1, ref2, GMP_RNDN); /* exact operation */
+ testfunc (res1, res1);
+ if (mpfr_compare (res1, ref1)) {
+ fprintf (stderr, "Error for %s(a, a) for a=%e\n", foo, mpfr_get_d (ref2));
+ fprintf (stderr, "expected %e, got %e\n", mpfr_get_d (ref1),
+ mpfr_get_d (res1));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref2);
+ mpfr_clear (res1);
+}
+
+/* one operand, two results */
+void test3a (char *foo, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpfr_t ref1, ref2, ref3;
+ mpfr_t res1, res2;
+ int i;
+
+#ifdef DEBUG
+ printf("checking %s\n", foo);
+#endif
+ mpfr_init2 (ref1, prec);
+ mpfr_init2 (ref2, prec);
+ mpfr_init2 (ref3, prec);
+ mpfr_init2 (res1, prec);
+ mpfr_init2 (res2, prec);
+ mpfr_random (ref3);
+
+ for (i=0; i<6; i++) {
+ if (i==0) MPFR_SET_NAN(ref3);
+ if (i==1) mpfr_set_d (ref3, 1.0/0.0, GMP_RNDN);
+ if (i==2) mpfr_set_d (ref3, -1.0/0.0, GMP_RNDN);
+ if (i==3) mpfr_set_d (ref3, 0.0, GMP_RNDN);
+ if (i==4) mpfr_set_d (ref3, -0.0, GMP_RNDN);
+ if (i==5) mpfr_random (ref3);
+
+ /* reference call: foo(a, b, c) */
+ testfunc (ref1, ref2, ref3, rnd);
+
+ /* foo(a, b, a) */
+ mpfr_set (res1, ref3, rnd); /* exact operation */
+ testfunc (res1, res2, res1, rnd);
+ if (mpfr_compare (res1, ref1) || mpfr_compare (res2, ref2)) {
+ fprintf (stderr, "Error for %s(a, b, a) for a=%e\n", foo, mpfr_get_d (ref3));
+ fprintf (stderr, "expected (%e,%e), got (%e,%e)\n", mpfr_get_d (ref1),
+ mpfr_get_d (ref2), mpfr_get_d (res1), mpfr_get_d (res2));
+ exit (1);
+ }
+
+ /* foo(a, b, b) */
+ mpfr_set (res2, ref3, rnd); /* exact operation */
+ testfunc (res1, res2, res2, rnd);
+ if (mpfr_compare (res1, ref1) || mpfr_compare (res2, ref2)) {
+ fprintf (stderr, "Error for %s(a, b, b) for b=%e\n", foo, mpfr_get_d (ref3));
+ fprintf (stderr, "expected (%e,%e), got (%e,%e)\n", mpfr_get_d (ref1),
+ mpfr_get_d (ref2), mpfr_get_d (res1), mpfr_get_d (res2));
+ exit (1);
+ }
+ }
+
+ mpfr_clear (ref1);
+ mpfr_clear (ref2);
+ mpfr_clear (ref3);
+ mpfr_clear (res1);
+ mpfr_clear (res2);
+}
+
+int main ()
+{
+ testfunc = mpfr_add; test3 ("mpfr_add", 53, GMP_RNDN);
+ testfunc = mpfr_add_ui; test2ui ("mpfr_add_ui", 53, GMP_RNDN);
+ testfunc = mpfr_agm; test3 ("mpfr_agm", 53, GMP_RNDN);
+ testfunc = mpfr_ceil; test2 ("mpfr_ceil", 53, GMP_RNDN);
+ testfunc = mpfr_div; test3 ("mpfr_div", 53, GMP_RNDN);
+ testfunc = mpfr_div_2exp; test2ui ("mpfr_div_2exp", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_div_ui; test2ui ("mpfr_div_ui", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_exp; test2 ("mpfr_exp", 53, GMP_RNDN);
+ testfunc = mpfr_floor; test2 ("mpfr_floor", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_log; test2 ("mpfr_log", 53, GMP_RNDN);
+ testfunc = mpfr_mul; test3 ("mpfr_mul", 53, GMP_RNDN);
+ testfunc = mpfr_mul_2exp; test2ui ("mpfr_mul_2exp", 53, GMP_RNDN);
+ testfunc = mpfr_mul_ui; test2ui ("mpfr_mul_ui", 53, GMP_RNDN);
+ testfunc = mpfr_neg; test2 ("mpfr_neg", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_pow_ui; test2ui ("mpfr_pow_ui", 53, GMP_RNDN);
+ testfunc = mpfr_reldiff; test3 ("mpfr_reldiff", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_sin_cos; test3a ("mpfr_sin_cos", 53, GMP_RNDN);
+ testfunc = mpfr_sub; test3 ("mpfr_sub", 53, GMP_RNDN);
+ testfunc = mpfr_sub_ui; test2ui ("mpfr_sub_ui", 53, GMP_RNDN);
+ testfunc = (void*) mpfr_sqrt; test2 ("mpfr_sqrt", 53, GMP_RNDN);
+ testfunc = mpfr_ui_div; testui2 ("mpfr_ui_div", 53, GMP_RNDN);
+ testfunc = mpfr_ui_sub; testui2 ("mpfr_ui_sub", 53, GMP_RNDN);
+ testfunc = mpfr_trunc; test2 ("mpfr_trunc", 53, GMP_RNDN);
+ return 0;
+}
diff --git a/mpfr/tests/tabs.c b/mpfr/tests/tabs.c
new file mode 100644
index 000000000..341ea516a
--- /dev/null
+++ b/mpfr/tests/tabs.c
@@ -0,0 +1,89 @@
+/* Test file for mpfr_abs.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 "mpfr-test.h"
+
+#define Infp (1/0.)
+
+extern int isnan();
+
+int main(int argc, char *argv[])
+{
+ mpfr_t x; int n, k, rnd; double d, dd;
+#ifdef __mips
+ /* to get denormalized numbers on IRIX64 */
+ union fpc_csr exp;
+ exp.fc_word = get_fpc_csr();
+ exp.fc_struct.flush = 0;
+ set_fpc_csr(exp.fc_word);
+#endif
+
+ mpfr_init2(x, 53);
+
+ mpfr_set_d(x, 1.0, GMP_RNDN);
+ mpfr_abs(x, x, GMP_RNDN);
+ if (mpfr_get_d(x) != 1.0) {
+ fprintf(stderr, "Error in mpfr_abs(1.0)\n"); exit(1);
+ }
+
+ mpfr_set_d(x, -1.0, GMP_RNDN);
+ mpfr_abs(x, x, GMP_RNDN);
+ if (mpfr_get_d(x) != 1.0) {
+ fprintf(stderr, "Error in mpfr_abs(-1.0)\n"); exit(1);
+ }
+
+ mpfr_set_d(x, -6/-0., GMP_RNDN);
+ mpfr_abs(x, x, GMP_RNDN);
+ if (mpfr_get_d(x) != Infp) {
+ fprintf(stderr, "Error in mpfr_abs(Inf).\n"); exit(1);
+ }
+
+ mpfr_set_d(x, 2/-0., GMP_RNDN);
+ mpfr_abs(x, x, GMP_RNDN);
+ if (mpfr_get_d(x) != Infp) {
+ fprintf(stderr, "Error in mpfr_abs(-Inf).\n"); exit(1);
+ }
+
+
+ n = (argc==1) ? 1000000 : atoi(argv[1]);
+ for (k = 1; k <= n; k++)
+ {
+ d = drand();
+ rnd = rand() % 4;
+ mpfr_set_d(x, d, 0);
+ mpfr_abs(x, x, rnd);
+ dd = mpfr_get_d(x);
+ if (dd != fabs(d) && !isnan(d))
+ {
+ fprintf(stderr,
+ "Mismatch on d = %1.18g\n", d);
+ mpfr_print_raw(x); putchar('\n');
+ exit(1);
+ }
+ }
+
+ mpfr_clear(x);
+ return 0;
+}
diff --git a/mpfr/tests/tadd.c b/mpfr/tests/tadd.c
index 517edb9e4..562d9b853 100644
--- a/mpfr/tests/tadd.c
+++ b/mpfr/tests/tadd.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_add and mpfr_sub.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -19,9 +19,6 @@ along with the MPFR 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. */
-/* #define DEBUG */
-/* #define VERBOSE */
-
#define N 100000
#include <math.h>
@@ -30,21 +27,26 @@ MA 02111-1307, USA. */
#include "gmp.h"
#include "mpfr.h"
#include "mpfr-impl.h"
-#ifdef IRIX64
-#include <sys/fpu.h>
-#endif
+#include "mpfr-test.h"
extern int isnan();
extern int getpid();
-
-#define ABS(x) (((x)>0) ? (x) : (-x))
+void check _PROTO((double, double, mp_rnd_t, unsigned int, unsigned int, unsigned int, double));
+void checknan _PROTO((double, double, mp_rnd_t, unsigned int, unsigned int, unsigned int));
+void check3 _PROTO((double, double, mp_rnd_t));
+void check4 _PROTO((double, double, mp_rnd_t));
+void check5 _PROTO((double, mp_rnd_t));
+void check2 _PROTO((double, int, double, int, int, int));
+void check2a _PROTO((double, int, double, int, int, int, char *));
+void check64 _PROTO((void));
+void check_same _PROTO((void));
/* checks that x+y gives the same results in double
and with mpfr with 53 bits of precision */
-void check(double x, double y, unsigned int rnd_mode, unsigned int px,
-unsigned int py, unsigned int pz, double res)
+void check (double x, double y, mp_rnd_t rnd_mode, unsigned int px,
+unsigned int py, unsigned int pz, double z1)
{
- double z1,z2; mpfr_t xx,yy,zz;
+ double z2; mpfr_t xx,yy,zz; int cert=0;
mpfr_init2(xx, px);
mpfr_init2(yy, py);
@@ -52,20 +54,46 @@ unsigned int py, unsigned int pz, double res)
mpfr_set_d(xx, x, rnd_mode);
mpfr_set_d(yy, y, rnd_mode);
mpfr_add(zz, xx, yy, rnd_mode);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd_mode);
- z1 = (res==0.0) ? x+y : res;
+ if (px==53 && py==53 && pz==53) cert=1;
+#endif
+ if (z1==0.0) z1=x+y; else cert=1;
z2 = mpfr_get_d(zz);
- if (px==53 && py==53 && pz==53) res=1.0;
- if (res!=0.0 && z1!=z2 && !(isnan(z1) && isnan(z2))) {
+ mpfr_set_d (yy, z2, GMP_RNDN);
+ if (!mpfr_cmp (xx, yy) && cert && z1!=z2 && !(isnan(z1) && isnan(z2))) {
printf("expected sum is %1.20e, got %1.20e\n",z1,z2);
- printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",x,y,rnd_mode);
+ printf("mpfr_add failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
+ x, y, mpfr_print_rnd_mode(rnd_mode));
exit(1);
}
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}
+void checknan (double x, double y, mp_rnd_t rnd_mode, unsigned int px,
+unsigned int py, unsigned int pz)
+{
+ double z2; mpfr_t xx, yy, zz;
+
+ mpfr_init2(xx, px);
+ mpfr_init2(yy, py);
+ mpfr_init2(zz, pz);
+ mpfr_set_d(xx, x, rnd_mode);
+ mpfr_set_d(yy, y, rnd_mode);
+ mpfr_add(zz, xx, yy, rnd_mode);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
+ if (MPFR_IS_NAN(zz) == 0) { printf("Error, not an MPFR_NAN for xx = %1.20e, y = %1.20e\n", x, y); exit(1); }
+ z2 = mpfr_get_d(zz);
+ if (!isnan(z2)) { printf("Error, not a NaN after conversion, xx = %1.20e yy = %1.20e\n", x, y); exit(1); }
+
+ mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+}
+
+#ifdef TEST
/* idem than check for mpfr_add(x, x, y) */
-void check3(double x, double y, unsigned int rnd_mode)
+void check3 (double x, double y, mp_rnd_t rnd_mode)
{
double z1,z2; mpfr_t xx,yy; int neg;
@@ -79,7 +107,8 @@ void check3(double x, double y, unsigned int rnd_mode)
mpfr_set_machine_rnd_mode(rnd_mode);
z1 = (neg) ? x-y : x+y;
z2 = mpfr_get_d(xx);
- if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
+ mpfr_set_d (yy, z2, GMP_RNDN);
+ if (!mpfr_cmp (xx, yy) && z1!=z2 && !(isnan(z1) && isnan(z2))) {
printf("expected result is %1.20e, got %1.20e\n",z1,z2);
printf("mpfr_%s(x,x,y) failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",
(neg) ? "sub" : "add",x,y,rnd_mode);
@@ -89,9 +118,11 @@ void check3(double x, double y, unsigned int rnd_mode)
}
/* idem than check for mpfr_add(x, y, x) */
-void check4(double x, double y, unsigned int rnd_mode)
+void check4 (double x, double y, mp_rnd_t rnd_mode)
{
- double z1,z2; mpfr_t xx,yy; int neg;
+ double z1, z2;
+ mpfr_t xx, yy;
+ int neg;
neg = rand() % 2;
mpfr_init2(xx, 53);
@@ -103,57 +134,56 @@ void check4(double x, double y, unsigned int rnd_mode)
mpfr_set_machine_rnd_mode(rnd_mode);
z1 = (neg) ? y-x : x+y;
z2 = mpfr_get_d(xx);
- if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
- printf("expected result is %1.20e, got %1.20e\n",z1,z2);
- printf("mpfr_%s(x,y,x) failed for x=%1.20e y=%1.20e with rnd_mode=%u\n",
- (neg) ? "sub" : "add",x,y,rnd_mode);
+ mpfr_set_d (yy, z2, GMP_RNDN);
+ /* check that xx is representable as a double and no overflow occurred */
+ if ((mpfr_cmp (xx, yy) == 0) && (z1 != z2)) {
+ printf("expected result is %1.20e, got %1.20e\n", z1, z2);
+ printf("mpfr_%s(x,y,x) failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
+ (neg) ? "sub" : "add", x, y, mpfr_print_rnd_mode(rnd_mode));
exit(1);
}
mpfr_clear(xx); mpfr_clear(yy);
}
/* idem than check for mpfr_add(x, x, x) */
-void check5(double x, unsigned int rnd_mode)
+void check5 (double x, mp_rnd_t rnd_mode)
{
- double z1,z2; mpfr_t xx; int neg;
+ double z1,z2; mpfr_t xx, yy; int neg;
- mpfr_init2(xx, 53); neg = rand() % 2;
+ mpfr_init2(xx, 53);
+ mpfr_init2(yy, 53);
+ neg = rand() % 2;
mpfr_set_d(xx, x, rnd_mode);
if (neg) mpfr_sub(xx, xx, xx, rnd_mode);
else mpfr_add(xx, xx, xx, rnd_mode);
mpfr_set_machine_rnd_mode(rnd_mode);
z1 = (neg) ? x-x : x+x;
z2 = mpfr_get_d(xx);
- if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
+ mpfr_set_d (yy, z2, GMP_RNDN);
+ if (!mpfr_cmp (xx, yy) && z1!=z2 && !(isnan(z1) && isnan(z2))) {
printf("expected result is %1.20e, got %1.20e\n",z1,z2);
- printf("mpfr_%s(x,x,x) failed for x=%1.20e with rnd_mode=%u\n",
- (neg) ? "sub" : "add",x,rnd_mode);
+ printf("mpfr_%s(x,x,x) failed for x=%1.20e with rnd_mode=%s\n",
+ (neg) ? "sub" : "add", x, mpfr_print_rnd_mode (rnd_mode));
exit(1);
}
mpfr_clear(xx);
+ mpfr_clear(yy);
}
-void check2(x,px,y,py,pz,rnd_mode) double x,y; int px,py,pz,rnd_mode;
+void check2 (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode)
{
mpfr_t xx, yy, zz; double z,z2; int u;
-#ifdef DEBUG
-printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%d\n",x,px,y,py,pz,rnd_mode);
-#endif
mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
mpfr_set_d(xx, x, rnd_mode);
mpfr_set_d(yy, y, rnd_mode);
mpfr_add(zz, xx, yy, rnd_mode);
mpfr_set_machine_rnd_mode(rnd_mode);
z = x+y; z2=mpfr_get_d(zz); u=ulp(z,z2);
-#ifdef DEBUG
- printf("x+y=%1.20e,%d (%d ulp) rnd_mode=%d\n",z2,pz,u,rnd_mode);
- mpfr_set_d(zz, z2, rnd_mode);
- printf("i.e."); mpfr_print_raw(zz); putchar('\n');
-#endif
- /* one ulp difference is possible due to composed rounding */
+ /* one ulp difference is possible due to composed rounding */
if (px>=53 && py>=53 && pz>=53 && ABS(u)>1) {
- printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%d\n",x,px,y,py,pz,rnd_mode);
+ printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
+ x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
printf("got %1.20e\n",z2);
printf("result should be %1.20e (diff=%d ulp)\n",z,u);
mpfr_set_d(zz, z, rnd_mode);
@@ -161,12 +191,53 @@ printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%d\n",x,px,y,py,pz,rnd_mode);
exit(1); }
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}
+#endif
-void check64()
+void check2a (double x, int px, double y, int py, int pz, mp_rnd_t rnd_mode,
+ char *res)
+{
+ mpfr_t xx, yy, zz;
+
+ mpfr_init2(xx,px); mpfr_init2(yy,py); mpfr_init2(zz,pz);
+ mpfr_set_d(xx, x, rnd_mode);
+ mpfr_set_d(yy, y, rnd_mode);
+ mpfr_add(zz, xx, yy, rnd_mode);
+ mpfr_set_prec(xx, pz);
+ mpfr_set_str(xx, res, 16, GMP_RNDN);
+ if (mpfr_cmp(xx, zz)) {
+ printf("x=%1.20e,%d y=%1.20e,%d pz=%d,rnd=%s\n",
+ x,px,y,py,pz,mpfr_print_rnd_mode(rnd_mode));
+ printf("got "); mpfr_print_raw(zz); putchar('\n');
+ printf("instead of "); mpfr_print_raw(xx); putchar('\n');
+ exit(1);
+ }
+ mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+}
+
+void check64 ()
{
mpfr_t x, t, u;
mpfr_init(x); mpfr_init(t); mpfr_init(u);
+
+ mpfr_set_prec (x, 112); mpfr_set_prec (t, 98); mpfr_set_prec (u, 54);
+ mpfr_set_str_raw (x, "-0.11111100100000000011000011100000101101010001000111E-401");
+ mpfr_set_str_raw (t, "0.10110000100100000101101100011111111011101000111000101E-464");
+ mpfr_add (u, x, t, GMP_RNDN);
+ if (mpfr_cmp (u, x)) {
+ fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=112, prec(t)=98\n");
+ exit (1);
+ }
+
+ mpfr_set_prec (x, 92); mpfr_set_prec (t, 86); mpfr_set_prec (u, 53);
+ mpfr_set_d (x, -5.03525136761487735093e-74, GMP_RNDN);
+ mpfr_set_d (t, 8.51539046314262304109e-91, GMP_RNDN);
+ mpfr_add (u, x, t, GMP_RNDN);
+ if (mpfr_get_d (u) != -5.0352513676148773509283672e-74) {
+ fprintf (stderr, "mpfr_add(u, x, t) failed for prec(x)=92, prec(t)=86\n");
+ exit (1);
+ }
+
mpfr_set_prec(x, 53); mpfr_set_prec(t, 76); mpfr_set_prec(u, 76);
mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
mpfr_set_str_raw(t, "-0.1011000101110010000101111111011111010001110011110111100110101011110010011111");
@@ -174,7 +245,7 @@ void check64()
mpfr_set_str_raw(t, "0.1011000101110010000101111111011100111111101010011011110110101011101000000100");
if (mpfr_cmp(u,t)) {
printf("expect "); mpfr_print_raw(t); putchar('\n');
- fprintf(stderr, "mpfr_add failed for precisions 53-76\n"); exit(1);
+ fprintf (stderr, "mpfr_add failed for precisions 53-76\n"); exit(1);
}
mpfr_set_prec(x, 53); mpfr_set_prec(t, 108); mpfr_set_prec(u, 108);
mpfr_set_str_raw(x, "-0.10010010001001011011110000000000001010011011011110001E-32");
@@ -210,7 +281,7 @@ void check64()
printf("Error in mpfr_sub: u=x-t and x=x-t give different results\n");
exit(1);
}
- if ((MANT(u)[(PREC(u)-1)/mp_bits_per_limb] &
+ if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
printf("Error in mpfr_sub: result is not msb-normalized\n"); exit(1);
}
@@ -221,22 +292,68 @@ void check64()
if (mpfr_get_d(u) != 9.4349060620538533806e167) { /* 2^558 */
printf("Error (1) in mpfr_sub\n"); exit(1);
}
+
mpfr_set_prec(x, 64); mpfr_set_prec(t, 64); mpfr_set_prec(u, 64);
mpfr_set_str_raw(x, "0.1000011110101111011110111111000011101011101111101101101100000100E-220");
mpfr_set_str_raw(t, "0.1000011110101111011110111111000011101011101111101101010011111101E-220");
mpfr_add(u, x, t, GMP_RNDU);
- if ((MANT(u)[0] & 1) != 1) {
+ if ((MPFR_MANT(u)[0] & 1) != 1) {
printf("error in mpfr_add with rnd_mode=GMP_RNDU\n");
printf("b= "); mpfr_print_raw(x); putchar('\n');
printf("c= "); mpfr_print_raw(t); putchar('\n');
printf("b+c="); mpfr_print_raw(u); putchar('\n');
exit(1);
}
+
+ /* bug found by Norbert Mueller, 14 Sep 2000 */
+ mpfr_set_prec(x, 56); mpfr_set_prec(t, 83); mpfr_set_prec(u, 10);
+ mpfr_set_str_raw(x, "0.10001001011011001111101100110100000101111010010111010111E-7");
+ mpfr_set_str_raw(t, "0.10001001011011001111101100110100000101111010010111010111000000000111110110110000100E-7");
+ mpfr_sub(u, x, t, GMP_RNDU);
+
+ /* array bound write found by Norbert Mueller, 26 Sep 2000 */
+ mpfr_set_prec(x, 109); mpfr_set_prec(t, 153); mpfr_set_prec(u, 95);
+ mpfr_set_str_raw(x,"0.1001010000101011101100111000110001111111111111111111111111111111111111111111111111111111111111100000000000000E33");
+ mpfr_set_str_raw(t,"-0.100101000010101110110011100011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011100101101000000100100001100110111E33");
+ mpfr_add(u, x, t, GMP_RNDN);
+
+ /* array bound writes found by Norbert Mueller, 27 Sep 2000 */
+ mpfr_set_prec(x, 106); mpfr_set_prec(t, 53); mpfr_set_prec(u, 23);
+ mpfr_set_str_raw(x, "-0.1000011110101111111001010001000100001011000000000000000000000000000000000000000000000000000000000000000000E-59");
+ mpfr_set_str_raw(t, "-0.10000111101011111110010100010001101100011100110100000E-59");
+ mpfr_sub(u, x, t, GMP_RNDN);
+ mpfr_set_prec(x, 177); mpfr_set_prec(t, 217); mpfr_set_prec(u, 160);
+ mpfr_set_str_raw(x, "-0.111010001011010000111001001010010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E35");
+ mpfr_set_str_raw(t, "0.1110100010110100001110010010100100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111011010011100001111001E35");
+ mpfr_add(u, x, t, GMP_RNDN);
+ mpfr_set_prec(x, 214); mpfr_set_prec(t, 278); mpfr_set_prec(u, 207);
+ mpfr_set_str_raw(x, "0.1000100110100110101101101101000000010000100111000001001110001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E66");
+ mpfr_set_str_raw(t, "-0.10001001101001101011011011010000000100001001110000010011100010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001111011111001001100011E66");
+ mpfr_add(u, x, t, GMP_RNDN);
+ mpfr_set_prec(x, 32); mpfr_set_prec(t, 247); mpfr_set_prec(u, 223);
+ mpfr_set_str_raw(x, "0.10000000000000000000000000000000E1");
+ mpfr_set_str_raw(t, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100000100011110000101110110011101110100110110111111011010111100100000000000000000000000000E0");
+ mpfr_sub(u, x, t, GMP_RNDN);
+ if ((MPFR_MANT(u)[(MPFR_PREC(u)-1)/mp_bits_per_limb] &
+ ((mp_limb_t)1<<(mp_bits_per_limb-1)))==0) {
+ printf("Error in mpfr_sub: result is not msb-normalized\n"); exit(1);
+ }
+
+ /* checks that NaN flag is correctly reset */
+ mpfr_set_d (t, 1.0, GMP_RNDN);
+ mpfr_set_d (u, 1.0, GMP_RNDN);
+ MPFR_SET_NAN(x);
+ mpfr_add (x, t, u, GMP_RNDN);
+ if (mpfr_cmp_ui (x, 2)) {
+ fprintf (stderr, "Error in mpfr_add: 1+1 gives %e\n", mpfr_get_d (x));
+ exit (1);
+ }
+
mpfr_clear(x); mpfr_clear(t); mpfr_clear(u);
}
/* checks when source and destination are equal */
-void check_same()
+void check_same ()
{
mpfr_t x;
@@ -248,10 +365,17 @@ void check_same()
mpfr_clear(x);
}
+#define check53(x, y, r, z) check(x, y, r, 53, 53, 53, z)
+#define check53nan(x, y, r) checknan(x, y, r, 53, 53, 53);
+
int main(argc,argv) int argc; char *argv[];
{
- double x,y; int i,prec,rnd_mode,px,py,pz,rnd;
-#ifdef IRIX64
+ int prec, rnd_mode;
+#ifdef TEST
+ int i, rnd;
+ double x, y;
+#endif
+#ifdef __mips
/* to get denormalized numbers on IRIX64 */
union fpc_csr exp;
exp.fc_word = get_fpc_csr();
@@ -260,155 +384,227 @@ int main(argc,argv) int argc; char *argv[];
#endif
check64();
+ check(293607738.0, 1.9967571564050541e-5, GMP_RNDU, 64, 53, 53,
+ 2.9360773800002003e8);
+ check(880524.0, -2.0769715792901673e-5, GMP_RNDN, 64, 53, 53,
+ 8.8052399997923023e5);
+ check(1196426492.0, -1.4218093058435347e-3, GMP_RNDN, 64, 53, 53,
+ 1.1964264919985781e9);
+ check(982013018.0, -8.941829477291838e-7, GMP_RNDN, 64, 53, 53,
+ 9.8201301799999905e8);
+ check(1092583421.0, 1.0880649218158844e9, GMP_RNDN, 64, 53, 53,
+ 2.1806483428158846e9);
+ check(1.8476886419022969e-6, 961494401.0, GMP_RNDN, 53, 64, 53,
+ 9.6149440100000179e8);
+ check(-2.3222118418069868e5, 1229318102.0, GMP_RNDN, 53, 64, 53,
+ 1.2290858808158193e9);
+ check(-3.0399171300395734e-6, 874924868.0, GMP_RNDN, 53, 64, 53,
+ 8.749248679999969e8);
+ check(9.064246624706179e1, 663787413.0, GMP_RNDN, 53, 64, 53,
+ 6.6378750364246619e8);
+ check(-1.0954322421551264e2, 281806592.0, GMP_RNDD, 53, 64, 53,
+ 2.8180648245677572e8);
+ check(5.9836930386056659e-8, 1016217213.0, GMP_RNDN, 53, 64, 53,
+ 1.0162172130000001e9);
+ check(-1.2772161928500301e-7, 1237734238.0, GMP_RNDN, 53, 64, 53,
+ 1.2377342379999998e9);
+ check(-4.567291988483277e8, 1262857194.0, GMP_RNDN, 53, 64, 53,
+ 8.0612799515167236e8);
+ check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 53, 53,
+ 2.4380935175292528e8);
+ check(4.7719471752925262e7, 196089880.0, GMP_RNDN, 53, 64, 53,
+ 2.4380935175292528e8);
+ check(-1.716113812768534e-140, 1271212614.0, GMP_RNDZ, 53, 64, 53,
+ 1.2712126139999998e9);
+ check(-1.2927455200185474e-50, 1675676122.0, GMP_RNDD, 53, 64, 53,
+ 1.6756761219999998e9);
+ check53(1.22191250737771397120e+20, 948002822.0, GMP_RNDN,
+ 122191250738719408128.0);
+ check53(9966027674114492.0, 1780341389094537.0, GMP_RNDN,
+ 11746369063209028.0);
+ check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
+ 3.5274425367757071711e272);
check_same();
- check(6.14384195492641560499e-02, -6.14384195401037683237e-02,
- GMP_RNDU, 53, 53, 53, 0.0);
- check(1.16809465359248765399e+196, 7.92883212101990665259e+196,
- GMP_RNDU, 53, 53, 53, 0.0);
- check(3.14553393112021279444e-67, 3.14553401015952024126e-67,
- GMP_RNDU, 53, 53, 53, 0.0);
-#ifdef VERBOSE
- printf("Checking random precisions\n");
-#endif
+ check53(6.14384195492641560499e-02, -6.14384195401037683237e-02, GMP_RNDU,
+ 9.1603877261370314499e-12);
+ check53(1.16809465359248765399e+196, 7.92883212101990665259e+196, GMP_RNDU,
+ 9.0969267746123943065e196);
+ check53(3.14553393112021279444e-67, 3.14553401015952024126e-67, GMP_RNDU,
+ 6.2910679412797336946e-67);
srand(getpid());
- check2(5.43885304644369509058e+185,53,-1.87427265794105342763e-57,53,53,0);
- check2(5.43885304644369509058e+185,53,-1.87427265794105342763e-57,53,53,1);
- check2(5.43885304644369509058e+185,53,-1.87427265794105342763e-57,53,53,2);
- check2(5.43885304644369509058e+185,53,-1.87427265794105342763e-57,53,53,3);
- check2(2.26531902208967707071e+168,99,-2.67795218510613988524e+168,67,94,2);
- check2(-1.21510626304662318398e+145,70,1.21367733647758957118e+145,65,61,3);
- check2(2.73028857032080744543e+155,83,-1.16446121423113355603e+163,59,125,1);
- check2(-4.38589520019641698848e+78,155,-1.09923643769309483415e+72,15,159,3);
- check2(-1.49963910666191123860e+265,76,-2.30915090591874527520e-191,8,75,1);
- check2(5.17945380930936917508e+112,119,1.11369077158813567738e+108,15,150,1);
- check2(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,1);
- check2(-1.87427265794105342764e-57,175,1.76570844587489516446e+190,2,115,1);
- check2(6.85523243386777784171e+107,187,-2.78148588123699111146e+48,87,178,3);
- check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,2);
- check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,3);
- check2(-3.31624349995221499866e-22,107,-8.20150212714204839621e+156,79,99,3);
- check2(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,3);
- check2(-1.08007920352320089721e+150,63,1.77607317509426332389e+73,64,64,0);
- check2(4.49465557237618783128e+53,108,-2.45103927353799477871e+48,60,105,0);
- check2(3.25471707846623300604e-160,81,-7.93846654265839958715e-274,58,54,0);
- check2(-8.88471912490080158206e+253,79,-7.84488427404526918825e+124,95,53,3);
- check2(-2.18548638152863831959e-125,61,-1.22788940592412363564e+226,71,54,0);
- check2(-7.94156823309993162569e+77,74,-5.26820160805275124882e+80,99,101,3);
- check2(-3.85170653452493859064e+189,62,2.18827389706660314090e+158,94,106,3);
- check2(1.07966151149311101950e+46,88,1.13198076934147457400e+46,67,53,0);
- check2(3.36768223223409657622e+209,55,-9.61624007357265441884e+219,113,53,0);
- check2(-6.47376909368979326475e+159,111,5.11127211034490340501e+159,99,62,3);
- check2(-4.95229483271607845549e+220,110,-6.06992115033276044773e+213,109,55,0);
- check2(-6.47376909368979326475e+159,74,5.11127211034490340501e+159,111,75,2);
- check2(2.26531902208967707070e+168,99,-2.67795218510613988525e+168,67,94,2);
- check2(-2.28886326552077479586e-188,67,3.41419438647157839320e-177,60,110,2);
- check2(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,1);
- check2(2.90983392714730768886e+50,101,2.31299792168440591870e+50,74,105,1);
- check2(2.72046257722708717791e+243,97,-1.62158447436486437113e+243,83,96,0);
- for (i=0;i<N;i++) {
-#ifdef DEBUG
-printf("\nTest i=%d\n",i);
-#endif
- px = 53 + (rand() % 64);
- py = 53 + (rand() % 64);
- pz = 53 + (rand() % 64);
- rnd_mode = rand() % 4;
- do { x = drand(); } while (isnan(x));
- do { y = drand(); } while (isnan(y));
- check2(x,px,y,py,pz,rnd_mode);
- }
-#ifdef VERBOSE
- printf("Checking double precision (53 bits)\n");
-#endif
+ check53(5.43885304644369509058e+185,-1.87427265794105342763e-57,GMP_RNDN,
+ 5.4388530464436950905e185);
+ check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDZ,
+ 5.4388530464436944867e185);
+ check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDU,
+ 5.4388530464436950905e185);
+ check53(5.43885304644369509058e+185,-1.87427265794105342763e-57, GMP_RNDD,
+ 5.4388530464436944867e185);
+ check2a(6.85523243386777784171e+107,187,-2.78148588123699111146e+48,87,178,
+ GMP_RNDD, "4.ab980a5cb9407ffffffffffffffffffffffffffffffe@89");
+ check2a(-1.21510626304662318398e+145,70,1.21367733647758957118e+145,65,61,
+ GMP_RNDD, "-1.2bfad031d94@118");
+ check2a(2.73028857032080744543e+155,83,-1.16446121423113355603e+163,59,125,
+ GMP_RNDZ, "-3.3c42dee09703d0639a6@135");
+ check2a(-4.38589520019641698848e+78,155,-1.09923643769309483415e+72,15,159,
+ GMP_RNDD, "-2.5e09955c663d@65");
+ check2a(-1.49963910666191123860e+265,76,-2.30915090591874527520e-191,8,75,
+ GMP_RNDZ, "-1.dc3ec027da54e@220");
+ check2a(3.25471707846623300604e-160,81,-7.93846654265839958715e-274,58,54,
+ GMP_RNDN, "4.936a52bc17254@-133");
+ check2a(5.17945380930936917508e+112,119,1.11369077158813567738e+108,15,150,
+ GMP_RNDZ, "5.62661692498ec@93");
+ check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
+ GMP_RNDZ, "-a.204acdd25d788@-44");
+ check2a(-1.87427265794105342764e-57,175,1.76570844587489516446e+190,2,115,
+ GMP_RNDZ, "b.fffffffffffffffffffffffffffe@157");
+ check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
+ GMP_RNDU, "-b.eae2643497ff6286b@-108");
+ check2a(-1.15706375390780417299e-135,94,-1.07455137477117851576e-129,66,111,
+ GMP_RNDD, "-b.eae2643497ff6286b@-108");
+ check2a(-3.31624349995221499866e-22,107,-8.20150212714204839621e+156,79,99,
+ GMP_RNDD, "-2.63b22b55697e8000000000008@130");
+ check2a(-1.08007920352320089721e+150,63,1.77607317509426332389e+73,64,64,
+ GMP_RNDN, "-5.4781549356e1c@124");
+ check2a(4.49465557237618783128e+53,108,-2.45103927353799477871e+48,60,105,
+ GMP_RNDN, "4.b14f230f909dc803e@44");
+ check2a(2.26531902208967707071e+168,99,-2.67795218510613988524e+168,67,94,
+ GMP_RNDU, "-1.bfd7ff2647098@139");
+ check2a(-8.88471912490080158206e+253,79,-7.84488427404526918825e+124,95,53,
+ GMP_RNDD, "-c.1e533b8d835@210");
+ check2a(-2.18548638152863831959e-125,61,-1.22788940592412363564e+226,71,54,
+ GMP_RNDN, "-8.4b0f99ffa3b58@187");
+ check2a(-7.94156823309993162569e+77,74,-5.26820160805275124882e+80,99,101,
+ GMP_RNDD, "-1.1cc90f11d6af26f4@67");
+ check2a(-3.85170653452493859064e+189,62,2.18827389706660314090e+158,94,106,
+ GMP_RNDD, "-3.753ac0935b701ffffffffffffd@157");
+ check2a(1.07966151149311101950e+46,88,1.13198076934147457400e+46,67,53,
+ GMP_RNDN, "3.dfbc152dd4368@38");
+ check2a(3.36768223223409657622e+209,55,-9.61624007357265441884e+219,113,53,
+ GMP_RNDN, "-6.cf7217a451388@182");
+ check2a(-6.47376909368979326475e+159,111,5.11127211034490340501e+159,99,62,
+ GMP_RNDD, "-1.8cf3aadf537c@132");
+ check2a(-4.95229483271607845549e+220,110,-6.06992115033276044773e+213,109,55,
+ GMP_RNDN, "-2.3129f1f63b31b@183");
+ check2a(-6.47376909368979326475e+159,74,5.11127211034490340501e+159,111,75,
+ GMP_RNDU, "-1.8cf3aadf537c@132");
+ check2a(2.26531902208967707070e+168,99,-2.67795218510613988525e+168,67,94,
+ GMP_RNDU, "-1.bfd7ff2647098@139");
+ check2a(-2.28886326552077479586e-188,67,3.41419438647157839320e-177,60,110,
+ GMP_RNDU, "3.75740b4fe8f17f90258907@-147");
+ check2a(-2.66910493504493276454e-52,117,1.61188644159592323415e-52,61,68,
+ GMP_RNDZ, "-a.204acdd25d788@-44");
+ check2a(2.90983392714730768886e+50,101,2.31299792168440591870e+50,74,105,
+ GMP_RNDZ, "1.655c53ff5719c8@42");
+ check2a(2.72046257722708717791e+243,97,-1.62158447436486437113e+243,83,96,
+ GMP_RNDN, "a.4cc63e002d2e8@201");
+ /* Checking double precision (53 bits) */
prec = (argc<2) ? 53 : atoi(argv[1]);
rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
- check(-8.22183238641455905806e-19, 7.42227178769761587878e-19,
- GMP_RNDD, 53, 53, 53, 0.0);
- check(5.82106394662028628236e+234, -5.21514064202368477230e+89,
- GMP_RNDD, 53, 53, 53, 0.0);
- check(5.72931679569871602371e+122, -5.72886070363264321230e+122,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(-5.09937369394650450820e+238, 2.70203299854862982387e+250,
- GMP_RNDD, 53, 53, 53, 0.0);
- check(-2.96695924472363684394e+27, 1.22842938251111500000e+16,
- GMP_RNDD, 53, 53, 53, 0.0);
- check(1.74693641655743793422e-227, -7.71776956366861843469e-229,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(-1.03432206392780011159e-125, 1.30127034799251347548e-133,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(1.05824655795525779205e+71, -1.06022698059744327881e+71,
- GMP_RNDZ, 53, 53, 53, 0.0);
- check(-5.84204911040921732219e+240, 7.26658169050749590763e+240,
- GMP_RNDD, 53, 53, 53, 0.0);
+ check53(-8.22183238641455905806e-19, 7.42227178769761587878e-19, GMP_RNDD,
+ -7.9956059871694317927e-20);
+ check53(5.82106394662028628236e+234, -5.21514064202368477230e+89, GMP_RNDD,
+ 5.8210639466202855763e234);
+ check53(5.72931679569871602371e+122, -5.72886070363264321230e+122, GMP_RNDN,
+ 4.5609206607281141508e118);
+ check53(-5.09937369394650450820e+238, 2.70203299854862982387e+250, GMP_RNDD,
+ 2.7020329985435301323e250);
+ check53(-2.96695924472363684394e+27, 1.22842938251111500000e+16, GMP_RNDD,
+ -2.96695924471135255027e27);
+ check53(1.74693641655743793422e-227, -7.71776956366861843469e-229, GMP_RNDN,
+ 1.669758720920751867e-227);
+ check53(-1.03432206392780011159e-125, 1.30127034799251347548e-133, GMP_RNDN,
+ -1.0343220509150965661e-125);
+ check53(1.05824655795525779205e+71, -1.06022698059744327881e+71, GMP_RNDZ,
+ -1.9804226421854867632e68);
+ check53(-5.84204911040921732219e+240, 7.26658169050749590763e+240, GMP_RNDD,
+ 1.4245325800982785854e240);
/* the following check double overflow */
- /* check(6.27557402141211962228e+307, 1.32141396570101687757e+308,
- GMP_RNDZ, 53, 53, 53, 0.0); */
- check(1.00944884131046636376e+221, 2.33809162651471520268e+215,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(4.29232078932667367325e-278, 1.07735250473897938332e-281,
- GMP_RNDU, 53, 53, 53, 0.0);
- check(5.27584773801377058681e-80, 8.91207657803547196421e-91,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(2.99280481918991653800e+272, 5.34637717585790933424e+271,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(4.67302514390488041733e-184, 2.18321376145645689945e-190,
- GMP_RNDN, 53, 53, 53, 0.0);
- check(5.57294120336300389254e+71, 2.60596167942024924040e+65,
- GMP_RNDZ, 53, 53, 53, 0.0);
- x=6151626677899716.0; for (i=0;i<30;i++) x = 2.0*x;
- check(x, 4938448004894539.0, GMP_RNDU, 53, 53, 53, 0.0);
- check(1.23056185051606761523e-190, 1.64589756643433857138e-181,
- GMP_RNDU, 53, 53, 53, 0.0);
- check(2.93231171510175981584e-280, 3.26266919161341483877e-273,
- GMP_RNDU, 53, 53, 53, 0.0);
- check(5.76707395945001907217e-58, 4.74752971449827687074e-51,
- GMP_RNDD, 53, 53, 53, 0.0);
- check(277363943109.0, 11.0, GMP_RNDN, 53, 53, 53, 0.0);
+ check53(6.27557402141211962228e+307, 1.32141396570101687757e+308,
+ GMP_RNDZ, 1.0/0.0);
+ check53(1.00944884131046636376e+221, 2.33809162651471520268e+215, GMP_RNDN,
+ 1.0094511794020929787e221);
+ check53(4.29232078932667367325e-278, 1.07735250473897938332e-281, GMP_RNDU,
+ 4.2933981418314132787e-278);
+ check53(5.27584773801377058681e-80, 8.91207657803547196421e-91, GMP_RNDN,
+ 5.2758477381028917269e-80);
+ check53(2.99280481918991653800e+272, 5.34637717585790933424e+271, GMP_RNDN,
+ 3.5274425367757071711e272);
+ check53(4.67302514390488041733e-184, 2.18321376145645689945e-190, GMP_RNDN,
+ 4.6730273271186420541e-184);
+ check53(5.57294120336300389254e+71, 2.60596167942024924040e+65, GMP_RNDZ,
+ 5.5729438093246831053e71);
+ check53(6.6052588496951015469e24, 4938448004894539.0, GMP_RNDU,
+ 6.6052588546335505068e24);
+ check53(1.23056185051606761523e-190, 1.64589756643433857138e-181, GMP_RNDU,
+ 1.6458975676649006598e-181);
+ check53(2.93231171510175981584e-280, 3.26266919161341483877e-273, GMP_RNDU,
+ 3.2626694848445867288e-273);
+ check53(5.76707395945001907217e-58, 4.74752971449827687074e-51, GMP_RNDD,
+ 4.747530291205672325e-51);
+ check53(277363943109.0, 11.0, GMP_RNDN, 277363943120.0);
/* test denormalized numbers too */
- check(8.06294740693074521573e-310, 6.95250701071929654575e-310,
- GMP_RNDU, 53, 53, 53, 0.0);
-#ifdef VERBOSE
- printf("Comparing to double precision using machine arithmetic\n");
-#endif
+ check53(8.06294740693074521573e-310, 6.95250701071929654575e-310, GMP_RNDU,
+ 1.5015454417650041761e-309);
+ check53(1/0., 6.95250701071929654575e-310, GMP_RNDU, 1/0.);
+ check53(-1/0., 6.95250701071929654575e-310, GMP_RNDU, -1/0.);
+ check53(6.95250701071929654575e-310, 1/0., GMP_RNDU, 1/0.);
+ check53(6.95250701071929654575e-310, -1/0., GMP_RNDU, -1/0.);
+ check53(1.44791789689198883921e-140, -1.90982880222349071284e-121,
+ GMP_RNDN, -1.90982880222349071e-121);
+
+ check53nan(1/0., -1/0., GMP_RNDN);
+
+#ifdef TEST
+ /* Comparing to double precision using machine arithmetic */
for (i=0;i<N;i++) {
x = drand();
y = drand();
- if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308)
+ if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
/* avoid denormalized numbers and overflows */
rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
check(x, y, rnd, prec, prec, prec, 0.0);
+ }
}
-#ifdef VERBOSE
- printf("Checking mpfr_add(x, x, y) with prec=53\n");
-#endif
+ /* tests with random precisions */
+ for (i=0;i<N;i++) {
+ int px, py, pz;
+ px = 53 + (rand() % 64);
+ py = 53 + (rand() % 64);
+ pz = 53 + (rand() % 64);
+ rnd_mode = rand() % 4;
+ do { x = drand(); } while (isnan(x));
+ do { y = drand(); } while (isnan(y));
+ check2(x,px,y,py,pz,rnd_mode);
+ }
+ /* Checking mpfr_add(x, x, y) with prec=53 */
for (i=0;i<N;i++) {
x = drand();
y = drand();
- if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308)
+ if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
/* avoid denormalized numbers and overflows */
rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
check3(x, y, rnd);
- }
-#ifdef VERBOSE
- printf("Checking mpfr_add(x, y, x) with prec=53\n");
-#endif
+ }
+ }
+ /* Checking mpfr_add(x, y, x) with prec=53 */
for (i=0;i<N;i++) {
x = drand();
y = drand();
- if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308)
+ if (ABS(x)>2.2e-307 && ABS(y)>2.2e-307 && x+y<1.7e+308 && x+y>-1.7e308) {
/* avoid denormalized numbers and overflows */
rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
check4(x, y, rnd);
- }
-#ifdef VERBOSE
- printf("Checking mpfr_add(x, x, x) with prec=53\n");
-#endif
+ }
+ }
+ /* Checking mpfr_add(x, x, x) with prec=53 */
for (i=0;i<N;i++) {
do { x = drand(); } while ((ABS(x)<2.2e-307) || (ABS(x)>0.8e308));
/* avoid denormalized numbers and overflows */
rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
check5(x, rnd);
- }
- exit (0);
+ }
+#endif
+ return 0;
}
diff --git a/mpfr/tests/tagm.c b/mpfr/tests/tagm.c
index 1782911f2..e6c411677 100644
--- a/mpfr/tests/tagm.c
+++ b/mpfr/tests/tagm.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_agm.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -26,9 +26,13 @@ MA 02111-1307, USA. */
#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
-extern int isnan();
+double drand_agm _PROTO((void));
+double dagm _PROTO((double, double));
+void check4 _PROTO((double, double, mp_rnd_t, double));
+void check_large _PROTO((void));
+void slave _PROTO((int, int));
double drand_agm()
{
@@ -46,27 +50,14 @@ double drand_agm()
return d;
}
-
-double max(double a,double b) {
- if (a>=b)
- return a;
- return b;
-}
-
-double min(double a,double b) {
- if (b>=a)
- return a;
- return b;
-}
-
-double dagm(double a, double b) {
- double u,v,tmpu,tmpv;
+double dagm (double a, double b) {
+ double u, v, tmpu, tmpv;
if ((isnan(a))||(isnan(b)))
return a+b;
- tmpv=max(a,b);
- tmpu=min(a,b);
+ tmpv=MAX(a,b);
+ tmpu=MIN(a,b);
do
{
@@ -83,7 +74,7 @@ double dagm(double a, double b) {
#define check(a,b,r) check4(a,b,r,0.0)
-void check4(double a, double b, unsigned char rnd_mode, double res1)
+void check4 (double a, double b, mp_rnd_t rnd_mode, double res1)
{
mpfr_t ta, tb, tres;
double res2;
@@ -97,9 +88,12 @@ void check4(double a, double b, unsigned char rnd_mode, double res1)
mpfr_set_d(tb, b, rnd_mode);
mpfr_agm(tres, ta, tb, rnd_mode);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
if (res1==0.0) res1=dagm(a,b); else ck=1;
+if (ck==0) printf("%1.20e\n", res1);
res2 = mpfr_get_d(tres);
if (ck && res1!=res2 && (!isnan(res1) || !isnan(res2))) {
@@ -111,7 +105,7 @@ void check4(double a, double b, unsigned char rnd_mode, double res1)
mpfr_clear(ta); mpfr_clear(tb); mpfr_clear(tres);
}
-void check_large()
+void check_large ()
{
mpfr_t a, b, agm;
@@ -126,7 +120,8 @@ void check_large()
mpfr_clear(a); mpfr_clear(b); mpfr_clear(agm);
}
-void slave(int N, int p) {
+void slave (int N, int p)
+{
int i;
double a,b;
mpfr_t ta, tb, tres;
@@ -153,7 +148,7 @@ int main(int argc, char* argv[]) {
if (argc==3) { /* tagm N p : N calculus with precision p*/
printf("Doing %d random tests in %d precision\n",atoi(argv[1]),atoi(argv[2]));
slave(atoi(argv[1]),atoi(argv[2]));
- exit (0);
+ return 0;
}
if (argc==2) { /* tagm N: N tests with random double's */
@@ -167,19 +162,21 @@ int main(int argc, char* argv[]) {
b = drand();
check(a, b, rand() % 4);
}
- exit (0);
+ return 0;
}
else {
check_large();
- check(2,1,GMP_RNDN);
- check(6,4,GMP_RNDN);
- check(62,61,GMP_RNDN);
- check(0.5,1,GMP_RNDN);
- check(1,2,GMP_RNDN);
- check4(234375765,234375000,GMP_RNDN,2.34375382499843955040e+08);
- check(8,1,GMP_RNDU);
- check(1,44,GMP_RNDU);
- check(1,3.725290298461914062500000e-9,GMP_RNDU);
+ check4(2.0, 1.0, GMP_RNDN, 1.45679103104690677029);
+ check4(6.0, 4.0, GMP_RNDN, 4.94936087247260925182);
+ check4(62.0, 61.0, GMP_RNDN, 6.14989837188450749750e+01);
+ check4(0.5, 1.0, GMP_RNDN, 7.28395515523453385143e-01);
+ check4(1.0, 2.0, GMP_RNDN, 1.45679103104690677029);
+ check4(234375765.0, 234375000.0, GMP_RNDN, 2.3437538249984395504e8);
+ check4(8.0, 1.0, GMP_RNDU, 3.615756177597362786);
+ check4(1.0, 44.0, GMP_RNDU, 1.33658354512981247808e1);
+ check4(1.0, 3.7252902984619140625e-9, GMP_RNDU, 7.55393356971199025907e-02);
}
- exit (0);
+
+ /* TODO : tests des infinis dans tagm.c */
+ return 0;
}
diff --git a/mpfr/tests/tcan_round.c b/mpfr/tests/tcan_round.c
index 7e452ab80..4d5fd0cb5 100644
--- a/mpfr/tests/tcan_round.c
+++ b/mpfr/tests/tcan_round.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_can_round.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -35,6 +35,14 @@ int main()
if (mpfr_can_round(x, 54, GMP_RNDZ, GMP_RNDZ, 53) != 0) {
fprintf(stderr, "Error in mpfr_can_round\n"); exit(1);
}
+
+ mpfr_set_str_raw(x, "-Inf");
+ if (mpfr_can_round(x, 2000, GMP_RNDZ, GMP_RNDZ, 2000) != 0) {
+ fprintf(stderr, "Error in mpfr_can_round\n"); exit(1);
+ }
mpfr_clear(x);
- exit (0);
+ return 0;
}
+
+
+
diff --git a/mpfr/tests/tcmp.c b/mpfr/tests/tcmp.c
index 0f063c235..3d033314d 100644
--- a/mpfr/tests/tcmp.c
+++ b/mpfr/tests/tcmp.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_cmp.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,17 +23,19 @@ MA 02111-1307, USA. */
#include <stdlib.h>
#include <math.h>
#include "gmp.h"
-#include "longlong.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
-int
-main()
-{
- double x,y; mpfr_t xx,yy; int i,c;
+#define Infp 1/0.
+#define Infm -1/0.
+
+extern int isnan();
- fprintf(stderr, "Test case 'tcmp' disabled\n");
- exit(0); /* THIS TEST CASE IS NOT WORKING */
+int main()
+{
+ double x, y;
+ mpfr_t xx, yy;
+ int i, c;
mpfr_init2(xx, 65); mpfr_init2(yy, 65);
mpfr_set_str_raw(xx, "0.10011010101000110101010000000011001001001110001011101011111011101E623");
@@ -41,8 +43,17 @@ main()
if (mpfr_cmp2(xx,yy)!=64) { printf("Error (1) in mpfr_cmp\n"); exit(1); }
mpfr_set_str_raw(xx, "0.10100010001110110111000010001000010011111101000100011101000011100");
mpfr_set_str_raw(yy, "0.10100010001110110111000010001000010011111101000100011101000011011");
- if (mpfr_cmp2(xx,yy)!=64) { printf("Error (1) in mpfr_cmp\n"); exit(1); }
- mpfr_set_prec(xx,53); mpfr_set_prec(yy,200);
+ if (mpfr_cmp2(xx,yy)!=64) { printf("Error (2) in mpfr_cmp\n"); exit(1); }
+
+ mpfr_set_prec (xx, 160); mpfr_set_prec (yy, 160);
+ mpfr_set_str_raw (xx, "0.1E1");
+ mpfr_set_str_raw (yy, "0.1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100000110001110100");
+ if (mpfr_cmp2 (xx, yy) != 144) {
+ printf("Error (3) in mpfr_cmp\n");
+ exit(1);
+ }
+
+ mpfr_set_prec(xx, 53); mpfr_set_prec(yy, 200);
mpfr_set_d(xx, 1.0, 0);
mpfr_set_d(yy, 1.0, 0);
if (mpfr_cmp(xx,yy)!=0) {
@@ -55,8 +66,76 @@ main()
printf("Error in mpfr_cmp: not 1.0000000002 > 1.0\n"); exit(1);
}
mpfr_set_prec(yy, 53);
- for (i=0;i<1000000;) {
- x=drand(); y=drand();
+
+ /* bug found by Gerardo Ballabio */
+ mpfr_set_d(xx, 0.0, GMP_RNDN);
+ mpfr_set_d(yy, 0.1, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) >= 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(0.0, 0.1), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infp, GMP_RNDN);
+ mpfr_set_d(yy, -23489745.0329, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) <= 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infp, 23489745.0329), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infp, GMP_RNDN);
+ mpfr_set_d(yy, Infm, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) <= 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infp, Infm), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infm, GMP_RNDN);
+ mpfr_set_d(yy, Infp, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) >= 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infm, Infp), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infp, GMP_RNDN);
+ mpfr_set_d(yy, Infp, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) != 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infp, Infp), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infm, GMP_RNDN);
+ mpfr_set_d(yy, Infm, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) != 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infm, Infm), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d(xx, Infm, GMP_RNDN);
+ mpfr_set_d(yy, 2346.09234, GMP_RNDN);
+ if (mpfr_cmp(xx, yy) >= 0) {
+ fprintf(stderr,
+ "Error in mpfr_cmp(Infm, 2346.09234), gives %d\n", mpfr_cmp(xx, yy));
+ exit(1);
+ }
+
+ mpfr_set_d (xx, 0.0, GMP_RNDN);
+ mpfr_set_d (yy, 1.0, GMP_RNDN);
+ if ((i = mpfr_cmp3 (xx, yy, 1)) >= 0) {
+ fprintf (stderr, "Error: mpfr_cmp3 (0, 1, 1) gives %d instead of a negative value\n", i);
+ exit (1);
+ }
+ if ((i = mpfr_cmp3 (xx, yy, -1)) <= 0) {
+ fprintf (stderr, "Error: mpfr_cmp3 (0, 1, -1) gives %d instead of a positive value\n", i);
+ exit (1);
+ }
+
+ for (i=0;i<1000000;) { x=drand(); y=drand();
if (!isnan(x) && !isnan(y)) {
i++;
mpfr_set_d(xx, x, 0);
@@ -68,6 +147,7 @@ main()
}
}
}
+
mpfr_clear(xx); mpfr_clear(yy);
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tcmp2.c b/mpfr/tests/tcmp2.c
index f7e3f3599..4ddb16bea 100644
--- a/mpfr/tests/tcmp2.c
+++ b/mpfr/tests/tcmp2.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_cmp2.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999-2000 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,18 +23,17 @@ MA 02111-1307, USA. */
#include <stdlib.h>
#include <math.h>
#include "gmp.h"
-#include "longlong.h"
#include "mpfr.h"
#include "mpfr-impl.h"
-#ifdef IRIX64
-#include <sys/fpu.h>
-#endif
+#include "mpfr-test.h"
-extern int isnan();
+void tcmp2 (double, double, int);
+void special (void);
-void tcmp2(x, y, i) double x, y; int i;
+void tcmp2 (double x, double y, int i)
{
- mpfr_t xx,yy; int j;
+ mpfr_t xx, yy;
+ int j;
if (i==-1) {
if (x==y) i=53;
@@ -43,27 +42,50 @@ void tcmp2(x, y, i) double x, y; int i;
mpfr_init2(xx, 53); mpfr_init2(yy, 53);
mpfr_set_d(xx, x, 0);
mpfr_set_d(yy, y, 0);
- if ((j=mpfr_cmp2(xx, yy)) != i) {
+ j = mpfr_cmp2(xx, yy);
+ if (j != i) {
printf("Error in mpfr_cmp2: x=%1.20e y=%1.20e mpfr_cmp2(x,y)=%d instead of %d\n",x,y,j,i);
exit(1);
}
- mpfr_set_prec(xx, 127); mpfr_set_prec(yy, 127);
- mpfr_set_str_raw(xx, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6");
- mpfr_set_str_raw(yy, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6");
- if ((j=mpfr_cmp2(xx, yy)) != 32) {
+ mpfr_clear(xx); mpfr_clear(yy);
+}
+
+void special ()
+{
+ mpfr_t x, y;
+ int j;
+
+ mpfr_init (x); mpfr_init (y);
+
+ mpfr_set_prec(x, 127); mpfr_set_prec(y, 127);
+ mpfr_set_str_raw(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6");
+ mpfr_set_str_raw(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6");
+ if ((j=mpfr_cmp2(x, y)) != 32) {
printf("Error in mpfr_cmp2:\n");
- printf("x="); mpfr_print_raw(xx); putchar('\n');
- printf("y="); mpfr_print_raw(yy); putchar('\n');
+ printf("x="); mpfr_print_raw(x); putchar('\n');
+ printf("y="); mpfr_print_raw(y); putchar('\n');
printf("got %d, expected 32\n", j);
exit(1);
}
- mpfr_clear(xx); mpfr_clear(yy);
+
+ mpfr_set_prec (x, 128); mpfr_set_prec (y, 239);
+ mpfr_set_str_raw (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167");
+ mpfr_set_str_raw (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167");
+ if ((j=mpfr_cmp2(x, y)) != 164) {
+ printf("Error in mpfr_cmp2:\n");
+ printf("x="); mpfr_print_raw(x); putchar('\n');
+ printf("y="); mpfr_print_raw(y); putchar('\n');
+ printf("got %d, expected 164\n", j);
+ exit(1);
+ }
+
+ mpfr_clear(x); mpfr_clear(y);
}
int main()
{
int i,j; double x=1.0, y, z;
-#ifdef IRIX64
+#ifdef __mips
/* to get denormalized numbers on IRIX64 */
union fpc_csr exp;
exp.fc_word = get_fpc_csr();
@@ -71,18 +93,20 @@ int main()
set_fpc_csr(exp.fc_word);
#endif
+ special ();
+ tcmp2(5.43885304644369510000e+185, -1.87427265794105340000e-57, 1);
tcmp2(1.06022698059744327881e+71, 1.05824655795525779205e+71, -1);
tcmp2(1.0, 1.0, 53);
for (i=0;i<54;i++) {
tcmp2(1.0, 1.0-x, i);
x /= 2.0;
}
- for (j=0;j<1000000;j++) {
+ for (j=0; j<100000; j++) {
x = drand48();
y = drand48();
if (x<y) { z=x; x=y; y=z; }
if (y != 0.0 && y != -0.0) tcmp2(x, y, -1);
}
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tcmp_ui.c b/mpfr/tests/tcmp_ui.c
index 61410edd0..ce0eb2ca5 100644
--- a/mpfr/tests/tcmp_ui.c
+++ b/mpfr/tests/tcmp_ui.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_cmp_ui.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,36 +23,48 @@ MA 02111-1307, USA. */
#include <stdlib.h>
#include <math.h>
#include "gmp.h"
-#include "longlong.h"
#include "mpfr.h"
-int
-main()
+int main()
{
mpfr_t x; unsigned long i; long s;
mpfr_init(x);
mpfr_set_ui(x, 3, GMP_RNDZ);
- if (mpfr_cmp_ui(x, i=3)!=0) {
- printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+ if (mpfr_cmp_ui(x, i=3)!=0) {
+ printf("Error in mpfr_cmp_ui(%1.20f,%lu)\n",mpfr_get_d(x), i);
+ mpfr_clear(x);
+ exit(1);
}
if (mpfr_cmp_ui(x, i=2)<=0) {
- printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+ printf("Error in mpfr_cmp_ui(%1.20f,%lu)\n",mpfr_get_d(x), i);
+ mpfr_clear(x);
+ exit(1);
}
if (mpfr_cmp_ui(x, i=4)>=0) {
- printf("Error in mpfr_cmp_ui(%1.20f,%d)\n",mpfr_get_d(x), i); exit(1);
+ printf("Error in mpfr_cmp_ui(%1.20f,%lu)\n",mpfr_get_d(x), i);
+ mpfr_clear(x);
+ exit(1);
}
mpfr_set_si(x, -3, GMP_RNDZ);
if (mpfr_cmp_si(x, s=-3)!=0) {
- printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+ printf("Error in mpfr_cmp_si(%1.20f,%ld)\n",mpfr_get_d(x), s);
+ mpfr_clear(x);
+ exit(1);
}
if (mpfr_cmp_si(x, s=-4)<=0) {
- printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+ printf("Error in mpfr_cmp_si(%1.20f,%ld)\n",mpfr_get_d(x), s);
+ mpfr_clear(x);
+ exit(1);
}
if (mpfr_cmp_si(x, s=1)>=0) {
- printf("Error in mpfr_cmp_si(%1.20f,%d)\n",mpfr_get_d(x), s); exit(1);
+ printf("Error in mpfr_cmp_si(%1.20f,%ld)\n",mpfr_get_d(x), s);
+ mpfr_clear(x);
+ exit(1);
}
- exit (0);
+
+ mpfr_clear(x);
+ return 0;
}
diff --git a/mpfr/tests/tdiv.c b/mpfr/tests/tdiv.c
index a7de5dde9..b97e5d1b4 100644
--- a/mpfr/tests/tdiv.c
+++ b/mpfr/tests/tdiv.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_div.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -22,62 +22,54 @@ MA 02111-1307, USA. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
+#include <unistd.h>
#include "gmp.h"
-#include "gmp-impl.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
-#ifdef IRIX64
-#include <sys/fpu.h>
-#endif
-
-extern int isnan();
+#include "mpfr-test.h"
-/* #define DEBUG */
+#define check53(n, d, rnd, res) check4(n, d, rnd, 53, res)
-#define check(n,d,r) check4(n,d,r,53)
+void check4 _PROTO((double, double, mp_rnd_t, int, double));
+void check24 _PROTO((float, float, mp_rnd_t, float));
+void check_float _PROTO((void));
+void check_convergence _PROTO((void));
-void check4(N, D, rnd_mode, p) double N, D; unsigned char rnd_mode; int p;
+void check4 (double N, double D, mp_rnd_t rnd_mode, int p, double Q)
{
- mpfr_t q, n, d; double Q,Q2;
+ mpfr_t q, n, d; double Q2;
-#ifdef DEBUG
- printf("N=%1.20e D=%1.20e rnd_mode=%d\n",N,D,rnd_mode);
-#endif
mpfr_init2(q, p); mpfr_init2(n, p); mpfr_init2(d, p);
mpfr_set_d(n, N, rnd_mode);
mpfr_set_d(d, D, rnd_mode);
mpfr_div(q, n, d, rnd_mode);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd_mode);
- Q = N/D;
- Q2 = mpfr_get_d(q);
-#ifdef DEBUG
- printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
- ulp(Q2,Q));
- mpfr_print_raw(q); putchar('\n');
#endif
- if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
- printf("mpfr_div failed for n=%1.20e, d=%1.20e, rnd_mode=%d\n",N,D,rnd_mode);
- printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,
- ulp(Q2,Q));
+ if (Q==0.0) Q = N/D;
+ Q2 = mpfr_get_d(q);
+ if (p==53 && Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
+ printf("mpfr_div failed for n=%1.20e, d=%1.20e, rnd_mode=%s\n",
+ N, D, mpfr_print_rnd_mode(rnd_mode));
+ printf("expected quotient is %1.20e, got %1.20e (%d ulp)\n", Q, Q2,
+ ulp(Q2, Q));
exit(1);
}
mpfr_clear(q); mpfr_clear(n); mpfr_clear(d);
}
-void check24(float N, float D, unsigned char rnd_mode)
+void check24 (float N, float D, mp_rnd_t rnd_mode, float Q)
{
- mpfr_t q, n, d; float Q,Q2;
+ mpfr_t q, n, d; float Q2;
mpfr_init2(q, 24); mpfr_init2(n, 24); mpfr_init2(d, 24);
mpfr_set_d(n, N, rnd_mode);
mpfr_set_d(d, D, rnd_mode);
mpfr_div(q, n, d, rnd_mode);
- mpfr_set_machine_rnd_mode(rnd_mode);
- Q = N/D;
Q2 = mpfr_get_d(q);
- if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
- printf("mpfr_div failed for n=%1.10e, d=%1.10e, prec=24, rnd_mode=%d\n",N,D,rnd_mode);
- printf("expected quotient is %1.10e, got %1.10e\n",Q,Q2);
+ if (Q!=Q2) {
+ printf("mpfr_div failed for n=%1.10e, d=%1.10e, prec=24, rnd_mode=%s\n",
+ N, D, mpfr_print_rnd_mode(rnd_mode));
+ printf("expected quotient is %1.10e, got %1.10e\n", Q, Q2);
exit(1);
}
mpfr_clear(q); mpfr_clear(n); mpfr_clear(d);
@@ -87,66 +79,144 @@ void check24(float N, float D, unsigned char rnd_mode)
Generation for Directed Rounding" from Michael Parks, Table 2 */
void check_float()
{
- int i; float b=8388608.0; /* 2^23 */
-
- for (i=0;i<4;i++) {
- check24(b*8388610.0, 8388609.0, i);
- check24(b*16777215.0, 16777213.0, i);
- check24(b*8388612.0, 8388611.0, i);
- check24(b*12582914.0, 12582911.0, i);
- check24(12582913.0, 12582910.0, i);
- check24(b*16777215.0, 8388609.0, i);
- check24(b*8388612.0, 8388609.0, i);
- check24(b*12582914.0, 8388610.0, i);
- check24(b*12582913.0, 8388610.0, i);
- }
+ float b=8388608.0; /* 2^23 */
+
+ check24(b*8388610.0, 8388609.0, GMP_RNDN, 8.388609e6);
+ check24(b*16777215.0, 16777213.0, GMP_RNDN, 8.388609e6);
+ check24(b*8388612.0, 8388611.0, GMP_RNDN, 8.388609e6);
+ check24(b*12582914.0, 12582911.0, GMP_RNDN, 8.38861e6);
+ check24(12582913.0, 12582910.0, GMP_RNDN, 1.000000238);
+ check24(b*16777215.0, 8388609.0, GMP_RNDN, 1.6777213e7);
+ check24(b*8388612.0, 8388609.0, GMP_RNDN, 8.388611e6);
+ check24(b*12582914.0, 8388610.0, GMP_RNDN, 1.2582911e7);
+ check24(b*12582913.0, 8388610.0, GMP_RNDN, 1.258291e7);
+
+ check24(b*8388610.0, 8388609.0, GMP_RNDZ, 8.388608e6);
+ check24(b*16777215.0, 16777213.0, GMP_RNDZ, 8.388609e6);
+ check24(b*8388612.0, 8388611.0, GMP_RNDZ, 8.388608e6);
+ check24(b*12582914.0, 12582911.0, GMP_RNDZ, 8.38861e6);
+ check24(12582913.0, 12582910.0, GMP_RNDZ, 1.000000238);
+ check24(b*16777215.0, 8388609.0, GMP_RNDZ, 1.6777213e7);
+ check24(b*8388612.0, 8388609.0, GMP_RNDZ, 8.38861e6);
+ check24(b*12582914.0, 8388610.0, GMP_RNDZ, 1.2582911e7);
+ check24(b*12582913.0, 8388610.0, GMP_RNDZ, 1.258291e7);
+
+ check24(b*8388610.0, 8388609.0, GMP_RNDU, 8.388609e6);
+ check24(b*16777215.0, 16777213.0, GMP_RNDU, 8.38861e6);
+ check24(b*8388612.0, 8388611.0, GMP_RNDU, 8.388609e6);
+ check24(b*12582914.0, 12582911.0, GMP_RNDU, 8.388611e6);
+ check24(12582913.0, 12582910.0, GMP_RNDU, 1.000000357);
+ check24(b*16777215.0, 8388609.0, GMP_RNDU, 1.6777214e7);
+ check24(b*8388612.0, 8388609.0, GMP_RNDU, 8.388611e6);
+ check24(b*12582914.0, 8388610.0, GMP_RNDU, 1.2582912e7);
+ check24(b*12582913.0, 8388610.0, GMP_RNDU, 1.2582911e7);
+
+ check24(b*8388610.0, 8388609.0, GMP_RNDD, 8.388608e6);
+ check24(b*16777215.0, 16777213.0, GMP_RNDD, 8.388609e6);
+ check24(b*8388612.0, 8388611.0, GMP_RNDD, 8.388608e6);
+ check24(b*12582914.0, 12582911.0, GMP_RNDD, 8.38861e6);
+ check24(12582913.0, 12582910.0, GMP_RNDD, 1.000000238);
+ check24(b*16777215.0, 8388609.0, GMP_RNDD, 1.6777213e7);
+ check24(b*8388612.0, 8388609.0, GMP_RNDD, 8.38861e6);
+ check24(b*12582914.0, 8388610.0, GMP_RNDD, 1.2582911e7);
+ check24(b*12582913.0, 8388610.0, GMP_RNDD, 1.258291e7);
}
-void check_convergence()
+void check_convergence ()
{
- mpfr_t x, y;
+ mpfr_t x, y; int i, j;
mpfr_init2(x, 130);
mpfr_set_str_raw(x, "0.1011111101011010101000001010011111101000011100011101010011111011000011001010000000111100100111110011001010110100100001001000111001E6944");
- mpfr_init_set_ui(y, 5, 130, GMP_RNDN);
+ mpfr_init2(y, 130);
+ mpfr_set_ui(y, 5, GMP_RNDN);
mpfr_div(x, x, y, GMP_RNDD); /* exact division */
+
+ mpfr_set_prec(x, 64);
+ mpfr_set_prec(y, 64);
+ mpfr_set_str_raw(x, "0.10010010011011010100101001010111100000101110010010101E55");
+ mpfr_set_str_raw(y, "0.1E585");
+ mpfr_div(x, x, y, GMP_RNDN);
+ mpfr_set_str_raw(y, "0.10010010011011010100101001010111100000101110010010101E-529");
+ if (mpfr_cmp(x, y)) {
+ fprintf(stderr, "Error in mpfr_div for prec=64, rnd=GMP_RNDN\n");
+ printf("got "); mpfr_print_raw(x); putchar('\n');
+ printf("instead of "); mpfr_print_raw(y); putchar('\n');
+ exit(1);
+ }
+
+ for (i=32; i<=64; i+=32) {
+ mpfr_set_prec(x, i);
+ mpfr_set_prec(y, i);
+ mpfr_set_d(x, 1.0, GMP_RNDN);
+ for (j=0;j<4;j++) {
+ mpfr_set_d(y, 1.0, GMP_RNDN);
+ mpfr_div(y, x, y, j);
+ if (mpfr_cmp_ui(y, 1)) {
+ fprintf(stderr, "mpfr_div failed for x=1.0, y=1.0, prec=%u rnd=%s\n",
+ i, mpfr_print_rnd_mode(j));
+ printf("got "); mpfr_print_raw(y); putchar('\n');
+ exit(1);
+ }
+ }
+ }
+
mpfr_clear(x); mpfr_clear(y);
}
int main(int argc, char *argv[])
{
- int i, N; double n, d, e;
-#ifdef IRIX64
+ int N;
+
+#ifdef TEST
+ int i; double n, d, e;
+#ifdef __mips
/* to get denormalized numbers on IRIX64 */
union fpc_csr exp;
exp.fc_word = get_fpc_csr();
exp.fc_struct.flush = 0;
set_fpc_csr(exp.fc_word);
#endif
+#endif
N = (argc>1) ? atoi(argv[1]) : 100000;
check_float(); /* checks single precision */
check_convergence();
- check(0.0, 1.0, 1);
- check(-7.49889692246885910000e+63, 4.88168664502887320000e+306, GMP_RNDD);
- check(-1.33225773037748601769e+199, 3.63449540676937123913e+79, GMP_RNDZ);
- d = 1.0; for (i=0;i<52;i++) d *= 2.0;
- check4(4.0, d, GMP_RNDZ, 62);
- check4(1.0, 2.10263340267725788209e+187, 2, 65);
- check4(2.44394909079968374564e-150, 2.10263340267725788209e+187, 2, 65);
- /* the following tests when d is an exact power of two */
- check(9.89438396044940256501e-134, 5.93472984109987421717e-67, 2);
- check(9.89438396044940256501e-134, -5.93472984109987421717e-67, 2);
- check(-4.53063926135729747564e-308, 7.02293374921793516813e-84, 3);
- check(6.25089225176473806123e-01, -2.35527154824420243364e-230, 3);
- check(6.52308934689126000000e+15, -1.62063546601505417497e+273, 0);
- check(1.04636807108079349236e-189, 3.72295730823253012954e-292, 1);
+ check53(0.0, 1.0, GMP_RNDZ, 0.0);
+ check53(-7.4988969224688591e63, 4.8816866450288732e306, GMP_RNDD,
+ -1.5361282826510687291e-243);
+ check53(-1.33225773037748601769e+199, 3.63449540676937123913e+79, GMP_RNDZ,
+ -3.6655920045905428978e119);
+ check4(4.0, 4.503599627370496e15, GMP_RNDZ, 62, 0.0);
+ check4(1.0, 2.10263340267725788209e+187, GMP_RNDU, 65, 0.0);
+ check4(2.44394909079968374564e-150, 2.10263340267725788209e+187, GMP_RNDU,
+ 65, 0.0);
+ check53(9.89438396044940256501e-134, 5.93472984109987421717e-67, GMP_RNDU,
+ 1.6672003992376663654e-67);
+ check53(1.0, sqrt(-1.0), GMP_RNDD, sqrt(-1.0));
+ check53(sqrt(-1.0), 1.0, GMP_RNDD, sqrt(-1.0));
+ check53(2.0/0.0, 1.0, GMP_RNDD, 1.0/0.0);
+ check53(1.0, 2.0/0.0, GMP_RNDD, 0.0);
+ check53(0.0, 0.0, GMP_RNDD, sqrt(-1.0));
+ check53(1.0/0.0, 1.0/0.0, GMP_RNDD, sqrt(-1.0));
+ check53(9.89438396044940256501e-134, -5.93472984109987421717e-67, GMP_RNDU,
+ -1.6672003992376663654e-67);
+ check53(-4.53063926135729747564e-308, 7.02293374921793516813e-84, GMP_RNDD,
+ -6.4512060388748850857e-225);
+ check53(6.25089225176473806123e-01, -2.35527154824420243364e-230, GMP_RNDD,
+ -2.6540006635008291192e229);
+ check53(6.52308934689126e15, -1.62063546601505417497e273, GMP_RNDN,
+ -4.0250194961676020848e-258);
+ check53(1.04636807108079349236e-189, 3.72295730823253012954e-292, GMP_RNDZ,
+ 2.810583051186143125e102);
+#ifdef TEST
srand48(getpid());
for (i=0;i<N;i++) {
do { n = drand(); d = drand(); e = fabs(n)/fabs(d); }
/* smallest normalized is 2^(-1022), largest is 2^(1023)*(2-2^(-52)) */
- while (e>=1.7976931348623157081e308 || e<2.225073858507201383e-308);
- check(n, d, rand() % 4);
+ while (e>=MAXNORM || e<MINNORM);
+ check4(n, d, rand() % 4, 53, 0.0);
}
- exit (0);
+#endif
+ return 0;
}
diff --git a/mpfr/tests/tdiv_ui.c b/mpfr/tests/tdiv_ui.c
index bb45287af..df6294793 100644
--- a/mpfr/tests/tdiv_ui.c
+++ b/mpfr/tests/tdiv_ui.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_div_ui.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -22,47 +22,65 @@ MA 02111-1307, USA. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
+#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
-#include "time.h"
+#include "mpfr-test.h"
-extern int isnan(), getpid();
+void check _PROTO((double, unsigned long, mp_rnd_t, double));
-void check(double d, unsigned long u, unsigned char rnd)
+void check (double d, unsigned long u, mp_rnd_t rnd, double e)
{
- mpfr_t x, y; double e, f;
+ mpfr_t x, y; double f;
mpfr_init2(x, 53); mpfr_init2(y, 53);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd);
- e = d / u;
+#endif
+ if (e==0.0) e = d / u;
mpfr_set_d(x, d, rnd);
mpfr_div_ui(y, x, u, rnd);
f = mpfr_get_d(y);
if (f != e && (!isnan(f) || !isnan(e))) {
- printf("mpfr_div_ui failed for x=%1.20e, u=%lu, rnd=%d\n",d,u,rnd);
- printf("expected result is %1.20e, got %1.20e, dif=%d ulp\n",e,f,
- ulp(e,f));
+ printf("mpfr_div_ui failed for x=%1.20e, u=%lu, rnd=%s\n", d, u,
+ mpfr_print_rnd_mode(rnd));
+ printf("expected result is %1.20e, got %1.20e, dif=%d ulp\n",e,f,ulp(e,f));
exit(1);
}
- mpfr_clear(x);
+ mpfr_clear(x); mpfr_clear(y);
}
int
main(int argc, char **argv)
{
+ mpfr_t x;
+#ifdef TEST
int i; unsigned long u; double d;
srand(getpid());
- check(1.0, 3, 0);
- check(1.0, 3, 1);
- check(1.0, 3, 2);
- check(1.0, 3, 3);
- check(1.0, 2116118, 0);
for (i=0;i<1000000;i++) {
do { u = lrand48(); } while (u==0);
do { d = drand(); } while (fabs(d/u)<2.2e-307);
- check(d, u, rand() % 4);
+ check(d, u, rand() % 4, 0.0);
}
- exit (0);
+#endif
+
+ check(1.0, 3, GMP_RNDN, 3.3333333333333331483e-1);
+ check(1.0, 3, GMP_RNDZ, 3.3333333333333331483e-1);
+ check(1.0, 3, GMP_RNDU, 3.3333333333333337034e-1);
+ check(1.0, 3, GMP_RNDD, 3.3333333333333331483e-1);
+ check(1.0, 2116118, GMP_RNDN, 4.7256343927890600483e-7);
+ check(1.098612288668109782, 5, GMP_RNDN, 0.21972245773362195087);
+
+ mpfr_init2(x, 100);
+ mpfr_set_prec(x, 53);
+ mpfr_set_ui(x, 3, GMP_RNDD);
+ mpfr_log(x, x, GMP_RNDD);
+ mpfr_div_ui(x, x, 5, GMP_RNDD);
+ if (mpfr_get_d(x) != 0.21972245773362189536) {
+ fprintf(stderr, "Error in mpfr_div_ui for x=ln(3), u=5\n"); exit(1);
+ }
+ mpfr_clear(x);
+
+ return 0;
}
diff --git a/mpfr/set_dfl_rnd.c b/mpfr/tests/tdump.c
index cf8fd0175..17b000496 100644
--- a/mpfr/set_dfl_rnd.c
+++ b/mpfr/tests/tdump.c
@@ -1,6 +1,6 @@
-/* mpfr_set_default_rounding_mode -- set the default rounding mode
+/* Test file for mpfr_dump.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 2000 Free Software Foundation.
This file is part of the MPFR Library.
@@ -19,19 +19,20 @@ along with the MPFR 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 <stdlib.h>
+#include <unistd.h>
#include "gmp.h"
-#include "gmp-impl.h"
+#include "mpfr.h"
-char __gmp_default_rounding_mode = 0;
-
-void
-#if __STDC__
-mpfr_set_default_rounding_mode (char rnd_mode)
-#else
-mpfr_set_default_rounding_mode (rnd_mode)
- char rnd_mode;
-#endif
+int main()
{
- __gmp_default_rounding_mode = rnd_mode;
+ mpfr_t z;
+
+ mpfr_init2(z, 100);
+ mpfr_set_ui(z, 0, GMP_RNDN);
+ mpfr_dump(z, GMP_RNDD);
+ printf(" ^--- 0.e1 printed above is ok\n");
+ mpfr_clear(z);
+ return 0;
}
diff --git a/mpfr/tests/teq.c b/mpfr/tests/teq.c
new file mode 100644
index 000000000..b5a3ff6f0
--- /dev/null
+++ b/mpfr/tests/teq.c
@@ -0,0 +1,76 @@
+/* Test file for mpfr_eq.
+
+Copyright (C) 1999-2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <math.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+void teq (mpfr_t);
+
+void teq (mpfr_t x)
+{
+ mpfr_t y; long k, px, mx;
+
+ mpfr_init2(y, MPFR_PREC(x));
+
+ mx = (MPFR_PREC(x) - 1)/mp_bits_per_limb;
+ px = mp_bits_per_limb - 2;
+
+ for (k = 2; k < MPFR_PREC(x); k++)
+ {
+ mpfr_set(y, x, GMP_RNDN);
+
+ MPFR_MANT(y) [mx] ^= (mp_limb_t) 1 << px;
+
+ if (mpfr_eq(y, x, k) ||
+ !mpfr_eq(y, x, k - 1))
+ {
+ fprintf(stderr, "Error in eq.\n");
+ printf("x = "); mpfr_print_raw(x); printf("\n");
+ printf("y = "); mpfr_print_raw(y); printf("\n");
+ printf("k = %ld\n", k);
+ printf("mpfr_eq(y, x, k) = %d\nmpfr_eq(y, x, k - 1) = %d\n", mpfr_eq(y, x, k),mpfr_eq(y, x, k - 1));
+ mpfr_clear(x); mpfr_clear(y);
+ exit(-1);
+ }
+
+ if (px) { --px; } else { --mx; px = mp_bits_per_limb - 1; }
+ }
+ mpfr_clear(y);
+}
+
+int main()
+{
+ int j; mpfr_t x;
+
+ mpfr_init2 (x, 1000);
+
+ for (j=0;j<1000;j++) {
+ mpfr_random (x);
+ teq (x);
+ }
+ mpfr_clear (x);
+ return 0;
+}
+
diff --git a/mpfr/tests/texp.c b/mpfr/tests/texp.c
index b23d28431..a441c2ea7 100644
--- a/mpfr/tests/texp.c
+++ b/mpfr/tests/texp.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_exp.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -21,45 +21,38 @@ MA 02111-1307, USA. */
#include <math.h>
#include <stdlib.h>
+#include <unistd.h>
#include "gmp.h"
-#include "gmp-impl.h"
#include "mpfr.h"
+#include "mpfr-test.h"
-/* #define DEBUG */
-
-extern int isnan();
-extern int getpid();
+int check3 _PROTO((double, mp_rnd_t, double));
+int check_large _PROTO((double, int, mp_rnd_t));
+int check_worst_case _PROTO((double, double));
+int check_worst_cases _PROTO((void));
+void compare_exp2_exp3 _PROTO((int));
int maxu=0;
-/* returns the number of ulp's between a and b */
-int ulp(a,b) double a,b;
-{
- double eps=1.1102230246251565404e-16; /* 2^(-53) */
- b = (a-b)/a; if (b<0) b = -b;
- return (int) floor(b/eps);
-}
-
#define check(d, r) check3(d, r, 0.0)
/* returns the number of ulp of error */
-int check3(double d, unsigned char rnd, double e)
+int check3 (double d, mp_rnd_t rnd, double e)
{
mpfr_t x, y; double f; int u=0, ck=0;
-#ifdef DEBUG
- printf("d=%1.20e rnd=%d\n",d,rnd);
-#endif
mpfr_init2(x, 53); mpfr_init2(y, 53);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd);
+#endif
if (e==0.0) e = exp(d); else ck=1; /* really check */
mpfr_set_d(x, d, rnd);
mpfr_exp(y, x, rnd);
f = mpfr_get_d(y);
if (f != e && (!isnan(f) || !isnan(e))) {
- u = ulp(e,f);
+ u = ulp(e, f);
if (u<0) {
- if (u == (mp_limb_t)1<<(BITS_PER_MP_LIMB-1)) u += 1;
+ if (u == (mp_limb_t) 1 << (mp_bits_per_limb-1)) u += 1;
u=-u;
}
if (u!=0) {
@@ -71,7 +64,8 @@ int check3(double d, unsigned char rnd, double e)
}
else if (u>maxu) {
maxu=u;
- printf("mpfr_exp differs from libm.a for x=%1.20e, rnd=%d\n",d,rnd);
+ printf("mpfr_exp differs from libm.a for x=%1.20e, rnd=%s\n",d,
+ mpfr_print_rnd_mode(rnd));
printf("libm.a gave %1.20e, mpfr_exp got %1.20e, dif=%d ulp\n",e,f,u);
}
}
@@ -81,7 +75,7 @@ int check3(double d, unsigned char rnd, double e)
}
/* computes n bits of exp(d) */
-int check_large (double d, int n, char rnd)
+int check_large (double d, int n, mp_rnd_t rnd)
{
mpfr_t x; mpfr_t y;
@@ -89,12 +83,12 @@ int check_large (double d, int n, char rnd)
if (d==0.0) { /* try exp(Pi*sqrt(163)/3)-640320 */
mpfr_set_d(x, 163.0, rnd);
mpfr_sqrt(x, x, rnd);
- mpfr_pi(y, rnd);
+ mpfr_const_pi(y, rnd);
mpfr_mul(x, x, y, rnd);
mpfr_div_ui(x, x, 3, rnd);
}
else mpfr_set_d(x, d, rnd);
- mpfr_exp(y, x, rnd);
+ mpfr_exp (y, x, rnd);
if (d==0.0) {
mpfr_set_d(x, 640320.0, rnd);
mpfr_sub(y, y, x, rnd);
@@ -105,36 +99,25 @@ int check_large (double d, int n, char rnd)
putchar('\n');
printf(" ="); mpfr_print_raw(y); putchar('\n');
if (n==53) printf(" =%1.20e\n", mpfr_get_d(y));
+
mpfr_clear(x); mpfr_clear(y);
return 0;
}
/* expx is the value of exp(X) rounded towards -infinity */
-int check_worst_case(double X, double expx)
+int check_worst_case (double X, double expx)
{
mpfr_t x, y;
mpfr_init2(x, 53); mpfr_init2(y, 53);
mpfr_set_d(x, X, GMP_RNDN);
-#ifdef DEBUG
- printf("x="); mpfr_print_raw(x); putchar('\n');
-#endif
mpfr_exp(y, x, GMP_RNDD);
-#ifdef DEBUG
- printf("D(exp(x))="); mpfr_print_raw(y); putchar('\n');
-#endif
if (mpfr_get_d(y) != expx) {
fprintf(stderr, "exp(x) rounded towards -infinity is wrong\n"); exit(1);
}
mpfr_exp(x, x, GMP_RNDN);
-#ifdef DEBUG
- printf("N(exp(x))="); mpfr_print_raw(x); putchar('\n');
-#endif
mpfr_set_d(x, X, GMP_RNDN);
mpfr_exp(x, x, GMP_RNDU);
-#ifdef DEBUG
- printf("U(exp(x))="); mpfr_print_raw(x); putchar('\n');
-#endif
mpfr_add_one_ulp(y);
if (mpfr_cmp(x,y)) {
fprintf(stderr, "exp(x) rounded towards +infinity is wrong\n"); exit(1);
@@ -144,11 +127,13 @@ int check_worst_case(double X, double expx)
}
/* worst cases communicated by Jean-Michel Muller and Vincent Lefevre */
-int check_worst_cases()
+int check_worst_cases ()
{
- mpfr_t x;
+ mpfr_t x; mpfr_t y;
+
+ mpfr_init(x);
+ mpfr_set_prec (x, 53);
- mpfr_init2(x, 53);
check_worst_case(4.44089209850062517562e-16, 1.00000000000000022204);
check_worst_case(6.39488462184069720009e-14, 1.00000000000006372680);
check_worst_case(1.84741111297455401935e-12, 1.00000000000184718907);
@@ -165,15 +150,54 @@ int check_worst_cases()
check3(-1.64310953745426656203e-10, GMP_RNDN, 0.9999999998356891017792);
check3(-1.38323574826034659172e-10, GMP_RNDZ, 0.9999999998616764251835);
check3(-1.23621668465115401498e-10, GMP_RNDZ, 0.9999999998763783315425);
- mpfr_clear(x);
+
+ mpfr_set_prec (x, 601);
+ mpfr_set_str_raw (x, "0.1000100010110110101110100101000100001110000100000100010100001110110111000010010110000111010010001011110010011101111111011101010001100110111100100001101101000111111011010010011001001100110111110010010010101010100011110110010010101111000111110011111110101101100111101100001000110000000111010100001111000000011101000011111101010011010010110101101010100010000000001001000111111111011011010011010100101101111101000101100011101111000110111010010100011001100000010001111011110110111101011011000100011000010100110101001101001111110110001111101000110010011101100100101000001010011011010010110100001101110100100E0");
+ mpfr_init2 (y, 601);
+ mpfr_exp2 (y, x, GMP_RNDD);
+ mpfr_exp3 (x, x, GMP_RNDD);
+ if (mpfr_cmp (x, y)) {
+ fprintf (stderr, "mpfr_exp2 and mpfr_exp3 for prec=601\n");
+ exit (1);
+ }
+
+ mpfr_clear (x);
+ mpfr_clear (y);
return 0;
}
+void compare_exp2_exp3 (int n)
+{
+ mpfr_t x, y, z; int prec; mp_rnd_t rnd;
+
+ mpfr_init(x); mpfr_init(y); mpfr_init(z);
+ for (prec=20;prec<=n;prec++) {
+ mpfr_set_prec(x, prec); mpfr_set_prec(y, prec); mpfr_set_prec(z, prec);
+ mpfr_random(x);
+ rnd = rand() % 4;
+ mpfr_exp2 (y, x, rnd);
+ mpfr_exp3 (z, x, rnd);
+ if (mpfr_cmp(y,z)) {
+ printf("mpfr_exp2 and mpfr_exp3 disagree for rnd=%s and\nx=",
+ mpfr_print_rnd_mode(rnd));
+ mpfr_print_raw(x); putchar('\n');
+ printf("mpfr_exp2 gives "); mpfr_print_raw(y); putchar('\n');
+ printf("mpfr_exp3 gives "); mpfr_print_raw(z); putchar('\n');
+ exit(1);
+ }
+ }
+ mpfr_clear(x); mpfr_clear(y); mpfr_clear(z);
+}
+
int
main(int argc, char **argv)
{
+#ifdef TEST
int i, N, s=0, e, maxe=0; double d, lo, hi;
+#endif
+ srand(getpid());
+ compare_exp2_exp3(1000);
if (argc==4) { check_large(atof(argv[1]), atoi(argv[2]), atoi(argv[3]));
exit(1); }
check_worst_cases();
@@ -188,6 +212,35 @@ main(int argc, char **argv)
check3(3.4657359027997265421, GMP_RNDN, 32.0);
check3(3.4657359027997265421, GMP_RNDU, 32.0);
check3(3.4657359027997265421, GMP_RNDD, 31.999999999999996447);
+ check3(2.26523754332090625496e+01, GMP_RNDD, 6.8833785261699581146e9);
+ check3(1.31478962104089092122e+01, GMP_RNDZ, 5.12930793917860137299e+05);
+ check3(4.25637507920002378103e-01, GMP_RNDU, 1.53056585656161181497e+00);
+ check3(6.26551618962329307459e-16, GMP_RNDU, 1.00000000000000066613e+00);
+ check3(-3.35589513871216568383e-03, GMP_RNDD, 9.96649729583626853291e-01);
+ check3(1.95151388850007272424e+01, GMP_RNDU, 2.98756340674767792225e+08);
+ check3(2.45045953503350730784e+01, GMP_RNDN, 4.38743344916128387451e+10);
+ check3(2.58165606081678085104e+01, GMP_RNDD, 1.62925781879432281494e+11);
+ check3(-2.36539020084338638128e+01, GMP_RNDZ, 5.33630792749924762447e-11);
+ check3(2.39211946135858077866e+01, GMP_RNDU, 2.44817704330214385986e+10);
+ check3(-2.78190533055889162029e+01, GMP_RNDZ, 8.2858803483596879512e-13);
+ check3(2.64028186174889789584e+01, GMP_RNDD, 2.9281844652878973388e11);
+ check3(2.92086338843268329413e+01, GMP_RNDZ, 4.8433797301907177734e12);
+ check3(-2.46355324071459982349e+01, GMP_RNDZ, 1.9995129297760994791e-11);
+ check3(-2.23509444608605427618e+01, GMP_RNDZ, 1.9638492867489702307e-10);
+ check3(-2.41175390197331687148e+01, GMP_RNDD, 3.3564940885530624592e-11);
+ check3(2.46363885231578088053e+01, GMP_RNDU, 5.0055014282693267822e10);
+ check3(1.111263531080090984914932e2, GMP_RNDN, 1.8262572323517295459e48);
+ check3(-3.56196340354684821250e+02, GMP_RNDN, 2.0225297096141478156e-155);
+ check3(6.59678273772710895173e+02, GMP_RNDU, 3.1234469273830195529e286);
+ check3(5.13772529701934331570e+02, GMP_RNDD, 1.3445427121297197752e223);
+ check3(3.57430211008718345056e+02, GMP_RNDD, 1.6981197246857298443e155);
+ check3(3.82001814471465536371e+02, GMP_RNDU, 7.9667300591087367805e165);
+ check3(5.92396038219384422518e+02, GMP_RNDD, 1.880747529554661989e257);
+ check3(-5.02678550462488090034e+02, GMP_RNDU, 4.8919201895446217839e-219);
+ check3(5.30015757134837031117e+02, GMP_RNDD, 1.5237672861171573939e230);
+ check3(5.16239362447650933063e+02, GMP_RNDZ, 1.5845518406744492105e224);
+ check3(6.00812634798592370977e-01, GMP_RNDN, 1.823600119339019443);
+#ifdef TEST
srand48(getpid());
N = (argc==1) ? 0 : atoi(argv[1]);
lo = (argc>=3) ? atof(argv[2]) : -7.083964185e2;
@@ -202,38 +255,6 @@ main(int argc, char **argv)
if (e>maxe) maxe=e;
}
if (N) printf("mean error=%1.2e max error=%d\n", (double)s/(double)N,maxe);
- check3(2.26523754332090625496e+01, 3, 6.8833785261699581146e9);
- /* errors found in libm.a on PC under Linux */
- check3(1.31478962104089092122e+01, GMP_RNDZ, 5.12930793917860137299e+05);
- check3(4.25637507920002378103e-01, GMP_RNDU, 1.53056585656161181497e+00);
- check3(6.26551618962329307459e-16, GMP_RNDU, 1.00000000000000066613e+00);
- check3(-3.35589513871216568383e-03, GMP_RNDD, 9.96649729583626853291e-01);
- check3(1.95151388850007272424e+01, GMP_RNDU, 2.98756340674767792225e+08);
- check3(2.45045953503350730784e+01, GMP_RNDN, 4.38743344916128387451e+10);
- check3(2.58165606081678085104e+01, GMP_RNDD, 1.62925781879432281494e+11);
- check3(-2.36539020084338638128e+01, GMP_RNDZ, 5.33630792749924762447e-11);
- check3(2.39211946135858077866e+01, GMP_RNDU, 2.44817704330214385986e+10);
- check3(-2.78190533055889162029e+01, 1, 8.2858803483596879512e-13);
- /* +45 ulp, wrong side */
- check3(2.64028186174889789584e+01, 3, 2.9281844652878973388e11); /* -45 ulp*/
- check3(2.92086338843268329413e+01, 1, 4.8433797301907177734e12); /* -45 ulp*/
- check3(-2.46355324071459982349e+01, 1, 1.9995129297760994791e-11);
- /* +45 ulp, wrong side */
- check3(-2.23509444608605427618e+01, 1, 1.9638492867489702307e-10);
- /* +45 ulp, wrong side */
- check3(-2.41175390197331687148e+01, 3, 3.3564940885530624592e-11);/*-45 ulp*/
- check3(2.46363885231578088053e+01, 2, 5.0055014282693267822e10); /* +45 ulp*/
- check3(1.111263531080090984914932e2, GMP_RNDN, 1.8262572323517295459e48);
- check3(-3.56196340354684821250e+02, 0, 2.0225297096141478156e-155); /*+352 */
- check3(6.59678273772710895173e+02, 2, 3.1234469273830195529e286); /* +459 */
- check3(5.13772529701934331570e+02, 3, 1.3445427121297197752e223); /* -469 */
- check3(3.57430211008718345056e+02, 3, 1.6981197246857298443e155); /* -610 */
- check3(3.82001814471465536371e+02, 2, 7.9667300591087367805e165); /* +705 */
- check3(5.92396038219384422518e+02, 3, 1.880747529554661989e257); /* -707 */
- check3(-5.02678550462488090034e+02, 2, 4.8919201895446217839e-219); /* +708*/
- check3(5.30015757134837031117e+02, 3, 1.5237672861171573939e230); /* -709 */
- check3(5.16239362447650933063e+02, 1, 1.5845518406744492105e224); /* -710 */
- /* between 1/2 and 1 */
- check3(6.00812634798592370977e-01, 0, 1.823600119339019443); /* +1 ulp */
- exit (0);
+#endif
+ return 0;
}
diff --git a/mpfr/tests/tget_str.c b/mpfr/tests/tget_str.c
index ef70086a0..b2b505b57 100644
--- a/mpfr/tests/tget_str.c
+++ b/mpfr/tests/tget_str.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_get_str.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,45 +23,89 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
-#include <time.h>
+#ifdef TEST
+#include "mpfr-test.h"
+#endif
-void check(d, rnd) double d; unsigned char rnd;
+void check _PROTO((double, mp_rnd_t));
+void check3 _PROTO((double, mp_rnd_t, char *));
+void check_small _PROTO((void));
+
+void check (double d, mp_rnd_t rnd)
{
- mpfr_t x; char *str, str2[30]; mp_exp_t e;
+ mpfr_t x; char *str; mp_exp_t e;
mpfr_init2(x, 53);
mpfr_set_d(x, d, rnd);
str = mpfr_get_str(NULL, &e, 10, 5, x, rnd);
- mpfr_set_machine_rnd_mode(rnd);
- sprintf(str2, "%1.4e", d);
mpfr_clear(x);
free(str);
}
+void check3 (double d, mp_rnd_t rnd, char *res)
+{
+ mpfr_t x; char *str; mp_exp_t e;
+
+ mpfr_init2(x, 53);
+ mpfr_set_d(x, d, rnd);
+ str = mpfr_get_str(NULL, &e, 10, 5, x, rnd);
+ if (strcmp(str, res)) {
+ fprintf(stderr, "Error in mpfr_get_str for x=%1.20e\n", d);
+ fprintf(stderr, "got %s instead of %s\n", str, res);
+ }
+ mpfr_clear(x);
+ free(str);
+}
+
+void check_small ()
+{
+ mpfr_t x; char *s, *t; mp_exp_t e;
+
+ mpfr_init(x);
+
+ /* problem found by Fabrice Rouillier */
+ mpfr_set_prec(x, 63);
+ mpfr_set_d(x, 5e14, GMP_RNDN);
+ s = mpfr_get_str(NULL, &e, 10, 18, x, GMP_RNDU);
+ free(s);
+
+ /* bug found by Johan Vervloet */
+ mpfr_set_prec(x, 6);
+ mpfr_set_d(x, 688.0, GMP_RNDN);
+ t = mpfr_get_str(NULL, &e, 2, 4, x, GMP_RNDU);
+ if (strcmp(t, "1011") || (e!=10)) {
+ fprintf(stderr, "Error in mpfr_get_str: 688 printed up to 4 bits should give 1.011e9\ninstead of ");
+ mpfr_out_str(stderr, 2, 4, x, GMP_RNDU); putchar('\n');
+ exit(1);
+ }
+ free(t);
+
+ mpfr_clear(x);
+}
+
int
main(int argc, char **argv)
{
+#ifdef TEST
int i; double d;
srand(getpid());
- /* printf seems to round towards nearest in all cases, at least with gcc */
- check(4.059650008e-83, 0);
- check(-6.606499965302424244461355e233, 0);
- check(-7.4, 0);
- check(0.997, 0);
- check(-4.53063926135729747564e-308, 0);
- check(2.14478198760196000000e+16, 0);
- check(7.02293374921793516813e-84, 0);
- check(-6.7274500420134077e-87,0);
for (i=0;i<100000;i++) {
do { d = drand(); } while (isnan(d));
- check(d, 0);
+ check(d, GMP_RNDN);
}
- exit (0);
+#endif
+ check_small();
+ check3(4.059650008e-83, GMP_RNDN, "40597");
+ check3(-6.606499965302424244461355e233, GMP_RNDN, "-66065");
+ check3(-7.4, GMP_RNDN, "-74000");
+ check3(0.997, GMP_RNDN, "99700");
+ check3(-4.53063926135729747564e-308, GMP_RNDN, "-45306");
+ check3(2.14478198760196000000e+16, GMP_RNDN, "21448");
+ check3(7.02293374921793516813e-84, GMP_RNDN, "70229");
+ check3(-6.7274500420134077e-87, GMP_RNDN, "-67275");
+ return 0;
}
-
-
-
diff --git a/mpfr/tests/tlog.c b/mpfr/tests/tlog.c
index 899c13135..95b27610a 100644
--- a/mpfr/tests/tlog.c
+++ b/mpfr/tests/tlog.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_log.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -21,13 +21,21 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
-#include <unistd.h>
#include <math.h>
+#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
+
+double drand_log _PROTO((void));
+int check1 _PROTO((double, mp_rnd_t, double, int));
+void check3 _PROTO((double, unsigned long, mp_rnd_t));
+void check4 _PROTO((int));
+void slave _PROTO((int, int));
+void check_worst_cases _PROTO((void));
+void special _PROTO((void));
-double drand_log()
+double drand_log ()
{
double d; long int *i;
@@ -35,23 +43,24 @@ double drand_log()
do {
i[0] = lrand48();
i[1] = lrand48();
- } while ((d<1e-153)||(d>1e153)); /* to avoid underflow or overflow
+ } while ((d<1e-153) || (d>1e153)); /* to avoid underflow or overflow
in double calculus in sqrt(u*v) */
return d;
}
+#define check2(a,rnd,res) check1(a,rnd,res,1)
#define check(a,r) check2(a,r,0.0)
-
-int check1(double a, unsigned char rnd_mode, double res1, int ck)
+int check1 (double a, mp_rnd_t rnd_mode, double res1, int ck)
{
mpfr_t ta, tres;
double res2;
/* ck=1 iff res1 is certified correct */
- mpfr_set_machine_rnd_mode(rnd_mode);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
if (ck==0 && res1==0.0) res1=log(a);
- /* printf("mpfr_log working on a=%1.20e, rnd_mode=%d\n",a,rnd_mode);*/
mpfr_init2(ta, 53);
mpfr_init2(tres, 53);
mpfr_set_d(ta, a, GMP_RNDN);
@@ -61,12 +70,14 @@ int check1(double a, unsigned char rnd_mode, double res1, int ck)
if (res1!=res2 && (!isnan(res1) || !isnan(res2))) {
if (ck) {
- printf("mpfr_log failed for a=%1.20e, rnd_mode=%d\n",a,rnd_mode);
+ printf("mpfr_log failed for a=%1.20e, rnd_mode=%s\n", a,
+ mpfr_print_rnd_mode(rnd_mode));
printf("correct result is %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
exit(1);
}
else {
- printf("mpfr_log differs from libm.a for a=%1.20e, rnd_mode=%d\n",a,rnd_mode);
+ printf("mpfr_log differs from libm.a for a=%1.20e, rnd_mode=%s\n", a,
+ mpfr_print_rnd_mode(rnd_mode));
printf(" double calculus gives %1.20e\n mpfr_log gives %1.20e (%d ulp)\n",res1,res2,ulp(res1,res2));
}
}
@@ -76,13 +87,7 @@ int check1(double a, unsigned char rnd_mode, double res1, int ck)
return 0;
}
-
-void check2(double a, unsigned char rnd_mode, double res1) {
- check1(a,rnd_mode, res1,1);
- return;
-}
-
-void check3(double d, unsigned long prec, unsigned char rnd)
+void check3 (double d, unsigned long prec, mp_rnd_t rnd)
{
mpfr_t x, y;
@@ -94,7 +99,8 @@ void check3(double d, unsigned long prec, unsigned char rnd)
mpfr_clear(x); mpfr_clear(y);
}
-void check4(int N) {
+void check4 (int N)
+{
int i, max=-1, sum=0, cur;
double d;
@@ -111,7 +117,8 @@ void check4(int N) {
printf("max error : %i \t mean error : %f (in ulps)\n",max,d);
}
-void slave(int N, int p) {
+void slave (int N, int p)
+{
int i;
double d;
mpfr_t ta, tres;
@@ -127,13 +134,11 @@ void slave(int N, int p) {
printf("fin\n");
}
-
-
/* examples from Jean-Michel Muller and Vincent Lefevre
Cf http://www.ens-lyon.fr/~jmmuller/Intro-to-TMD.htm
*/
-void check_worst_cases()
+void check_worst_cases ()
{
check2(1.00089971802309629645, GMP_RNDD, 8.99313519443722736088e-04);
check2(1.00089971802309629645, GMP_RNDN, 8.99313519443722844508e-04);
@@ -216,28 +221,44 @@ void check_worst_cases()
check2(428.315247165198229595, GMP_RNDU, 6.05985948325268353187);
}
+void special ()
+{
+ mpfr_t x, y;
+
+ mpfr_init2 (x, 53);
+ mpfr_init2 (y, 53);
+ mpfr_set_ui (x, 3, GMP_RNDD);
+ mpfr_log (y, x, GMP_RNDD);
+ if (mpfr_get_d (y) != 1.09861228866810956) {
+ fprintf (stderr, "Error in mpfr_log(3) for GMP_RNDD\n");
+ exit (1);
+ }
+ mpfr_clear (x);
+ mpfr_clear (y);
+}
+
int main(int argc, char *argv[]) {
int N=0;
srand48(getpid());
if (argc==4) { /* tlog x prec rnd */
check3(atof(argv[1]), atoi(argv[2]), atoi(argv[3]));
- exit (0);
+ return 0;
}
if (argc==3) { /* tlog N p : N calculus with precision p*/
printf("Doing %d random tests in %d precision\n",atoi(argv[1]),atoi(argv[2]));
slave(atoi(argv[1]),atoi(argv[2]));
- exit (0);
+ return 0;
}
if (argc==2) { /* tlog N: N tests with random double's */
N=atoi(argv[1]);
printf("Doing %d random tests in double precision\n", N);
- /*printf("GMP_RNDN : %i, GMP_RNDZ : %i,GMP_RNDU : %i,GMP_RNDD : %i\n",GMP_RNDN, GMP_RNDZ,GMP_RNDU, GMP_RNDD); */
check4(N);
}
else {
+ special ();
check_worst_cases();
check2(1.01979300812244555452, GMP_RNDN, 1.95996734891603664741e-02);
@@ -299,5 +320,5 @@ int main(int argc, char *argv[]) {
check2(7.34302197248998461006e+43,GMP_RNDZ,1.01004909469513179942e+02);
check2(6.09969788341579732815e+00,GMP_RNDD,1.80823924264386204363e+00);
}
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tlog2.c b/mpfr/tests/tlog2.c
index 0ba296a83..afa730748 100644
--- a/mpfr/tests/tlog2.c
+++ b/mpfr/tests/tlog2.c
@@ -1,6 +1,6 @@
-/* Test file for mpfr_log2.
+/* Test file for mpfr_const_log2.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -32,13 +32,13 @@ int main(argc, argv) int argc; char *argv[];
p = (argc>1) ? atoi(argv[1]) : 53;
rnd = (argc>2) ? atoi(argv[2]) : GMP_RNDZ;
mpfr_init2(x, p);
- mpfr_log2(x, rnd);
+ mpfr_const_log2(x, rnd);
if (argc>=2) {
printf("log(2)="); mpfr_out_str(stdout, 10, 0, x, rnd); putchar('\n');
}
else if (mpfr_get_d(x) != 6.9314718055994530941e-1) {
- fprintf(stderr, "mpfr_log2 failed for prec=53\n"); exit(1);
+ fprintf(stderr, "mpfr_const_log2 failed for prec=53\n"); exit(1);
}
mpfr_clear(x);
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tmul.c b/mpfr/tests/tmul.c
index 881d0eb64..9af97a998 100644
--- a/mpfr/tests/tmul.c
+++ b/mpfr/tests/tmul.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_mul.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -21,18 +21,24 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
+#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
-#define MINNORM 2.2250738585072013831e-308 /* 2^(-1022), smallest normalized */
+void check _PROTO((double, double, mp_rnd_t, unsigned int,
+ unsigned int, unsigned int, double));
+void check53 _PROTO((double, double, mp_rnd_t, double));
+void check24 _PROTO((float, float, mp_rnd_t, float));
+void check_float _PROTO((void));
+void check_sign _PROTO((void));
/* checks that x*y gives the same results in double
and with mpfr with 53 bits of precision */
-void check(double x, double y, unsigned int rnd_mode, unsigned int px,
-unsigned int py, unsigned int pz, double res)
+void check (double x, double y, mp_rnd_t rnd_mode, unsigned int px,
+ unsigned int py, unsigned int pz, double res)
{
- double z1,z2; mpfr_t xx,yy,zz;
+ double z1, z2; mpfr_t xx, yy, zz;
mpfr_init2(xx, px);
mpfr_init2(yy, py);
@@ -40,15 +46,37 @@ unsigned int py, unsigned int pz, double res)
mpfr_set_d(xx, x, rnd_mode);
mpfr_set_d(yy, y, rnd_mode);
mpfr_mul(zz, xx, yy, rnd_mode);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
z1 = (res==0.0) ? x*y : res;
z2 = mpfr_get_d(zz);
if (z1!=z2 && (z1>=MINNORM || z1<=-MINNORM)) {
printf("mpfr_mul ");
if (res==0.0) printf("differs from libm.a"); else printf("failed");
- printf(" for x=%1.20e y=%1.20e with rnd_mode=%u\n",x,y,rnd_mode);
- mpfr_print_raw(zz); putchar('\n');
- printf("libm.a gives %1.20e, mpfr_mul gives %1.20e\n",z1,z2);
+ printf(" for x=%1.20e y=%1.20e with rnd_mode=%s\n", x, y,
+ mpfr_print_rnd_mode(rnd_mode));
+ printf("libm.a gives %1.20e, mpfr_mul gives %1.20e\n", z1, z2);
+ if (res!=0.0) exit(1);
+ }
+ mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
+}
+
+void check53 (double x, double y, mp_rnd_t rnd_mode, double z1)
+{
+ double z2; mpfr_t xx, yy, zz;
+
+ mpfr_init2(xx, 53);
+ mpfr_init2(yy, 53);
+ mpfr_init2(zz, 53);
+ mpfr_set_d(xx, x, rnd_mode);
+ mpfr_set_d(yy, y, rnd_mode);
+ mpfr_mul(zz, xx, yy, rnd_mode);
+ z2 = mpfr_get_d(zz);
+ if (z1!=z2 && (!isnan(z1) || !isnan(z2))) {
+ printf("mpfr_mul failed for x=%1.20e y=%1.20e with rnd_mode=%s\n",
+ x, y, mpfr_print_rnd_mode(rnd_mode));
+ printf("libm.a gives %1.20e, mpfr_mul gives %1.20e\n", z1, z2);
exit(1);
}
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
@@ -56,9 +84,9 @@ unsigned int py, unsigned int pz, double res)
/* checks that x*y gives the same results in double
and with mpfr with 24 bits of precision */
-void check24(float x, float y, unsigned int rnd_mode, float res)
+void check24 (float x, float y, mp_rnd_t rnd_mode, float z1)
{
- float z1,z2; mpfr_t xx,yy,zz;
+ float z2; mpfr_t xx, yy, zz;
mpfr_init2(xx, 24);
mpfr_init2(yy, 24);
@@ -66,38 +94,39 @@ void check24(float x, float y, unsigned int rnd_mode, float res)
mpfr_set_d(xx, x, rnd_mode);
mpfr_set_d(yy, y, rnd_mode);
mpfr_mul(zz, xx, yy, rnd_mode);
- mpfr_set_machine_rnd_mode(rnd_mode);
- z1 = (res==0.0) ? x*y : res;
z2 = (float) mpfr_get_d(zz);
if (z1!=z2) {
- printf("mpfr_mul ");
- if (res==0.0) printf("differs from libm.a"); else printf("failed");
- printf(" for x=%1.10e y=%1.10e with prec=24 and rnd_mode=%u\n",x,y,rnd_mode);
- printf("libm.a gives %1.10e, mpfr_mul gives %1.10e\n",z1,z2);
- if (res!=0.0) exit(1);
+ printf("mpfr_mul failed for x=%1.0f y=%1.0f with prec=24 and rnd_mode=%s\n", x, y, mpfr_print_rnd_mode(rnd_mode));
+ printf("libm.a gives %1.0f, mpfr_mul gives %1.0f\n", z1, z2);
+ exit(1);
}
mpfr_clear(xx); mpfr_clear(yy); mpfr_clear(zz);
}
/* the following examples come from the paper "Number-theoretic Test
Generation for Directed Rounding" from Michael Parks, Table 1 */
-void check_float()
+void check_float ()
{
- int i;
- for (i=0;i<4;i++) {
- if (i!=2) {
- check24(8388609.0, 8388609.0, i, 0.0);
- check24(16777213.0, 8388609.0, i, 0.0);
- check24(8388611.0, 8388609.0, i, 0.0);
- check24(12582911.0, 8388610.0, i, 0.0);
- check24(12582914.0, 8388610.0, i, 0.0);
- check24(13981013.0, 8388611.0, i, 0.0);
- check24(11184811.0, 8388611.0, i, 0.0);
- check24(11184810.0, 8388611.0, i, 0.0);
- check24(13981014.0, 8388611.0, i, 0.0);
- }
- }
- i=GMP_RNDU;
+ check24(8388609.0, 8388609.0, GMP_RNDN, 70368760954880.0);
+ check24(16777213.0, 8388609.0, GMP_RNDN, 140737479966720.0);
+ check24(8388611.0, 8388609.0, GMP_RNDN, 70368777732096.0);
+ check24(12582911.0, 8388610.0, GMP_RNDN, 105553133043712.0);
+ check24(12582914.0, 8388610.0, GMP_RNDN, 105553158209536.0);
+ check24(13981013.0, 8388611.0, GMP_RNDN, 117281279442944.0);
+ check24(11184811.0, 8388611.0, GMP_RNDN, 93825028587520.0);
+ check24(11184810.0, 8388611.0, GMP_RNDN, 93825020198912.0);
+ check24(13981014.0, 8388611.0, GMP_RNDN, 117281287831552.0);
+
+ check24(8388609.0, 8388609.0, GMP_RNDZ, 70368760954880.0);
+ check24(16777213.0, 8388609.0, GMP_RNDZ, 140737471578112.0);
+ check24(8388611.0, 8388609.0, GMP_RNDZ, 70368777732096.0);
+ check24(12582911.0, 8388610.0, GMP_RNDZ, 105553124655104.0);
+ check24(12582914.0, 8388610.0, GMP_RNDZ, 105553158209536.0);
+ check24(13981013.0, 8388611.0, GMP_RNDZ, 117281271054336.0);
+ check24(11184811.0, 8388611.0, GMP_RNDZ, 93825028587520.0);
+ check24(11184810.0, 8388611.0, GMP_RNDZ, 93825011810304.0);
+ check24(13981014.0, 8388611.0, GMP_RNDZ, 117281287831552.0);
+
check24(8388609.0, 8388609.0, GMP_RNDU, 70368769343488.0);
check24(16777213.0, 8388609.0, GMP_RNDU, 140737479966720.0);
check24(8388611.0, 8388609.0, GMP_RNDU, 70368786120704.0);
@@ -107,10 +136,20 @@ void check_float()
check24(11184811.0, 8388611.0, GMP_RNDU, 93825036976128.0);
check24(11184810.0, 8388611.0, GMP_RNDU, 93825020198912.0);
check24(13981014.0, 8388611.0, GMP_RNDU, 117281296220160.0);
+
+ check24(8388609.0, 8388609.0, GMP_RNDD, 70368760954880.0);
+ check24(16777213.0, 8388609.0, GMP_RNDD, 140737471578112.0);
+ check24(8388611.0, 8388609.0, GMP_RNDD, 70368777732096.0);
+ check24(12582911.0, 8388610.0, GMP_RNDD, 105553124655104.0);
+ check24(12582914.0, 8388610.0, GMP_RNDD, 105553158209536.0);
+ check24(13981013.0, 8388611.0, GMP_RNDD, 117281271054336.0);
+ check24(11184811.0, 8388611.0, GMP_RNDD, 93825028587520.0);
+ check24(11184810.0, 8388611.0, GMP_RNDD, 93825011810304.0);
+ check24(13981014.0, 8388611.0, GMP_RNDD, 117281287831552.0);
}
/* check sign of result */
-void check_sign()
+void check_sign ()
{
mpfr_t a, b;
@@ -124,42 +163,44 @@ void check_sign()
mpfr_clear(a); mpfr_clear(b);
}
-int main(argc,argv) int argc; char *argv[];
+int main (int argc, char *argv[])
{
- double x,y,z; int i,prec,rnd_mode;
+#ifdef TEST
+ double x, y, z; int i, prec, rnd_mode;
+#endif
check_float();
- check(6.9314718055994530941514e-1, 0.0, GMP_RNDZ, 53, 53, 53, 0.0);
- check(0.0, 6.9314718055994530941514e-1, GMP_RNDZ, 53, 53, 53, 0.0);
+ check53(0.0, 1.0/0.0, GMP_RNDN, 0.0/0.0);
+ check53(1.0, 1.0/0.0, GMP_RNDN, 1.0/0.0);
+ check53(-1.0, 1.0/0.0, GMP_RNDN, -1.0/0.0);
+ check53(0.0/0.0, 0.0, GMP_RNDN, 0.0/0.0);
+ check53(1.0, 0.0/0.0, GMP_RNDN, 0.0/0.0);
+ check53(6.9314718055994530941514e-1, 0.0, GMP_RNDZ, 0.0);
+ check53(0.0, 6.9314718055994530941514e-1, GMP_RNDZ, 0.0);
+ check_sign();
+ check53(-4.165000000e4, -0.00004801920768307322868063274915, GMP_RNDN, 2.0);
+ check53(2.71331408349172961467e-08, -6.72658901114033715233e-165, GMP_RNDZ,
+ -1.8251348697787782844e-172);
+ check53(0.31869277231188065, 0.88642843322303122, GMP_RNDZ,
+ 2.8249833483992453642e-1);
+ check(8.47622108205396074254e-01, 3.24039313247872939883e-01, GMP_RNDU,
+ 28, 45, 1, 0.5);
+ check(2.63978122803639081440e-01, 6.8378615379333496093e-1, GMP_RNDN,
+ 34, 23, 31, 0.180504585267044603);
+ check(1.0, 0.11835170935876249132, GMP_RNDU, 6, 41, 36, 0.1183517093595583);
+ check53(67108865.0, 134217729.0, GMP_RNDN, 9.007199456067584e15);
+ check(1.37399642157394197284e-01, 2.28877275604219221350e-01, GMP_RNDN,
+ 49, 15, 32, 0.0314472340833162888);
+ check(4.03160720978664954828e-01, 5.85483042917246621073e-01, GMP_RNDZ,
+ 51, 22, 32, 0.2360436821472831);
+ check(3.90798504668055102229e-14, 9.85394674650308388664e-04, GMP_RNDN,
+ 46, 22, 12, 0.385027296503914762e-16);
+ check(4.58687081072827851358e-01, 2.20543551472118792844e-01, GMP_RNDN,
+ 49, 3, 1, 0.125);
+#ifdef TEST
+ srand48(getpid());
prec = (argc<2) ? 53 : atoi(argv[1]);
rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
- check_sign();
- check(-4.165000000e4, -0.00004801920768307322868063274915, GMP_RNDN,
- 53, 53, 53, 2.0);
- check(2.71331408349172961467e-08, -6.72658901114033715233e-165,
- GMP_RNDZ, 53, 53, 53, 0.0);
- x=0.31869277231188065; y=0.88642843322303122;
- check(x, y, GMP_RNDZ, 53, 53, 53, 0.0);
- x=8.47622108205396074254e-01; y=3.24039313247872939883e-01;
- check(x, y, GMP_RNDU, 28, 45, 1, 0.5);
- x=2.63978122803639081440e-01;
- y=5736014.0/8388608.0; /* 6.83786096444222835089e-01; */
- check(x, y, GMP_RNDN, 34, 23, 31, 0.180504585267044603);
- x=9.84891017624509146344e-01; /* rounded to 1.0 with prec=6 */
- x=1.0;
- y=1.18351709358762491320e-01;
- check(x, y, GMP_RNDU, 6, 41, 36, 0.1183517093595583);
- /* the following checks that rounding to nearest sets the last
- bit to zero in case of equal distance */
- check(67108865.0, 134217729.0, GMP_RNDN, 53, 53, 53, 0.0);
- x=1.37399642157394197284e-01; y=2.28877275604219221350e-01;
- check(x, y, GMP_RNDN, 49, 15, 32, 0.0314472340833162888);
- x=4.03160720978664954828e-01; y=5.85483042917246621073e-01;
- check(x, y, GMP_RNDZ, 51, 22, 32, 0.2360436821472831);
- x=3.90798504668055102229e-14; y=9.85394674650308388664e-04;
- check(x, y, GMP_RNDN, 46, 22, 12, 0.385027296503914762e-16);
- x=4.58687081072827851358e-01; y=2.20543551472118792844e-01;
- check(x, y, GMP_RNDN, 49, 3, 1, 0.125);
for (i=0;i<1000000;) {
x = drand();
y = drand();
@@ -170,6 +211,7 @@ int main(argc,argv) int argc; char *argv[];
prec, prec, prec, 0.0);
}
}
- exit (0);
+#endif
+ return 0;
}
diff --git a/mpfr/tests/tmul_2exp.c b/mpfr/tests/tmul_2exp.c
index c47031313..f5b9490c6 100644
--- a/mpfr/tests/tmul_2exp.c
+++ b/mpfr/tests/tmul_2exp.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_mul_2exp.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -24,6 +24,7 @@ MA 02111-1307, USA. */
#include <time.h>
#include "gmp.h"
#include "mpfr.h"
+#include "mpfr-impl.h"
/* checks that x*y gives the same results in double
and with mpfr with 53 bits of precision */
@@ -31,19 +32,42 @@ MA 02111-1307, USA. */
int
main(argc,argv) int argc; char *argv[];
{
- double x, z; mpfr_t w;
+ double x, z; mpfr_t w; unsigned long k;
mpfr_init2(w, 53);
- srand48(time(NULL));
- x = drand48();
- mpfr_set_d(w, x, 0);
+ mpfr_set_d(w, 1.0/0.0, 0);
mpfr_mul_2exp(w, w, 10, GMP_RNDZ);
- if (x != (z = mpfr_get_d(w)/1024))
- {
- fprintf(stderr, "%f != %f\n", x, z);
- exit (1);
- }
- exit (0);
+ if (!MPFR_IS_INF(w)) { fprintf(stderr, "Inf != Inf"); exit(-1); }
+
+ mpfr_set_d(w, 0.0/0.0, 0);
+ mpfr_mul_2exp(w, w, 10, GMP_RNDZ);
+ if (!MPFR_IS_NAN(w)) { fprintf(stderr, "NaN != NaN"); exit(-1); }
+
+ for (k = 0; k < 100000; k++) {
+ srand48(time(NULL));
+ x = drand48();
+ mpfr_set_d(w, x, 0);
+ mpfr_mul_2exp(w, w, 10, GMP_RNDZ);
+ if (x != (z = mpfr_get_d(w)/1024))
+ {
+ fprintf(stderr, "%f != %f\n", x, z);
+ return (-1);
+ }
+
+ mpfr_set_d(w, x, 0);
+ mpfr_div_2exp(w, w, 10, GMP_RNDZ);
+ if (x != (z = mpfr_get_d(w)*1024))
+ {
+ fprintf(stderr, "%f != %f\n", x, z);
+ mpfr_clear(w);
+ return (-1);
+ }
+ }
+
+
+
+ mpfr_clear(w);
+ return (0);
}
diff --git a/mpfr/tests/tmul_ui.c b/mpfr/tests/tmul_ui.c
index fab7d50fd..e226b35ba 100644
--- a/mpfr/tests/tmul_ui.c
+++ b/mpfr/tests/tmul_ui.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_mul_ui.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -22,9 +22,8 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"
-#include "gmp-impl.h"
#include "mpfr.h"
-#include "time.h"
+#include "mpfr-impl.h"
int
main(int argc, char **argv)
@@ -36,7 +35,7 @@ main(int argc, char **argv)
/* checks that result is normalized */
mpfr_set_d(y, 6.93147180559945286227e-01, GMP_RNDZ);
mpfr_mul_ui(x, y, 1, GMP_RNDZ);
- if (MANT(x)[PREC(x)/BITS_PER_MP_LIMB] >> (BITS_PER_MP_LIMB-1) == 0) {
+ if (MPFR_MANT(x)[MPFR_PREC(x)/mp_bits_per_limb] >> (mp_bits_per_limb-1) == 0) {
fprintf(stderr, "Error in mpfr_mul_ui: result not normalized\n");
exit(1);
}
@@ -47,6 +46,25 @@ main(int argc, char **argv)
exit(1);
}
+
+ mpfr_set_d(x, 1.0/0.0, GMP_RNDZ);
+ mpfr_mul_ui(x, x, 3, GMP_RNDU);
+ if (mpfr_get_d(x) != 1.0/0.0) {
+ fprintf(stderr, "Error in mpfr_mul_ui: Inf*3 does not give Inf\n"); exit(1);
+ }
+
+ mpfr_set_d(x, -1.0/0.0, GMP_RNDZ);
+ mpfr_mul_ui(x, x, 3, GMP_RNDU);
+ if (mpfr_get_d(x) != -1.0/0.0) {
+ fprintf(stderr, "Error in mpfr_mul_ui: -Inf*3 does not give -Inf\n"); exit(1);
+ }
+
+ mpfr_set_d(x, 0.0/0.0, GMP_RNDZ);
+ mpfr_mul_ui(x, x, 3, GMP_RNDU);
+ if (!isnan(mpfr_get_d(x))) {
+ fprintf(stderr, "Error in mpfr_mul_ui: NaN*3 does not give NaN\n"); exit(1);
+ }
+
mpfr_set_d(x, 1.0/3.0, GMP_RNDZ);
mpfr_mul_ui(x, x, 3, GMP_RNDU);
if (mpfr_get_d(x) != 1.0) {
@@ -57,10 +75,34 @@ main(int argc, char **argv)
mpfr_set_d(x, -2.0, GMP_RNDZ);
mpfr_set_d(y, 3.0, GMP_RNDZ);
mpfr_mul_ui(x, y, 4, GMP_RNDZ);
- if (SIGN(x) != 1) {
- fprintf(stderr, "Error in mpfr_mul_ui: 4*3.0 does not give a positive result\n"); exit(1);
+ if (mpfr_cmp_ui(x, 0) <= 0) {
+ fprintf(stderr, "Error in mpfr_mul_ui: 4*3.0 does not give a positive result:\n");
+ mpfr_print_raw(x); putchar('\n');
+ printf("mpfr_cmp_ui(x, 0) = %d\n", mpfr_cmp_ui(x, 0));
+ exit(1);
+ }
+
+ mpfr_set_prec(x, 9);
+ mpfr_set_prec(y, 9);
+ mpfr_set_str_raw(y,"0.100001111E9");
+ mpfr_mul_ui(x, y, 1335, GMP_RNDN);
+ mpfr_set_str_raw(y,"0.100111001E19");
+ if (mpfr_cmp(x, y)) {
+ fprintf(stderr, "Error in mul_ui for 1335*(0.100001111E9)\n"); exit(1);
+ }
+
+ mpfr_set_prec(y, 100);
+ mpfr_set_prec(x, 100);
+ /* y = 1199781142214086656 */
+ mpfr_set_str_raw(y, "0.1000010100110011110101001011110010101111000100001E61");
+ mpfr_mul_ui(x, y, 121, GMP_RNDD);
+ /* 121*y = 145173518207904485376, representable exactly */
+ mpfr_set_str_raw(y, "0.1111101111010101111111100011010010111010111110110011001E67");
+ if (mpfr_cmp(x, y)) {
+ printf("Error for 121*y: expected result is:\n");
+ mpfr_print_raw(y); putchar('\n');
}
mpfr_clear(x); mpfr_clear(y);
- exit (0);
+ return(0);
}
diff --git a/mpfr/tests/tout_str.c b/mpfr/tests/tout_str.c
index cfb8a822e..2b5cc4c66 100644
--- a/mpfr/tests/tout_str.c
+++ b/mpfr/tests/tout_str.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_out_str.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -24,36 +24,90 @@ MA 02111-1307, USA. */
#include <stdlib.h>
#include <string.h>
#include <time.h>
+#include <unistd.h>
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
FILE *fout;
#define check(d,r,b) check4(d,r,b,53)
-void check4(d, rnd, base, prec) double d; unsigned char rnd; int base, prec;
+void check4 _PROTO((double, mp_rnd_t, int, int));
+void check_large _PROTO((void));
+
+void check4(double d, mp_rnd_t rnd, int base, int prec)
{
mpfr_t x;
mpfr_init2(x, prec);
mpfr_set_d(x, d, rnd);
- mpfr_set_machine_rnd_mode(rnd);
fprintf(fout, "%1.19e base %d rnd %d:\n ", d, base, rnd);
mpfr_out_str(fout, base, (base==2) ? prec : 0, x, rnd);
fputc('\n', fout);
mpfr_clear(x);
}
+void check_large ()
+{
+ mpfr_t x; mp_exp_t e; char *s;
+
+ mpfr_init(x);
+
+ mpfr_set_prec(x, 7);
+ mpfr_set_str_raw(x, "0.1010101E10");
+ s = mpfr_get_str(NULL, &e, 10, 2, x, GMP_RNDU);
+ free(s);
+
+ /* checks rounding of negative numbers */
+ mpfr_set_d(x, -1.5, GMP_RNDN);
+ s = mpfr_get_str(NULL, &e, 10, 1, x, GMP_RNDD);
+ if (strcmp(s, "-2")) {
+ fprintf(stderr, "Error in mpfr_get_str for x=-1.5 and rnd=GMP_RNDD\n");
+ free(s); mpfr_clear(x);
+ exit(1);
+ }
+ free(s);
+
+ s = mpfr_get_str(NULL, &e, 10, 1, x, GMP_RNDU);
+ if (strcmp(s, "-1")) {
+ fprintf(stderr, "Error in mpfr_get_str for x=-1.5 and rnd=GMP_RNDU\n");
+ free(s);
+ mpfr_clear(x);
+ exit(1);
+ }
+
+ free(s);
+
+ /* bug found by Jean-Pierre Merlet, produced error in mpfr_get_str */
+ mpfr_set_prec(x, 128);
+ mpfr_set_str_raw(x, "0.10111001100110011001100110011001100110011001100110011001100110011001100110011001100110011001100110011001100110011001100110011010E3");
+ s = mpfr_get_str(NULL, &e, 10, 0, x, GMP_RNDU);
+ free(s);
+
+ mpfr_set_prec(x, 381);
+ mpfr_set_str_raw(x, "0.111111111111111111111111111111111111111111111111111111111111111111101110110000100110011101101101001010111000101111000100100011110101010110101110100000010100001000110100000100011111001000010010000010001010111001011110000001110010111101100001111000101101100000010110000101100100000101010110010110001010100111001111100011100101100000100100111001100010010011110011011010110000001000010");
+ s = mpfr_get_str (NULL, &e, 10, 0, x, GMP_RNDD);
+ if (e != 0) {
+ fprintf(stderr, "Error in mpfr_get_str for x=0.999999..., exponent is %d instead of 0\n", (int) e);
+ exit(1);
+ }
+ free(s);
+
+ mpfr_clear(x);
+}
+
int
main(int argc, char **argv)
{
- int i,N=100,r,p; double d;
+ int i,N=10000,r,p; double d;
+ check_large();
/* with no argument: prints to /dev/null,
tout_str N: prints N tests to stdout */
if (argc==1) fout=fopen("/dev/null", "w");
else { fout=stdout; N=atoi(argv[1]); }
+ check(-1.37247529013405550000e+15, GMP_RNDN, 7);
check(-1.5674376729569697500e+15, GMP_RNDN, 19);
check(-5.71262771772792640000e-79, GMP_RNDU, 16);
check(-0.0, GMP_RNDU, 7);
@@ -78,8 +132,5 @@ main(int argc, char **argv)
p = 2 + rand()%35;
check(d, r, p);
}
- exit (0);
+ return 0;
}
-
-
-
diff --git a/mpfr/tests/tpi.c b/mpfr/tests/tpi.c
index 487983c5f..692008ebf 100644
--- a/mpfr/tests/tpi.c
+++ b/mpfr/tests/tpi.c
@@ -1,6 +1,6 @@
-/* Test file for mpfr_pi.
+/* Test file for mpfr_const_pi.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -32,13 +32,13 @@ int main(argc, argv) int argc; char *argv[];
p = (argc>1) ? atoi(argv[1]) : 53;
rnd = (argc>2) ? atoi(argv[2]) : GMP_RNDZ;
mpfr_init2(x, p);
- mpfr_pi(x, rnd);
+ mpfr_const_pi(x, rnd);
if (argc>=2) {
printf("Pi="); mpfr_out_str(stdout, 10, 0, x, rnd); putchar('\n');
}
else if (mpfr_get_d(x) != 3.141592653589793116) {
- fprintf(stderr, "mpfr_pi failed for prec=53\n"); exit(1);
+ fprintf(stderr, "mpfr_const_pi failed for prec=53\n"); exit(1);
}
mpfr_clear(x);
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tpow.c b/mpfr/tests/tpow.c
new file mode 100644
index 000000000..4bad9b73f
--- /dev/null
+++ b/mpfr/tests/tpow.c
@@ -0,0 +1,66 @@
+/* Test file for mpfr_pow and mpfr_pow_ui.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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"
+
+void check_pow_ui (void);
+
+void check_pow_ui ()
+{
+ mpfr_t a, b;
+
+ mpfr_init2 (a, 53);
+ mpfr_init2 (b, 53);
+
+ /* check in-place operations */
+ mpfr_set_d (b, 0.6926773, GMP_RNDN);
+ mpfr_pow_ui (a, b, 10, GMP_RNDN);
+ mpfr_pow_ui (b, b, 10, GMP_RNDN);
+ if (mpfr_cmp (a, b)) {
+ fprintf (stderr, "Error for mpfr_pow_ui (b, b, ...)\n"); exit (1);
+ }
+
+ /* check large exponents */
+ mpfr_set_d (b, 1, GMP_RNDN);
+ mpfr_pow_ui (a, b, (unsigned long) 4294967295UL, GMP_RNDN);
+
+ mpfr_set_d (a, -1.0/0.0, GMP_RNDN);
+ mpfr_pow_ui (a, a, 4049053855UL, GMP_RNDN);
+ if (mpfr_get_d (a) != -1.0/0.0) {
+ fprintf (stderr, "Error for (-Inf)^4049053855\n"); exit (1);
+ }
+
+ mpfr_set_d (a, -1.0/0.0, GMP_RNDN);
+ mpfr_pow_ui (a, a, 30002752UL, GMP_RNDN);
+ if (mpfr_get_d (a) != 1.0/0.0) {
+ fprintf (stderr, "Error for (-Inf)^30002752\n"); exit (1);
+ }
+
+ mpfr_clear (a);
+ mpfr_clear (b);
+}
+
+int main ()
+{
+ check_pow_ui ();
+ return 0;
+}
diff --git a/mpfr/tests/trandom.c b/mpfr/tests/trandom.c
new file mode 100644
index 000000000..7229bd46f
--- /dev/null
+++ b/mpfr/tests/trandom.c
@@ -0,0 +1,182 @@
+/* Test file for the various mpfr_random fonctions.
+
+Copyright (C) 1999-2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <math.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+void test_random (unsigned long, unsigned long, int);
+void test_random2 (unsigned long, unsigned long, int);
+void test_urandomb (unsigned long, unsigned long, int);
+
+void test_random (unsigned long nbtests, unsigned long prec, int verbose)
+{
+ mpfr_t x;
+ int *tab, size_tab, k;
+ double d, av = 0, var = 0, chi2 = 0, th;
+
+ mpfr_init2(x, prec);
+
+ size_tab = (nbtests < 1000 ? nbtests / 50 : 20);
+ tab = (int *) malloc (size_tab * sizeof(int));
+ for (k = 0; k < size_tab; ++k) tab[k] = 0;
+
+ for (k = 0; k < nbtests; k++) {
+ mpfr_random(x);
+ d = mpfr_get_d(x); av += d; var += d*d;
+ tab[(int)(size_tab * d)]++;
+ }
+
+ mpfr_clear(x);
+ if (!verbose) { free(tab); return; }
+
+ av /= nbtests;
+ var = (var /nbtests) - av*av;
+
+ th = (double)nbtests / size_tab;
+
+ printf("Average = %.5f\nVariance = %.5f\n", av, var);
+ printf("Repartition for random. Each integer should be close to %d.\n",
+ (int)th);
+
+ for (k = 0; k < size_tab; k++) {
+ chi2 += (tab[k] - th) * (tab[k] - th) / th;
+ printf("%d ", tab[k]);
+ if (((k+1) & 7) == 0) printf("\n");
+ }
+
+ printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n\n",
+ size_tab - 1, chi2);
+
+ printf("\n");
+
+ free(tab);
+ return;
+}
+
+void test_random2 (unsigned long nbtests, unsigned long prec, int verbose)
+{
+ mpfr_t x;
+ int *tab, size_tab, k;
+ double d, av = 0, var = 0, chi2 = 0, th;
+
+ mpfr_init2(x, prec);
+
+ size_tab = (nbtests < 1000 ? nbtests / 50 : 20);
+ tab = (int *) malloc (size_tab * sizeof(int));
+ for (k = 0; k < size_tab; ++k) tab[k] = 0;
+
+ for (k = 0; k < nbtests; k++) {
+ mpfr_random2 (x, MPFR_ABSSIZE(x), 0);
+ d = mpfr_get_d(x); av += d; var += d*d;
+ if (d < 1)
+ tab[(int)(size_tab * d)]++;
+ }
+
+ mpfr_clear(x);
+ if (!verbose) { free(tab); return; }
+
+ av /= nbtests;
+ var = (var /nbtests) - av*av;
+
+ th = (double)nbtests / size_tab;
+ printf("Average = %.5f\nVariance = %.5f\n", av, var);
+ printf("Repartition for random2 (taking only values < 1 into account.\n");
+
+ for (k = 0; k < size_tab; k++) {
+ chi2 += (tab[k] - th) * (tab[k] - th) / th;
+ printf("%d ", tab[k]);
+ if (((k+1) & 7) == 0) printf("\n");
+ }
+
+ printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n\n",
+ size_tab - 1, chi2);
+
+ free(tab);
+ return;
+}
+
+void test_urandomb (unsigned long nbtests, unsigned long prec, int verbose)
+{
+ mpfr_t x;
+ int *tab, size_tab, k;
+ gmp_randstate_t state;
+ double d, av = 0, var = 0, chi2 = 0, th;
+
+ mpfr_init2(x, prec);
+
+ size_tab = (nbtests < 1000 ? nbtests / 50 : 20);
+ tab = (int *) malloc (size_tab * sizeof(int));
+ for (k = 0; k < size_tab; ++k) tab[k] = 0;
+
+ gmp_randinit(state, GMP_RAND_ALG_LC, 128);
+ gmp_randseed_ui(state, (unsigned long int)time(NULL));
+
+ for (k = 0; k < nbtests; k++) {
+ mpfr_urandomb(x, state);
+ d = mpfr_get_d(x); av += d; var += d*d;
+ tab[(int)(size_tab * d)]++;
+ }
+
+ mpfr_clear(x);
+ gmp_randclear(state);
+ if (!verbose) { free(tab); return; }
+
+ av /= nbtests;
+ var = (var /nbtests) - av*av;
+
+ th = (double)nbtests / size_tab;
+ printf("Average = %.5f\nVariance = %.5f\n", av, var);
+ printf("Repartition for urandomb. Each integer should be close to %d.\n",
+ (int)th);
+
+ for (k = 0; k < size_tab; k++) {
+ chi2 += (tab[k] - th) * (tab[k] - th) / th;
+ printf("%d ", tab[k]);
+ if (((k+1) & 7) == 0) printf("\n");
+ }
+
+ printf("\nChi2 statistics value (with %d degrees of freedom) : %.5f\n\n",
+ size_tab - 1, chi2);
+
+ free(tab);
+ return;
+}
+
+int main (int argc, char **argv)
+{
+ unsigned long nbtests, prec; int verbose = 0;
+
+ if (argc > 1) verbose = 1;
+
+ if (argc == 1) { nbtests = 10000; } else nbtests = atoi(argv[1]);
+ if (argc <= 2) { prec = 1000; } else prec = atoi(argv[2]);
+
+ test_random(nbtests, prec, verbose);
+ test_random2(nbtests, prec, verbose);
+ test_urandomb(nbtests, prec, verbose);
+
+ return 0;
+}
+
diff --git a/mpfr/tests/tround.c b/mpfr/tests/tround.c
index 3ccd6ff7f..5aff9413e 100644
--- a/mpfr/tests/tround.c
+++ b/mpfr/tests/tround.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_round.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -24,22 +24,21 @@ MA 02111-1307, USA. */
#include "gmp.h"
#include "mpfr.h"
-int
-main()
+int main()
{
mpfr_t x;
/* checks that rounds to nearest sets the last
bit to zero in case of equal distance */
mpfr_init2(x, 2);
- mpfr_set_d(x, 5.0, 0);
+ mpfr_set_d(x, 5.0, GMP_RNDN);
if (mpfr_get_d(x) != 4.0) { printf("Error in tround: got %1.1f instead of 4.0\n",mpfr_get_d(x)); }
- mpfr_set_d(x, 0.00098539467465030839, 0);
+ mpfr_set_d(x, 0.00098539467465030839, GMP_RNDN);
mpfr_set_d(x, 9.84891017624509146344e-01, GMP_RNDU);
if (mpfr_get_d(x) != 1.0) { printf("Error in tround: got %f instead of 1.0\n",mpfr_get_d(x)); exit(1); }
mpfr_clear(x);
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tset_d.c b/mpfr/tests/tset_d.c
index c542dc1e2..2ab6eb1b2 100644
--- a/mpfr/tests/tset_d.c
+++ b/mpfr/tests/tset_d.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_set_d and mpfr_get_d.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -22,24 +22,17 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
-#ifdef IRIX64
-#include <sys/fpu.h>
-#endif
#include "gmp.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
+#include "mpfr-test.h"
extern int isnan();
int
-main(int argc, char **argv)
+main (int argc, char **argv)
{
- fprintf(stderr, "Test case tset_d disabled\n");
- exit(0); /* THIS TEST CASE IS NOT WORKING */
-
-#if 0
mpfr_t x,y,z; unsigned long k,n; double d, dd;
-#ifdef IRIX64
+#ifdef __mips
/* to get denormalized numbers on IRIX64 */
union fpc_csr exp;
exp.fc_word = get_fpc_csr();
@@ -47,13 +40,15 @@ main(int argc, char **argv)
set_fpc_csr(exp.fc_word);
#endif
+ mpfr_init(x);
+
mpfr_init2(z, 32);
mpfr_set_d(z, 1.0, 0);
if (mpfr_get_d(z) != 1.0) {
mpfr_print_raw(z); putchar('\n');
printf("Error: 1.0 != 1.0\n"); exit(1);
}
- mpfr_init2(x, 53); mpfr_init2(y, 53);
+ mpfr_set_prec(x, 53); mpfr_init2(y, 53);
mpfr_set_d(x, d=-1.08007920352320089721e+150, 0);
if (mpfr_get_d(x) != d) {
mpfr_print_raw(x); putchar('\n');
@@ -80,11 +75,11 @@ main(int argc, char **argv)
fprintf(stderr,
"Mismatch on : %1.18g != %1.18g\n", d, mpfr_get_d(x));
mpfr_print_raw(x); putchar('\n');
+ exit(1);
}
}
mpfr_clear(x); mpfr_clear(y); mpfr_clear(z);
- exit (0);
-#endif /* 0 */
+ return 0;
}
diff --git a/mpfr/tests/tset_f.c b/mpfr/tests/tset_f.c
index 4a82f0a40..82b09ad76 100644
--- a/mpfr/tests/tset_f.c
+++ b/mpfr/tests/tset_f.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_set_f.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -32,43 +32,50 @@ MA 02111-1307, USA. */
int
main()
{
- mpfr_t x; mpf_t y; mpf_t z; unsigned long k, pr; double f;
+ mpfr_t x, u; mpf_t y, z; unsigned long k, pr;
mpfr_init2(x, 100);
mpf_init(y);
+ mpf_init(z);
mpf_set_d(y, 0.0);
mpfr_set_f(x, y, GMP_RNDN);
- srandom(time(NULL));
+ srandom((int)time(NULL));
mpf_random2(y, 10, 0);
mpfr_set_f(x, y, rand() & 3);
- /* bug found by Jean-Pierre Merlet on February 3, 2000 */
- mpfr_set_prec(x, 256); mpf_init2(y, 256);
- mpfr_set_machine_rnd_mode(GMP_RNDD);
+ /* bug found by Jean-Pierre Merlet */
+ mpfr_set_prec(x, 256); mpf_set_prec(y, 256);
+ mpfr_init2(u, 256);
+ mpfr_set_str(u,
+ "7.f10872b020c49ba5e353f7ced916872b020c49ba5e353f7ced916872b020c498@2",
+ 16, GMP_RNDN);
mpf_set_str(y, "2033.033", 10);
mpfr_set_f(x, y, GMP_RNDN);
- f = mpfr_get_d(x);
- if (f != 2033.0329999999999017745722085) {
+ if (mpfr_cmp(x, u)) {
fprintf(stderr, "mpfr_set_f failed for y=2033.033\n"); exit(1);
}
mpf_set_str(y, "-2033.033", 10);
mpfr_set_f(x, y, GMP_RNDN);
- f = mpfr_get_d(x);
- if (f != -2033.0330000000001291482476518) {
+ mpfr_neg(u, u, GMP_RNDN);
+ if (mpfr_cmp(x, u)) {
fprintf(stderr, "mpfr_set_f failed for y=-2033.033\n"); exit(1);
}
- mpf_clear(y); mpfr_clear(x);
+ mpfr_clear(u);
+ mpfr_clear(x);
for (k = 1; k <= 100000; k++)
{
- pr = 1 + (rand()&255);
- mpf_init2(z, pr);
+ pr = 1 + (rand()&255);
+ mpf_set_prec(z, pr);
mpf_random2(z, z->_mp_prec, 0);
mpfr_init2(x, pr);
mpfr_set_f(x, z, 0);
+ mpfr_clear(x);
}
- exit (0);
+ mpf_clear(y);
+ mpf_clear(z);
+ return(0);
}
diff --git a/mpfr/tests/tset_q.c b/mpfr/tests/tset_q.c
new file mode 100644
index 000000000..d833378cd
--- /dev/null
+++ b/mpfr/tests/tset_q.c
@@ -0,0 +1,78 @@
+/* Test file for mpfr_set_q.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <stdlib.h>
+#include <unistd.h>
+#include "gmp.h"
+#include "mpfr.h"
+
+void check _PROTO((long int, long int, mp_rnd_t, double));
+
+void check(long int n, long int d, mp_rnd_t rnd, double y)
+{
+ mpq_t q; mpfr_t x; double z;
+
+ mpfr_init2(x, 53); mpq_init(q);
+ mpq_set_si(q, n, d);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd);
+ y = (double) n / d;
+#endif
+ mpfr_set_q(x, q, rnd);
+ z = mpfr_get_d(x);
+ if (y != z) {
+ fprintf(stderr, "Error for q=%ld/%lu and rnd=%s\n", n, d,
+ mpfr_print_rnd_mode(rnd));
+ fprintf(stderr, "libm.a gives %1.20e, mpfr_set_q gives %1.20e\n",
+ y, z);
+ exit(1);
+ }
+ mpfr_clear(x); mpq_clear(q);
+}
+
+int main()
+{
+#ifdef TEST
+ long int i, n;
+ unsigned long int d;
+ double y;
+ unsigned char rnd;
+
+ srand48(getpid());
+ for (i=0;i<1000000;i++) {
+ n = lrand48();
+ d = lrand48();
+ if (lrand48()%2) n = -n;
+ rnd = lrand48() % 4;
+ y = (double) n / d;
+ check(n, d, rnd, y);
+ }
+#endif
+ check(-1647229822, 40619231, GMP_RNDZ, -4.055295438754120596e1);
+ check(-148939696, 1673285490, GMP_RNDZ, -8.9010331404953485501e-2);
+ check(-441322590, 273662545, GMP_RNDZ, -1.6126525096812205362);
+ check(-1631156895, 1677687197, GMP_RNDU, -9.722652100563177191e-1);
+ check(2141332571, 3117601, GMP_RNDZ, 6.8685267004982347316e2);
+ check(75504803, 400207282, GMP_RNDU, 1.8866424074712365155e-1);
+ check(643562308, 23100894, GMP_RNDD, 2.7858762002890447462e1);
+ check(632549085, 1831935802, GMP_RNDN, 3.4528998467600230393e-1);
+ return 0;
+}
diff --git a/mpfr/tests/tset_si.c b/mpfr/tests/tset_si.c
index db4e73ac7..0d20cbe4c 100644
--- a/mpfr/tests/tset_si.c
+++ b/mpfr/tests/tset_si.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_set_si and mpfr_set_ui.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,6 +23,7 @@ MA 02111-1307, USA. */
#include <stdlib.h>
#include "gmp.h"
#include "mpfr.h"
+#include "mpfr-impl.h"
#include "time.h"
int
@@ -41,11 +42,9 @@ main(int argc, char **argv)
z = random() - (1 << 30);
mpfr_set_si(x, z, GMP_RNDZ);
d = (int)mpfr_get_d(x);
- if (d != z)
- {
- fprintf(stderr, "Expected %ld got %ld\n", z, d);
- exit(1);
- }
+ if (d != z) {
+ fprintf(stderr, "Error in mpfr_set_si: expected %ld got %ld\n", z, d); exit(1);
+ }
}
for (k = 1; k <= N; k++)
@@ -53,13 +52,39 @@ main(int argc, char **argv)
zl = random();
mpfr_set_ui(x, zl, GMP_RNDZ);
dl = (unsigned int) mpfr_get_d(x);
- if (dl != zl)
- {
- fprintf(stderr, "Expected %lu got %lu\n", zl, dl);
- exit(1);
- }
+ if (dl != zl) {
+ fprintf(stderr, "Error in mpfr_set_ui: expected %lu got %lu\n", zl, dl); exit(1);
+ }
}
+ mpfr_set_prec(x, 3);
+ mpfr_set_si(x, 77617, GMP_RNDD); /* should be 65536 */
+ if (MPFR_MANT(x)[0] != ((mp_limb_t)1 << (mp_bits_per_limb-1))) {
+ fprintf(stderr, "Error in mpfr_set_si(x:3, 77617, GMP_RNDD)\n");
+ mpfr_print_raw(x); putchar('\n');
+ exit(1);
+ }
+ mpfr_set_ui(x, 77617, GMP_RNDD); /* should be 65536 */
+ if (MPFR_MANT(x)[0] != ((mp_limb_t)1 << (mp_bits_per_limb-1))) {
+ fprintf(stderr, "Error in mpfr_set_ui(x:3, 77617, GMP_RNDD)\n");
+ mpfr_print_raw(x); putchar('\n');
+ exit(1);
+ }
+
+ mpfr_set_prec(x, 1);
+ mpfr_set_si(x, 33096, GMP_RNDU);
+ if (mpfr_get_d(x) != 65536.0) {
+ fprintf(stderr, "Error in mpfr_set_si, expected 65536, got %lu\n",
+ (unsigned long) mpfr_get_d(x));
+ exit(1);
+ }
+ mpfr_set_ui(x, 33096, GMP_RNDU);
+ if (mpfr_get_d(x) != 65536.0) {
+ fprintf(stderr, "Error in mpfr_set_ui, expected 65536, got %lu\n",
+ (unsigned long) mpfr_get_d(x));
+ exit(1);
+ }
+
mpfr_clear(x);
- exit (0);
+ return(0);
}
diff --git a/mpfr/tests/tset_str.c b/mpfr/tests/tset_str.c
index a1cd3000d..cda222d30 100644
--- a/mpfr/tests/tset_str.c
+++ b/mpfr/tests/tset_str.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_set_str.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -21,32 +21,30 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
#include "gmp.h"
#include "mpfr.h"
+#include "mpfr-impl.h"
#include <time.h>
-#if defined (hpux)
-#define srandom srand48
-#define random mrand48
-#endif
-
int
main(int argc, char **argv)
{
- mpfr_t x; unsigned long k, bd, nc; char *str, *str2;
+ mpfr_t x, y; unsigned long k, bd, nc, i; char *str, *str2; mp_exp_t e;
+ int base, logbase, prec, baseprec;
if (argc==2) { /* tset_str <string> */
mpfr_init2(x, 53);
mpfr_set_str_raw(x, argv[1]);
printf("%1.20e\n", mpfr_get_d(x));
mpfr_clear(x);
- exit (0);
+ return 0;
}
srandom(time(NULL));
if (argc > 1) { nc = atoi(argv[1]); } else { nc = 53; }
- if (nc < 24) { nc = 24; }
+ if (nc < 100) { nc = 100; }
bd = random()&8;
@@ -67,21 +65,71 @@ main(int argc, char **argv)
}
*(str2++) = 'e';
- sprintf(str2, "%d", random() - (1 << 30));
+ sprintf(str2, "%d", (int) random() - (1 << 30));
- /* printf("%s\n", str); */
mpfr_init2(x, nc + 10);
- mpfr_set_str_raw(x, str);
- /* mpfr_print_raw(x); printf("\n"); */
+ mpfr_set_str_raw(x, str);
+
+ mpfr_set_prec(x, 54);
+ mpfr_set_str_raw(x, "0.100100100110110101001010010101111000001011100100101010E-529");
+ mpfr_init2(y, 54);
+ mpfr_set_str(y, "4.936a52bc17254@-133", 16, GMP_RNDN);
+ if (mpfr_cmp(x, y)) {
+ fprintf(stderr, "Error in mpfr_set_str\n");
+ mpfr_print_raw(x); putchar('\n');
+ mpfr_print_raw(y); putchar('\n');
+ mpfr_clear(x); mpfr_clear(y);
+ exit(1);
+ }
+
+ free(str);
mpfr_set_prec(x, 53);
mpfr_set_str_raw(x, "+110101100.01010000101101000000100111001000101011101110E00");
mpfr_set_str_raw(x, "1.0");
if (mpfr_get_d(x) != 1.0) {
- fprintf(stderr, "Error in mpfr_set_str_raw for s=1.0\n"); exit(1);
+ fprintf(stderr, "Error in mpfr_set_str_raw for s=1.0\n");
+ mpfr_clear(x); mpfr_clear(y);
+ exit(1);
+ }
+
+ mpfr_set_str(x, "+243495834958.53452345E1", 10, GMP_RNDN);
+ mpfr_set_str(x, "9007199254740993", 10, GMP_RNDN);
+ mpfr_set_str(x, "9007199254740992", 10, GMP_RNDU);
+ mpfr_set_str(x, "9007199254740992", 10, GMP_RNDD);
+ mpfr_set_str(x, "9007199254740992", 10, GMP_RNDZ);
+
+ /* check a random number printed and read is not modified */
+ prec = 53;
+ mpfr_set_prec(x, prec);
+ mpfr_set_prec(y, prec);
+ for (i=0;i<100000;i++) {
+ mpfr_random(x);
+ k = rand() % 4;
+ logbase = (rand() % 5) + 1;
+ base = 1 << logbase;
+ /* Warning: the number of bits needed to print exactly a number of
+ 'prec' bits in base 2^logbase may be greater than ceil(prec/logbase),
+ for example 0.11E-1 in base 2 cannot be written exactly with only
+ one digit in base 4 */
+ if (base==2) baseprec=prec;
+ else baseprec=1+(prec-2+logbase)/logbase;
+ str = mpfr_get_str(NULL, &e, base, baseprec, x, k);
+ mpfr_set_str(y, str, base, k);
+ MPFR_EXP(y) += logbase*(e-strlen(str));
+ if (mpfr_cmp(x, y)) {
+ fprintf(stderr, "mpfr_set_str o mpfr_get_str <> id for rnd_mode=%s\n",
+ mpfr_print_rnd_mode(k));
+ printf("x="); mpfr_print_raw(x); putchar('\n');
+ printf("s=%s, exp=%d, base=%d\n", str, (int) e, base);
+ printf("y="); mpfr_print_raw(y); putchar('\n');
+ mpfr_clear(x); mpfr_clear(y);
+ exit(1);
+ }
+ free(str);
}
- mpfr_clear(x); free(str);
- exit (0);
+ mpfr_clear(x); mpfr_clear(y);
+ return 0;
}
diff --git a/mpfr/tests/tset_z.c b/mpfr/tests/tset_z.c
index 61fe654f3..ae069f38f 100644
--- a/mpfr/tests/tset_z.c
+++ b/mpfr/tests/tset_z.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_set_z.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -24,7 +24,8 @@ MA 02111-1307, USA. */
#include "gmp.h"
#include "mpfr.h"
-/* tset_z z rnd prec */
+void check _PROTO((long, unsigned char));
+void check_large _PROTO((void));
void check(long i, unsigned char rnd) {
mpfr_t f; mpz_t z;
@@ -39,7 +40,7 @@ void check(long i, unsigned char rnd) {
mpfr_clear(f); mpz_clear(z);
}
-void check_large()
+void check_large ()
{
mpz_t z; mpfr_t x,y;
@@ -53,14 +54,17 @@ void check_large()
mpz_clear(z); mpfr_clear(x); mpfr_clear(y);
}
-int main(argc,argv) int argc; char *argv[];
+/* tset_z z rnd prec */
+
+int main(int argc, char *argv[])
{
long j;
check_large();
srand(getpid());
+ check(0, 0);
for (j=0; j<1000000; j++)
check(lrand48(), rand()%4);
- exit (0);
+ return 0;
}
diff --git a/mpfr/tests/tsin_cos.c b/mpfr/tests/tsin_cos.c
new file mode 100644
index 000000000..817d10e32
--- /dev/null
+++ b/mpfr/tests/tsin_cos.c
@@ -0,0 +1,101 @@
+/* Test file for mpfr_sin_cos.
+
+Copyright (C) 2000 Free Software Foundation, Inc.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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"
+
+void large_test (int, int);
+void check53 (double, double, double, mp_rnd_t);
+
+void large_test (int prec, int N)
+{
+ int i;
+ mpfr_t x, s, c;
+
+ mpfr_init2 (x, prec);
+ mpfr_init2 (s, prec);
+ mpfr_init2 (c, prec);
+ mpfr_set_d (x, 3.0, GMP_RNDN);
+ mpfr_sqrt (x, x, GMP_RNDN);
+ for (i=0; i<N; i++) mpfr_sin_cos (s, c, x, GMP_RNDN);
+ mpfr_out_str (stdout, 10, 0, s, GMP_RNDN); putchar('\n');
+ mpfr_clear (x);
+ mpfr_clear (s);
+ mpfr_clear (c);
+}
+
+void check53 (double x, double sin_x, double cos_x, mp_rnd_t rnd_mode)
+{
+ mpfr_t xx, s, c;
+
+ mpfr_init2 (xx, 53);
+ mpfr_init2 (s, 53);
+ mpfr_init2 (c, 53);
+ mpfr_set_d (xx, x, rnd_mode); /* should be exact */
+ mpfr_sin_cos (s, c, xx, rnd_mode);
+ if (mpfr_get_d (s) != sin_x && (!isnan(sin_x) || !isnan(mpfr_get_d(s)))) {
+ fprintf (stderr, "mpfr_sin_cos failed for x=%1.20e, rnd=%s\n", x,
+ mpfr_print_rnd_mode (rnd_mode));
+ fprintf (stderr, "mpfr_sin_cos gives sin(x)=%1.20e, expected %1.20e\n",
+ mpfr_get_d (s), sin_x);
+ exit(1);
+ }
+ if (mpfr_get_d (c) != cos_x && (!isnan(cos_x) || !isnan(mpfr_get_d(c)))) {
+ fprintf (stderr, "mpfr_sin_cos failed for x=%1.20e, rnd=%s\n", x,
+ mpfr_print_rnd_mode (rnd_mode));
+ fprintf (stderr, "mpfr_sin_cos gives cos(x)=%1.20e, expected %1.20e\n",
+ mpfr_get_d (c), cos_x);
+ exit(1);
+ }
+ mpfr_clear (xx);
+ mpfr_clear (s);
+ mpfr_clear (c);
+}
+
+/* tsin_cos prec [N] performs N tests with prec bits */
+int main(int argc, char *argv[])
+{
+ if (argc > 1) {
+ large_test (atoi (argv[1]), (argc > 2) ? atoi (argv[2]) : 1);
+ }
+
+ check53(0.0/0.0, 0.0/0.0, 0.0/0.0, GMP_RNDN);
+ check53(1.0/0.0, 0.0/0.0, 0.0/0.0, GMP_RNDN);
+ check53(-1.0/0.0, 0.0/0.0, 0.0/0.0, GMP_RNDN);
+ /* worst case from PhD thesis of Vincent Lefe`vre: x=8980155785351021/2^54 */
+ check53 (4.984987858808754279e-1, 4.781075595393330379e-1,
+ 8.783012931285841817e-1, GMP_RNDN);
+ check53 (4.984987858808754279e-1, 4.781075595393329824e-1,
+ 8.783012931285840707e-1, GMP_RNDD);
+ check53 (4.984987858808754279e-1, 4.781075595393329824e-1,
+ 8.783012931285840707e-1, GMP_RNDZ);
+ check53 (4.984987858808754279e-1, 4.781075595393330379e-1,
+ 8.783012931285841817e-1, GMP_RNDU);
+ check53 (1.00031274099908640274, 8.416399183372403892e-1,
+ 0.540039116973283217504, GMP_RNDN);
+ check53 (1.00229256850978698523, 8.427074524447979442e-1,
+ 0.538371757797526551137, GMP_RNDZ);
+ check53 (1.00288304857059840103, 8.430252033025980029e-1,
+ 0.537874062022526966409, GMP_RNDZ);
+ check53 (1.00591265847407274059, 8.446508805292128885e-1,
+ 0.53531755997839769456, GMP_RNDN);
+ return 0;
+}
diff --git a/mpfr/tests/tsqrt.c b/mpfr/tests/tsqrt.c
index ccffad0c5..940ffbd44 100644
--- a/mpfr/tests/tsqrt.c
+++ b/mpfr/tests/tsqrt.c
@@ -1,6 +1,6 @@
/* Test file for mpfr_sqrt.
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
+Copyright (C) 1999 Free Software Foundation.
This file is part of the MPFR Library.
@@ -23,58 +23,90 @@ MA 02111-1307, USA. */
#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"
-#include "gmp-impl.h"
#include "mpfr.h"
-#include "mpfr-impl.h"
-#ifdef IRIX64
-#include <sys/fpu.h>
-#endif
-
-extern int isnan(), getpid();
+#include "mpfr-test.h"
#define check(a,r) check3(a,r,-1.0)
int maxulp=0;
-void check3(a, rnd_mode, Q) double a; unsigned char rnd_mode; double Q;
+void check3 _PROTO((double, mp_rnd_t, double));
+void check4 _PROTO((double, mp_rnd_t, char *));
+void check24 _PROTO((float, mp_rnd_t, float));
+void check_float _PROTO((void));
+void special _PROTO((void));
+
+void check3(double a, mp_rnd_t rnd_mode, double Q)
{
mpfr_t q; double Q2; int ck,u;
ck = (Q!=-1.0); /* if ck=1, then Q is certified correct */
mpfr_init2(q, 53);
mpfr_set_d(q, a, rnd_mode);
+#ifdef TEST
mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
mpfr_sqrt(q, q, rnd_mode);
if (ck==0) Q = sqrt(a);
+ else {
+ if (Q != sqrt(a) && (!isnan(Q) || !isnan(sqrt(a)))) {
+ fprintf(stderr, "you've found a bug in your machine's sqrt for x=%1.20e\n", a);
+ mpfr_clear(q);
+ exit(1);
+
+ }
+ }
Q2 = mpfr_get_d(q);
if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
u = ulp(Q2,Q);
if (ck) {
- printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%d\n",a,rnd_mode);
+ printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
+ a, mpfr_print_rnd_mode(rnd_mode));
printf("expected sqrt is %1.20e, got %1.20e (%d ulp)\n",Q,Q2,u);
+ mpfr_clear(q);
exit(1);
}
else if (u>maxulp || u<-maxulp) {
maxulp = (u>maxulp) ? u : -u;
- printf("libm.a differs from mpfr_sqrt for a=%1.20e, rnd_mode=%d\n",a,rnd_mode);
+ printf("libm.a differs from mpfr_sqrt for a=%1.20e, rnd_mode=%s\n",
+ a, mpfr_print_rnd_mode(rnd_mode));
printf("libm.a gives %1.20e, mpfr_sqrt gives %1.20e (%d ulp)\n",Q,Q2,u);
}
}
mpfr_clear(q);
}
-void check24(a, rnd_mode) float a; unsigned char rnd_mode;
+void check4 (double a, mp_rnd_t rnd_mode, char *Q)
{
- mpfr_t q; float Q,Q2;
+ mpfr_t q, res;
+
+ mpfr_init2(q, 53); mpfr_init2(res, 53);
+ mpfr_set_d(q, a, rnd_mode);
+ mpfr_sqrt(q, q, rnd_mode);
+ mpfr_set_str(res, Q, 16, GMP_RNDN);
+ if (mpfr_cmp(q, res)) {
+ printf("mpfr_sqrt failed for a=%1.20e, rnd_mode=%s\n",
+ a, mpfr_print_rnd_mode(rnd_mode));
+ printf("expected "); mpfr_print_raw(res); putchar('\n');
+ printf("got "); mpfr_print_raw(q); putchar('\n');
+ mpfr_clear(q); mpfr_clear(res);
+ exit(1);
+ }
+ mpfr_clear(res);
+ mpfr_clear(q);
+}
+
+void check24 (float a, mp_rnd_t rnd_mode, float Q)
+{
+ mpfr_t q; float Q2;
mpfr_init2(q, 24);
mpfr_set_d(q, a, rnd_mode);
- mpfr_set_machine_rnd_mode(rnd_mode);
mpfr_sqrt(q, q, rnd_mode);
- Q = sqrt(a);
Q2 = mpfr_get_d(q);
- if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
- printf("mpfr_sqrt failed for a=%1.10e, rnd_mode=%d\n",a,rnd_mode);
+ if (Q!=Q2) {
+ printf("mpfr_sqrt failed for a=%1.10e, prec=24, rnd_mode=%s\n",
+ a, mpfr_print_rnd_mode(rnd_mode));
printf("expected sqrt is %1.10e, got %1.10e\n",Q,Q2);
exit(1);
}
@@ -83,28 +115,89 @@ void check24(a, rnd_mode) float a; unsigned char rnd_mode;
/* the following examples come from the paper "Number-theoretic Test
Generation for Directed Rounding" from Michael Parks, Table 3 */
-void check_float()
+void check_float ()
{
- int i; float b = 8388608.0; /* 2^23 */
-
- for (i=0;i<4;i++) {
- check24(b*8388610.0, i);
- check24(b*2.0*16777214.0, i);
- check24(b*8388612.0, i);
- check24(b*2.0*16777212.0, i);
- check24(b*11946704.0, i);
- check24(b*14321479.0, i);
- check24(b*2.0*13689673.0, i);
- check24(b*8388614.0, i);
- check24(b*2.0*16777210.0, i);
- check24(b*10873622.0, i);
+ float b = 8388608.0; /* 2^23 */
+
+ check24(b*8388610.0, GMP_RNDN, 8.388609e6);
+ check24(b*2.0*16777214.0, GMP_RNDN, 1.6777215e7);
+ check24(b*8388612.0, GMP_RNDN, 8.388610e6);
+ check24(b*2.0*16777212.0, GMP_RNDN, 1.6777214e7);
+ check24(b*11946704.0, GMP_RNDN, 1.0010805e7);
+ check24(b*14321479.0, GMP_RNDN, 1.0960715e7);
+ check24(b*2.0*13689673.0, GMP_RNDN, 1.5155019e7);
+ check24(b*8388614.0, GMP_RNDN, 8.388611e6);
+ check24(b*2.0*16777210.0, GMP_RNDN, 1.6777213e7);
+ check24(b*10873622.0, GMP_RNDN, 9.550631e6);
+
+ check24(b*8388610.0, GMP_RNDZ, 8.388608e6);
+ check24(b*2.0*16777214.0, GMP_RNDZ, 1.6777214e7);
+ check24(b*8388612.0, GMP_RNDZ, 8.388609e6);
+ check24(b*2.0*16777212.0, GMP_RNDZ, 1.6777213e7);
+ check24(b*11946704.0, GMP_RNDZ, 1.0010805e7);
+ check24(b*14321479.0, GMP_RNDZ, 1.0960715e7);
+ check24(b*2.0*13689673.0, GMP_RNDZ, 1.5155019e7);
+ check24(b*8388614.0, GMP_RNDZ, 8.38861e6);
+ check24(b*2.0*16777210.0, GMP_RNDZ, 1.6777212e7);
+ check24(b*10873622.0, GMP_RNDZ, 9.550631e6);
+
+ check24(b*8388610.0, GMP_RNDU, 8.388609e6);
+ check24(b*2.0*16777214.0, GMP_RNDU, 1.6777215e7);
+ check24(b*8388612.0, GMP_RNDU, 8.388610e6);
+ check24(b*2.0*16777212.0, GMP_RNDU, 1.6777214e7);
+ check24(b*11946704.0, GMP_RNDU, 1.0010806e7);
+ check24(b*14321479.0, GMP_RNDU, 1.0960716e7);
+ check24(b*2.0*13689673.0, GMP_RNDU, 1.515502e7);
+ check24(b*8388614.0, GMP_RNDU, 8.388611e6);
+ check24(b*2.0*16777210.0, GMP_RNDU, 1.6777213e7);
+ check24(b*10873622.0, GMP_RNDU, 9.550632e6);
+
+ check24(b*8388610.0, GMP_RNDD, 8.388608e6);
+ check24(b*2.0*16777214.0, GMP_RNDD, 1.6777214e7);
+ check24(b*8388612.0, GMP_RNDD, 8.388609e6);
+ check24(b*2.0*16777212.0, GMP_RNDD, 1.6777213e7);
+ check24(b*11946704.0, GMP_RNDD, 1.0010805e7);
+ check24(b*14321479.0, GMP_RNDD, 1.0960715e7);
+ check24(b*2.0*13689673.0, GMP_RNDD, 1.5155019e7);
+ check24(b*8388614.0, GMP_RNDD, 8.38861e6);
+ check24(b*2.0*16777210.0, GMP_RNDD, 1.6777212e7);
+ check24(b*10873622.0, GMP_RNDD, 9.550631e6);
+}
+
+void special ()
+{
+ mpfr_t x, z;
+
+ mpfr_init2 (x, 1);
+ mpfr_init2 (z, 1);
+
+ /* checks the sign is correctly set */
+ mpfr_set_d (x, 1.0, GMP_RNDN);
+ mpfr_set_d (z, -1.0, GMP_RNDN);
+ mpfr_sqrt (z, x, GMP_RNDN);
+ if (mpfr_cmp_ui (z, 0) < 0) {
+ fprintf (stderr, "Error: square root of %e gives %e\n",
+ mpfr_get_d (x), mpfr_get_d (z));
+ exit (1);
}
+
+ mpfr_set_prec (x, 192);
+ mpfr_set_prec (z, 160);
+ mpfr_set_str_raw (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1");
+ mpfr_set_prec(x, 160);
+ mpfr_sqrt(x, z, GMP_RNDN);
+ mpfr_sqrt(z, x, GMP_RNDN);
+
+ mpfr_clear (x);
+ mpfr_clear (z);
}
int main()
{
- int i; double a;
-#ifdef IRIX64
+ double a;
+#ifdef TEST
+ int i;
+#ifdef __mips
/* to get denormalized numbers on IRIX64 */
union fpc_csr exp;
exp.fc_word = get_fpc_csr();
@@ -112,38 +205,80 @@ int main()
set_fpc_csr(exp.fc_word);
#endif
- check_float();
srand(getpid());
-/* the following examples come from the paper "Number-theoretic Test
- Generation for Directed Rounding" from Michael Parks, Table 4 */
- a = 4503599627370496.0; /* 2^52 */
- for (i=0;i<4;i++) {
- check(a*2.0*8732221479794286.0, i);
- check(a*8550954388695124.0, i);
- check(a*7842344481681754.0, i);
- check(a*5935035262218600.0, i);
- check(a*5039650445085418.0, i);
- check(a*5039721545366078.0, i);
- check(a*8005963117781324.0, i);
- check(a*6703494707970582.0, i);
- check(a*8010323124937260.0, i);
- check(a*2.0*8010776873384260.0, i);
- }
- check(6.37983013646045901440e+32, GMP_RNDN);
- check(1.0, GMP_RNDN);
- check(1.0, GMP_RNDZ);
- check(3.725290298461914062500000e-9, GMP_RNDN);
- check(3.725290298461914062500000e-9, GMP_RNDZ);
- a=1190456976439861.0;
- check3(a, GMP_RNDZ, dbl(4630914205854029.0,-27));
- check3(1024.0*a, GMP_RNDZ, dbl(4630914205854029.0,-22));
- /* the following examples are bugs in Cygnus compiler/system, found by
- Fabrice Rouillier while porting mpfr to Windows */
- check3(9.89438396044940256501e-134, GMP_RNDU, 3.14553397063986684729e-67);
- check3(7.86528588050363751914e+31, GMP_RNDZ, 8.86864469944739400000e+15);
for (i=0;i<100000;i++) {
a = drand();
+ if (a < 0.0) a = -a; /* ensures a is positive */
check(a, rand() % 4);
}
- exit (0);
+#endif
+ special ();
+ check_float();
+ check3(0.0/0.0, GMP_RNDN, 0.0/0.0);
+ check3(-1.0, GMP_RNDN, 0.0/0.0);
+ check3(1.0/0.0, GMP_RNDN, 1.0/0.0);
+ check3(-1.0/0.0, GMP_RNDN, 0.0/0.0);
+ check3(-0.0, GMP_RNDN, 0.0);
+ check4(6.37983013646045901440e+32, GMP_RNDN, "5.9bc5036d09e0c@13");
+ check4(1.0, GMP_RNDN, "1");
+ check4(1.0, GMP_RNDZ, "1");
+ check4(3.725290298461914062500000e-9, GMP_RNDN, "4@-4");
+ check4(3.725290298461914062500000e-9, GMP_RNDZ, "4@-4");
+ a=1190456976439861.0;
+ check4(a, GMP_RNDZ, "2.0e7957873529a@6");
+ check4(1024.0*a, GMP_RNDZ, "4.1cf2af0e6a534@7");
+ /* the following examples are bugs in Cygnus compiler/system, found by
+ Fabrice Rouillier while porting mpfr to Windows */
+ check4(9.89438396044940256501e-134, GMP_RNDU, "8.7af7bf0ebbee@-56");
+ check4(7.86528588050363751914e+31, GMP_RNDZ, "1.f81fc40f32062@13");
+ check4(0.99999999999999988897, GMP_RNDN, "f.ffffffffffff8@-1");
+ check4(1.00000000000000022204, GMP_RNDN, "1");
+ /* the following examples come from the paper "Number-theoretic Test
+ Generation for Directed Rounding" from Michael Parks, Table 4 */
+ a = 4503599627370496.0; /* 2^52 */
+
+ check4(a*2.0*8732221479794286.0, GMP_RNDN, "1.f81fc40f32063@13");
+ check4(a*8550954388695124.0, GMP_RNDN, "1.60c012a92fc65@13");
+ check4(a*7842344481681754.0, GMP_RNDN, "1.51d17526c7161@13");
+ check4(a*5935035262218600.0, GMP_RNDN, "1.25e19302f7e51@13");
+ check4(a*5039650445085418.0, GMP_RNDN, "1.0ecea7dd2ec3d@13");
+ check4(a*5039721545366078.0, GMP_RNDN, "1.0ecf250e8e921@13");
+ check4(a*8005963117781324.0, GMP_RNDN, "1.5552f3eedcf33@13");
+ check4(a*6703494707970582.0, GMP_RNDN, "1.3853ee10c9c99@13");
+ check4(a*8010323124937260.0, GMP_RNDN, "1.556abe212b56f@13");
+ check4(a*2.0*8010776873384260.0, GMP_RNDN, "1.e2d9a51977e6e@13");
+
+ check4(a*2.0*8732221479794286.0, GMP_RNDZ, "1.f81fc40f32062@13");
+ check4(a*8550954388695124.0, GMP_RNDZ, "1.60c012a92fc64@13");
+ check4(a*7842344481681754.0, GMP_RNDZ, "1.51d17526c716@13");
+ check4(a*5935035262218600.0, GMP_RNDZ, "1.25e19302f7e5@13");
+ check4(a*5039650445085418.0, GMP_RNDZ, "1.0ecea7dd2ec3c@13");
+ check4(a*5039721545366078.0, GMP_RNDZ, "1.0ecf250e8e92@13");
+ check4(a*8005963117781324.0, GMP_RNDZ, "1.5552f3eedcf32@13");
+ check4(a*6703494707970582.0, GMP_RNDZ, "1.3853ee10c9c98@13");
+ check4(a*8010323124937260.0, GMP_RNDZ, "1.556abe212b56e@13");
+ check4(a*2.0*8010776873384260.0, GMP_RNDZ, "1.e2d9a51977e6d@13");
+
+ check4(a*2.0*8732221479794286.0, GMP_RNDU, "1.f81fc40f32063@13");
+ check4(a*8550954388695124.0, GMP_RNDU, "1.60c012a92fc65@13");
+ check4(a*7842344481681754.0, GMP_RNDU, "1.51d17526c7161@13");
+ check4(a*5935035262218600.0, GMP_RNDU, "1.25e19302f7e51@13");
+ check4(a*5039650445085418.0, GMP_RNDU, "1.0ecea7dd2ec3d@13");
+ check4(a*5039721545366078.0, GMP_RNDU, "1.0ecf250e8e921@13");
+ check4(a*8005963117781324.0, GMP_RNDU, "1.5552f3eedcf33@13");
+ check4(a*6703494707970582.0, GMP_RNDU, "1.3853ee10c9c99@13");
+ check4(a*8010323124937260.0, GMP_RNDU, "1.556abe212b56f@13");
+ check4(a*2.0*8010776873384260.0, GMP_RNDU, "1.e2d9a51977e6e@13");
+
+ check4(a*2.0*8732221479794286.0, GMP_RNDD, "1.f81fc40f32062@13");
+ check4(a*8550954388695124.0, GMP_RNDD, "1.60c012a92fc64@13");
+ check4(a*7842344481681754.0, GMP_RNDD, "1.51d17526c716@13");
+ check4(a*5935035262218600.0, GMP_RNDD, "1.25e19302f7e5@13");
+ check4(a*5039650445085418.0, GMP_RNDD, "1.0ecea7dd2ec3c@13");
+ check4(a*5039721545366078.0, GMP_RNDD, "1.0ecf250e8e92@13");
+ check4(a*8005963117781324.0, GMP_RNDD, "1.5552f3eedcf32@13");
+ check4(a*6703494707970582.0, GMP_RNDD, "1.3853ee10c9c98@13");
+ check4(a*8010323124937260.0, GMP_RNDD, "1.556abe212b56e@13");
+ check4(a*2.0*8010776873384260.0, GMP_RNDD, "1.e2d9a51977e6d@13");
+ return 0;
}
diff --git a/mpfr/tests/tsqrt_ui.c b/mpfr/tests/tsqrt_ui.c
new file mode 100644
index 000000000..f803e5317
--- /dev/null
+++ b/mpfr/tests/tsqrt_ui.c
@@ -0,0 +1,79 @@
+/* Test file for mpfr_sqrt_ui.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "mpfr-test.h"
+
+void check _PROTO((unsigned long, mp_rnd_t, double));
+
+int maxulp=0;
+
+void check (unsigned long a, mp_rnd_t rnd_mode, double Q)
+{
+ mpfr_t q; double Q2; int u, ck;
+
+ mpfr_init2(q, 53);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
+ mpfr_sqrt_ui(q, a, rnd_mode);
+ ck = (Q >= 0.0);
+ if (!ck) Q = sqrt(1.0 * a);
+ Q2 = mpfr_get_d(q);
+ if (Q!=Q2 && (!isnan(Q) || !isnan(Q2))) {
+ u = ulp(Q2,Q);
+ if (ck) printf("mpfr_sqrt_ui failed");
+ else printf("mpfr_sqrt_ui differs from sqrt");
+ printf(" for a=%lu, rnd_mode=%s\n",
+ a, mpfr_print_rnd_mode(rnd_mode));
+ printf("sqrt gives %1.20e, mpfr_sqrt_ui gives %1.20e (%d ulp)\n",Q,Q2,u);
+ exit(1);
+ }
+ mpfr_clear(q);
+}
+
+int main()
+{
+#ifdef TEST
+ int i; unsigned long a;
+#ifdef __mips
+ /* to get denormalized numbers on IRIX64 */
+ union fpc_csr exp;
+ exp.fc_word = get_fpc_csr();
+ exp.fc_struct.flush = 0;
+ set_fpc_csr(exp.fc_word);
+#endif
+
+ srand(getpid());
+ for (i=0;i<1000000;i++) {
+ a = lrand48();
+ /* machine arithmetic must agree if a <= 2.0^53 */
+ if (1.0*a < 9007199254872064.0) check(a, rand() % 4, -1.0);
+ }
+#endif
+ check(0, GMP_RNDN, 0.0);
+ check(2116118, GMP_RNDU, 1.45468828276026215e3);
+ return 0;
+}
diff --git a/mpfr/tests/tswap.c b/mpfr/tests/tswap.c
new file mode 100644
index 000000000..1cdbcd3ac
--- /dev/null
+++ b/mpfr/tests/tswap.c
@@ -0,0 +1,42 @@
+/* Test file for mpfr_swap.
+
+Copyright (C) 2000 Free Software Foundation, Inc.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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"
+
+int main()
+{
+ mpfr_t u, v;
+
+ mpfr_init2 (u, 24);
+ mpfr_init2 (v, 53);
+ mpfr_set_ui (u, 16777215, GMP_RNDN); /* 2^24 - 1 */
+ mpfr_set_d (v, 9007199254740991.0, GMP_RNDN); /* 2^53 - 1 */
+ mpfr_swap (u, v);
+ mpfr_swap (u, v);
+ if (mpfr_cmp_ui (u, 16777215) || (mpfr_get_d (v) != 9007199254740991.0)) {
+ fprintf (stderr, "Error in mpfr_swap\n");
+ exit (1);
+ }
+ mpfr_clear (u);
+ mpfr_clear (v);
+ return 0;
+}
diff --git a/mpfr/tests/ttrunc.c b/mpfr/tests/ttrunc.c
new file mode 100644
index 000000000..41987170d
--- /dev/null
+++ b/mpfr/tests/ttrunc.c
@@ -0,0 +1,144 @@
+/* Test file for mpfr_trunc, mpfr_ceil, mpfr_floor.
+
+Copyright (C) 1999-2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 "mpfr-impl.h"
+
+#define SIZEX 100
+
+int main()
+{
+ int j, k; mpfr_t x, y, z, t, y2, z2, t2;
+
+ mpfr_init2(x, SIZEX);
+ mpfr_init2(y, SIZEX);
+ mpfr_init2(z, SIZEX);
+ mpfr_init2(t, SIZEX);
+ mpfr_init2(y2, SIZEX);
+ mpfr_init2(z2, SIZEX);
+ mpfr_init2(t2, SIZEX);
+
+ mpfr_set_d(x, 0.5, GMP_RNDN);
+ mpfr_ceil(y, x);
+ if (mpfr_get_d(y) != 1.0) {
+ fprintf(stderr, "Error in mpfr_ceil for x=0.5: expected 1.0, got %f\n",
+ mpfr_get_d(y)); exit(1);
+ }
+
+ mpfr_set_d(x, 0.0, GMP_RNDN);
+ mpfr_ceil(y, x);
+ if (mpfr_get_d(y) != 0.0) {
+ fprintf(stderr, "Error in mpfr_ceil for x=0.0: expected 0.0, got %f\n",
+ mpfr_get_d(y)); exit(1);
+ }
+
+ mpfr_set_d(x, 1.0, GMP_RNDN);
+ mpfr_ceil(y, x);
+ if (mpfr_get_d(y) != 1.0) {
+ fprintf(stderr, "Error in mpfr_ceil for x=1.0: expected 1.0, got %f\n",
+ mpfr_get_d(y)); exit(1);
+ }
+
+ for (j=0;j<1000;j++) {
+
+ mpfr_random(x);
+ MPFR_EXP(x) = 1;
+
+ for (k = 1; k <= SIZEX; k++)
+ {
+ mpfr_set_prec(y, k);
+ mpfr_set_prec(y2, k);
+ mpfr_set_prec(z, k);
+ mpfr_set_prec(z2, k);
+ mpfr_set_prec(t, k);
+ mpfr_set_prec(t2, k);
+
+ mpfr_floor(y, x);
+ mpfr_set(y2, x, GMP_RNDD);
+
+ mpfr_trunc(z, x);
+ mpfr_set(z2, x, GMP_RNDZ);
+
+ mpfr_ceil(t, x);
+ mpfr_set(t2, x, GMP_RNDU);
+
+ if (!mpfr_eq(y, y2, k))
+ {
+ printf("Error in floor, x = "); mpfr_print_raw(x); printf("\n");
+ printf("floor(x) = "); mpfr_print_raw(y); printf("\n");
+ printf("round(x, RNDD) = "); mpfr_print_raw(y2); printf("\n");
+ mpfr_clear(x);
+ mpfr_clear(y);
+ mpfr_clear(y2);
+ mpfr_clear(z);
+ mpfr_clear(z2);
+ mpfr_clear(t);
+ mpfr_clear(t2);
+ exit(-1);
+ }
+
+ if (!mpfr_eq(z, z2, k))
+ {
+ printf("Error in trunc, x = "); mpfr_print_raw(x); printf("\n");
+ printf("trunc(x) = "); mpfr_print_raw(z); printf("\n");
+ printf("round(x, RNDZ) = "); mpfr_print_raw(z2); printf("\n");
+ mpfr_clear(x);
+ mpfr_clear(y);
+ mpfr_clear(y2);
+ mpfr_clear(z);
+ mpfr_clear(z2);
+ mpfr_clear(t);
+ mpfr_clear(t2);
+ exit(-1);
+ }
+
+ if (!mpfr_eq(y, y2, k))
+ {
+ printf("Error in ceil, x = "); mpfr_print_raw(x); printf("\n");
+ printf("ceil(x) = "); mpfr_print_raw(t); printf("\n");
+ printf("round(x, RNDU) = "); mpfr_print_raw(t2); printf("\n");
+ mpfr_clear(x);
+ mpfr_clear(y);
+ mpfr_clear(y2);
+ mpfr_clear(z);
+ mpfr_clear(z2);
+ mpfr_clear(t);
+ mpfr_clear(t2);
+ exit(-1);
+ }
+ MPFR_EXP(x)++;
+ }
+ }
+
+ mpfr_clear(x);
+ mpfr_clear(y);
+ mpfr_clear(y2);
+ mpfr_clear(z);
+ mpfr_clear(z2);
+ mpfr_clear(t);
+ mpfr_clear(t2);
+
+ return 0;
+}
+
diff --git a/mpfr/tests/tui_div.c b/mpfr/tests/tui_div.c
new file mode 100644
index 000000000..b3e812378
--- /dev/null
+++ b/mpfr/tests/tui_div.c
@@ -0,0 +1,98 @@
+/* Test file for mpfr_ui_div.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#ifdef TEST
+#include "mpfr-test.h"
+#endif
+
+void check _PROTO((unsigned long, double, mp_rnd_t, double));
+
+/* checks that y/x gives the same results in double
+ and with mpfr with 53 bits of precision */
+void check (unsigned long y, double x, mp_rnd_t rnd_mode, double z1)
+{
+ double z2; mpfr_t xx, zz;
+
+ mpfr_init2(xx, 53);
+ mpfr_init2(zz, 53);
+ mpfr_set_d(xx, x, rnd_mode);
+ mpfr_ui_div(zz, y, xx, rnd_mode);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
+ if (z1==0.0) z1 = y/x;
+ z2 = mpfr_get_d(zz);
+ if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
+ printf("expected quotient is %1.20e, got %1.20e\n",z1,z2);
+ printf("mpfr_ui_div failed for y=%lu x=%1.20e with rnd_mode=%s\n",
+ y, x, mpfr_print_rnd_mode(rnd_mode));
+ exit(1);
+ }
+ mpfr_clear(xx); mpfr_clear(zz);
+}
+
+int main(argc,argv) int argc; char *argv[];
+{
+#ifdef TEST
+ double x; unsigned long y, N; int i,rnd_mode,rnd;
+#ifdef __mips
+ /* to get denormalized numbers on IRIX64 */
+ union fpc_csr exp;
+ exp.fc_word = get_fpc_csr();
+ exp.fc_struct.flush = 0;
+ set_fpc_csr(exp.fc_word);
+#endif
+
+ srand48(getpid());
+ N = (argc<2) ? 1000000 : atoi(argv[1]);
+ rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
+ for (i=0;i<1000000;i++) {
+ x = drand();
+ y = lrand48();
+ if (ABS(x)>4e-286) {
+ /* avoid denormalized numbers and overflows */
+ rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
+ check(y, x, rnd, 0.0);
+ }
+ }
+#endif
+ check(1, 1.0/0.0, GMP_RNDN, 0.0);
+ check(1, -1.0/0.0, GMP_RNDN, -0.0);
+ check(1, 0.0/0.0, GMP_RNDN, 0.0/0.0);
+ check(0, 0.0, GMP_RNDN, 0.0/0.0);
+ check(948002822, 1.22191250737771397120e+20, GMP_RNDN,
+ 7.758352715731357946e-12);
+ check(1976245324, 1.25296395864546893357e+232, GMP_RNDZ,
+ 1.5772563211925444801e-223);
+ check(740454110, 2.11496253355831863313e+183, GMP_RNDZ,
+ 3.5010270784996976041e-175);
+ check(1690540942, 1.28278599852446657468e-276, GMP_RNDU,
+ 1.3178666932321966062e285);
+ check(1476599377, -2.14191393656148625995e+305, GMP_RNDD,
+ -6.8938315017943889615e-297);
+ return 0;
+}
+
diff --git a/mpfr/tests/tui_sub.c b/mpfr/tests/tui_sub.c
new file mode 100644
index 000000000..9456e7d10
--- /dev/null
+++ b/mpfr/tests/tui_sub.c
@@ -0,0 +1,99 @@
+/* Test file for mpfr_ui_sub.
+
+Copyright (C) 2000 Free Software Foundation.
+
+This file is part of the MPFR Library.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR 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 Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR 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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "gmp.h"
+#include "mpfr.h"
+#ifdef TEST
+#include "mpfr-test.h"
+#endif
+
+void check _PROTO((unsigned long, double, mp_rnd_t, double));
+
+/* checks that y/x gives the same results in double
+ and with mpfr with 53 bits of precision */
+void check (unsigned long y, double x, mp_rnd_t rnd_mode, double z1)
+{
+ double z2; mpfr_t xx, zz;
+
+ mpfr_init2(xx, 53);
+ mpfr_init2(zz, 53);
+ mpfr_set_d(xx, x, rnd_mode);
+ mpfr_ui_sub(zz, y, xx, rnd_mode);
+#ifdef TEST
+ mpfr_set_machine_rnd_mode(rnd_mode);
+#endif
+ if (z1==0.0) z1 = y-x;
+ z2 = mpfr_get_d(zz);
+ if (z1!=z2 && !(isnan(z1) && isnan(z2))) {
+ printf("expected difference is %1.20e, got %1.20e\n",z1,z2);
+ printf("mpfr_ui_sub failed for y=%lu x=%1.20e with rnd_mode=%s\n",
+ y, x, mpfr_print_rnd_mode(rnd_mode));
+ exit(1);
+ }
+ mpfr_clear(xx); mpfr_clear(zz);
+}
+
+int main(argc,argv) int argc; char *argv[];
+{
+#ifdef TEST
+ double x; unsigned long y, N; int i,rnd_mode,rnd;
+#ifdef __mips
+ /* to get denormalized numbers on IRIX64 */
+ union fpc_csr exp;
+ exp.fc_word = get_fpc_csr();
+ exp.fc_struct.flush = 0;
+ set_fpc_csr(exp.fc_word);
+#endif
+
+ srand48(getpid());
+ N = (argc<2) ? 1000000 : atoi(argv[1]);
+ rnd_mode = (argc<3) ? -1 : atoi(argv[2]);
+ for (i=0;i<1000000;i++) {
+ x = drand();
+ y = lrand48();
+ if (ABS(x)>2.2e-307) {
+ /* avoid denormalized numbers and overflows */
+ rnd = (rnd_mode==-1) ? lrand48()%4 : rnd_mode;
+ check(y, x, rnd, 0.0);
+ }
+ }
+#endif
+ check(1, 1.0/0.0, GMP_RNDN, -1.0/0.0);
+ check(1, -1.0/0.0, GMP_RNDN, 1.0/0.0);
+ check(1, 0.0/0.0, GMP_RNDN, 0.0/0.0);
+ check(1196426492, 1.4218093058435347e-3, GMP_RNDN, 1.1964264919985781e9);
+ check(1092583421, -1.0880649218158844e9, GMP_RNDN, 2.1806483428158845901e9);
+ check(948002822, 1.22191250737771397120e+20, GMP_RNDN,
+ -1.2219125073682338611e20);
+ check(832100416, 4.68311314939691330000e-215, GMP_RNDD,
+ 8.3210041599999988079e8);
+ check(1976245324, 1.25296395864546893357e+232, GMP_RNDZ,
+ -1.2529639586454686577e232);
+ check(2128997392, -1.08496826129284207724e+187, GMP_RNDU,
+ 1.0849682612928422704e187);
+ check(293607738, -1.9967571564050541e-5, GMP_RNDU, 2.9360773800002003e8);
+ check(354270183, 2.9469161763489528e3, GMP_RNDN, 3.5426723608382362e8);
+ return 0;
+}
+
diff --git a/mpfr/tests/tzeta.c b/mpfr/tests/tzeta.c
deleted file mode 100644
index d9495e57b..000000000
--- a/mpfr/tests/tzeta.c
+++ /dev/null
@@ -1,70 +0,0 @@
-/* Test file for mpfr_zeta.
-
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
-
-This file is part of the MPFR Library.
-
-The MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Library General Public License as published by
-the Free Software Foundation; either version 2 of the License, or (at your
-option) any later version.
-
-The MPFR 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 Library General Public
-License for more details.
-
-You should have received a copy of the GNU Library General Public License
-along with the MPFR 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 <math.h>
-#include "gmp.h"
-#include "gmp-impl.h"
-#include "longlong.h"
-#include "mpfr.h"
-
-/* #define DEBUG */
-
-int main()
-{mpfr_t p,result,res_p; int succes; size_t t;
-mpfr_init2(result, 53);
-mpfr_init2(res_p, 53);
-mpfr_init2(p, 53); mpfr_set_d(p, 2.0, GMP_RNDN);
-
-mpfr_zeta(result,p,GMP_RNDN);
-/*printf("%d\n",succes);*/
-#ifdef DEBUG
-printf("Valeur de zeta(2) avec prec=53 et arrondi au plus pres:\n");
-mpfr_print_raw(result);printf("\n");
-t=mpfr_out_str(stdout,10,0,result,GMP_RNDN);printf("\n");
-#endif
- if (mpfr_get_d(result) != 1.64493406684822640607e+00) {
- fprintf(stderr, "mpfr_zeta fails for s=2 and rnd_mode=GMP_RNDN\n");
- exit(1);
- }
-
-/*Test de comparaison avec pi^2/6*/
-mpfr_set_prec(p, 67);
-mpfr_pi(p,GMP_RNDN);
-mpfr_mul(p,p,p,GMP_RNDN);
-mpfr_set_ui(res_p,6,GMP_RNDN);
-mpfr_div(p,p,res_p,GMP_RNDN);
-/*mpfr_print_raw(p);printf("\n");
-t=mpfr_out_str(stdout,10,0,p,GMP_RNDN);printf("\n");*/
-succes=mpfr_can_round(p,63,GMP_RNDN,GMP_RNDN,53);
-mpfr_set(res_p,p,GMP_RNDN);
-/*printf("%d\n",succes);*/
-#ifdef DEBUG
-printf("Valeur de pi^2/6 avec prec=53 et arrondi au plus pres:\n");
-mpfr_print_raw(res_p);printf("\n");
-t=mpfr_out_str(stdout,10,0,res_p,GMP_RNDN);printf("\n");
-#endif
-
-mpfr_clear(p);
-mpfr_clear(result);
-mpfr_clear(res_p);
-exit (0);
-}
diff --git a/mpfr/zeta.c b/mpfr/zeta.c
deleted file mode 100644
index 5250a9dd7..000000000
--- a/mpfr/zeta.c
+++ /dev/null
@@ -1,116 +0,0 @@
-/* mpfr_zeta -- Riemann Zeta function at a floating-point number
-
-Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
-
-This file is part of the MPFR Library.
-
-The MPFR Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Library General Public License as published by
-the Free Software Foundation; either version 2 of the License, or (at your
-option) any later version.
-
-The MPFR 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 Library General Public
-License for more details.
-
-You should have received a copy of the GNU Library General Public License
-along with the MPFR 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 <math.h>
-#include "gmp.h"
-#include "gmp-impl.h"
-#include "longlong.h"
-#include "mpfr.h"
-
-int
-#if __STDC__
-mpfr_zeta(mpfr_ptr result, mpfr_srcptr op, unsigned char rnd_mode)
-#else
-mpfr_zeta(result, op, rnd_mode)
- mpfr_ptr result;
- mpfr_srcptr op;
- unsigned char rnd_mode;
-#endif
-{
- mpfr_t s,s2,x,y,u,b,v,nn,z,z2;
- int i,n,succes;
-
- /* first version */
- if (mpfr_get_d(op) != 2.0 || rnd_mode != GMP_RNDN
- || PREC(result) != 53) {
- fprintf(stderr, "not yet implemented\n"); exit(1);
- }
-
-mpfr_set_default_prec(67);
-mpfr_init(x);
-mpfr_init(y);
-mpfr_init(s);
-mpfr_init(s2);
-mpfr_init(u);
-mpfr_init(b);
-mpfr_init(v);
-mpfr_init(nn);
-mpfr_init(z);
-mpfr_init(z2);
-mpfr_set_ui(u,1,GMP_RNDN);
-mpfr_set_ui(s,0,GMP_RNDN);
-/*s=Somme des 1/i^2 (i=100...2)*/
-n=100;
-for (i=n; i>1; i--)
-{
-mpfr_set_ui(x,i*i,GMP_RNDN);
-mpfr_div(y,u,x,GMP_RNDN);
-mpfr_add(s,s,y,GMP_RNDN);
-};
-/*mpfr_print_raw(s);printf("\n");
-t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
-/*Application d'Euler-Maclaurin, jusqu'au terme 1/n^7 - n=100)*/
-mpfr_set_ui(nn,n,GMP_RNDN);
-mpfr_div(z,u,nn,GMP_RNDN);
-mpfr_set(s2,z,GMP_RNDN);
-mpfr_mul(z2,z,z,GMP_RNDN);
-mpfr_div_2exp(v,z2,1,GMP_RNDN);
-mpfr_sub(s2,s2,v,GMP_RNDN);
-mpfr_set_ui(b,6,GMP_RNDN);
-mpfr_mul(z,z,z2,GMP_RNDN);
-mpfr_div(v,z,b,GMP_RNDN);
-mpfr_add(s2,s2,v,GMP_RNDN);
-mpfr_set_si(b,-30,GMP_RNDN);
-mpfr_mul(z,z,z2,GMP_RNDN);
-mpfr_div(v,z,b,GMP_RNDN);
-mpfr_add(s2,s2,v,GMP_RNDN);
-mpfr_set_si(b,42,GMP_RNDN);
-mpfr_mul(z,z,z2,GMP_RNDN);
-mpfr_div(v,z,b,GMP_RNDN);
-mpfr_add(s2,s2,v,GMP_RNDN);
-/*mpfr_print_raw(s2);printf("\n");
-t=mpfr_out_str(stdout,10,0,s2,GMP_RNDN);printf("\n");*/
-mpfr_add(s,s,s2,GMP_RNDN);
-/*mpfr_print_raw(s);printf("\n");
-t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
-mpfr_add(s,s,u,GMP_RNDN);
-/*mpfr_print_raw(s);printf("\n");
-t=mpfr_out_str(stdout,10,0,s,GMP_RNDN);printf("\n");*/
-/*Peut-on arrondir ? La reponse est oui*/
-succes=mpfr_can_round(s,57,GMP_RNDN,GMP_RNDN,53);
-if (succes) mpfr_set(result,s,GMP_RNDN);
-else {
- fprintf(stderr, "can't round in mpfr_zeta\n"); exit(1);
-}
-
-mpfr_clear(x);
-mpfr_clear(y);
-mpfr_clear(s);
-mpfr_clear(s2);
-mpfr_clear(u);
-mpfr_clear(b);
-mpfr_clear(v);
-mpfr_clear(nn);
-mpfr_clear(z);
-mpfr_clear(z2);
-return 1; /* result is inexact */
-}