diff options
author | Kevin Ryde <user42@zip.com.au> | 2001-02-01 23:31:48 +0100 |
---|---|---|
committer | Kevin Ryde <user42@zip.com.au> | 2001-02-01 23:31:48 +0100 |
commit | 4b133dd08ea286095823032802a0fe78b15e2715 (patch) | |
tree | 5a645599aa66bd9f76c9c2cec26cefb63c3c1f3f /mpfr | |
parent | 618be102b4c29f1e02e4a7d2244ae57b4bb544c9 (diff) | |
download | gmp-4b133dd08ea286095823032802a0fe78b15e2715.tar.gz |
* mpfr/*: Update to mpfr 2001.
Diffstat (limited to 'mpfr')
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 */ -} |