system/freebsd-src
system/freebsd-src
1
0
Fork 0
mirror of https://github.com/freebsd/freebsd-src synced 2024-09-29 21:26:51 +00:00
Code Issues Packages Projects Releases Wiki Activity

GNU MP (Multiprecision) library. This is needed by secure RPC (being

Browse source 
done by Bill Paul) and various other BSD programs.
Obtained from:FSF
This commit is contained in:
Mark Murray 1995-11-12 14:40:41 +00:00
parent 4540b99acc
commit ae82e96f8c
Notes: svn2git 2020-12-20 02:59:44 +00:00 
svn path=/vendor/misc-GNU/dist2/; revision=12234
131 changed files with 16601 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

339
gnu/lib/libgmp/COPYING Normal file
View file

@ -0,0 +1,339 @@
GNU GENERAL PUBLIC LICENSE
Version 2, June 1991
Copyright (C) 1989, 1991 Free Software Foundation, Inc.
675 Mass Ave, Cambridge, MA 02139, USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
Preamble
The licenses for most software are designed to take away your
freedom to share and change it. By contrast, the GNU General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users. This
General Public License applies to most of the Free Software
Foundation's software and to any other program whose authors commit to
using it. (Some other Free Software Foundation software is covered by
the GNU Library General Public License instead.) You can apply it to
your programs, 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 software, or if you modify it.
For example, if you distribute copies of such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have. You must make sure that they, too, receive or can get the
source code. And you must show them these terms so they know their
rights.
We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.
Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software. If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, 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 redistributors of a free
program will individually obtain patent licenses, in effect making the
program proprietary. To prevent this, we have made it clear that any
patent must be licensed for everyone's free use or not licensed at all.
The precise terms and conditions for copying, distribution and
modification follow.
GNU GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License applies to any program or other work which contains
a notice placed by the copyright holder saying it may be distributed
under the terms of this General Public License. The "Program", below,
refers to any such program or work, and a "work based on the Program"
means either the Program or any derivative work under copyright law:
that is to say, a work containing the Program or a portion of it,
either verbatim or with modifications and/or translated into another
language. (Hereinafter, translation is included without limitation in
the term "modification".) Each licensee is addressed as "you".
Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope. The act of
running the Program is not restricted, and the output from the Program
is covered only if its contents constitute a work based on the
Program (independent of having been made by running the Program).
Whether that is true depends on what the Program does.
1. You may copy and distribute verbatim copies of the Program's
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 give any other recipients of the Program a copy of this License
along with the Program.
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 Program or any portion
of it, thus forming a work based on the Program, 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) You must cause the modified files to carry prominent notices
stating that you changed the files and the date of any change.
b) You must cause any work that you distribute or publish, that in
whole or in part contains or is derived from the Program or any
part thereof, to be licensed as a whole at no charge to all third
parties under the terms of this License.
c) If the modified program normally reads commands interactively
when run, you must cause it, when started running for such
interactive use in the most ordinary way, to print or display an
announcement including an appropriate copyright notice and a
notice that there is no warranty (or else, saying that you provide
a warranty) and that users may redistribute the program under
these conditions, and telling the user how to view a copy of this
License. (Exception: if the Program itself is interactive but
does not normally print such an announcement, your work based on
the Program is not required to print an announcement.)
These requirements apply to the modified work as a whole. If
identifiable sections of that work are not derived from the Program,
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 Program, 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 Program.
In addition, mere aggregation of another work not based on the Program
with the Program (or with a work based on the Program) on a volume of
a storage or distribution medium does not bring the other work under
the scope of this License.
3. You may copy and distribute the Program (or a work based on it,
under Section 2) in object code or executable form under the terms of
Sections 1 and 2 above provided that you also do one of the following:
a) 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; or,
b) Accompany it with a written offer, valid for at least three
years, to give any third party, for a charge no more than your
cost of physically performing source distribution, a complete
machine-readable copy of the corresponding source code, to be
distributed under the terms of Sections 1 and 2 above on a medium
customarily used for software interchange; or,
c) Accompany it with the information you received as to the offer
to distribute corresponding source code. (This alternative is
allowed only for noncommercial distribution and only if you
received the program in object code or executable form with such
an offer, in accord with Subsection b above.)
The source code for a work means the preferred form of the work for
making modifications to it. For an executable work, 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 executable. 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.
If distribution of executable or 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 counts as
distribution of the source code, even though third parties are not
compelled to copy the source along with the object code.
4. You may not copy, modify, sublicense, or distribute the Program
except as expressly provided under this License. Any attempt
otherwise to copy, modify, sublicense or distribute the Program 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.
5. 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 Program or its derivative works. These actions are
prohibited by law if you do not accept this License. Therefore, by
modifying or distributing the Program (or any work based on the
Program), you indicate your acceptance of this License to do so, and
all its terms and conditions for copying, distributing or modifying
the Program or works based on it.
6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the
original licensor to copy, distribute or modify the Program 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.
7. 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 Program at all. For example, if a patent
license would not permit royalty-free redistribution of the Program 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 Program.
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.
8. If the distribution and/or use of the Program is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Program 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.
9. The Free Software Foundation may publish revised and/or new versions
of the 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 Program
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 Program does not specify a version number of
this License, you may choose any version ever published by the Free Software
Foundation.
10. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, 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
11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM "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 PROGRAM IS WITH YOU. SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.
12. 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 PROGRAM 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 PROGRAM (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 PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
END OF TERMS AND CONDITIONS
Appendix: How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
free software which everyone can redistribute and change under these terms.
To do so, attach the following notices to the program. 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 program's name and a brief idea of what it does.>
Copyright (C) 19yy <name of author>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
Gnomovision version 69, Copyright (C) 19yy name of author
Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
The hypothetical commands `show w' and `show c' should show the appropriate
parts of the General Public License. Of course, the commands you use may
be called something other than `show w' and `show c'; they could even be
mouse-clicks or menu items--whatever suits your program.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the program, if
necessary. Here is a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the program
`Gnomovision' (which makes passes at compilers) written by James Hacker.
<signature of Ty Coon>, 1 April 1989
Ty Coon, President of Vice
This General Public License does not permit incorporating your program into
proprietary programs. If your program is a subroutine library, you may
consider it more useful to permit linking proprietary applications with the
library. If this is what you want to do, use the GNU Library General
Public License instead of this License.

1347
gnu/lib/libgmp/ChangeLog Normal file
View file

File diff suppressed because it is too large Load diff

34
gnu/lib/libgmp/INSTALL Normal file
View file

@ -0,0 +1,34 @@
Here is how to compile GNU MP.
You probably want to use the GNU C compiler to build this library.
With other compilers the speed of the library will be 3-10 times
slower for many CPU:s. The reason for this is that the GNU C compiler
will use inline assembler for some important operations, while other C
compilers will have to stick to plain C code.
This is how to build the library:
Type "make" to build libgmp.a and libmp.a. The former is the main
GNU MP library. The latter is the Berkeley MP compatible library.
If you don't have GCC, type "make CC=cc". The compilation should, at
least with GCC, proceed without any kind of warnings from the compiler
programs. On the DEC Alpha, you have to use GCC because of bugs in DEC's
own compiler. GCC 2.3.3 for x86, Alpha, and HP-PA has bugs that make
several functions be mis-optimized. Later version of GCC does not have
this problem.
To build and run the tests, do "make check".
The documentation is an a texinfo file, gmp.texi.
To create the documentation from the texinfo source, type "make doc".
This requires the "tex" and "makeinfo" commands to be available in
your search path. If you have only one of them, you can create the
dvi file (for the paper manual) with "make gmp.dvi", and the info file
(for the GNU online manual facility) with "make gmp.info".
You need version 2.06 or later of texinfo in order to build the
documentation.
Please report problems to tege@gnu.ai.mit.edu.

61
gnu/lib/libgmp/README Normal file
View file

@ -0,0 +1,61 @@
THE GNU MP LIBRARY
GNU MP is a library for arbitrary precision arithmetic, operating on
signed integers and rational numbers. It has a rich set of functions,
and the functions have a regular interface.
I have tried to make these functions as fast as possible, both for small
operands and for huge operands. The speed is achieved by using fullwords
as the basic arithmetic type, by using fast algorithms, by defining inline
assembler for mixed sized multiplication and division (i.e 32*32->64 bit
multiplication and 64/32->32,32 bit division), and by hacking the code
with emphasis on speed (and not simplicity and elegance).
The speed of GNU MP is about 5 to 100 times that of Berkeley MP for
small operands. The speed-up increases with the operand sizes for
certain operations, for which GNU MP has asymptotically faster algorithms.
There are four classes of functions in GNU MP.
1. Signed integer arithmetic functions, mpz_*. The set of functions are
intended to be easy to use, being rich and regular.
To use these functions, include the file "gmp.h".
2. Rational arithmetic functions, mpq_*. For now, just a small set of
functions necessary for basic rational arithmetics.
To use these functions, include the file "gmp.h".
3. Positive-integer, low-level, harder-to-use, but for small operands
about twice as fast than the mpz_* functions are the functions in the
mpn_* class. No memory management is performed. The caller must
ensure enough space is available for the results. The set of
functions is not quite regular, nor is the calling interface. These
functions accept input arguments in the form of pairs consisting of a
pointer to the least significant word, and a integral size telling how
many limbs (= words) the pointer points to.
Almost all calculations, in the entire package, are made in these
low-level functions.
These functions are not fully documented in this release. They will
probably be so in a future release.
4. Berkeley MP compatible functions.
To use these functions, include the file "mp.h". You can test if you
are using the GNU version by testing if the symbol __GNU_MP__ is
defined.
REPORTING BUGS
If you find a bug in the library, please make sure to tell us about it!
You can report bugs, and propose modifications and enhancements to
tege@gnu.ai.mit.edu. How to report a bug is further described in
the texinfo documentation, see the file gmp.texi.

184
gnu/lib/libgmp/TODO Normal file
View file

@ -0,0 +1,184 @@
THINGS TO WORK ON
Note that many of these things mentioned here are already fixed in GMP 2.0.
* Improve speed for non-gcc compilers by defining umul_ppmm, udiv_qrnnd,
etc, to call __umul_ppmm, __udiv_qrnnd. A typical definition for
umul_ppmm would be
#define umul_ppmm(ph,pl,m0,m1) \
{unsigned long __ph; (pl) = __umul_ppmm (&__ph, (m0), (m1)); (ph) = __ph;}
In order to maintain just one version of longlong.h (gmp and gcc), this
has to be done outside of longlong.h.
* Change mpn-routines to not deal with normalisation?
mpn_add: Unchanged.
mpn_sub: Remove normalization loop. Does it assume normalised input?
mpn_mul: Make it return most sign limb, to simplify normalisation.
Karatsubas algorith will be greatly simplified if mpn_add and
mpn_sub doesn't normalise their results.
mpn_div: Still requires strict normalisation.
Beware of problems with mpn_cmp (and similar), a larger size does not
ensure that an operand is larger, since it may be "less normalised".
Normalization has to be moved into mpz-functions.
Bennet Yee at CMU proposes:
* mpz_{put,get}_raw for memory oriented I/O like other *_raw functions.
* A function mpfatal that is called for exceptions. The user may override
the default definition.
* mout should group in 10-digit groups.
* ASCII dependence?
* Error reporting from I/O functions (linkoping)?
* Make all computation mpz_* functions return a signed int indicating if
the result was zero, positive, or negative?
* Implement mpz_cmpabs, mpz_xor, mpz_to_double, mpz_to_si, mpz_lcm,
mpz_dpb, mpz_ldb, various bit string operations like mpz_cntbits. Also
mpz_@_si for most @??
Brian Beuning proposes:
1. An array of small primes
3. A function to factor an MINT
4. A routine to look for "small" divisors of an MINT
5. A 'multiply mod n' routine based on Montgomery's algorithm.
Doug Lea proposes:
1. A way to find out if an integer fits into a signed int, and if so, a
way to convert it out.
2. Similarly for double precision float conversion.
3. A function to convert the ratio of two integers to a double. This
can be useful for mixed mode operations with integers, rationals, and
doubles.
5. Bit-setting, clearing, and testing operations, as in
mpz_setbit(MP_INT* dest, MP_INT* src, unsigned long bit_number),
and used, for example in
mpz_setbit(x, x, 123)
to directly set the 123rd bit of x.
If these are supported, you don't first have to set up
an otherwise unnecessary mpz holding a shifted value, then
do an "or" operation.
Elliptic curve method descrition in the Chapter `Algorithms in Number
Theory' in the Handbook of Theoretical Computer Science, Elsevier,
Amsterdam, 1990. Also in Carl Pomerance's lecture notes on Cryptology and
Computational Number Theory, 1990.
* New function: mpq_get_ifstr (int_str, frac_str, base,
precision_in_som_way, rational_number). Convert RATIONAL_NUMBER to a
string in BASE and put the integer part in INT_STR and the fraction part
in FRAC_STR. (This function would do a division of the numerator and the
denominator.)
* Should mpz_powm* handle negative exponents?
* udiv_qrnnd: If the denominator is normalized, the n0 argument has very
little effect on the quotient. Maybe we can assume it is 0, and
compensate at a later stage?
* Better sqrt: First calculate the reciprocal square root, then multiply by
the operand to get the square root. The reciprocal square root can be
obtained through Newton-Raphson without division. The iteration is x :=
x*(3-a*x^2)/2, where a is the operand.
* Newton-Raphson using multiplication: We get twice as many correct digits
in each iteration. So if we square x(k) as part of the iteration, the
result will have the leading digits in common with the entire result from
iteration k-1. A _mpn_mul_lowpart could implement this.
* Peter Montgomery: If 0 <= a, b < p < 2^31 and I want a modular product
a*b modulo p and the long long type is unavailable, then I can write
typedef signed long slong;
typedef unsigned long ulong;
slong a, b, p, quot, rem;
quot = (slong) (0.5 + (double)a * (double)b / (double)p);
rem = (slong)((ulong)a * (ulong)b - (ulong)p * (ulong)q);
if (rem < 0} {rem += p; quot--;}
FFT:
{
* Multiplication could be done with Montgomery's method combined with
the "three primes" method described in Lipson. Maybe this would be
faster than to Nussbaumer's method with 3 (simple) moduli?
* Maybe the modular tricks below are not needed: We are using very
special numbers, Fermat numbers with a small base and a large exponent,
and maybe it's possible to just subtract and add?
* Modify Nussbaumer's convolution algorithm, to use 3 words for each
coefficient, calculating in 3 relatively prime moduli (e.g.
0xffffffff, 0x100000000, and 0x7fff on a 32-bit computer). Both all
operations and CRR would be very fast with such numbers.
* Optimize the Shoenhage-Stassen multiplication algorithm. Take
advantage of the real valued input to save half of the operations and
half of the memory. Try recursive variants with large, optimized base
cases. Use recursive FFT with large base cases, since recursive FFT
has better memory locality. A normal FFT get 100% cache miss.
}
* Speed modulo arithmetic, using Montgomery's method or my pre-invertion
method. In either case, special arithmetic calls would be needed,
mpz_mmmul, mpz_mmadd, mpz_mmsub, plus some kind of initialization
functions.
* mpz_powm* should not use division to reduce the result in the loop, but
instead pre-compute the reciprocal of the MOD argument and do reduced_val
= val-val*reciprocal(MOD)*MOD, or use Montgomery's method.
* mpz_mod_2expplussi -- to reduce a bignum modulo (2**n)+s
* It would be a quite important feature never to allocate more memory than
really necessary for a result. Sometimes we can achieve this cheaply, by
deferring reallocation until the result size is known.
* New macro in longlong.h: shift_rhl that extracts a word by shifting two
words as a unit. (Supported by i386, i860, HP-PA, RS6000, 29k.) Useful
for shifting multiple precision numbers.
* The installation procedure should make a test run of multiplication to
decide the threshold values for algorithm switching between the available
methods.
* The gcd algorithm could probably be improved with a divide-and-conquer
(DAC) approach. At least the bulk of the operations should be done with
single precision.
* Fast output conversion of x to base B:
1. Find n, such that (B^n > x).
2. Set y to (x*2^m)/(B^n), where m large enough to make 2^n ~~ B^n
3. Multiply the low half of y by B^(n/2), and recursively convert the
result. Truncate the low half of y and convert that recursively.
Complexity: O(M(n)log(n))+O(D(n))!
* Extensions for floating-point arithmetic.
* Improve special cases for division.
1. When the divisor is just one word, normalization is not needed for
most CPUs, and can be done in the division loop for CPUs that need
normalization.
2. Even when the result is going to be very small, (i.e. nsize-dsize is
small) normalization should also be done in the division loop.
To fix this, a new routine mpn_div_unnormalized is needed.
* Never allocate temporary space for a source param that overlaps with a
destination param needing reallocation. Instead malloc a new block for
the destination (and free the source before returning to the caller).
* When any of the source operands overlap with the destination, mult (and
other routines) slow down. This is so because the need of temporary
allocation (with alloca) and copying. If a new destination were
malloc'ed instead (and the overlapping source free'd before return) no
copying would be needed. Is GNU malloc quick enough to make this faster
even for reasonably small operands?
Local Variables:
mode: text
fill-column: 75
version-control: never
End:

1
gnu/lib/libgmp/VERSION Normal file
View file

@ -0,0 +1 @@
GNU MP version 1.3.2

309
gnu/lib/libgmp/_mpz_get_str.c Normal file
View file

@ -0,0 +1,309 @@
/* _mpz_get_str (string, base, mp_src) -- Convert the multiple precision
number MP_SRC to a string STRING of base BASE. If STRING is NULL
allocate space for the result. In any case, return a pointer to the
result. If STRING is not NULL, the caller must ensure enough space is
available to store the result.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#ifndef UMUL_TIME
#define UMUL_TIME 1
#endif
#ifndef UDIV_TIME
#define UDIV_TIME UMUL_TIME
#endif
#define udiv_qrnndx(q, r, nh, nl, d, di) \
do { \
unsigned long int _q, _ql, _r; \
unsigned long int _xh, _xl; \
umul_ppmm (_q, _ql, (nh), (di)); \
_q += (nh); /* DI is 2**32 too small. Compensate */\
if (_q < (nh)) \
{ \
/* Got carry. Propagate it in the multiplication. */ \
umul_ppmm (_xh, _xl, (d), _q); \
_xh += (d); \
} \
else \
umul_ppmm (_xh, _xl, (d), _q); \
sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl); \
if (_xh != 0) \
{ \
sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
_q += 1; \
if (_xh != 0) \
{ \
sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
_q += 1; \
} \
} \
if (_r >= (d)) \
{ \
_r -= (d); \
_q += 1; \
} \
(r) = _r; \
(q) = _q; \
} while (0)
char *
#ifdef __STDC__
_mpz_get_str (char *str, int base, const MP_INT *m)
#else
_mpz_get_str (str, base, m)
char *str;
int base;
const MP_INT *m;
#endif
{
mp_ptr tp;
mp_size msize;
mp_limb big_base;
#if UDIV_NEEDS_NORMALIZATION || UDIV_TIME > 2 * UMUL_TIME
int normalization_steps;
#if UDIV_TIME > 2 * UMUL_TIME
mp_limb big_base_inverted;
#endif
#endif
unsigned int dig_per_u;
mp_size out_len;
char *s;
char *num_to_ascii;
if (base >= 0)
{
if (base == 0)
base = 10;
num_to_ascii = "0123456789abcdefghijklmnopqrstuvwxyz";
}
else
{
base = -base;
num_to_ascii = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
}
dig_per_u = __mp_bases[base].chars_per_limb;
out_len = mpz_sizeinbase (m, base) + 1;
big_base = __mp_bases[base].big_base;
msize = m->size;
if (str == NULL)
str = (char *) (*_mp_allocate_func) (out_len + (msize < 0));
if (msize < 0)
*str++ = '-';
s = str;
msize = ABS (msize);
/* Special case zero, as the code below doesn't handle it. */
if (msize == 0)
{
s[0] = '0';
s[1] = 0;
return str;
}
if ((base & (base - 1)) == 0)
{
/* The base is a power of 2. Make conversion from most
significant side. */
mp_limb n1, n0;
int bits_per_digit = big_base;
int x;
int bit_pos;
int i;
unsigned mask = (1 << bits_per_digit) - 1;
tp = m->d;
n1 = tp[msize - 1];
count_leading_zeros (x, n1);
/* BIT_POS should be R when input ends in least sign. nibble,
R + bits_per_digit * n when input ends in n:th least significant
nibble. */
{
int bits;
bits = BITS_PER_MP_LIMB * msize - x;
x = bits % bits_per_digit;
if (x != 0)
bits += bits_per_digit - x;
bit_pos = bits - (msize - 1) * BITS_PER_MP_LIMB;
}
/* Fast loop for bit output. */
i = msize - 1;
for (;;)
{
bit_pos -= bits_per_digit;
while (bit_pos >= 0)
{
*s++ = num_to_ascii[(n1 >> bit_pos) & mask];
bit_pos -= bits_per_digit;
}
i--;
if (i < 0)
break;
n0 = (n1 << -bit_pos) & mask;
n1 = tp[i];
bit_pos += BITS_PER_MP_LIMB;
*s++ = num_to_ascii[n0 | (n1 >> bit_pos)];
}
*s = 0;
}
else
{
/* General case. The base is not a power of 2. Make conversion
from least significant end. */
/* If udiv_qrnnd only handles divisors with the most significant bit
set, prepare BIG_BASE for being a divisor by shifting it to the
left exactly enough to set the most significant bit. */
#if UDIV_NEEDS_NORMALIZATION || UDIV_TIME > 2 * UMUL_TIME
count_leading_zeros (normalization_steps, big_base);
big_base <<= normalization_steps;
#if UDIV_TIME > 2 * UMUL_TIME
/* Get the fixed-point approximation to 1/BIG_BASE. */
big_base_inverted = __mp_bases[base].big_base_inverted;
#endif
#endif
out_len--; /* now not include terminating \0 */
s += out_len;
/* Allocate temporary space and move the multi prec number to
convert there, as we need to overwrite it below, while
computing the successive remainders. */
tp = (mp_ptr) alloca ((msize + 1) * BYTES_PER_MP_LIMB);
MPN_COPY (tp, m->d, msize);
while (msize != 0)
{
int i;
mp_limb n0, n1;
#if UDIV_NEEDS_NORMALIZATION || UDIV_TIME > 2 * UMUL_TIME
/* If we shifted BIG_BASE above, shift the dividend too, to get
the right quotient. We need to do this every loop,
as the intermediate quotients are OK, but the quotient from
one turn in the loop is going to be the dividend in the
next turn, and the dividend needs to be up-shifted. */
if (normalization_steps != 0)
{
n0 = mpn_lshift (tp, tp, msize, normalization_steps);
/* If the shifting gave a carry out limb, store it and
increase the length. */
if (n0 != 0)
{
tp[msize] = n0;
msize++;
}
}
#endif
/* Divide the number at TP with BIG_BASE to get a quotient and a
remainder. The remainder is our new digit in base BIG_BASE. */
i = msize - 1;
n1 = tp[i];
if (n1 >= big_base)
n1 = 0;
else
{
msize--;
i--;
}
for (; i >= 0; i--)
{
n0 = tp[i];
#if UDIV_TIME > 2 * UMUL_TIME
udiv_qrnndx (tp[i], n1, n1, n0, big_base, big_base_inverted);
#else
udiv_qrnnd (tp[i], n1, n1, n0, big_base);
#endif
}
#if UDIV_NEEDS_NORMALIZATION || UDIV_TIME > 2 * UMUL_TIME
/* If we shifted above (at previous UDIV_NEEDS_NORMALIZATION tests)
the remainder will be up-shifted here. Compensate. */
n1 >>= normalization_steps;
#endif
/* Convert N1 from BIG_BASE to a string of digits in BASE
using single precision operations. */
for (i = dig_per_u - 1; i >= 0; i--)
{
*--s = num_to_ascii[n1 % base];
n1 /= base;
/* Break from the loop as soon as we would only write zeros. */
if (n1 == 0 && msize == 0)
break;
}
}
/* There should be no leading zeros. */
if (*s == '0')
abort ();
if (s == str)
{
/* This should be the common case. */
s[out_len] = 0;
}
else if (s == str + 1)
{
/* The string became 1 digit shorter than its maximum. */
/* Need to copy it back one char pos. */
out_len--;
#ifndef HAS_MEMMOVE
{
size_t i;
for (i = 0; i < out_len; i++)
str[i] = s[i];
}
#else
memmove (str, s, out_len);
#endif
str[out_len] = 0;
}
else
{
/* Hopefully never. */
abort ();
}
}
alloca (0);
/* Ugly, we incremented str for negative numbers. Fix that here. */
return str - (m->size < 0);
}

258
gnu/lib/libgmp/_mpz_set_str.c Normal file
View file

@ -0,0 +1,258 @@
/* _mpz_set_str(mp_dest, string, base) -- Convert the \0-terminated
string STRING in base BASE to multiple precision integer in
MP_DEST. Allow white space in the string. If BASE == 0 determine
the base in the C standard way, i.e. 0xhh...h means base 16,
0oo...o means base 8, otherwise assume base 10.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
enum char_type
{
XX = -3,
SPC = -2,
EOF = -1
};
static signed char ascii_to_num[256] =
{
EOF,XX, XX, XX, XX, XX, XX, XX, XX, SPC,SPC,XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
SPC,XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, XX, XX, XX, XX, XX, XX,
XX, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, XX, XX, XX, XX, XX,
XX, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX,
XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX, XX
};
int
#ifdef __STDC__
_mpz_set_str (MP_INT *x, const char *str, int base)
#else
_mpz_set_str (x, str, base)
MP_INT *x;
const char *str;
int base;
#endif
{
mp_ptr xp;
mp_size size;
mp_limb big_base;
int indigits_per_limb;
int negative = 0;
int inp_rawchar;
mp_limb inp_digit;
mp_limb res_digit;
size_t str_len;
mp_size i;
if (str[0] == '-')
{
negative = 1;
str++;
}
if (base == 0)
{
if (str[0] == '0')
{
if (str[1] == 'x' || str[1] == 'X')
base = 16;
else
base = 8;
}
else
base = 10;
}
big_base = __mp_bases[base].big_base;
indigits_per_limb = __mp_bases[base].chars_per_limb;
str_len = strlen (str);
size = str_len / indigits_per_limb + 1;
if (x->alloc < size)
_mpz_realloc (x, size);
xp = x->d;
size = 0;
if ((base & (base - 1)) == 0)
{
/* The base is a power of 2. Read the input string from
least to most significant character/digit. */
const char *s;
int next_bitpos;
int bits_per_indigit = big_base;
/* Accept and ignore 0x or 0X before hexadecimal numbers. */
if (base == 16 && str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
{
str += 2;
str_len -= 2;
}
res_digit = 0;
next_bitpos = 0;
for (s = str + str_len - 1; s >= str; s--)
{
inp_rawchar = *s;
inp_digit = ascii_to_num[inp_rawchar];
if (inp_digit >= base)
{
/* Was it white space? Just ignore it. */
if ((char) inp_digit == (char) SPC)
continue;
/* We found rubbish in the string. Return -1 to indicate
the error. */
return -1;
}
res_digit |= inp_digit << next_bitpos;
next_bitpos += bits_per_indigit;
if (next_bitpos >= BITS_PER_MP_LIMB)
{
xp[size] = res_digit;
size++;
next_bitpos -= BITS_PER_MP_LIMB;
res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
}
}
xp[size] = res_digit;
size++;
for (i = size - 1; i >= 0; i--)
{
if (xp[i] != 0)
break;
}
size = i + 1;
}
else
{
/* General case. The base is not a power of 2. */
mp_size i;
int j;
mp_limb cy;
for (;;)
{
res_digit = 0;
for (j = 0; j < indigits_per_limb; )
{
inp_rawchar = (unsigned char) *str++;
inp_digit = ascii_to_num[inp_rawchar];
/* Negative means that the character was not a proper digit. */
if (inp_digit >= base)
{
/* Was it white space? Just ignore it. */
if ((char) inp_digit == (char) SPC)
continue;
goto end_or_error;
}
res_digit = res_digit * base + inp_digit;
/* Increment the loop counter here, since it mustn't be
incremented when we do "continue" above. */
j++;
}
cy = res_digit;
/* Insert RES_DIGIT into the result multi prec integer. */
for (i = 0; i < size; i++)
{
mp_limb p1, p0;
umul_ppmm (p1, p0, big_base, xp[i]);
p0 += cy;
cy = p1 + (p0 < cy);
xp[i] = p0;
}
if (cy != 0)
{
xp[size] = cy;
size++;
}
}
end_or_error:
/* We probably have some digits in RES_DIGIT (J tells how many). */
if ((char) inp_digit != (char) EOF)
{
/* Error return. */
return -1;
}
/* J contains number of digits (in base BASE) remaining in
RES_DIGIT. */
if (j > 0)
{
big_base = 1;
do
{
big_base *= base;
j--;
}
while (j > 0);
cy = res_digit;
/* Insert ultimate RES_DIGIT into the result multi prec integer. */
for (i = 0; i < size; i++)
{
mp_limb p1, p0;
umul_ppmm (p1, p0, big_base, xp[i]);
p0 += cy;
cy = p1 + (p0 < cy);
xp[i] = p0;
}
if (cy != 0)
{
xp[size] = cy;
size++;
}
}
}
if (negative)
size = -size;
x->size = size;
return 0;
}

466
gnu/lib/libgmp/alloca.c Normal file
View file

@ -0,0 +1,466 @@
/* alloca.c -- allocate automatically reclaimed memory
(Mostly) portable public-domain implementation -- D A Gwyn
This implementation of the PWB library alloca function,
which is used to allocate space off the run-time stack so
that it is automatically reclaimed upon procedure exit,
was inspired by discussions with J. Q. Johnson of Cornell.
J.Otto Tennant <jot@cray.com> contributed the Cray support.
There are some preprocessor constants that can
be defined when compiling for your specific system, for
improved efficiency; however, the defaults should be okay.
The general concept of this implementation is to keep
track of all alloca-allocated blocks, and reclaim any
that are found to be deeper in the stack than the current
invocation. This heuristic does not reclaim storage as
soon as it becomes invalid, but it will do so eventually.
As a special case, alloca(0) reclaims storage without
allocating any. It is a good idea to use alloca(0) in
your main control loop, etc. to force garbage collection. */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
/* If compiling with GCC, this file's not needed. */
#ifndef alloca
#ifdef emacs
#ifdef static
/* actually, only want this if static is defined as ""
-- this is for usg, in which emacs must undefine static
in order to make unexec workable
*/
#ifndef STACK_DIRECTION
you
lose
-- must know STACK_DIRECTION at compile-time
#endif /* STACK_DIRECTION undefined */
#endif /* static */
#endif /* emacs */
#ifdef emacs
#define free xfree
#endif
/* If your stack is a linked list of frames, you have to
provide an "address metric" ADDRESS_FUNCTION macro. */
#ifdef CRAY
long i00afunc ();
#define ADDRESS_FUNCTION(arg) (char *) i00afunc (&(arg))
#else
#define ADDRESS_FUNCTION(arg) &(arg)
#endif
#if __STDC__
typedef void *pointer;
#else
typedef char *pointer;
#endif
#define NULL 0
extern pointer (*_mp_allocate_func) ();
/* Define STACK_DIRECTION if you know the direction of stack
growth for your system; otherwise it will be automatically
deduced at run-time.
STACK_DIRECTION > 0 => grows toward higher addresses
STACK_DIRECTION < 0 => grows toward lower addresses
STACK_DIRECTION = 0 => direction of growth unknown */
#ifndef STACK_DIRECTION
#define STACK_DIRECTION 0 /* Direction unknown. */
#endif
#if STACK_DIRECTION != 0
#define STACK_DIR STACK_DIRECTION /* Known at compile-time. */
#else /* STACK_DIRECTION == 0; need run-time code. */
static int stack_dir; /* 1 or -1 once known. */
#define STACK_DIR stack_dir
static void
find_stack_direction ()
{
static char *addr = NULL; /* Address of first `dummy', once known. */
auto char dummy; /* To get stack address. */
if (addr == NULL)
{ /* Initial entry. */
addr = ADDRESS_FUNCTION (dummy);
find_stack_direction (); /* Recurse once. */
}
else
{
/* Second entry. */
if (ADDRESS_FUNCTION (dummy) > addr)
stack_dir = 1; /* Stack grew upward. */
else
stack_dir = -1; /* Stack grew downward. */
}
}
#endif /* STACK_DIRECTION == 0 */
/* An "alloca header" is used to:
(a) chain together all alloca'ed blocks;
(b) keep track of stack depth.
It is very important that sizeof(header) agree with malloc
alignment chunk size. The following default should work okay. */
#ifndef ALIGN_SIZE
#define ALIGN_SIZE sizeof(double)
#endif
typedef union hdr
{
char align[ALIGN_SIZE]; /* To force sizeof(header). */
struct
{
union hdr *next; /* For chaining headers. */
char *deep; /* For stack depth measure. */
} h;
} header;
static header *last_alloca_header = NULL; /* -> last alloca header. */
/* Return a pointer to at least SIZE bytes of storage,
which will be automatically reclaimed upon exit from
the procedure that called alloca. Originally, this space
was supposed to be taken from the current stack frame of the
caller, but that method cannot be made to work for some
implementations of C, for example under Gould's UTX/32. */
pointer
alloca (size)
unsigned size;
{
auto char probe; /* Probes stack depth: */
register char *depth = ADDRESS_FUNCTION (probe);
#if STACK_DIRECTION == 0
if (STACK_DIR == 0) /* Unknown growth direction. */
find_stack_direction ();
#endif
/* Reclaim garbage, defined as all alloca'd storage that
was allocated from deeper in the stack than currently. */
{
register header *hp; /* Traverses linked list. */
for (hp = last_alloca_header; hp != NULL;)
if ((STACK_DIR > 0 && hp->h.deep > depth)
|| (STACK_DIR < 0 && hp->h.deep < depth))
{
register header *np = hp->h.next;
free ((pointer) hp); /* Collect garbage. */
hp = np; /* -> next header. */
}
else
break; /* Rest are not deeper. */
last_alloca_header = hp; /* -> last valid storage. */
}
if (size == 0)
return NULL; /* No allocation required. */
/* Allocate combined header + user data storage. */
{
register pointer new = (*_mp_allocate_func) (sizeof (header) + size);
/* Address of header. */
((header *) new)->h.next = last_alloca_header;
((header *) new)->h.deep = depth;
last_alloca_header = (header *) new;
/* User storage begins just after header. */
return (pointer) ((char *) new + sizeof (header));
}
}
#ifdef CRAY
#ifdef DEBUG_I00AFUNC
#include <stdio.h>
#endif
#ifndef CRAY_STACK
#define CRAY_STACK
#ifndef CRAY2
/* Stack structures for CRAY-1, CRAY X-MP, and CRAY Y-MP */
struct stack_control_header
{
long shgrow:32; /* Number of times stack has grown. */
long shaseg:32; /* Size of increments to stack. */
long shhwm:32; /* High water mark of stack. */
long shsize:32; /* Current size of stack (all segments). */
};
/* The stack segment linkage control information occurs at
the high-address end of a stack segment. (The stack
grows from low addresses to high addresses.) The initial
part of the stack segment linkage control information is
0200 (octal) words. This provides for register storage
for the routine which overflows the stack. */
struct stack_segment_linkage
{
long ss[0200]; /* 0200 overflow words. */
long sssize:32; /* Number of words in this segment. */
long ssbase:32; /* Offset to stack base. */
long:32;
long sspseg:32; /* Offset to linkage control of previous
segment of stack. */
long:32;
long sstcpt:32; /* Pointer to task common address block. */
long sscsnm; /* Private control structure number for
microtasking. */
long ssusr1; /* Reserved for user. */
long ssusr2; /* Reserved for user. */
long sstpid; /* Process ID for pid based multi-tasking. */
long ssgvup; /* Pointer to multitasking thread giveup. */
long sscray[7]; /* Reserved for Cray Research. */
long ssa0;
long ssa1;
long ssa2;
long ssa3;
long ssa4;
long ssa5;
long ssa6;
long ssa7;
long sss0;
long sss1;
long sss2;
long sss3;
long sss4;
long sss5;
long sss6;
long sss7;
};
#else /* CRAY2 */
/* The following structure defines the vector of words
returned by the STKSTAT library routine. */
struct stk_stat
{
long now; /* Current total stack size. */
long maxc; /* Amount of contiguous space which would
be required to satisfy the maximum
stack demand to date. */
long high_water; /* Stack high-water mark. */
long overflows; /* Number of stack overflow ($STKOFEN) calls. */
long hits; /* Number of internal buffer hits. */
long extends; /* Number of block extensions. */
long stko_mallocs; /* Block allocations by $STKOFEN. */
long underflows; /* Number of stack underflow calls ($STKRETN). */
long stko_free; /* Number of deallocations by $STKRETN. */
long stkm_free; /* Number of deallocations by $STKMRET. */
long segments; /* Current number of stack segments. */
long maxs; /* Maximum number of stack segments so far. */
long pad_size; /* Stack pad size. */
long current_address; /* Current stack segment address. */
long current_size; /* Current stack segment size. This
number is actually corrupted by STKSTAT to
include the fifteen word trailer area. */
long initial_address; /* Address of initial segment. */
long initial_size; /* Size of initial segment. */
};
/* The following structure describes the data structure which trails
any stack segment. I think that the description in 'asdef' is
out of date. I only describe the parts that I am sure about. */
struct stk_trailer
{
long this_address; /* Address of this block. */
long this_size; /* Size of this block (does not include
this trailer). */
long unknown2;
long unknown3;
long link; /* Address of trailer block of previous
segment. */
long unknown5;
long unknown6;
long unknown7;
long unknown8;
long unknown9;
long unknown10;
long unknown11;
long unknown12;
long unknown13;
long unknown14;
};
#endif /* CRAY2 */
#endif /* not CRAY_STACK */
#ifdef CRAY2
/* Determine a "stack measure" for an arbitrary ADDRESS.
I doubt that "lint" will like this much. */
static long
i00afunc (long *address)
{
struct stk_stat status;
struct stk_trailer *trailer;
long *block, size;
long result = 0;
/* We want to iterate through all of the segments. The first
step is to get the stack status structure. We could do this
more quickly and more directly, perhaps, by referencing the
$LM00 common block, but I know that this works. */
STKSTAT (&status);
/* Set up the iteration. */
trailer = (struct stk_trailer *) (status.current_address
+ status.current_size
- 15);
/* There must be at least one stack segment. Therefore it is
a fatal error if "trailer" is null. */
if (trailer == 0)
abort ();
/* Discard segments that do not contain our argument address. */
while (trailer != 0)
{
block = (long *) trailer->this_address;
size = trailer->this_size;
if (block == 0 || size == 0)
abort ();
trailer = (struct stk_trailer *) trailer->link;
if ((block <= address) && (address < (block + size)))
break;
}
/* Set the result to the offset in this segment and add the sizes
of all predecessor segments. */
result = address - block;
if (trailer == 0)
{
return result;
}
do
{
if (trailer->this_size <= 0)
abort ();
result += trailer->this_size;
trailer = (struct stk_trailer *) trailer->link;
}
while (trailer != 0);
/* We are done. Note that if you present a bogus address (one
not in any segment), you will get a different number back, formed
from subtracting the address of the first block. This is probably
not what you want. */
return (result);
}
#else /* not CRAY2 */
/* Stack address function for a CRAY-1, CRAY X-MP, or CRAY Y-MP.
Determine the number of the cell within the stack,
given the address of the cell. The purpose of this
routine is to linearize, in some sense, stack addresses
for alloca. */
static long
i00afunc (long address)
{
long stkl = 0;
long size, pseg, this_segment, stack;
long result = 0;
struct stack_segment_linkage *ssptr;
/* Register B67 contains the address of the end of the
current stack segment. If you (as a subprogram) store
your registers on the stack and find that you are past
the contents of B67, you have overflowed the segment.
B67 also points to the stack segment linkage control
area, which is what we are really interested in. */
stkl = CRAY_STACKSEG_END ();
ssptr = (struct stack_segment_linkage *) stkl;
/* If one subtracts 'size' from the end of the segment,
one has the address of the first word of the segment.
If this is not the first segment, 'pseg' will be
nonzero. */
pseg = ssptr->sspseg;
size = ssptr->sssize;
this_segment = stkl - size;
/* It is possible that calling this routine itself caused
a stack overflow. Discard stack segments which do not
contain the target address. */
while (!(this_segment <= address && address <= stkl))
{
#ifdef DEBUG_I00AFUNC
fprintf (stderr, "%011o %011o %011o\n", this_segment, address, stkl);
#endif
if (pseg == 0)
break;
stkl = stkl - pseg;
ssptr = (struct stack_segment_linkage *) stkl;
size = ssptr->sssize;
pseg = ssptr->sspseg;
this_segment = stkl - size;
}
result = address - this_segment;
/* If you subtract pseg from the current end of the stack,
you get the address of the previous stack segment's end.
This seems a little convoluted to me, but I'll bet you save
a cycle somewhere. */
while (pseg != 0)
{
#ifdef DEBUG_I00AFUNC
fprintf (stderr, "%011o %011o\n", pseg, size);
#endif
stkl = stkl - pseg;
ssptr = (struct stack_segment_linkage *) stkl;
size = ssptr->sssize;
pseg = ssptr->sspseg;
result += size;
}
return (result);
}
#endif /* not CRAY2 */
#endif /* CRAY */
#endif /* no alloca */

141
gnu/lib/libgmp/cre-conv-tab.c Normal file
View file

@ -0,0 +1,141 @@
/* cre-conv-tab.c -- Create conversion table in a wordsize-dependent way.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
extern double floor ();
extern double log ();
static unsigned long int
upow (b, e)
unsigned long int b;
unsigned int e;
{
unsigned long int y = 1;
while (e != 0)
{
while ((e & 1) == 0)
{
b = b * b;
e >>= 1;
}
y = y * b;
e -= 1;
}
return y;
}
unsigned int
ulog2 (x)
unsigned long int x;
{
unsigned int i;
for (i = 0; x != 0; i++)
x >>= 1;
return i;
}
main ()
{
int i;
unsigned long idig;
unsigned long big_base, big_base_inverted;
double fdig;
int dummy;
int normalization_steps;
unsigned long int max_uli;
int bits_uli;
max_uli = 1;
for (i = 1; ; i++)
{
if ((max_uli << 1) == 0)
break;
max_uli <<= 1;
}
bits_uli = i;
puts ("/* __mp_bases -- Structure for conversion between internal binary");
puts (" format and strings in base 2..36. The fields are explained in");
puts (" gmp-impl.h.");
puts ("");
puts (" ***** THIS FILE WAS CREATED BY A PROGRAM. DON'T EDIT IT! *****");
puts ("");
puts ("Copyright (C) 1991 Free Software Foundation, Inc.");
puts ("");
puts ("This file is part of the GNU MP Library.");
puts ("");
puts ("The GNU MP Library is free software; you can redistribute it and/or");
puts ("modify it under the terms of the GNU General Public License as");
puts ("published by the Free Software Foundation; either version 2, or");
puts ("(at your option) any later version.");
puts ("");
puts ("The GNU MP Library is distributed in the hope that it will be");
puts ("useful, but WITHOUT ANY WARRANTY; without even the implied warranty");
puts ("of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the");
puts ("GNU General Public License for more details.");
puts ("");
puts ("You should have received a copy of the GNU General Public License");
puts ("along with the GNU MP Library; see the file COPYING. If not, write");
puts ("to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,");
puts ("USA. */");
puts ("");
puts ("#include \"gmp.h\"");
puts ("#include \"gmp-impl.h\"");
puts ("");
puts ("const struct bases __mp_bases[37] =\n{");
puts (" /* 0 */ {0, 0, 0, 0.0},");
puts (" /* 1 */ {0, 0, 0, 0.0},");
for (i = 2; i <= 36; i++)
{
/* The weird expression here is because many /bin/cc compilers
generate incorrect code for conversions from large unsigned
integers to double. */
fdig = log(2.0)/log((double) i);
idig = floor(bits_uli * fdig);
if ((i & (i - 1)) == 0)
{
big_base = ulog2 (i) - 1;
big_base_inverted = 0;
}
else
{
big_base = upow (i, idig);
for (normalization_steps = 0;
(long int) (big_base << normalization_steps) >= 0;
normalization_steps++)
;
udiv_qrnnd (big_base_inverted, dummy,
-(big_base << normalization_steps), 0,
big_base << normalization_steps);
}
printf (" /* %2u */ {%lu, 0x%lX, 0x%lX, %.8f},\n",
i, idig, big_base, big_base_inverted, fdig);
}
puts ("};");
exit (0);
}

118
gnu/lib/libgmp/cre-mparam.c Normal file
View file

@ -0,0 +1,118 @@
/* cre-mparam.c -- Create machine-depedent parameter file.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
unsigned int
ulog2 (x)
unsigned long int x;
{
unsigned int i;
for (i = 0; x != 0; i++)
x >>= 1;
return i;
}
main ()
{
int i;
unsigned long int max_uli;
int bits_uli;
unsigned long int max_ui;
int bits_ui;
unsigned long int max_usi;
int bits_usi;
unsigned long int max_uc;
int bits_uc;
max_uli = 1;
for (i = 0; ; i++)
{
if (max_uli == 0)
break;
max_uli <<= 1;
}
bits_uli = i;
max_ui = 1;
for (i = 0; ; i++)
{
if ((unsigned int) max_ui == 0)
break;
max_ui <<= 1;
}
bits_ui = i;
max_usi = 1;
for (i = 0; ; i++)
{
if ((unsigned short int) max_usi == 0)
break;
max_usi <<= 1;
}
bits_usi = i;
max_uc = 1;
for (i = 0; ; i++)
{
if ((unsigned char) max_uc == 0)
break;
max_uc <<= 1;
}
bits_uc = i;
puts ("/* gmp-mparam.h -- Compiler/machine parameter header file.");
puts ("");
puts (" ***** THIS FILE WAS CREATED BY A PROGRAM. DON'T EDIT IT! *****");
puts ("");
puts ("Copyright (C) 1991 Free Software Foundation, Inc.");
puts ("");
puts ("This file is part of the GNU MP Library.");
puts ("");
puts ("The GNU MP Library is free software; you can redistribute it and/or");
puts ("modify it under the terms of the GNU General Public License as");
puts ("published by the Free Software Foundation; either version 2, or");
puts ("(at your option) any later version.");
puts ("");
puts ("The GNU MP Library is distributed in the hope that it will be");
puts ("useful, but WITHOUT ANY WARRANTY; without even the implied warranty");
puts ("of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the");
puts ("GNU General Public License for more details.");
puts ("");
puts ("You should have received a copy of the GNU General Public License");
puts ("along with the GNU MP Library; see the file COPYING. If not, write");
puts ("to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,");
puts ("USA. */");
puts ("");
printf ("#define BITS_PER_MP_LIMB %d\n", bits_uli);
printf ("#define BYTES_PER_MP_LIMB %d\n", sizeof(mp_limb));
printf ("#define BITS_PER_LONGINT %d\n", bits_uli);
printf ("#define BITS_PER_INT %d\n", bits_ui);
printf ("#define BITS_PER_SHORTINT %d\n", bits_usi);
printf ("#define BITS_PER_CHAR %d\n", bits_uc);
exit (0);
}

42
gnu/lib/libgmp/cre-stddefh.c Normal file
View file

@ -0,0 +1,42 @@
/* cre-stddefh.c -- Check the size of a pointer and output an
appropriate size_t declaration.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
main (argc, argv)
int argc;
char **argv;
{
if (sizeof (int *) == sizeof (unsigned long int))
puts ("typedef unsigned long int size_t;");
else
if (sizeof (int *) == sizeof (unsigned int))
puts ("typedef unsigned int size_t;");
else
{
fprintf (stderr,
"%s: Can't find a reasonable definition for \"size_t\".\n",
argv[0]);
exit (1);
}
exit (0);
}

126
gnu/lib/libgmp/gmp-impl.h Normal file
View file

@ -0,0 +1,126 @@
/* Include file for internal GNU MP types and definitions.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#if defined (__GNUC__) || defined (__sparc__) || defined (sparc)
#define alloca __builtin_alloca
#endif
#ifndef NULL
#define NULL 0L
#endif
#if defined (__GNUC__)
volatile void abort (void);
#else
#define inline /* Empty */
void *alloca();
#endif
#define ABS(x) (x >= 0 ? x : -x)
#include "gmp-mparam.h"
#ifdef __STDC__
void *malloc (size_t);
void *realloc (void *, size_t);
void free (void *);
extern void * (*_mp_allocate_func) (size_t);
extern void * (*_mp_reallocate_func) (void *, size_t, size_t);
extern void (*_mp_free_func) (void *, size_t);
void *_mp_default_allocate (size_t);
void *_mp_default_reallocate (void *, size_t, size_t);
void _mp_default_free (void *, size_t);
char *_mpz_get_str (char *, int, const MP_INT *);
int _mpz_set_str (MP_INT *, const char *, int);
void _mpz_impl_sqrt (MP_INT *, MP_INT *, const MP_INT *);
#else
#define const /* Empty */
#define signed /* Empty */
void *malloc ();
void *realloc ();
void free ();
extern void * (*_mp_allocate_func) ();
extern void * (*_mp_reallocate_func) ();
extern void (*_mp_free_func) ();
void *_mp_default_allocate ();
void *_mp_default_reallocate ();
void _mp_default_free ();
char *_mpz_get_str ();
int _mpz_set_str ();
void _mpz_impl_sqrt ();
#endif
/* Copy NLIMBS *limbs* from SRC to DST. */
#define MPN_COPY(DST, SRC, NLIMBS) \
do { \
mp_size i; \
for (i = 0; i < (NLIMBS); i++) \
(DST)[i] = (SRC)[i]; \
} while (0)
/* Zero NLIMBS *limbs* AT DST. */
#define MPN_ZERO(DST, NLIMBS) \
do { \
mp_size i; \
for (i = 0; i < (NLIMBS); i++) \
(DST)[i] = 0; \
} while (0)
/* Initialize the MP_INT X with space for NLIMBS limbs.
X should be a temporary variable, and it will be automatically
cleared out when the running function returns. */
#define MPZ_TMP_INIT(X, NLIMBS) \
do { \
(X)->alloc = (NLIMBS); \
(X)->d = (mp_ptr) alloca ((NLIMBS) * BYTES_PER_MP_LIMB); \
} while (0)
/* Structure for conversion between internal binary format and
strings in base 2..36. */
struct bases
{
/* Number of digits in the conversion base that always fits in
an mp_limb. For example, for base 10 this is 10, since
2**32 = 4294967296 has ten digits. */
int chars_per_limb;
/* big_base is conversion_base**chars_per_limb, i.e. the biggest
number that fits a word, built by factors of conversion_base.
Exception: For 2, 4, 8, etc, big_base is log2(base), i.e. the
number of bits used to represent each digit in the base. */
mp_limb big_base;
/* big_base_inverted is a BITS_PER_MP_LIMB bit approximation to
1/big_base, represented as a fixed-point number. Instead of
dividing by big_base an application can choose to multiply
by big_base_inverted. */
mp_limb big_base_inverted;
/* log(2)/log(conversion_base) */
float chars_per_bit_exactly;
};
extern const struct bases __mp_bases[37];

302
gnu/lib/libgmp/gmp.h Normal file
View file

@ -0,0 +1,302 @@
/* gmp.h -- Definitions for GNU multiple precision functions.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#ifndef __GMP_H__
#define __GMP_H__
#define __GNU_MP__
#ifndef __MP_H__
#define __need_size_t
#include <stddef.h>
#endif
#ifndef MINT
#ifndef __MP_SMALL__
typedef struct
{
long int alloc; /* Number of *limbs* allocated and pointed
to by the D field. */
long int size; /* abs(SIZE) is the number of limbs
the last field points to. If SIZE
is negative this is a negative
number. */
unsigned long int *d; /* Pointer to the limbs. */
} __MP_INT;
#else
typedef struct
{
short int alloc; /* Number of *limbs* allocated and pointed
to by the D field. */
short int size; /* abs(SIZE) is the number of limbs
the last field points to. If SIZE
is negative this is a negative
number. */
unsigned long int *d; /* Pointer to the limbs. */
} __MP_INT;
#endif
#endif
#define MP_INT __MP_INT
typedef unsigned long int mp_limb;
typedef long int mp_limb_signed;
typedef mp_limb * mp_ptr;
#ifdef __STDC__
typedef const mp_limb * mp_srcptr;
#else
typedef mp_limb * mp_srcptr;
#endif
typedef long int mp_size;
/* Structure for rational numbers. Zero is represented as 0/any, i.e.
the denominator is ignored. Negative numbers have the sign in
the numerator. */
typedef struct
{
MP_INT num;
MP_INT den;
#if 0
long int num_alloc; /* Number of limbs allocated
for the numerator. */
long int num_size; /* The absolute value of this field is the
length of the numerator; the sign is the
sign of the entire rational number. */
mp_ptr num; /* Pointer to the numerator limbs. */
long int den_alloc; /* Number of limbs allocated
for the denominator. */
long int den_size; /* Length of the denominator. (This field
should always be positive.) */
mp_ptr den; /* Pointer to the denominator limbs. */
#endif
} MP_RAT;
#ifdef __STDC__
void mp_set_memory_functions (void *(*) (size_t),
void *(*) (void *, size_t, size_t),
void (*) (void *, size_t));
/**************** Integer (i.e. Z) routines. ****************/
void mpz_init (MP_INT *);
void mpz_set (MP_INT *, const MP_INT *);
void mpz_set_ui (MP_INT *, unsigned long int);
void mpz_set_si (MP_INT *, signed long int);
int mpz_set_str (MP_INT *, const char *, int);
void mpz_init_set (MP_INT *, const MP_INT *);
void mpz_init_set_ui (MP_INT *, unsigned long int);
void mpz_init_set_si (MP_INT *, signed long int);
int mpz_init_set_str (MP_INT *, const char *, int);
unsigned long int mpz_get_ui (const MP_INT *);
signed long int mpz_get_si (const MP_INT *);
char * mpz_get_str (char *, int, const MP_INT *);
void mpz_clear (MP_INT *);
void * _mpz_realloc (MP_INT *, mp_size);
void mpz_add (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_add_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_sub (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_sub_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_mul (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_mul_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_div (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_div_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_mod (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_mod_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_divmod (MP_INT *, MP_INT *, const MP_INT *, const MP_INT *);
void mpz_divmod_ui (MP_INT *, MP_INT *, const MP_INT *, unsigned long int);
void mpz_mdiv (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_mdiv_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_mmod (MP_INT *, const MP_INT *, const MP_INT *);
unsigned long int mpz_mmod_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_mdivmod (MP_INT *, MP_INT *, const MP_INT *, const MP_INT *);
unsigned long int mpz_mdivmod_ui (MP_INT *, MP_INT *, const MP_INT *,
unsigned long int);
void mpz_sqrt (MP_INT *, const MP_INT *);
void mpz_sqrtrem (MP_INT *, MP_INT *, const MP_INT *);
int mpz_perfect_square_p (const MP_INT *);
int mpz_probab_prime_p (const MP_INT *, int);
void mpz_powm (MP_INT *, const MP_INT *, const MP_INT *, const MP_INT *);
void mpz_powm_ui (MP_INT *, const MP_INT *, unsigned long int, const MP_INT *);
void mpz_pow_ui (MP_INT *, const MP_INT *, unsigned long int);
void mpz_fac_ui (MP_INT *, unsigned long int);
void mpz_gcd (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_gcdext (MP_INT *, MP_INT *, MP_INT *, const MP_INT *, const MP_INT *);
void mpz_neg (MP_INT *, const MP_INT *);
void mpz_com (MP_INT *, const MP_INT *);
void mpz_abs (MP_INT *, const MP_INT *);
int mpz_cmp (const MP_INT *, const MP_INT *);
int mpz_cmp_ui (const MP_INT *, unsigned long int);
int mpz_cmp_si (const MP_INT *, signed long int);
void mpz_mul_2exp (MP_INT *, const MP_INT *, unsigned long int);
void mpz_div_2exp (MP_INT *, const MP_INT *, unsigned long int);
void mpz_mod_2exp (MP_INT *, const MP_INT *, unsigned long int);
void mpz_and (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_ior (MP_INT *, const MP_INT *, const MP_INT *);
void mpz_xor (MP_INT *, const MP_INT *, const MP_INT *);
#if defined (FILE) || defined (_STDIO_H) || defined (__STDIO_H__)
void mpz_inp_raw (MP_INT *, FILE *);
void mpz_inp_str (MP_INT *, FILE *, int);
void mpz_out_raw (FILE *, const MP_INT *);
void mpz_out_str (FILE *, int, const MP_INT *);
#endif
void mpz_array_init (MP_INT [], size_t, mp_size);
void mpz_random (MP_INT *, mp_size);
void mpz_random2 (MP_INT *, mp_size);
size_t mpz_size (const MP_INT *);
size_t mpz_sizeinbase (const MP_INT *, int);
/**************** Rational (i.e. Q) routines. ****************/
void mpq_init (MP_RAT *);
void mpq_clear (MP_RAT *);
void mpq_set (MP_RAT *, const MP_RAT *);
void mpq_set_ui (MP_RAT *, unsigned long int, unsigned long int);
void mpq_set_si (MP_RAT *, signed long int, unsigned long int);
void mpq_add (MP_RAT *, const MP_RAT *, const MP_RAT *);
void mpq_sub (MP_RAT *, const MP_RAT *, const MP_RAT *);
void mpq_mul (MP_RAT *, const MP_RAT *, const MP_RAT *);
void mpq_div (MP_RAT *, const MP_RAT *, const MP_RAT *);
void mpq_neg (MP_RAT *, const MP_RAT *);
int mpq_cmp (const MP_RAT *, const MP_RAT *);
void mpq_inv (MP_RAT *, const MP_RAT *);
void mpq_set_num (MP_RAT *, const MP_INT *);
void mpq_set_den (MP_RAT *, const MP_INT *);
void mpq_get_num (MP_INT *, const MP_RAT *);
void mpq_get_den (MP_INT *, const MP_RAT *);
/************ Low level positive-integer (i.e. N) routines. ************/
mp_limb mpn_add (mp_ptr, mp_srcptr, mp_size, mp_srcptr, mp_size);
mp_size mpn_sub (mp_ptr, mp_srcptr, mp_size, mp_srcptr, mp_size);
mp_size mpn_mul (mp_ptr, mp_srcptr, mp_size, mp_srcptr, mp_size);
mp_size mpn_div (mp_ptr, mp_ptr, mp_size, mp_srcptr, mp_size);
mp_limb mpn_divmod_1 (mp_ptr, mp_srcptr, mp_size, mp_limb);
mp_limb mpn_mod_1 (mp_srcptr, mp_size, mp_limb);
mp_limb mpn_lshift (mp_ptr, mp_srcptr, mp_size, unsigned int);
mp_size mpn_rshift (mp_ptr, mp_srcptr, mp_size, unsigned int);
mp_size mpn_rshiftci (mp_ptr, mp_srcptr, mp_size, unsigned int, mp_limb);
mp_size mpn_sqrt (mp_ptr, mp_ptr, mp_srcptr, mp_size);
int mpn_cmp (mp_srcptr, mp_srcptr, mp_size);
#else /* ! __STDC__ */
void mp_set_memory_functions ();
/**************** Integer (i.e. Z) routines. ****************/
void mpz_init ();
void mpz_set ();
void mpz_set_ui ();
void mpz_set_si ();
int mpz_set_str ();
void mpz_init_set ();
void mpz_init_set_ui ();
void mpz_init_set_si ();
int mpz_init_set_str ();
unsigned long int mpz_get_ui ();
long int mpz_get_si ();
char * mpz_get_str ();
void mpz_clear ();
void * _mpz_realloc ();
void mpz_add ();
void mpz_add_ui ();
void mpz_sub ();
void mpz_sub_ui ();
void mpz_mul ();
void mpz_mul_ui ();
void mpz_div ();
void mpz_div_ui ();
void mpz_mod ();
void mpz_mod_ui ();
void mpz_divmod ();
void mpz_divmod_ui ();
void mpz_mdiv ();
void mpz_mdiv_ui ();
void mpz_mmod ();
unsigned long int mpz_mmod_ui ();
void mpz_mdivmod ();
unsigned long int mpz_mdivmod_ui ();
void mpz_sqrt ();
void mpz_sqrtrem ();
int mpz_perfect_square_p ();
int mpz_probab_prime_p ();
void mpz_powm ();
void mpz_powm_ui ();
void mpz_pow_ui ();
void mpz_fac_ui ();
void mpz_gcd ();
void mpz_gcdext ();
void mpz_neg ();
void mpz_com ();
void mpz_abs ();
int mpz_cmp ();
int mpz_cmp_ui ();
int mpz_cmp_si ();
void mpz_mul_2exp ();
void mpz_div_2exp ();
void mpz_mod_2exp ();
void mpz_and ();
void mpz_ior ();
void mpz_xor ();
void mpz_inp_raw ();
void mpz_inp_str ();
void mpz_out_raw ();
void mpz_out_str ();
void mpz_array_init ();
void mpz_random ();
void mpz_random2 ();
size_t mpz_size ();
size_t mpz_sizeinbase ();
/**************** Rational (i.e. Q) routines. ****************/
void mpq_init ();
void mpq_clear ();
void mpq_set ();
void mpq_set_ui ();
void mpq_set_si ();
void mpq_add ();
void mpq_sub ();
void mpq_mul ();
void mpq_div ();
void mpq_neg ();
int mpq_cmp ();
void mpq_inv ();
void mpq_set_num ();
void mpq_set_den ();
void mpq_get_num ();
void mpq_get_den ();
/************ Low level positive-integer (i.e. N) routines. ************/
mp_limb mpn_add ();
mp_size mpn_sub ();
mp_size mpn_mul ();
mp_size mpn_div ();
mp_limb mpn_lshift ();
mp_size mpn_rshift ();
mp_size mpn_rshiftci ();
int mpn_cmp ();
#endif /* __STDC__ */
#endif /* __GMP_H__ */

1291
gnu/lib/libgmp/gmp.texi Normal file
View file

File diff suppressed because it is too large Load diff

53
gnu/lib/libgmp/itom.c Normal file
View file

@ -0,0 +1,53 @@
/* itom -- BSD compatible allocate and initiate a MINT.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
MINT *
#ifdef __STDC__
itom (signed short int n)
#else
itom (n)
short int n;
#endif
{
MINT *x;
mp_ptr xp;
x = (MINT *) (*_mp_allocate_func) (sizeof (MINT));
x->alloc = 1;
x->d = xp = (mp_ptr) (*_mp_allocate_func) (x->alloc * BYTES_PER_MP_LIMB);
if (n > 0)
{
x->size = 1;
xp[0] = n;
}
else if (n < 0)
{
x->size = -1;
xp[0] = -n;
}
else
x->size = 0;
return x;
}

1027
gnu/lib/libgmp/longlong.h Normal file
View file

File diff suppressed because it is too large Load diff

38
gnu/lib/libgmp/mdiv.c Normal file
View file

@ -0,0 +1,38 @@
/* mdiv -- BSD compatible divide producing both remainder and quotient.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mdiv (const MINT *num, const MINT *den, MINT *quot, MINT *rem)
#else
mdiv (num, den, quot, rem)
const MINT *num;
const MINT *den;
MINT *quot;
MINT *rem;
#endif
#define COMPUTE_QUOTIENT
#include "mpz_dmincl.c"

96
gnu/lib/libgmp/memory.c Normal file
View file

@ -0,0 +1,96 @@
/* Memory allocation routines.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include "gmp.h"
#include "gmp-impl.h"
#ifdef __NeXT__
#define static
#endif
#ifdef __STDC__
void * (*_mp_allocate_func) (size_t) = _mp_default_allocate;
void * (*_mp_reallocate_func) (void *, size_t, size_t)
= _mp_default_reallocate;
void (*_mp_free_func) (void *, size_t) = _mp_default_free;
#else
void * (*_mp_allocate_func) () = _mp_default_allocate;
void * (*_mp_reallocate_func) () = _mp_default_reallocate;
void (*_mp_free_func) () = _mp_default_free;
#endif
/* Default allocation functions. In case of failure to allocate/reallocate
an error message is written to stderr and the program aborts. */
void *
#ifdef __STDC__
_mp_default_allocate (size_t size)
#else
_mp_default_allocate (size)
size_t size;
#endif
{
void *ret;
ret = malloc (size);
if (ret == 0)
{
perror ("cannot allocate in libmp");
abort ();
}
return ret;
}
void *
#ifdef __STDC__
_mp_default_reallocate (void *oldptr, size_t old_size, size_t new_size)
#else
_mp_default_reallocate (oldptr, old_size, new_size)
void *oldptr;
size_t old_size;
size_t new_size;
#endif
{
void *ret;
ret = realloc (oldptr, new_size);
if (ret == 0)
{
perror ("cannot allocate in libmp");
abort ();
}
return ret;
}
void
#ifdef __STDC__
_mp_default_free (void *blk_ptr, size_t blk_size)
#else
_mp_default_free (blk_ptr, blk_size)
void *blk_ptr;
size_t blk_size;
#endif
{
free (blk_ptr);
}

35
gnu/lib/libgmp/mfree.c Normal file
View file

@ -0,0 +1,35 @@
/* mfree -- BSD compatible mfree.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mfree (MINT *m)
#else
mfree (m)
MINT *m;
#endif
{
(*_mp_free_func) (m->d, m->alloc * BYTES_PER_MP_LIMB);
(*_mp_free_func) (m, sizeof (MINT));
}

64
gnu/lib/libgmp/min.c Normal file
View file

@ -0,0 +1,64 @@
/* min(MINT) -- Do decimal input from standard input and store result in
MINT.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include <ctype.h>
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
min (MINT *x)
#else
min (x)
MINT *x;
#endif
{
char *str;
size_t str_size;
size_t i;
int c;
str_size = 100;
str = (char *) (*_mp_allocate_func) (str_size);
for (i = 0; ; i++)
{
if (i >= str_size)
{
size_t old_str_size = str_size;
str_size = str_size * 3 / 2;
str = (char *) (*_mp_reallocate_func) (str, old_str_size, str_size);
}
c = getc (stdin);
if (!(isdigit(c) || c == ' ' || c == '\t'))
break;
str[i] = c;
}
ungetc (c, stdin);
str[i] = 0;
_mpz_set_str (x, str, 10);
(*_mp_free_func) (str, str_size);
}

42
gnu/lib/libgmp/mout.c Normal file
View file

@ -0,0 +1,42 @@
/* mout(MINT) -- Do decimal output of MINT to standard output.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mout (const MINT *x)
#else
mout (x)
const MINT *x;
#endif
{
char *str;
size_t str_size;
str_size = ((size_t) (ABS (x->size) * BITS_PER_MP_LIMB
* __mp_bases[10].chars_per_bit_exactly)) + 3;
str = (char *) alloca (str_size);
_mpz_get_str (str, 10, x);
puts (str);
alloca (0);
}

45
gnu/lib/libgmp/move.c Normal file
View file

@ -0,0 +1,45 @@
/* move -- BSD compatible assignment.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "mp.h"
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
move (const MINT *u, MINT *w)
#else
move (u, w)
const MINT *u;
MINT *w;
#endif
{
mp_size usize;
mp_size abs_usize;
usize = u->size;
abs_usize = ABS (usize);
if (w->alloc < abs_usize)
_mpz_realloc (w, abs_usize);
w->size = usize;
MPN_COPY (w->d, u->d, abs_usize);
}

103
gnu/lib/libgmp/mp.h Normal file
View file

@ -0,0 +1,103 @@
/* mp.h -- Definitions for Berkeley compatible multiple precision functions.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#ifndef __MP_H__
#define __MP_H__
#define __GNU_MP__
#ifndef __GMP_H__
#define __need_size_t
#include <stddef.h>
#endif
#ifndef MP_INT
#ifndef __MP_SMALL__
typedef struct
{
long int alloc; /* Number of *limbs* allocated and pointed
to by the D field. */
long int size; /* abs(SIZE) is the number of limbs
the last field points to. If SIZE
is negative this is a negative
number. */
unsigned long int *d; /* Pointer to the limbs. */
} __MP_INT;
#else
typedef struct
{
short int alloc; /* Number of *limbs* allocated and pointed
to by the D field. */
short int size; /* abs(SIZE) is the number of limbs
the last field points to. If SIZE
is negative this is a negative
number. */
unsigned long int *d; /* Pointer to the limbs. */
} __MP_INT;
#endif
#endif
#define MINT __MP_INT
#ifdef __STDC__
void mp_set_memory_functions (void *(*) (size_t),
void *(*) (void *, size_t, size_t),
void (*) (void *, size_t));
MINT *itom (signed short int);
MINT *xtom (const char *);
void move (const MINT *, MINT *);
void madd (const MINT *, const MINT *, MINT *);
void msub (const MINT *, const MINT *, MINT *);
void mult (const MINT *, const MINT *, MINT *);
void mdiv (const MINT *, const MINT *, MINT *, MINT *);
void sdiv (const MINT *, signed short int, MINT *, signed short int *);
void msqrt (const MINT *, MINT *, MINT *);
void pow (const MINT *, const MINT *, const MINT *, MINT *);
void rpow (const MINT *, signed short int, MINT *);
void gcd (const MINT *, const MINT *, MINT *);
int mcmp (const MINT *, const MINT *);
void min (MINT *);
void mout (const MINT *);
char *mtox (const MINT *);
void mfree (MINT *);
#else
void mp_set_memory_functions ();
MINT *itom ();
MINT *xtom ();
void move ();
void madd ();
void msub ();
void mult ();
void mdiv ();
void sdiv ();
void msqrt ();
void pow ();
void rpow ();
void gcd ();
int mcmp ();
void min ();
void mout ();
char *mtox ();
void mfree ();
#endif
#endif /* __MP_H__ */

36
gnu/lib/libgmp/mp_clz_tab.c Normal file
View file

@ -0,0 +1,36 @@
/* __clz_tab -- support for longlong.h
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* BOTCH: This ought to be made machine-independent. */
const unsigned char __clz_tab[] =
{
0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,
};

47
gnu/lib/libgmp/mp_set_fns.c Normal file
View file

@ -0,0 +1,47 @@
/* mp_set_memory_functions -- Set the allocate, reallocate, and free functions
for use by the mp package.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mp_set_memory_functions (void *(*alloc_func) (size_t),
void *(*realloc_func) (void *, size_t, size_t),
void (*free_func) (void *, size_t))
#else
mp_set_memory_functions (alloc_func, realloc_func, free_func)
void *(*alloc_func) ();
void *(*realloc_func) ();
void (*free_func) ();
#endif
{
if (alloc_func == 0)
alloc_func = _mp_default_allocate;
if (realloc_func == 0)
realloc_func = _mp_default_reallocate;
if (free_func == 0)
free_func = _mp_default_free;
_mp_allocate_func = alloc_func;
_mp_reallocate_func = realloc_func;
_mp_free_func = free_func;
}

141
gnu/lib/libgmp/mpn_add.c Normal file
View file

@ -0,0 +1,141 @@
/* mpn_add -- Add two low-level integers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Add ADD1_PTR/ADD1_SIZE and ADD2_PTR/ADD2_SIZE and store the first
ADD1_SIZE words of the result at SUM_PTR.
Return 1 if carry out was generated, return 0 otherwise.
Argument constraint: ADD1_SIZE >= ADD2_SIZE.
The size of SUM can be calculated as ADD1_SIZE + the return value. */
mp_limb
#ifdef __STDC__
mpn_add (mp_ptr sum_ptr,
mp_srcptr add1_ptr, mp_size add1_size,
mp_srcptr add2_ptr, mp_size add2_size)
#else
mpn_add (sum_ptr, add1_ptr, add1_size, add2_ptr, add2_size)
mp_ptr sum_ptr;
mp_srcptr add1_ptr;
mp_size add1_size;
mp_srcptr add2_ptr;
mp_size add2_size;
#endif
{
mp_limb a1, a2, sum;
mp_size j;
/* The loop counter and index J goes from some negative value to zero.
This way the loops become faster. Need to offset the base pointers
to take care of the negative indices. */
j = -add2_size;
if (j == 0)
goto add2_finished;
add1_ptr -= j;
add2_ptr -= j;
sum_ptr -= j;
/* There are two do-loops, marked NON-CARRY LOOP and CARRY LOOP that
jump between each other. The first loop is for when the previous
addition didn't produce a carry-out; the second is for the
complementary case. */
/* NON-CARRY LOOP */
do
{
a1 = add1_ptr[j];
a2 = add2_ptr[j];
sum = a1 + a2;
sum_ptr[j] = sum;
if (sum < a2)
goto cy_loop;
ncy_loop:
j++;
}
while (j < 0);
/* We have exhausted ADD2. Just copy ADD1 to SUM, and return
0 as an indication of no carry-out. */
add2_finished:
/* Immediate return if the copy would be a no-op. */
if (sum_ptr == add1_ptr)
return 0;
j = add2_size - add1_size;
add1_ptr -= j;
sum_ptr -= j;
while (j < 0)
{
sum_ptr[j] = add1_ptr[j];
j++;
}
return 0;
/* CARRY LOOP */
do
{
a1 = add1_ptr[j];
a2 = add2_ptr[j];
sum = a1 + a2 + 1;
sum_ptr[j] = sum;
if (sum > a2)
goto ncy_loop;
cy_loop:
j++;
}
while (j < 0);
j = add2_size - add1_size;
add1_ptr -= j;
sum_ptr -= j;
while (j < 0)
{
a1 = add1_ptr[j];
sum = a1 + 1;
sum_ptr[j] = sum;
if (sum > 0)
goto copy_add1;
j++;
}
return 1;
copy_add1:
if (sum_ptr == add1_ptr)
return 0;
j++;
while (j < 0)
{
sum_ptr[j] = add1_ptr[j];
j++;
}
return 0;
}

52
gnu/lib/libgmp/mpn_cmp.c Normal file
View file

@ -0,0 +1,52 @@
/* mpn_cmp -- Compare two low-level natural-number integers.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Compare OP1_PTR/OP1_SIZE with OP2_PTR/OP2_SIZE.
There are no restrictions on the relative sizes of
the two arguments.
Return 1 if OP1 > OP2, 0 if they are equal, and -1 if OP1 < OP2. */
int
#ifdef __STDC__
mpn_cmp (mp_srcptr op1_ptr, mp_srcptr op2_ptr, mp_size size)
#else
mpn_cmp (op1_ptr, op2_ptr, size)
mp_srcptr op1_ptr;
mp_srcptr op2_ptr;
mp_size size;
#endif
{
mp_size i;
mp_limb op1_word, op2_word;
for (i = size - 1; i >= 0; i--)
{
op1_word = op1_ptr[i];
op2_word = op2_ptr[i];
if (op1_word != op2_word)
goto diff;
}
return 0;
diff:
return (op1_word > op2_word) ? 1 : -1;
}

321
gnu/lib/libgmp/mpn_div.c Normal file
View file

@ -0,0 +1,321 @@
/* mpn_div -- Divide natural numbers, producing both remainder and
quotient.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
/* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
the quotient at QUOT_PTR and the remainder at NUM_PTR.
Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE)
or (NSIZE - DSIZE + 1).
Argument constraints:
1. The most significant bit of d must be set.
2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage
area must be distinct from either input operands.
The exact sizes of the quotient and remainder must be determined
by the caller, in spite of the return value. The return value just
informs the caller about if the highest digit is written or not, and
it may very well be 0. */
/* THIS WILL BE IMPROVED SOON. MORE COMMENTS AND FASTER CODE. */
mp_size
#ifdef __STDC__
mpn_div (mp_ptr quot_ptr,
mp_ptr num_ptr, mp_size num_size,
mp_srcptr den_ptr, mp_size den_size)
#else
mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size)
mp_ptr quot_ptr;
mp_ptr num_ptr;
mp_size num_size;
mp_srcptr den_ptr;
mp_size den_size;
#endif
{
mp_size q_is_long = 0;
switch (den_size)
{
case 0:
/* We are asked to divide by zero, so go ahead and do it!
(To make the compiler not remove this statement, assign NUM_SIZE
and fall through.) */
num_size = 1 / den_size;
case 1:
{
mp_size i;
mp_limb n1, n0;
mp_limb d;
d = den_ptr[0];
i = num_size - 1;
n1 = num_ptr[i];
i--;
if (n1 >= d)
{
q_is_long = 1;
n1 = 0;
i++;
}
for (; i >= 0; i--)
{
n0 = num_ptr[i];
udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
}
num_ptr[0] = n1;
}
break;
case 2:
{
mp_size i;
mp_limb n0, n1, n2;
mp_limb d0, d1;
num_ptr += num_size - 2;
d0 = den_ptr[1];
d1 = den_ptr[0];
n0 = num_ptr[1];
n1 = num_ptr[0];
if (n0 >= d0)
{
q_is_long = 1;
n1 = n0;
n0 = 0;
num_ptr++;
num_size++;
}
for (i = num_size - den_size - 1; i >= 0; i--)
{
mp_limb q;
mp_limb r;
num_ptr--;
if (n0 == d0)
{
/* Q should be either 111..111 or 111..110. Need special
treatment of this rare case as normal division would
give overflow. */
q = ~0;
r = n1 + d0;
if (r < d0) /* Carry in the addition? */
{
n2 = num_ptr[0];
add_ssaaaa (n0, n1, r - d1, n2, 0, d1);
quot_ptr[i] = q;
continue;
}
n0 = d1 - (d1 != 0);
n1 = -d1;
}
else
{
udiv_qrnnd (q, r, n0, n1, d0);
umul_ppmm (n0, n1, d1, q);
}
n2 = num_ptr[0];
q_test:
if (n0 > r || (n0 == r && n1 > n2))
{
/* The estimated Q was too large. */
q--;
sub_ddmmss (n0, n1, n0, n1, 0, d1);
r += d0;
if (r >= d0) /* If not carry, test q again. */
goto q_test;
}
quot_ptr[i] = q;
sub_ddmmss (n0, n1, r, n2, n0, n1);
}
num_ptr[1] = n0;
num_ptr[0] = n1;
}
break;
default:
{
mp_size i;
mp_limb d0 = den_ptr[den_size - 1];
mp_limb d1 = den_ptr[den_size - 2];
mp_limb n0 = num_ptr[num_size - 1];
int ugly_hack_flag = 0;
if (n0 >= d0)
{
/* There's a problem with this case, which shows up later in the
code. q becomes 1 (and sometimes 0) the first time when
we've been here, and w_cy == 0 after the main do-loops below.
But c = num_ptr[j] reads rubbish outside the num_ptr vector!
Maybe I can solve this cleanly when I fix the early-end
optimization here in the default case. For now, I change the
add_back entering condition, to kludge. Leaving the stray
memref behind!
HACK: Added ugly_hack_flag to make it work. */
q_is_long = 1;
n0 = 0;
num_size++;
ugly_hack_flag = 1;
}
num_ptr += num_size;
den_ptr += den_size;
for (i = num_size - den_size - 1; i >= 0; i--)
{
mp_limb q;
mp_limb n1;
mp_limb w_cy;
mp_limb d, c;
mp_size j;
num_ptr--;
if (n0 == d0)
/* This might over-estimate q, but it's probably not worth
the extra code here to find out. */
q = ~0;
else
{
mp_limb r;
udiv_qrnnd (q, r, n0, num_ptr[-1], d0);
umul_ppmm (n1, n0, d1, q);
while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
{
q--;
r += d0;
if (r < d0) /* I.e. "carry in previous addition?" */
break;
n1 -= n0 < d1;
n0 -= d1;
}
}
w_cy = 0;
j = -den_size;
do
{
d = den_ptr[j];
c = num_ptr[j];
umul_ppmm (n1, n0, d, q);
n0 += w_cy;
w_cy = (n0 < w_cy) + n1;
n0 = c - n0;
num_ptr[j] = n0;
if (n0 > c)
goto cy_loop;
ncy_loop:
j++;
}
while (j < 0);
if (ugly_hack_flag)
{
c = 0;
ugly_hack_flag = 0;
}
else
c = num_ptr[j];
if (c >= w_cy)
goto store_q;
goto add_back;
do
{
d = den_ptr[j];
c = num_ptr[j];
umul_ppmm (n1, n0, d, q);
n0 += w_cy;
w_cy = (n0 < w_cy) + n1;
n0 = c - n0 - 1;
num_ptr[j] = n0;
if (n0 < c)
goto ncy_loop;
cy_loop:
j++;
}
while (j < 0);
if (ugly_hack_flag)
{
c = 0;
ugly_hack_flag = 0;
}
else
c = num_ptr[j];
w_cy++;
if (c >= w_cy)
goto store_q;
add_back:
j = -den_size;
do
{
d = den_ptr[j];
n0 = num_ptr[j] + d;
num_ptr[j] = n0;
if (n0 < d)
goto ab_cy_loop;
ab_ncy_loop:
j++;
}
while (j < 0);
abort (); /* We should always have a carry out! */
do
{
d = den_ptr[j];
n0 = num_ptr[j] + d + 1;
num_ptr[j] = n0;
if (n0 > d)
goto ab_ncy_loop;
ab_cy_loop:
j++;
}
while (j < 0);
q--;
store_q:
quot_ptr[i] = q;
}
}
}
return q_is_long;
}

185
gnu/lib/libgmp/mpn_dm_1.c Normal file
View file

@ -0,0 +1,185 @@
/* mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
Return the single-limb remainder.
There are no constraints on the value of the divisor.
QUOT_PTR and DIVIDEND_PTR might point to the same limb.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#ifndef UMUL_TIME
#define UMUL_TIME 1
#endif
#ifndef UDIV_TIME
#define UDIV_TIME UMUL_TIME
#endif
#if UDIV_TIME > 2 * UMUL_TIME
#undef UDIV_NEEDS_NORMALIZATION
#define UDIV_NEEDS_NORMALIZATION 1
#endif
#define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
do { \
unsigned long int _q, _ql, _r; \
unsigned long int _xh, _xl; \
umul_ppmm (_q, _ql, (nh), (di)); \
_q += (nh); /* DI is 2**BITS_PER_MP_LIMB too small. */\
umul_ppmm (_xh, _xl, _q, (d)); \
sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl); \
if (_xh != 0) \
{ \
sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
_q += 1; \
if (_xh != 0) \
{ \
sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
_q += 1; \
} \
} \
if (_r >= (d)) \
{ \
_r -= (d); \
_q += 1; \
} \
(r) = _r; \
(q) = _q; \
} while (0)
mp_limb
#ifdef __STDC__
mpn_divmod_1 (mp_ptr quot_ptr,
mp_srcptr dividend_ptr, mp_size dividend_size,
unsigned long int divisor_limb)
#else
mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
mp_ptr quot_ptr;
mp_srcptr dividend_ptr;
mp_size dividend_size;
unsigned long int divisor_limb;
#endif
{
mp_size i;
mp_limb n1, n0, r;
/* Botch: Should this be handled at all? Rely on callers? */
if (dividend_size == 0)
return 0;
if (UDIV_NEEDS_NORMALIZATION)
{
int normalization_steps;
count_leading_zeros (normalization_steps, divisor_limb);
if (normalization_steps != 0)
{
divisor_limb <<= normalization_steps;
n1 = dividend_ptr[dividend_size - 1];
r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
/* Possible optimization:
if (r == 0
&& divisor_limb > ((n1 << normalization_steps)
| (dividend_ptr[dividend_size - 2] >> ...)))
...one division less...
[Don't forget to zero most sign. quotient limb!] */
/* If multiplication is much faster than division, and the
dividend is large, pre-invert the divisor, and use
only multiplications in the inner loop. */
if (UDIV_TIME > 2 * UMUL_TIME && dividend_size >= 4)
{
mp_limb divisor_limb_inverted;
int dummy;
/* Compute (2**64 - 2**32 * DIVISOR_LIMB) / DIVISOR_LIMB.
The result is an 33-bit approximation to 1/DIVISOR_LIMB,
with the most significant bit (weight 2**32) implicit. */
/* Special case for DIVISOR_LIMB == 100...000. */
if (divisor_limb << 1 == 0)
divisor_limb_inverted = ~0;
else
udiv_qrnnd (divisor_limb_inverted, dummy,
-divisor_limb, 0, divisor_limb);
for (i = dividend_size - 2; i >= 0; i--)
{
n0 = dividend_ptr[i];
udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
((n1 << normalization_steps)
| (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
divisor_limb, divisor_limb_inverted);
n1 = n0;
}
udiv_qrnnd_preinv (quot_ptr[0], r, r,
n1 << normalization_steps,
divisor_limb, divisor_limb_inverted);
return r >> normalization_steps;
}
else
{
for (i = dividend_size - 2; i >= 0; i--)
{
n0 = dividend_ptr[i];
udiv_qrnnd (quot_ptr[i + 1], r, r,
((n1 << normalization_steps)
| (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
divisor_limb);
n1 = n0;
}
udiv_qrnnd (quot_ptr[0], r, r,
n1 << normalization_steps,
divisor_limb);
return r >> normalization_steps;
}
}
}
/* No normalization needed, either because udiv_qrnnd doesn't require
it, or because DIVISOR_LIMB is already normalized. */
i = dividend_size - 1;
r = dividend_ptr[i];
if (r >= divisor_limb)
{
r = 0;
}
else
{
/* Callers expect the quotient to be DIVIDEND_SIZE limbs. Store
a leading zero to make that expectation come true. */
quot_ptr[i] = 0;
i--;
}
for (; i >= 0; i--)
{
n0 = dividend_ptr[i];
udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
}
return r;
}

83
gnu/lib/libgmp/mpn_lshift.c Normal file
View file

@ -0,0 +1,83 @@
/* mpn_lshift -- Shift left low level.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the left
and store the USIZE least significant digits of the result at WP.
Return the bits shifted out from the most significant digit.
Argument constraints:
0. U must be normalized (i.e. it's most significant digit != 0).
1. 0 <= CNT < BITS_PER_MP_LIMB
2. If the result is to be written over the input, WP must be >= UP.
*/
mp_limb
#ifdef __STDC__
mpn_lshift (mp_ptr wp,
mp_srcptr up, mp_size usize,
unsigned cnt)
#else
mpn_lshift (wp, up, usize, cnt)
mp_ptr wp;
mp_srcptr up;
mp_size usize;
unsigned cnt;
#endif
{
mp_limb high_limb, low_limb;
unsigned sh_1, sh_2;
mp_size i;
mp_limb retval;
if (usize == 0)
return 0;
sh_1 = cnt;
if (sh_1 == 0)
{
if (wp != up)
{
/* Copy from high end to low end, to allow specified input/output
overlapping. */
for (i = usize - 1; i >= 0; i--)
wp[i] = up[i];
}
return 0;
}
wp += 1;
sh_2 = BITS_PER_MP_LIMB - sh_1;
i = usize - 1;
low_limb = up[i];
retval = low_limb >> sh_2;
high_limb = low_limb;
while (--i >= 0)
{
low_limb = up[i];
wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
high_limb = low_limb;
}
wp[i] = high_limb << sh_1;
return retval;
}

104
gnu/lib/libgmp/mpn_mod_1.c Normal file
View file

@ -0,0 +1,104 @@
/* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
Return the single-limb remainder.
There are no constraints on the value of the divisor.
QUOT_PTR and DIVIDEND_PTR might point to the same limb.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
mp_limb
#ifdef __STDC__
mpn_mod_1 (mp_srcptr dividend_ptr, mp_size dividend_size,
unsigned long int divisor_limb)
#else
mpn_mod_1 (dividend_ptr, dividend_size, divisor_limb)
mp_srcptr dividend_ptr;
mp_size dividend_size;
unsigned long int divisor_limb;
#endif
{
int normalization_steps;
mp_size i;
mp_limb n1, n0, r;
int dummy;
/* Botch: Should this be handled at all? Rely on callers? */
if (dividend_size == 0)
return 0;
if (UDIV_NEEDS_NORMALIZATION)
{
count_leading_zeros (normalization_steps, divisor_limb);
if (normalization_steps != 0)
{
divisor_limb <<= normalization_steps;
n1 = dividend_ptr[dividend_size - 1];
r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
/* Possible optimization:
if (r == 0
&& divisor_limb > ((n1 << normalization_steps)
| (dividend_ptr[dividend_size - 2] >> ...)))
...one division less...
*/
for (i = dividend_size - 2; i >= 0; i--)
{
n0 = dividend_ptr[i];
udiv_qrnnd (dummy, r, r,
((n1 << normalization_steps)
| (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
divisor_limb);
n1 = n0;
}
udiv_qrnnd (dummy, r, r,
n1 << normalization_steps,
divisor_limb);
return r >> normalization_steps;
}
}
/* No normalization needed, either because udiv_qrnnd doesn't require
it, or because DIVISOR_LIMB is already normalized. */
i = dividend_size - 1;
r = dividend_ptr[i];
if (r >= divisor_limb)
{
r = 0;
}
else
{
i--;
}
for (; i >= 0; i--)
{
n0 = dividend_ptr[i];
udiv_qrnnd (dummy, r, r, n0, divisor_limb);
}
return r;
}

414
gnu/lib/libgmp/mpn_mul.c Normal file
View file

@ -0,0 +1,414 @@
/* mpn_mul -- Multiply two natural numbers.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#ifdef DEBUG
#define MPN_MUL_VERIFY(res_ptr,res_size,op1_ptr,op1_size,op2_ptr,op2_size) \
mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size)
#include <stdio.h>
static void
mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size)
mp_ptr res_ptr, op1_ptr, op2_ptr;
mp_size res_size, op1_size, op2_size;
{
mp_ptr tmp_ptr;
mp_size tmp_size;
tmp_ptr = alloca ((op1_size + op2_size) * BYTES_PER_MP_LIMB);
if (op1_size >= op2_size)
tmp_size = mpn_mul_classic (tmp_ptr,
op1_ptr, op1_size, op2_ptr, op2_size);
else
tmp_size = mpn_mul_classic (tmp_ptr,
op2_ptr, op2_size, op1_ptr, op1_size);
if (tmp_size != res_size
|| mpn_cmp (tmp_ptr, res_ptr, tmp_size) != 0)
{
fprintf (stderr, "GNU MP internal error: Wrong result in mpn_mul.\n");
fprintf (stderr, "op1{%d} = ", op1_size); mpn_dump (op1_ptr, op1_size);
fprintf (stderr, "op2{%d} = ", op2_size); mpn_dump (op2_ptr, op2_size);
abort ();
}
}
#else
#define MPN_MUL_VERIFY(a,b,c,d,e,f)
#endif
/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
and v (pointed to by VP, with VSIZE limbs), and store the result at
PRODP. USIZE + VSIZE limbs are always stored, but if the input
operands are normalized, the return value will reflect the true
result size (which is either USIZE + VSIZE, or USIZE + VSIZE -1).
NOTE: The space pointed to by PRODP is overwritten before finished
with U and V, so overlap is an error.
Argument constraints:
1. USIZE >= VSIZE.
2. PRODP != UP and PRODP != VP, i.e. the destination
must be distinct from the multiplier and the multiplicand. */
/* If KARATSUBA_THRESHOLD is not already defined, define it to a
value which is good on most machines. */
#ifndef KARATSUBA_THRESHOLD
#define KARATSUBA_THRESHOLD 8
#endif
/* The code can't handle KARATSUBA_THRESHOLD smaller than 4. */
#if KARATSUBA_THRESHOLD < 4
#undef KARATSUBA_THRESHOLD
#define KARATSUBA_THRESHOLD 4
#endif
mp_size
#ifdef __STDC__
mpn_mul (mp_ptr prodp,
mp_srcptr up, mp_size usize,
mp_srcptr vp, mp_size vsize)
#else
mpn_mul (prodp, up, usize, vp, vsize)
mp_ptr prodp;
mp_srcptr up;
mp_size usize;
mp_srcptr vp;
mp_size vsize;
#endif
{
mp_size n;
mp_size prod_size;
mp_limb cy;
if (vsize < KARATSUBA_THRESHOLD)
{
/* Handle simple cases with traditional multiplication.
This is the most critical code of the entire function. All
multiplies rely on this, both small and huge. Small ones arrive
here immediately. Huge ones arrive here as this is the base case
for the recursive algorithm below. */
mp_size i, j;
mp_limb prod_low, prod_high;
mp_limb cy_limb;
mp_limb v_limb;
if (vsize == 0)
return 0;
/* Offset UP and PRODP so that the inner loop can be faster. */
up += usize;
prodp += usize;
/* Multiply by the first limb in V separately, as the result can
be stored (not added) to PROD. We also avoid a loop for zeroing. */
v_limb = vp[0];
if (v_limb <= 1)
{
if (v_limb == 1)
MPN_COPY (prodp - usize, up - usize, usize);
else
MPN_ZERO (prodp - usize, usize);
cy_limb = 0;
}
else
{
cy_limb = 0;
j = -usize;
do
{
umul_ppmm (prod_high, prod_low, up[j], v_limb);
add_ssaaaa (cy_limb, prodp[j], prod_high, prod_low, 0, cy_limb);
j++;
}
while (j < 0);
}
prodp[0] = cy_limb;
prodp++;
/* For each iteration in the outer loop, multiply one limb from
U with one limb from V, and add it to PROD. */
for (i = 1; i < vsize; i++)
{
v_limb = vp[i];
if (v_limb <= 1)
{
cy_limb = 0;
if (v_limb == 1)
cy_limb = mpn_add (prodp - usize,
prodp - usize, usize, up - usize, usize);
}
else
{
cy_limb = 0;
j = -usize;
do
{
umul_ppmm (prod_high, prod_low, up[j], v_limb);
add_ssaaaa (cy_limb, prod_low,
prod_high, prod_low, 0, cy_limb);
add_ssaaaa (cy_limb, prodp[j],
cy_limb, prod_low, 0, prodp[j]);
j++;
}
while (j < 0);
}
prodp[0] = cy_limb;
prodp++;
}
return usize + vsize - (cy_limb == 0);
}
n = (usize + 1) / 2;
/* Is USIZE larger than 1.5 times VSIZE? Avoid Karatsuba's algorithm. */
if (2 * usize > 3 * vsize)
{
/* If U has at least twice as many limbs as V. Split U in two
pieces, U1 and U0, such that U = U0 + U1*(2**BITS_PER_MP_LIMB)**N,
and recursively multiply the two pieces separately with V. */
mp_size u0_size;
mp_ptr tmp;
mp_size tmp_size;
/* V1 (the high part of V) is zero. */
/* Calculate the length of U0. It is normally equal to n, but
of course not for sure. */
for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--)
;
/* Perform (U0 * V). */
if (u0_size >= vsize)
prod_size = mpn_mul (prodp, up, u0_size, vp, vsize);
else
prod_size = mpn_mul (prodp, vp, vsize, up, u0_size);
MPN_MUL_VERIFY (prodp, prod_size, up, u0_size, vp, vsize);
/* We have to zero-extend the lower partial product to n limbs,
since the mpn_add some lines below expect the first n limbs
to be well defined. (This is normally a no-op. It may
do something when U1 has many leading 0 limbs.) */
while (prod_size < n)
prodp[prod_size++] = 0;
tmp = (mp_ptr) alloca ((usize + vsize - n) * BYTES_PER_MP_LIMB);
/* Perform (U1 * V). Make sure the first source argument to mpn_mul
is not less than the second source argument. */
if (vsize <= usize - n)
tmp_size = mpn_mul (tmp, up + n, usize - n, vp, vsize);
else
tmp_size = mpn_mul (tmp, vp, vsize, up + n, usize - n);
MPN_MUL_VERIFY (tmp, tmp_size, up + n, usize - n, vp, vsize);
/* In this addition hides a potentially large copying of TMP. */
if (prod_size - n >= tmp_size)
cy = mpn_add (prodp + n, prodp + n, prod_size - n, tmp, tmp_size);
else
cy = mpn_add (prodp + n, tmp, tmp_size, prodp + n, prod_size - n);
if (cy)
abort (); /* prodp[prod_size] = cy; */
alloca (0);
return tmp_size + n;
}
else
{
/* Karatsuba's divide-and-conquer algorithm.
Split U in two pieces, U1 and U0, such that
U = U0 + U1*(B**n),
and V in V1 and V0, such that
V = V0 + V1*(B**n).
UV is then computed recursively using the identity
2n n n n
UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
1 1 1 0 0 1 0 0
Where B = 2**BITS_PER_MP_LIMB.
*/
/* It's possible to decrease the temporary allocation by using the
prodp area for temporary storage of the middle term, and doing
that recursive multiplication first. (Do this later.) */
mp_size u0_size;
mp_size v0_size;
mp_size u0v0_size;
mp_size u1v1_size;
mp_ptr temp;
mp_size temp_size;
mp_size utem_size;
mp_size vtem_size;
mp_ptr ptem;
mp_size ptem_size;
int negflg;
mp_ptr pp;
pp = (mp_ptr) alloca (4 * n * BYTES_PER_MP_LIMB);
/* Calculate the lengths of U0 and V0. They are normally equal
to n, but of course not for sure. */
for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--)
;
for (v0_size = n; v0_size > 0 && vp[v0_size - 1] == 0; v0_size--)
;
/*** 1. PROD]2n..0] := U0 x V0
(Recursive call to mpn_mul may NOT overwrite input operands.)
________________ ________________
|________________||____U0 x V0_____| */
if (u0_size >= v0_size)
u0v0_size = mpn_mul (pp, up, u0_size, vp, v0_size);
else
u0v0_size = mpn_mul (pp, vp, v0_size, up, u0_size);
MPN_MUL_VERIFY (pp, u0v0_size, up, u0_size, vp, v0_size);
/* Zero-extend to 2n limbs. */
while (u0v0_size < 2 * n)
pp[u0v0_size++] = 0;
/*** 2. PROD]4n..2n] := U1 x V1
(Recursive call to mpn_mul may NOT overwrite input operands.)
________________ ________________
|_____U1 x V1____||____U0 x V0_____| */
u1v1_size = mpn_mul (pp + 2*n,
up + n, usize - n,
vp + n, vsize - n);
MPN_MUL_VERIFY (pp + 2*n, u1v1_size,
up + n, usize - n, vp + n, vsize - n);
prod_size = 2 * n + u1v1_size;
/*** 3. PTEM]2n..0] := (U1-U0) x (V0-V1)
(Recursive call to mpn_mul may overwrite input operands.)
________________
|_(U1-U0)(V0-V1)_| */
temp = (mp_ptr) alloca ((2 * n + 1) * BYTES_PER_MP_LIMB);
if (usize - n > u0_size
|| (usize - n == u0_size
&& mpn_cmp (up + n, up, u0_size) >= 0))
{
utem_size = usize - n
+ mpn_sub (temp, up + n, usize - n, up, u0_size);
negflg = 0;
}
else
{
utem_size = u0_size
+ mpn_sub (temp, up, u0_size, up + n, usize - n);
negflg = 1;
}
if (vsize - n > v0_size
|| (vsize - n == v0_size
&& mpn_cmp (vp + n, vp, v0_size) >= 0))
{
vtem_size = vsize - n
+ mpn_sub (temp + n, vp + n, vsize - n, vp, v0_size);
negflg ^= 1;
}
else
{
vtem_size = v0_size
+ mpn_sub (temp + n, vp, v0_size, vp + n, vsize - n);
/* No change of NEGFLG. */
}
ptem = (mp_ptr) alloca (2 * n * BYTES_PER_MP_LIMB);
if (utem_size >= vtem_size)
ptem_size = mpn_mul (ptem, temp, utem_size, temp + n, vtem_size);
else
ptem_size = mpn_mul (ptem, temp + n, vtem_size, temp, utem_size);
MPN_MUL_VERIFY (ptem, ptem_size, temp, utem_size, temp + n, vtem_size);
/*** 4. TEMP]2n..0] := PROD]2n..0] + PROD]4n..2n]
________________
|_____U1 x V1____|
________________
|_____U0_x_V0____| */
cy = mpn_add (temp, pp, 2*n, pp + 2*n, u1v1_size);
if (cy != 0)
{
temp[2*n] = cy;
temp_size = 2*n + 1;
}
else
{
/* Normalize temp. pp[2*n-1] might have been zero in the
mpn_add call above, and thus temp might be unnormalized. */
for (temp_size = 2*n; temp_size > 0 && temp[temp_size - 1] == 0;
temp_size--)
;
}
if (prod_size - n >= temp_size)
cy = mpn_add (pp + n, pp + n, prod_size - n, temp, temp_size);
else
{
/* This is a weird special case that should not happen (often)! */
cy = mpn_add (pp + n, temp, temp_size, pp + n, prod_size - n);
prod_size = temp_size + n;
}
if (cy != 0)
{
pp[prod_size] = cy;
prod_size++;
}
#ifdef DEBUG
if (prod_size > 4 * n)
abort();
#endif
if (negflg)
prod_size = prod_size
+ mpn_sub (pp + n, pp + n, prod_size - n, ptem, ptem_size);
else
{
if (prod_size - n < ptem_size)
abort();
cy = mpn_add (pp + n, pp + n, prod_size - n, ptem, ptem_size);
if (cy != 0)
{
pp[prod_size] = cy;
prod_size++;
#ifdef DEBUG
if (prod_size > 4 * n)
abort();
#endif
}
}
MPN_COPY (prodp, pp, prod_size);
alloca (0);
return prod_size;
}
}

90
gnu/lib/libgmp/mpn_rshift.c Normal file
View file

@ -0,0 +1,90 @@
/* mpn_rshift -- Shift right a low-level natural-number integer.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
and store the USIZE least significant limbs of the result at WP.
Return the size of the result.
Argument constraints:
0. U must be normalized (i.e. it's most significant limb != 0).
1. 0 <= CNT < BITS_PER_MP_LIMB
2. If the result is to be written over the input, WP must be <= UP.
*/
mp_size
#ifdef __STDC__
mpn_rshift (mp_ptr wp,
mp_srcptr up, mp_size usize,
unsigned cnt)
#else
mpn_rshift (wp, up, usize, cnt)
mp_ptr wp;
mp_srcptr up;
mp_size usize;
unsigned cnt;
#endif
{
mp_limb high_limb, low_limb;
unsigned sh_1, sh_2;
mp_size i;
if (usize == 0)
return 0;
sh_1 = cnt;
if (sh_1 == 0)
{
if (wp != up)
{
/* Copy from low end to high end, to allow specified input/output
overlapping. */
for (i = 0; i < usize; i++)
wp[i] = up[i];
}
return usize;
}
wp -= 1;
sh_2 = BITS_PER_MP_LIMB - sh_1;
high_limb = up[0];
#if 0
if (cy_limb != NULL)
*cy_limb = high_limb << sh_2;
#endif
low_limb = high_limb;
for (i = 1; i < usize; i++)
{
high_limb = up[i];
wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
low_limb = high_limb;
}
low_limb >>= sh_1;
if (low_limb != 0)
{
wp[i] = low_limb;
return usize;
}
return usize - 1;
}

86
gnu/lib/libgmp/mpn_rshiftci.c Normal file
View file

@ -0,0 +1,86 @@
/* mpn_rshiftci -- Shift a low level natural-number integer with carry in.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the right
and store the USIZE least significant digits of the result at WP.
Return the size of the result.
Argument constraints:
0. U must be normalized (i.e. it's most significant digit != 0).
1. 0 <= CNT < BITS_PER_MP_LIMB
2. If the result is to be written over the input, WP must be <= UP.
*/
mp_size
#ifdef __STDC__
mpn_rshiftci (mp_ptr wp,
mp_srcptr up, mp_size usize,
unsigned cnt,
mp_limb carry_in)
#else
mpn_rshiftci (wp, up, usize, cnt, carry_in)
mp_ptr wp;
mp_srcptr up;
mp_size usize;
unsigned cnt;
mp_limb carry_in;
#endif
{
mp_limb high_limb, low_limb;
unsigned sh_1, sh_2;
mp_size i;
if (usize <= 0)
return 0;
sh_1 = cnt;
if (sh_1 == 0)
{
if (wp != up)
{
/* Copy from low end to high end, to allow specified input/output
overlapping. */
for (i = 0; i < usize; i++)
wp[i] = up[i];
}
return usize;
}
wp -= 1;
sh_2 = BITS_PER_MP_LIMB - sh_1;
low_limb = up[0];
for (i = 1; i < usize; i++)
{
high_limb = up[i];
wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
low_limb = high_limb;
}
low_limb = (low_limb >> sh_1) | (carry_in << sh_2);
if (low_limb != 0)
{
wp[i] = low_limb;
return usize;
}
return usize - 1;
}

479
gnu/lib/libgmp/mpn_sqrt.c Normal file
View file

@ -0,0 +1,479 @@
/* mpn_sqrt(root_ptr, rem_ptr, op_ptr, op_size)
Write the square root of {OP_PTR, OP_SIZE} at ROOT_PTR.
Write the remainder at REM_PTR, if REM_PTR != NULL.
Return the size of the remainder.
(The size of the root is always half of the size of the operand.)
OP_PTR and ROOT_PTR may not point to the same object.
OP_PTR and REM_PTR may point to the same object.
If REM_PTR is NULL, only the root is computed and the return value of
the function is 0 if OP is a perfect square, and *any* non-zero number
otherwise.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
/* This code is just correct if "unsigned char" has at least 8 bits. It
doesn't help to use CHAR_BIT from limits.h, as the real problem is
the static arrays. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
/* Square root algorithm:
1. Shift OP (the input) to the left an even number of bits s.t. there
are an even number of words and either (or both) of the most
significant bits are set. This way, sqrt(OP) has exactly half as
many words as OP, and has its most significant bit set.
2. Get a 9-bit approximation to sqrt(OP) using the pre-computed tables.
This approximation is used for the first single-precision
iterations of Newton's method, yielding a full-word approximation
to sqrt(OP).
3. Perform multiple-precision Newton iteration until we have the
exact result. Only about half of the input operand is used in
this calculation, as the square root is perfectly determinable
from just the higher half of a number. */
/* Define this macro for IEEE P854 machines with a fast sqrt instruction. */
#if defined __GNUC__
#if defined __sparc__
#define SQRT(a) \
({ \
double __sqrt_res; \
asm ("fsqrtd %1,%0" : "=f" (__sqrt_res) : "f" (a)); \
__sqrt_res; \
})
#endif
#if defined __HAVE_68881__
#define SQRT(a) \
({ \
double __sqrt_res; \
asm ("fsqrtx %1,%0" : "=f" (__sqrt_res) : "f" (a)); \
__sqrt_res; \
})
#endif
#if defined __hppa
#define SQRT(a) \
({ \
double __sqrt_res; \
asm ("fsqrt,dbl %1,%0" : "=fx" (__sqrt_res) : "fx" (a)); \
__sqrt_res; \
})
#endif
#endif
#ifndef SQRT
/* Tables for initial approximation of the square root. These are
indexed with bits 1-8 of the operand for which the square root is
calculated, where bit 0 is the most significant non-zero bit. I.e.
the most significant one-bit is not used, since that per definition
is one. Likewise, the tables don't return the highest bit of the
result. That bit must be inserted by or:ing the returned value with
0x100. This way, we get a 9-bit approximation from 8-bit tables! */
/* Table to be used for operands with an even total number of bits.
(Exactly as in the decimal system there are similarities between the
square root of numbers with the same initial digits and an even
difference in the total number of digits. Consider the square root
of 1, 10, 100, 1000, ...) */
static unsigned char even_approx_tab[256] =
{
0x6a, 0x6a, 0x6b, 0x6c, 0x6c, 0x6d, 0x6e, 0x6e,
0x6f, 0x70, 0x71, 0x71, 0x72, 0x73, 0x73, 0x74,
0x75, 0x75, 0x76, 0x77, 0x77, 0x78, 0x79, 0x79,
0x7a, 0x7b, 0x7b, 0x7c, 0x7d, 0x7d, 0x7e, 0x7f,
0x80, 0x80, 0x81, 0x81, 0x82, 0x83, 0x83, 0x84,
0x85, 0x85, 0x86, 0x87, 0x87, 0x88, 0x89, 0x89,
0x8a, 0x8b, 0x8b, 0x8c, 0x8d, 0x8d, 0x8e, 0x8f,
0x8f, 0x90, 0x90, 0x91, 0x92, 0x92, 0x93, 0x94,
0x94, 0x95, 0x96, 0x96, 0x97, 0x97, 0x98, 0x99,
0x99, 0x9a, 0x9b, 0x9b, 0x9c, 0x9c, 0x9d, 0x9e,
0x9e, 0x9f, 0xa0, 0xa0, 0xa1, 0xa1, 0xa2, 0xa3,
0xa3, 0xa4, 0xa4, 0xa5, 0xa6, 0xa6, 0xa7, 0xa7,
0xa8, 0xa9, 0xa9, 0xaa, 0xaa, 0xab, 0xac, 0xac,
0xad, 0xad, 0xae, 0xaf, 0xaf, 0xb0, 0xb0, 0xb1,
0xb2, 0xb2, 0xb3, 0xb3, 0xb4, 0xb5, 0xb5, 0xb6,
0xb6, 0xb7, 0xb7, 0xb8, 0xb9, 0xb9, 0xba, 0xba,
0xbb, 0xbb, 0xbc, 0xbd, 0xbd, 0xbe, 0xbe, 0xbf,
0xc0, 0xc0, 0xc1, 0xc1, 0xc2, 0xc2, 0xc3, 0xc3,
0xc4, 0xc5, 0xc5, 0xc6, 0xc6, 0xc7, 0xc7, 0xc8,
0xc9, 0xc9, 0xca, 0xca, 0xcb, 0xcb, 0xcc, 0xcc,
0xcd, 0xce, 0xce, 0xcf, 0xcf, 0xd0, 0xd0, 0xd1,
0xd1, 0xd2, 0xd3, 0xd3, 0xd4, 0xd4, 0xd5, 0xd5,
0xd6, 0xd6, 0xd7, 0xd7, 0xd8, 0xd9, 0xd9, 0xda,
0xda, 0xdb, 0xdb, 0xdc, 0xdc, 0xdd, 0xdd, 0xde,
0xde, 0xdf, 0xe0, 0xe0, 0xe1, 0xe1, 0xe2, 0xe2,
0xe3, 0xe3, 0xe4, 0xe4, 0xe5, 0xe5, 0xe6, 0xe6,
0xe7, 0xe7, 0xe8, 0xe8, 0xe9, 0xea, 0xea, 0xeb,
0xeb, 0xec, 0xec, 0xed, 0xed, 0xee, 0xee, 0xef,
0xef, 0xf0, 0xf0, 0xf1, 0xf1, 0xf2, 0xf2, 0xf3,
0xf3, 0xf4, 0xf4, 0xf5, 0xf5, 0xf6, 0xf6, 0xf7,
0xf7, 0xf8, 0xf8, 0xf9, 0xf9, 0xfa, 0xfa, 0xfb,
0xfb, 0xfc, 0xfc, 0xfd, 0xfd, 0xfe, 0xfe, 0xff,
};
/* Table to be used for operands with an odd total number of bits.
(Further comments before previous table.) */
static unsigned char odd_approx_tab[256] =
{
0x00, 0x00, 0x00, 0x01, 0x01, 0x02, 0x02, 0x03,
0x03, 0x04, 0x04, 0x05, 0x05, 0x06, 0x06, 0x07,
0x07, 0x08, 0x08, 0x09, 0x09, 0x0a, 0x0a, 0x0b,
0x0b, 0x0c, 0x0c, 0x0d, 0x0d, 0x0e, 0x0e, 0x0f,
0x0f, 0x10, 0x10, 0x10, 0x11, 0x11, 0x12, 0x12,
0x13, 0x13, 0x14, 0x14, 0x15, 0x15, 0x16, 0x16,
0x16, 0x17, 0x17, 0x18, 0x18, 0x19, 0x19, 0x1a,
0x1a, 0x1b, 0x1b, 0x1b, 0x1c, 0x1c, 0x1d, 0x1d,
0x1e, 0x1e, 0x1f, 0x1f, 0x20, 0x20, 0x20, 0x21,
0x21, 0x22, 0x22, 0x23, 0x23, 0x23, 0x24, 0x24,
0x25, 0x25, 0x26, 0x26, 0x27, 0x27, 0x27, 0x28,
0x28, 0x29, 0x29, 0x2a, 0x2a, 0x2a, 0x2b, 0x2b,
0x2c, 0x2c, 0x2d, 0x2d, 0x2d, 0x2e, 0x2e, 0x2f,
0x2f, 0x30, 0x30, 0x30, 0x31, 0x31, 0x32, 0x32,
0x32, 0x33, 0x33, 0x34, 0x34, 0x35, 0x35, 0x35,
0x36, 0x36, 0x37, 0x37, 0x37, 0x38, 0x38, 0x39,
0x39, 0x39, 0x3a, 0x3a, 0x3b, 0x3b, 0x3b, 0x3c,
0x3c, 0x3d, 0x3d, 0x3d, 0x3e, 0x3e, 0x3f, 0x3f,
0x40, 0x40, 0x40, 0x41, 0x41, 0x41, 0x42, 0x42,
0x43, 0x43, 0x43, 0x44, 0x44, 0x45, 0x45, 0x45,
0x46, 0x46, 0x47, 0x47, 0x47, 0x48, 0x48, 0x49,
0x49, 0x49, 0x4a, 0x4a, 0x4b, 0x4b, 0x4b, 0x4c,
0x4c, 0x4c, 0x4d, 0x4d, 0x4e, 0x4e, 0x4e, 0x4f,
0x4f, 0x50, 0x50, 0x50, 0x51, 0x51, 0x51, 0x52,
0x52, 0x53, 0x53, 0x53, 0x54, 0x54, 0x54, 0x55,
0x55, 0x56, 0x56, 0x56, 0x57, 0x57, 0x57, 0x58,
0x58, 0x59, 0x59, 0x59, 0x5a, 0x5a, 0x5a, 0x5b,
0x5b, 0x5b, 0x5c, 0x5c, 0x5d, 0x5d, 0x5d, 0x5e,
0x5e, 0x5e, 0x5f, 0x5f, 0x60, 0x60, 0x60, 0x61,
0x61, 0x61, 0x62, 0x62, 0x62, 0x63, 0x63, 0x63,
0x64, 0x64, 0x65, 0x65, 0x65, 0x66, 0x66, 0x66,
0x67, 0x67, 0x67, 0x68, 0x68, 0x68, 0x69, 0x69,
};
#endif
mp_size
#ifdef __STDC__
mpn_sqrt (mp_ptr root_ptr, mp_ptr rem_ptr, mp_srcptr op_ptr, mp_size op_size)
#else
mpn_sqrt (root_ptr, rem_ptr, op_ptr, op_size)
mp_ptr root_ptr;
mp_ptr rem_ptr;
mp_srcptr op_ptr;
mp_size op_size;
#endif
{
/* R (root result) */
mp_ptr rp; /* Pointer to least significant word */
mp_size rsize; /* The size in words */
/* T (OP shifted to the left a.k.a. normalized) */
mp_ptr tp; /* Pointer to least significant word */
mp_size tsize; /* The size in words */
mp_ptr t_end_ptr; /* Pointer right beyond most sign. word */
mp_limb t_high0, t_high1; /* The two most significant words */
/* TT (temporary for numerator/remainder) */
mp_ptr ttp; /* Pointer to least significant word */
/* X (temporary for quotient in main loop) */
mp_ptr xp; /* Pointer to least significant word */
mp_size xsize; /* The size in words */
unsigned cnt;
mp_limb initial_approx; /* Initially made approximation */
mp_size tsizes[BITS_PER_MP_LIMB]; /* Successive calculation precisions */
mp_size tmp;
mp_size i;
/* If OP is zero, both results are zero. */
if (op_size == 0)
return 0;
count_leading_zeros (cnt, op_ptr[op_size - 1]);
tsize = op_size;
if ((tsize & 1) != 0)
{
cnt += BITS_PER_MP_LIMB;
tsize++;
}
rsize = tsize / 2;
rp = root_ptr;
/* Shift OP an even number of bits into T, such that either the most or
the second most significant bit is set, and such that the number of
words in T becomes even. This way, the number of words in R=sqrt(OP)
is exactly half as many as in OP, and the most significant bit of R
is set.
Also, the initial approximation is simplified by this up-shifted OP.
Finally, the Newtonian iteration which is the main part of this
program performs division by R. The fast division routine expects
the divisor to be "normalized" in exactly the sense of having the
most significant bit set. */
tp = (mp_ptr) alloca (tsize * BYTES_PER_MP_LIMB);
t_high0 = mpn_lshift (tp + cnt / BITS_PER_MP_LIMB, op_ptr, op_size,
(cnt & ~1) % BITS_PER_MP_LIMB);
if (cnt >= BITS_PER_MP_LIMB)
tp[0] = 0;
t_high0 = tp[tsize - 1];
t_high1 = tp[tsize - 2]; /* Never stray. TSIZE is >= 2. */
/* Is there a fast sqrt instruction defined for this machine? */
#ifdef SQRT
{
initial_approx = SQRT (t_high0 * 2.0
* ((mp_limb) 1 << (BITS_PER_MP_LIMB - 1))
+ t_high1);
/* If t_high0,,t_high1 is big, the result in INITIAL_APPROX might have
become incorrect due to overflow in the conversion from double to
mp_limb above. It will typically be zero in that case, but might be
a small number on some machines. The most significant bit of
INITIAL_APPROX should be set, so that bit is a good overflow
indication. */
if ((mp_limb_signed) initial_approx >= 0)
initial_approx = ~0;
}
#else
/* Get a 9 bit approximation from the tables. The tables expect to
be indexed with the 8 high bits right below the highest bit.
Also, the highest result bit is not returned by the tables, and
must be or:ed into the result. The scheme gives 9 bits of start
approximation with just 256-entry 8 bit tables. */
if ((cnt & 1) == 0)
{
/* The most sign bit of t_high0 is set. */
initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 1);
initial_approx &= 0xff;
initial_approx = even_approx_tab[initial_approx];
}
else
{
/* The most significant bit of T_HIGH0 is unset,
the second most significant is set. */
initial_approx = t_high0 >> (BITS_PER_MP_LIMB - 8 - 2);
initial_approx &= 0xff;
initial_approx = odd_approx_tab[initial_approx];
}
initial_approx |= 0x100;
initial_approx <<= BITS_PER_MP_LIMB - 8 - 1;
/* Perform small precision Newtonian iterations to get a full word
approximation. For small operands, these iteration will make the
entire job. */
if (t_high0 == ~0)
initial_approx = t_high0;
else
{
mp_limb quot;
if (t_high0 >= initial_approx)
initial_approx = t_high0 + 1;
/* First get about 18 bits with pure C arithmetics. */
quot = t_high0 / (initial_approx >> BITS_PER_MP_LIMB/2) << BITS_PER_MP_LIMB/2;
initial_approx = (initial_approx + quot) / 2;
initial_approx |= (mp_limb) 1 << (BITS_PER_MP_LIMB - 1);
/* Now get a full word by one (or for > 36 bit machines) several
iterations. */
for (i = 16; i < BITS_PER_MP_LIMB; i <<= 1)
{
mp_limb ignored_remainder;
udiv_qrnnd (quot, ignored_remainder,
t_high0, t_high1, initial_approx);
initial_approx = (initial_approx + quot) / 2;
initial_approx |= (mp_limb) 1 << (BITS_PER_MP_LIMB - 1);
}
}
#endif
rp[0] = initial_approx;
rsize = 1;
xp = (mp_ptr) alloca (tsize * BYTES_PER_MP_LIMB);
ttp = (mp_ptr) alloca (tsize * BYTES_PER_MP_LIMB);
t_end_ptr = tp + tsize;
#ifdef DEBUG
printf ("\n\nT = ");
_mp_mout (tp, tsize);
#endif
if (tsize > 2)
{
/* Determine the successive precisions to use in the iteration. We
minimize the precisions, beginning with the highest (i.e. last
iteration) to the lowest (i.e. first iteration). */
tmp = tsize / 2;
for (i = 0;;i++)
{
tsize = (tmp + 1) / 2;
if (tmp == tsize)
break;
tsizes[i] = tsize + tmp;
tmp = tsize;
}
/* Main Newton iteration loop. For big arguments, most of the
time is spent here. */
/* It is possible to do a great optimization here. The successive
divisors in the mpn_div call below has more and more leading
words equal to its predecessor. Therefore the beginning of
each division will repeat the same work as did the last
division. If we could guarantee that the leading words of two
consecutive divisors are the same (i.e. in this case, a later
divisor has just more digits at the end) it would be a simple
matter of just using the old remainder of the last division in
a subsequent division, to take care of this optimization. This
idea would surely make a difference even for small arguments. */
/* Loop invariants:
R <= shiftdown_to_same_size(floor(sqrt(OP))) < R + 1.
X - 1 < shiftdown_to_same_size(floor(sqrt(OP))) <= X.
R <= shiftdown_to_same_size(X). */
while (--i >= 0)
{
mp_limb cy;
#ifdef DEBUG
mp_limb old_least_sign_r = rp[0];
mp_size old_rsize = rsize;
printf ("R = ");
_mp_mout (rp, rsize);
#endif
tsize = tsizes[i];
/* Need to copy the numerator into temporary space, as
mpn_div overwrites its numerator argument with the
remainder (which we currently ignore). */
MPN_COPY (ttp, t_end_ptr - tsize, tsize);
cy = mpn_div (xp, ttp, tsize, rp, rsize);
xsize = tsize - rsize;
cy = cy ? xp[xsize] : 0;
#ifdef DEBUG
printf ("X =%d", cy);
_mp_mout (xp, xsize);
#endif
/* Add X and R with the most significant limbs aligned,
temporarily ignoring at least one limb at the low end of X. */
tmp = xsize - rsize;
cy += mpn_add (xp + tmp, rp, rsize, xp + tmp, rsize);
/* If T begins with more than 2 x BITS_PER_MP_LIMB of ones, we get
intermediate roots that'd need an extra bit. We don't want to
handle that since it would make the subsequent divisor
non-normalized, so round such roots down to be only ones in the
current precision. */
if (cy == 2)
{
mp_size j;
for (j = xsize; j >= 0; j--)
xp[j] = ~(mp_limb)0;
}
/* Divide X by 2 and put the result in R. This is the new
approximation. Shift in the carry from the addition. */
rsize = mpn_rshiftci (rp, xp, xsize, 1, (mp_limb) 1);
#ifdef DEBUG
if (old_least_sign_r != rp[rsize - old_rsize])
printf (">>>>>>>> %d: %08x, %08x <<<<<<<<\n",
i, old_least_sign_r, rp[rsize - old_rsize]);
#endif
}
}
#ifdef DEBUG
printf ("(final) R = ");
_mp_mout (rp, rsize);
#endif
/* We computed the square root of OP * 2**(2*floor(cnt/2)).
This has resulted in R being 2**floor(cnt/2) to large.
Shift it down here to fix that. */
rsize = mpn_rshift (rp, rp, rsize, cnt/2);
/* Calculate the remainder. */
tsize = mpn_mul (tp, rp, rsize, rp, rsize);
if (op_size < tsize
|| (op_size == tsize && mpn_cmp (op_ptr, tp, op_size) < 0))
{
/* R is too large. Decrement it. */
mp_limb one = 1;
tsize = tsize + mpn_sub (tp, tp, tsize, rp, rsize);
tsize = tsize + mpn_sub (tp, tp, tsize, rp, rsize);
tsize = tsize + mpn_add (tp, tp, tsize, &one, 1);
(void) mpn_sub (rp, rp, rsize, &one, 1);
#ifdef DEBUG
printf ("(adjusted) R = ");
_mp_mout (rp, rsize);
#endif
}
if (rem_ptr != NULL)
{
mp_size retval = op_size + mpn_sub (rem_ptr, op_ptr, op_size, tp, tsize);
alloca (0);
return retval;
}
else
{
mp_size retval = (op_size != tsize || mpn_cmp (op_ptr, tp, op_size));
alloca (0);
return retval;
}
}
#ifdef DEBUG
_mp_mout (mp_srcptr p, mp_size size)
{
mp_size ii;
for (ii = size - 1; ii >= 0; ii--)
printf ("%08X", p[ii]);
puts ("");
}
#endif

162
gnu/lib/libgmp/mpn_sub.c Normal file
View file

@ -0,0 +1,162 @@
/* mpn_sub -- Subtract two low-level natural-number integers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Subtract SUB_PTR/SUB_SIZE from MIN_PTR/MIN_SIZE and store the
result (MIN_SIZE words) at DIF_PTR.
Return 1 if min < sub (result is negative). Otherwise, return the
negative difference between the number of words in dif and min.
(I.e. return 0 if the result has MIN_SIZE words, -1 if it has
MIN_SIZE - 1 words, etc.)
Argument constraint: MIN_SIZE >= SUB_SIZE.
The size of DIF can be calculated as MIN_SIZE + the return value. */
mp_size
#ifdef __STDC__
mpn_sub (mp_ptr dif_ptr,
mp_srcptr min_ptr, mp_size min_size,
mp_srcptr sub_ptr, mp_size sub_size)
#else
mpn_sub (dif_ptr, min_ptr, min_size, sub_ptr, sub_size)
mp_ptr dif_ptr;
mp_srcptr min_ptr;
mp_size min_size;
mp_srcptr sub_ptr;
mp_size sub_size;
#endif
{
mp_limb m, s, dif;
mp_size j;
/* The loop counter and index J goes from some negative value to zero.
This way the loops are faster. Need to offset the base pointers
to take care of the negative indices. */
j = -sub_size;
if (j == 0)
goto sub_finished;
min_ptr -= j;
sub_ptr -= j;
dif_ptr -= j;
/* There are two do-loops, marked NON-CARRY LOOP and CARRY LOOP that
jump between each other. The first loop is for when the previous
subtraction didn't produce a carry-out; the second is for the
complementary case. */
/* NON-CARRY LOOP */
do
{
m = min_ptr[j];
s = sub_ptr[j];
dif = m - s;
dif_ptr[j] = dif;
if (dif > m)
goto cy_loop;
ncy_loop:
j++;
}
while (j < 0);
/* We have exhausted SUB, with no carry out. Copy remaining part of
MIN to DIF. */
sub_finished:
j = sub_size - min_size;
/* If there's no difference between the length of the operands, the
last words might have become zero, and re-normalization is needed. */
if (j == 0)
goto normalize;
min_ptr -= j;
dif_ptr -= j;
goto copy;
/* CARRY LOOP */
do
{
m = min_ptr[j];
s = sub_ptr[j];
dif = m - s - 1;
dif_ptr[j] = dif;
if (dif < m)
goto ncy_loop;
cy_loop:
j++;
}
while (j < 0);
/* We have exhausted SUB, but need to propagate carry. */
j = sub_size - min_size;
if (j == 0)
return 1; /* min < sub. Flag it to the caller */
min_ptr -= j;
dif_ptr -= j;
/* Propagate carry. Sooner or later the carry will cancel with a
non-zero word, because the minuend is normalized. Considering this,
there's no need to test the index J. */
for (;;)
{
m = min_ptr[j];
dif = m - 1;
dif_ptr[j] = dif;
j++;
if (dif < m)
break;
}
if (j == 0)
goto normalize;
copy:
/* Don't copy the remaining words of MIN to DIF if MIN_PTR and DIF_PTR
are equal. It would just be a no-op copying. Return 0, as the length
of the result equals that of the minuend. */
if (dif_ptr == min_ptr)
return 0;
do
{
dif_ptr[j] = min_ptr[j];
j++;
}
while (j < 0);
return 0;
normalize:
for (j = -1; j >= -min_size; j--)
{
if (dif_ptr[j] != 0)
return j + 1;
}
return -min_size;
}

85
gnu/lib/libgmp/mpq_add.c Normal file
View file

@ -0,0 +1,85 @@
/* mpq_add -- add two rational numbers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_add (MP_RAT *sum, const MP_RAT *a1, const MP_RAT *a2)
#else
mpq_add (sum, a1, a2)
MP_RAT *sum;
const MP_RAT *a1;
const MP_RAT *a2;
#endif
{
MP_INT gcd1, gcd2;
MP_INT tmp1, tmp2;
mpz_init (&gcd1);
mpz_init (&gcd2);
mpz_init (&tmp1);
mpz_init (&tmp2);
/* SUM might be identical to either operand, so don't store the
result there until we are finished with the input operands. We
dare to overwrite the numerator of SUM when we are finished
with the numerators of A1 and A2. */
mpz_gcd (&gcd1, &(a1->den), &(a2->den));
if (gcd1.size > 1 || gcd1.d[0] != 1)
{
MP_INT t;
mpz_init (&t);
mpz_div (&tmp1, &(a2->den), &gcd1);
mpz_mul (&tmp1, &(a1->num), &tmp1);
mpz_div (&tmp2, &(a1->den), &gcd1);
mpz_mul (&tmp2, &(a2->num), &tmp2);
mpz_add (&t, &tmp1, &tmp2);
mpz_gcd (&gcd2, &t, &gcd1);
mpz_div (&(sum->num), &t, &gcd2);
mpz_div (&tmp1, &(a1->den), &gcd1);
mpz_div (&tmp2, &(a2->den), &gcd2);
mpz_mul (&(sum->den), &tmp1, &tmp2);
mpz_clear (&t);
}
else
{
/* The common divisior is 1. This is the case (for random input) with
probability 6/(pi**2). */
mpz_mul (&tmp1, &(a1->num), &(a2->den));
mpz_mul (&tmp2, &(a2->num), &(a1->den));
mpz_add (&(sum->num), &tmp1, &tmp2);
mpz_mul (&(sum->den), &(a1->den), &(a2->den));
}
mpz_clear (&tmp2);
mpz_clear (&tmp1);
mpz_clear (&gcd2);
mpz_clear (&gcd1);
}

34
gnu/lib/libgmp/mpq_clear.c Normal file
View file

@ -0,0 +1,34 @@
/* mpq_clear -- free the space occupied by a MP_RAT.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_clear (MP_RAT *m)
#else
mpq_clear (m)
MP_RAT *m;
#endif
{
(*_mp_free_func) (m->num.d, m->num.alloc * BYTES_PER_MP_LIMB);
(*_mp_free_func) (m->den.d, m->den.alloc * BYTES_PER_MP_LIMB);
}

76
gnu/lib/libgmp/mpq_cmp.c Normal file
View file

@ -0,0 +1,76 @@
/* mpq_cmp(u,v) -- Compare U, V. Return positive, zero, or negative
based on if U > V, U == V, or U < V.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
int
#ifdef __STDC__
mpq_cmp (const MP_RAT *op1, const MP_RAT *op2)
#else
mpq_cmp (op1, op2)
const MP_RAT *op1;
const MP_RAT *op2;
#endif
{
mp_size num1_size = op1->num.size;
mp_size den1_size = op1->den.size;
mp_size num2_size = op2->num.size;
mp_size den2_size = op2->den.size;
mp_size tmp1_size, tmp2_size;
mp_ptr tmp1_ptr, tmp2_ptr;
mp_size num1_sign;
int cc;
if (num1_size == 0)
return -num2_size;
if (num2_size == 0)
return num1_size;
if ((num1_size ^ num2_size) < 0) /* I.e. are the signs different? */
return num1_size;
num1_sign = num1_size;
num1_size = ABS (num1_size);
num2_size = ABS (num2_size);
tmp1_size = num1_size + den2_size;
tmp2_size = num2_size + den1_size;
if (tmp1_size != tmp2_size)
return (tmp1_size - tmp2_size) ^ num1_sign;
tmp1_ptr = (mp_ptr) alloca (tmp1_size * BYTES_PER_MP_LIMB);
tmp2_ptr = (mp_ptr) alloca (tmp2_size * BYTES_PER_MP_LIMB);
tmp1_size = (num1_size >= den2_size)
? mpn_mul (tmp1_ptr, op1->num.d, num1_size, op2->den.d, den2_size)
: mpn_mul (tmp1_ptr, op2->den.d, den2_size, op1->num.d, num1_size);
tmp2_size = (num2_size >= den1_size)
? mpn_mul (tmp2_ptr, op2->num.d, num2_size, op1->den.d, den1_size)
: mpn_mul (tmp2_ptr, op1->den.d, den1_size, op2->num.d, num2_size);
cc = tmp1_size - tmp2_size != 0
? tmp1_size - tmp2_size : mpn_cmp (tmp1_ptr, tmp2_ptr, tmp1_size);
alloca (0);
return (num1_sign < 0) ? -cc : cc;
}

92
gnu/lib/libgmp/mpq_div.c Normal file
View file

@ -0,0 +1,92 @@
/* mpq_div -- divide two rational numbers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_div (MP_RAT *quot, const MP_RAT *dividend, const MP_RAT *divisor)
#else
mpq_div (quot, dividend, divisor)
MP_RAT *quot;
const MP_RAT *dividend;
const MP_RAT *divisor;
#endif
{
MP_INT gcd1, gcd2;
MP_INT tmp1, tmp2;
MP_INT numtmp;
mpz_init (&gcd1);
mpz_init (&gcd2);
mpz_init (&tmp1);
mpz_init (&tmp2);
mpz_init (&numtmp);
/* QUOT might be identical to either operand, so don't store the
result there until we are finished with the input operands. We
dare to overwrite the numerator of QUOT when we are finished
with the numerators of DIVIDEND and DIVISOR. */
mpz_gcd (&gcd1, &(dividend->num), &(divisor->num));
mpz_gcd (&gcd2, &(divisor->den), &(dividend->den));
if (gcd1.size > 1 || gcd1.d[0] != 1)
mpz_div (&tmp1, &(dividend->num), &gcd1);
else
mpz_set (&tmp1, &(dividend->num));
if (gcd2.size > 1 || gcd2.d[0] != 1)
mpz_div (&tmp2, &(divisor->den), &gcd2);
else
mpz_set (&tmp2, &(divisor->den));
mpz_mul (&numtmp, &tmp1, &tmp2);
if (gcd1.size > 1 || gcd1.d[0] != 1)
mpz_div (&tmp1, &(divisor->num), &gcd1);
else
mpz_set (&tmp1, &(divisor->num));
if (gcd2.size > 1 || gcd2.d[0] != 1)
mpz_div (&tmp2, &(dividend->den), &gcd2);
else
mpz_set (&tmp2, &(dividend->den));
mpz_mul (&(quot->den), &tmp1, &tmp2);
/* We needed to go via NUMTMP to take care of QUOT being the same
as either input operands. Now move NUMTMP to QUOT->NUM. */
mpz_set (&(quot->num), &numtmp);
/* Keep the denominator positive. */
if (quot->den.size < 0)
{
quot->den.size = -quot->den.size;
quot->num.size = -quot->num.size;
}
mpz_clear (&numtmp);
mpz_clear (&tmp2);
mpz_clear (&tmp1);
mpz_clear (&gcd2);
mpz_clear (&gcd1);
}

40
gnu/lib/libgmp/mpq_get_den.c Normal file
View file

@ -0,0 +1,40 @@
/* mpq_get_den(den,rat_src) -- Set DEN to the denominator of RAT_SRC.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_get_den (MP_INT *den, const MP_RAT *src)
#else
mpq_get_den (den, src)
MP_INT *den;
const MP_RAT *src;
#endif
{
mp_size size = src->den.size;
if (den->alloc < size)
_mpz_realloc (den, size);
MPN_COPY (den->d, src->den.d, size);
den->size = size;
}

41
gnu/lib/libgmp/mpq_get_num.c Normal file
View file

@ -0,0 +1,41 @@
/* mpq_get_num(num,rat_src) -- Set NUM to the numerator of RAT_SRC.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_get_num (MP_INT *num, const MP_RAT *src)
#else
mpq_get_num (num, src)
MP_INT *num;
const MP_RAT *src;
#endif
{
mp_size size = src->num.size;
mp_size abs_size = ABS (size);
if (num->alloc < abs_size)
_mpz_realloc (num, abs_size);
MPN_COPY (num->d, src->num.d, abs_size);
num->size = size;
}

39
gnu/lib/libgmp/mpq_init.c Normal file
View file

@ -0,0 +1,39 @@
/* mpq_init -- Make a new rational number with value 0/1.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_init (MP_RAT *x)
#else
mpq_init (x)
MP_RAT *x;
#endif
{
x->num.alloc = 1;
x->num.d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->num.alloc);
x->num.size = 0;
x->den.alloc = 1;
x->den.d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->den.alloc);
x->den.d[0] = 1;
x->den.size = 1;
}

74
gnu/lib/libgmp/mpq_inv.c Normal file
View file

@ -0,0 +1,74 @@
/* mpq_inv(dest,src) -- invert a rational number, i.e. set DEST to SRC
with the numerator and denominator swapped.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_inv (MP_RAT *dest, const MP_RAT *src)
#else
mpq_inv (dest, src)
MP_RAT *dest;
const MP_RAT *src;
#endif
{
mp_size num_size = src->num.size;
mp_size den_size = src->den.size;
if (num_size == 0)
num_size = 1 / num_size; /* Divide by zero! */
if (num_size < 0)
{
num_size = -num_size;
den_size = -den_size;
}
dest->den.size = num_size;
dest->num.size = den_size;
/* If dest == src we may just swap the numerator and denominator, but
we have to ensure the new denominator is positive. */
if (dest == src)
{
mp_size alloc = dest->num.alloc;
mp_ptr limb_ptr = dest->num.d;
dest->num.alloc = dest->den.alloc;
dest->num.d = dest->den.d;
dest->den.alloc = alloc;
dest->den.d = limb_ptr;
}
else
{
den_size = ABS (den_size);
if (dest->num.alloc < den_size)
_mpz_realloc (&(dest->num), den_size);
if (dest->den.alloc < num_size)
_mpz_realloc (&(dest->den), num_size);
MPN_COPY (dest->num.d, src->den.d, den_size);
MPN_COPY (dest->den.d, src->num.d, num_size);
}
}

78
gnu/lib/libgmp/mpq_mul.c Normal file
View file

@ -0,0 +1,78 @@
/* mpq_mul -- mutiply two rational numbers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_mul (MP_RAT *prod, const MP_RAT *m1, const MP_RAT *m2)
#else
mpq_mul (prod, m1, m2)
MP_RAT *prod;
const MP_RAT *m1;
const MP_RAT *m2;
#endif
{
MP_INT gcd1, gcd2;
MP_INT tmp1, tmp2;
mpz_init (&gcd1);
mpz_init (&gcd2);
mpz_init (&tmp1);
mpz_init (&tmp2);
/* PROD might be identical to either operand, so don't store the
result there until we are finished with the input operands. We
dare to overwrite the numerator of PROD when we are finished
with the numerators of M1 and M1. */
mpz_gcd (&gcd1, &(m1->num), &(m2->den));
mpz_gcd (&gcd2, &(m2->num), &(m1->den));
if (gcd1.size > 1 || gcd1.d[0] != 1)
mpz_div (&tmp1, &(m1->num), &gcd1);
else
mpz_set (&tmp1, &(m1->num));
if (gcd2.size > 1 || gcd2.d[0] != 1)
mpz_div (&tmp2, &(m2->num), &gcd2);
else
mpz_set (&tmp2, &(m2->num));
mpz_mul (&(prod->num), &tmp1, &tmp2);
if (gcd1.size > 1 || gcd1.d[0] != 1)
mpz_div (&tmp1, &(m2->den), &gcd1);
else
mpz_set (&tmp1, &(m2->den));
if (gcd2.size > 1 || gcd2.d[0] != 1)
mpz_div (&tmp2, &(m1->den), &gcd2);
else
mpz_set (&tmp2, &(m1->den));
mpz_mul (&(prod->den), &tmp1, &tmp2);
mpz_clear (&tmp2);
mpz_clear (&tmp1);
mpz_clear (&gcd2);
mpz_clear (&gcd1);
}

35
gnu/lib/libgmp/mpq_neg.c Normal file
View file

@ -0,0 +1,35 @@
/* mpq_neg(dst, src) -- Assign the negated value of SRC to DST.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_neg (MP_RAT *dst, const MP_RAT *src)
#else
mpq_neg (dst, src)
MP_RAT *dst;
const MP_RAT *src;
#endif
{
mpz_neg (&dst->num, &src->num);
mpz_set (&dst->den, &src->den);
}

48
gnu/lib/libgmp/mpq_set.c Normal file
View file

@ -0,0 +1,48 @@
/* mpq_set(dest,src) -- Set DEST to SRC.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_set (MP_RAT *dest, const MP_RAT *src)
#else
mpq_set (dest, src)
MP_RAT *dest;
const MP_RAT *src;
#endif
{
mp_size num_size, den_size;
mp_size abs_num_size;
num_size = src->num.size;
abs_num_size = ABS (num_size);
if (dest->num.alloc < abs_num_size)
_mpz_realloc (&(dest->num), abs_num_size);
MPN_COPY (dest->num.d, src->num.d, abs_num_size);
dest->num.size = num_size;
den_size = src->den.size;
if (dest->den.alloc < den_size)
_mpz_realloc (&(dest->den), den_size);
MPN_COPY (dest->den.d, src->den.d, den_size);
dest->den.size = den_size;
}

46
gnu/lib/libgmp/mpq_set_den.c Normal file
View file

@ -0,0 +1,46 @@
/* mpq_set_den(dest,den) -- Set the denominator of DEST from DEN.
If DEN < 0 change the sign of the numerator of DEST.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_set_den (MP_RAT *dest, const MP_INT *den)
#else
mpq_set_den (dest, den)
MP_RAT *dest;
const MP_INT *den;
#endif
{
mp_size size = den->size;
mp_size abs_size = ABS (size);
if (dest->den.alloc < abs_size)
_mpz_realloc (&(dest->den), abs_size);
MPN_COPY (dest->den.d, den->d, abs_size);
dest->den.size = abs_size;
/* The denominator is always positive; move the sign to the numerator. */
if (size < 0)
dest->num.size = -dest->num.size;
}

41
gnu/lib/libgmp/mpq_set_num.c Normal file
View file

@ -0,0 +1,41 @@
/* mpq_set_num(dest,num) -- Set the numerator of DEST from NUM.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_set_num (MP_RAT *dest, const MP_INT *num)
#else
mpq_set_num (dest, num)
MP_RAT *dest;
const MP_INT *num;
#endif
{
mp_size size = num->size;
mp_size abs_size = ABS (size);
if (dest->num.alloc < abs_size)
_mpz_realloc (&(dest->num), abs_size);
MPN_COPY (dest->num.d, num->d, abs_size);
dest->num.size = size;
}

76
gnu/lib/libgmp/mpq_set_si.c Normal file
View file

@ -0,0 +1,76 @@
/* mpq_set_si(dest,ulong_num,ulong_den) -- Set DEST to the retional number
ULONG_NUM/ULONG_DEN.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
static unsigned long int
gcd (x, y)
unsigned long int x, y;
{
for (;;)
{
x = x % y;
if (x == 0)
return y;
y = y % x;
if (y == 0)
return x;
}
}
void
#ifdef __STDC__
mpq_set_si (MP_RAT *dest, signed long int num, unsigned long int den)
#else
mpq_set_si (dest, num, den)
MP_RAT *dest;
signed long int num;
unsigned long int den;
#endif
{
unsigned long int g;
unsigned long int abs_num;
abs_num = ABS (num);
if (num == 0)
{
/* Canonicalize 0/d to 0/1. */
den = 1;
dest->num.size = 0;
}
else
{
/* Remove any common factor in NUM and DEN. */
/* Pass DEN as the second argument to gcd, in order to make the
gcd function divide by zero if DEN is zero. */
g = gcd (abs_num, den);
abs_num /= g;
den /= g;
dest->num.d[0] = abs_num;
dest->num.size = num > 0 ? 1 : -1;
}
dest->den.d[0] = den;
dest->den.size = 1;
}

73
gnu/lib/libgmp/mpq_set_ui.c Normal file
View file

@ -0,0 +1,73 @@
/* mpq_set_ui(dest,ulong_num,ulong_den) -- Set DEST to the retional number
ULONG_NUM/ULONG_DEN.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
static unsigned long int
gcd (x, y)
unsigned long int x, y;
{
for (;;)
{
x = x % y;
if (x == 0)
return y;
y = y % x;
if (y == 0)
return x;
}
}
void
#ifdef __STDC__
mpq_set_ui (MP_RAT *dest, unsigned long int num, unsigned long int den)
#else
mpq_set_ui (dest, num, den)
MP_RAT *dest;
unsigned long int num;
unsigned long int den;
#endif
{
unsigned long int g;
if (num == 0)
{
/* Canonicalize 0/n to 0/1. */
den = 1;
dest->num.size = 0;
}
else
{
/* Remove any common factor in NUM and DEN. */
/* Pass DEN as the second argument to gcd, in order to make the
gcd function divide by zero if DEN is zero. */
g = gcd (num, den);
num /= g;
den /= g;
dest->num.d[0] = num;
dest->num.size = 1;
}
dest->den.d[0] = den;
dest->den.size = 1;
}

85
gnu/lib/libgmp/mpq_sub.c Normal file
View file

@ -0,0 +1,85 @@
/* mpq_sub -- subtract two rational numbers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpq_sub (MP_RAT *dif, const MP_RAT *min, const MP_RAT *sub)
#else
mpq_sub (dif, min, sub)
MP_RAT *dif;
const MP_RAT *min;
const MP_RAT *sub;
#endif
{
MP_INT gcd1, gcd2;
MP_INT tmp1, tmp2;
mpz_init (&gcd1);
mpz_init (&gcd2);
mpz_init (&tmp1);
mpz_init (&tmp2);
/* DIF might be identical to either operand, so don't store the
result there until we are finished with the input operands. We
dare to overwrite the numerator of DIF when we are finished
with the numerators of MIN and SUB. */
mpz_gcd (&gcd1, &(min->den), &(sub->den));
if (gcd1.size > 1 || gcd1.d[0] != 1)
{
MP_INT t;
mpz_init (&t);
mpz_div (&tmp1, &(sub->den), &gcd1);
mpz_mul (&tmp1, &(min->num), &tmp1);
mpz_div (&tmp2, &(min->den), &gcd1);
mpz_mul (&tmp2, &(sub->num), &tmp2);
mpz_sub (&t, &tmp1, &tmp2);
mpz_gcd (&gcd2, &t, &gcd1);
mpz_div (&(dif->num), &t, &gcd2);
mpz_div (&tmp1, &(min->den), &gcd1);
mpz_div (&tmp2, &(sub->den), &gcd2);
mpz_mul (&(dif->den), &tmp1, &tmp2);
mpz_clear (&t);
}
else
{
/* The common divisior is 1. This is the case (for random input) with
probability 6/(pi**2). */
mpz_mul (&tmp1, &(min->num), &(sub->den));
mpz_mul (&tmp2, &(sub->num), &(min->den));
mpz_sub (&(dif->num), &tmp1, &tmp2);
mpz_mul (&(dif->den), &(min->den), &(sub->den));
}
mpz_clear (&tmp2);
mpz_clear (&tmp1);
mpz_clear (&gcd2);
mpz_clear (&gcd1);
}

44
gnu/lib/libgmp/mpz_abs.c Normal file
View file

@ -0,0 +1,44 @@
/* mpz_abs(MP_INT *dst, MP_INT *src) -- Assign the absolute value of SRC to DST.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_abs (MP_INT *dst, const MP_INT *src)
#else
mpz_abs (dst, src)
MP_INT *dst;
const MP_INT *src;
#endif
{
mp_size src_size = ABS (src->size);
if (src != dst)
{
if (dst->alloc < src_size)
_mpz_realloc (dst, src_size);
MPN_COPY (dst->d, src->d, src_size);
}
dst->size = src_size;
}

121
gnu/lib/libgmp/mpz_add.c Normal file
View file

@ -0,0 +1,121 @@
/* mpz_add -- Add two integers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#ifndef BERKELEY_MP
void
#ifdef __STDC__
mpz_add (MP_INT *sum, const MP_INT *u, const MP_INT *v)
#else
mpz_add (sum, u, v)
MP_INT *sum;
const MP_INT *u;
const MP_INT *v;
#endif
#else /* BERKELEY_MP */
void
#ifdef __STDC__
madd (const MP_INT *u, const MP_INT *v, MP_INT *sum)
#else
madd (u, v, sum)
const MP_INT *u;
const MP_INT *v;
MP_INT *sum;
#endif
#endif /* BERKELEY_MP */
{
mp_srcptr up, vp;
mp_ptr sump;
mp_size usize, vsize, sumsize;
mp_size abs_usize;
mp_size abs_vsize;
usize = u->size;
vsize = v->size;
abs_usize = ABS (usize);
abs_vsize = ABS (vsize);
if (abs_usize < abs_vsize)
{
/* Swap U and V. */
{const MP_INT *t = u; u = v; v = t;}
{mp_size t = usize; usize = vsize; vsize = t;}
{mp_size t = abs_usize; abs_usize = abs_vsize; abs_vsize = t;}
}
/* True: abs(USIZE) >= abs(VSIZE) */
/* If not space for sum (and possible carry), increase space. */
sumsize = abs_usize + 1;
if (sum->alloc < sumsize)
_mpz_realloc (sum, sumsize);
/* These must be after realloc (u or v may be the same as sum). */
up = u->d;
vp = v->d;
sump = sum->d;
if (usize >= 0)
{
if (vsize >= 0)
{
sumsize = mpn_add (sump, up, abs_usize, vp, abs_vsize);
if (sumsize != 0)
sump[abs_usize] = 1;
sumsize = sumsize + abs_usize;
}
else
{
/* The signs are different. Need exact comparision to determine
which operand to subtract from which. */
if (abs_usize == abs_vsize && mpn_cmp (up, vp, abs_usize) < 0)
sumsize = -(abs_usize
+ mpn_sub (sump, vp, abs_usize, up, abs_usize));
else
sumsize = (abs_usize
+ mpn_sub (sump, up, abs_usize, vp, abs_vsize));
}
}
else
{
if (vsize >= 0)
{
/* The signs are different. Need exact comparision to determine
which operand to subtract from which. */
if (abs_usize == abs_vsize && mpn_cmp (up, vp, abs_usize) < 0)
sumsize = (abs_usize
+ mpn_sub (sump, vp, abs_usize, up, abs_usize));
else
sumsize = -(abs_usize
+ mpn_sub (sump, up, abs_usize, vp, abs_vsize));
}
else
{
sumsize = mpn_add (sump, up, abs_usize, vp, abs_vsize);
if (sumsize != 0)
sump[abs_usize] = 1;
sumsize = -(sumsize + abs_usize);
}
}
sum->size = sumsize;
}

84
gnu/lib/libgmp/mpz_add_ui.c Normal file
View file

@ -0,0 +1,84 @@
/* mpz_add_ui -- Add an MP_INT and an unsigned one-word integer.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_add_ui (MP_INT *sum, const MP_INT *add1, mp_limb add2)
#else
mpz_add_ui (sum, add1, add2)
MP_INT *sum;
const MP_INT *add1;
mp_limb add2;
#endif
{
mp_srcptr add1p;
mp_ptr sump;
mp_size add1size, sumsize;
mp_size abs_add1size;
add1size = add1->size;
abs_add1size = ABS (add1size);
/* If not space for SUM (and possible carry), increase space. */
sumsize = abs_add1size + 1;
if (sum->alloc < sumsize)
_mpz_realloc (sum, sumsize);
/* These must be after realloc (ADD1 may be the same as SUM). */
add1p = add1->d;
sump = sum->d;
if (add2 == 0)
{
MPN_COPY (sump, add1p, abs_add1size);
sum->size = add1size;
return;
}
if (abs_add1size == 0)
{
sump[0] = add2;
sum->size = 1;
return;
}
if (add1size >= 0)
{
sumsize = mpn_add (sump, add1p, abs_add1size, &add2, 1);
if (sumsize != 0)
sump[abs_add1size] = 1;
sumsize = sumsize + abs_add1size;
}
else
{
/* The signs are different. Need exact comparision to determine
which operand to subtract from which. */
if (abs_add1size == 1 && add1p[0] < add2)
sumsize = (abs_add1size
+ mpn_sub (sump, &add2, 1, add1p, 1));
else
sumsize = -(abs_add1size
+ mpn_sub (sump, add1p, abs_add1size, &add2, 1));
}
sum->size = sumsize;
}

267
gnu/lib/libgmp/mpz_and.c Normal file
View file

@ -0,0 +1,267 @@
/* mpz_and -- Logical and.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#define min(l,o) ((l) < (o) ? (l) : (o))
#define max(h,i) ((h) > (i) ? (h) : (i))
void
#ifdef __STDC__
mpz_and (MP_INT *res, const MP_INT *op1, const MP_INT *op2)
#else
mpz_and (res, op1, op2)
MP_INT *res;
const MP_INT *op1;
const MP_INT *op2;
#endif
{
mp_srcptr op1_ptr, op2_ptr;
mp_size op1_size, op2_size;
mp_ptr res_ptr;
mp_size res_size;
mp_size i;
op1_size = op1->size;
op2_size = op2->size;
op1_ptr = op1->d;
op2_ptr = op2->d;
res_ptr = res->d;
if (op1_size >= 0)
{
if (op2_size >= 0)
{
res_size = min (op1_size, op2_size);
/* First loop finds the size of the result. */
for (i = res_size - 1; i >= 0; i--)
if ((op1_ptr[i] & op2_ptr[i]) != 0)
break;
res_size = i + 1;
/* Handle allocation, now when we know exactly how much space is
needed for the result. */
if (res->alloc < res_size)
{
_mpz_realloc (res, res_size);
op1_ptr = op1->d;
op2_ptr = op2->d;
res_ptr = res->d;
}
/* Second loop computes the real result. */
for (i = res_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] & op2_ptr[i];
res->size = res_size;
return;
}
else /* op2_size < 0 */
/* Fall through to the code at the end of the function. */
;
}
else
{
if (op2_size < 0)
{
mp_ptr opx;
mp_limb cy;
mp_limb one = 1;
mp_size res_alloc;
/* Both operands are negative, so will be the result.
-((-OP1) & (-OP2)) = -(~(OP1 - 1) & ~(OP2 - 1)) =
= ~(~(OP1 - 1) & ~(OP2 - 1)) + 1 =
= ((OP1 - 1) | (OP2 - 1)) + 1 */
op1_size = -op1_size;
op2_size = -op2_size;
res_alloc = 1 + max (op1_size, op2_size);
opx = (mp_ptr) alloca (op1_size * BYTES_PER_MP_LIMB);
op1_size += mpn_sub (opx, op1_ptr, op1_size, &one, 1);
op1_ptr = opx;
opx = (mp_ptr) alloca (op2_size * BYTES_PER_MP_LIMB);
op2_size += mpn_sub (opx, op2_ptr, op2_size, &one, 1);
op2_ptr = opx;
if (res->alloc < res_alloc)
{
_mpz_realloc (res, res_alloc);
res_ptr = res->d;
/* Don't re-read OP1_PTR and OP2_PTR. They point to
temporary space--never to the space RES->D used
to point to before reallocation. */
}
if (op1_size >= op2_size)
{
MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size,
op1_size - op2_size);
for (i = op2_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] | op2_ptr[i];
res_size = op1_size;
}
else
{
MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size,
op2_size - op1_size);
for (i = op1_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] | op2_ptr[i];
res_size = op2_size;
}
if (res_size != 0)
{
cy = mpn_add (res_ptr, res_ptr, res_size, &one, 1);
if (cy)
{
res_ptr[res_size] = cy;
res_size++;
}
}
else
{
res_ptr[0] = 1;
res_size = 1;
}
res->size = -res_size;
return;
}
else
{
/* We should compute -OP1 & OP2. Swap OP1 and OP2 and fall
through to the code that handles OP1 & -OP2. */
{const MP_INT *t = op1; op1 = op2; op2 = t;}
{mp_srcptr t = op1_ptr; op1_ptr = op2_ptr; op2_ptr = t;}
{mp_size t = op1_size; op1_size = op2_size; op2_size = t;}
}
}
{
#if 0
mp_size op2_lim;
/* OP2 must be negated as with infinite precision.
Scan from the low end for a non-zero limb. The first non-zero
limb is simply negated (two's complement). Any subsequent
limbs are one's complemented. Of course, we don't need to
handle more limbs than there are limbs in the other, positive
operand as the result for those limbs is going to become zero
anyway. */
/* Scan for the least significant. non-zero OP2 limb, and zero the
result meanwhile for those limb positions. (We will surely
find a non-zero limb, so we can write the loop with one
termination condition only.) */
for (i = 0; op2_ptr[i] == 0; i++)
res_ptr[i] = 0;
op2_lim = i;
op2_size = -op2_size;
if (op1_size <= op2_size)
{
/* The ones-extended OP2 is >= than the zero-extended OP1.
RES_SIZE <= OP1_SIZE. Find the exact size. */
for (i = op1_size - 1; i > op2_lim; i--)
if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
break;
res_size = i + 1;
}
else
{
/* The ones-extended OP2 is < than the zero-extended OP1.
RES_SIZE == OP1_SIZE, since OP1 is normalized. */
res_size = op1_size;
}
#endif
/* OP1 is positive and zero-extended,
OP2 is negative and ones-extended.
The result will be positive.
OP1 & -OP2 = OP1 & ~(OP2 - 1). */
mp_ptr opx;
const mp_limb one = 1;
op2_size = -op2_size;
opx = (mp_ptr) alloca (op2_size * BYTES_PER_MP_LIMB);
op2_size += mpn_sub (opx, op2_ptr, op2_size, &one, 1);
op2_ptr = opx;
if (op1_size > op2_size)
{
/* The result has the same size as OP1, since OP1 is normalized
and longer than the ones-extended OP2. */
res_size = op1_size;
/* Handle allocation, now when we know exactly how much space is
needed for the result. */
if (res->alloc < res_size)
{
_mpz_realloc (res, res_size);
res_ptr = res->d;
op1_ptr = op1->d;
/* Don't re-read OP2_PTR. It points to temporary space--never
to the space RES->D used to point to before reallocation. */
}
MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size,
res_size - op2_size);
for (i = op2_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
res->size = res_size;
}
else
{
/* Find out the exact result size. Ignore the high limbs of OP2,
OP1 is zero-extended and would make the result zero. */
for (i = op1_size - 1; i >= 0; i--)
if ((op1_ptr[i] & ~op2_ptr[i]) != 0)
break;
res_size = i + 1;
/* Handle allocation, now when we know exactly how much space is
needed for the result. */
if (res->alloc < res_size)
{
_mpz_realloc (res, res_size);
res_ptr = res->d;
op1_ptr = op1->d;
/* Don't re-read OP2_PTR. It points to temporary space--never
to the space RES->D used to point to before reallocation. */
}
for (i = res_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] & ~op2_ptr[i];
res->size = res_size;
}
}
}

34
gnu/lib/libgmp/mpz_clear.c Normal file
View file

@ -0,0 +1,34 @@
/* mpz_clear -- de-allocate the space occupied by the dynamic digit space of
an integer.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_clear (MP_INT *m)
#else
mpz_clear (m)
MP_INT *m;
#endif
{
(*_mp_free_func) (m->d, m->alloc * BYTES_PER_MP_LIMB);
}

124
gnu/lib/libgmp/mpz_clrbit.c Normal file
View file

@ -0,0 +1,124 @@
/* mpz_clrbit -- clear a specified bit.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#define MPN_NORMALIZE(p, size) \
do { \
mp_size i; \
for (i = (size) - 1; i >= 0; i--) \
if ((p)[i] != 0) \
break; \
(size) = i + 1; \
} while (0)
void
#ifdef __STDC__
mpz_clrbit (MP_INT *d, unsigned long int bit_index)
#else
mpz_clrbit (d, bit_index)
MP_INT *d;
unsigned long int bit_index;
#endif
{
mp_size dsize = d->size;
mp_ptr dp = d->d;
mp_size limb_index;
limb_index = bit_index / BITS_PER_MP_LIMB;
if (dsize >= 0)
{
if (limb_index < dsize)
{
dp[limb_index] &= ~((mp_limb) 1 << (bit_index % BITS_PER_MP_LIMB));
MPN_NORMALIZE (dp, dsize);
d->size = dsize;
}
else
;
}
else
{
mp_size zero_bound;
/* Simulate two's complement arithmetic, i.e. simulate
1. Set OP = ~(OP - 1) [with infinitely many leading ones].
2. clear the bit.
3. Set OP = ~OP + 1. */
dsize = -dsize;
/* No upper bound on this loop, we're sure there's a non-zero limb
sooner ot later. */
for (zero_bound = 0; ; zero_bound++)
if (dp[zero_bound] != 0)
break;
if (limb_index > zero_bound)
{
if (limb_index < dsize)
{
dp[limb_index] |= ((mp_limb) 1 << (bit_index % BITS_PER_MP_LIMB));
}
else
{
/* Ugh. The bit should be cleared outside of the end of the
number. We have to increase the size of the number. */
if (d->alloc < limb_index + 1)
{
_mpz_realloc (d, limb_index + 1);
dp = d->d;
}
MPN_ZERO (dp + dsize, limb_index - dsize);
dp[limb_index] = ((mp_limb) 1 << (bit_index % BITS_PER_MP_LIMB));
d->size = -(limb_index + 1);
}
}
else if (limb_index == zero_bound)
{
dp[limb_index] = ((dp[limb_index] - 1)
| ((mp_limb) 1 << (bit_index % BITS_PER_MP_LIMB))) + 1;
if (dp[limb_index] == 0)
{
mp_size i;
for (i = limb_index + 1; i < dsize; i++)
{
dp[i] += 1;
if (dp[i] != 0)
goto fin;
}
/* We got carry all way out beyond the end of D. Increase
its size (and allocation if necessary). */
dsize++;
if (d->alloc < dsize)
{
_mpz_realloc (d, dsize);
dp = d->d;
}
dp[i] = 1;
d->size = -dsize;
fin:;
}
}
else
;
}
}

84
gnu/lib/libgmp/mpz_cmp.c Normal file
View file

@ -0,0 +1,84 @@
/* mpz_cmp(u,v) -- Compare U, V. Return postive, zero, or negative
based on if U > V, U == V, or U < V.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#ifdef BERKELEY_MP
#include "mp.h"
#endif
#include "gmp.h"
#include "gmp-impl.h"
#ifndef BERKELEY_MP
int
#ifdef __STDC__
mpz_cmp (const MP_INT *u, const MP_INT *v)
#else
mpz_cmp (u, v)
const MP_INT *u;
const MP_INT *v;
#endif
#else /* BERKELEY_MP */
int
#ifdef __STDC__
mcmp (const MP_INT *u, const MP_INT *v)
#else
mcmp (u, v)
const MP_INT *u;
const MP_INT *v;
#endif
#endif /* BERKELEY_MP */
{
mp_size usize = u->size;
mp_size vsize = v->size;
mp_size size;
mp_size i;
mp_limb a, b;
mp_srcptr up, vp;
if (usize != vsize)
return usize - vsize;
if (usize == 0)
return 0;
size = ABS (usize);
up = u->d;
vp = v->d;
i = size - 1;
do
{
a = up[i];
b = vp[i];
i--;
if (i < 0)
break;
}
while (a == b);
if (a == b)
return 0;
if ((a < b) == (usize < 0))
return 1;
else
return -1;
}

62
gnu/lib/libgmp/mpz_cmp_si.c Normal file
View file

@ -0,0 +1,62 @@
/* mpz_cmp_si(u,v) -- Compare an integer U with a single-word int V.
Return positive, zero, or negative based on if U > V, U == V, or U < V.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
int
#ifdef __STDC__
mpz_cmp_si (const MP_INT *u, signed long int v_digit)
#else
mpz_cmp_si (u, v_digit)
const MP_INT *u;
signed long int v_digit;
#endif
{
mp_size usize = u->size;
mp_size vsize;
mp_limb u_digit;
vsize = 0;
if (v_digit > 0)
vsize = 1;
else if (v_digit < 0)
{
vsize = -1;
v_digit = -v_digit;
}
if (usize != vsize)
return usize - vsize;
if (usize == 0)
return 0;
u_digit = u->d[0];
if (u_digit == v_digit)
return 0;
if ((u_digit < v_digit) == (usize < 0))
return 1;
else
return -1;
}

52
gnu/lib/libgmp/mpz_cmp_ui.c Normal file
View file

@ -0,0 +1,52 @@
/* mpz_cmp_ui.c -- Compare a MP_INT a with an mp_limb b. Return positive,
zero, or negative based on if a > b, a == b, or a < b.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
int
#ifdef __STDC__
mpz_cmp_ui (const MP_INT *u, mp_limb v_digit)
#else
mpz_cmp_ui (u, v_digit)
const MP_INT *u;
mp_limb v_digit;
#endif
{
mp_size usize = u->size;
if (usize == 0)
return -(v_digit != 0);
if (usize == 1)
{
mp_limb u_digit;
u_digit = u->d[0];
if (u_digit > v_digit)
return 1;
if (u_digit < v_digit)
return -1;
return 0;
}
return (usize > 0) ? 1 : -1;
}

96
gnu/lib/libgmp/mpz_com.c Normal file
View file

@ -0,0 +1,96 @@
/* mpz_com(MP_INT *dst, MP_INT *src) -- Assign the bit-complemented value of
SRC to DST.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_com (MP_INT *dst, const MP_INT *src)
#else
mpz_com (dst, src)
MP_INT *dst;
const MP_INT *src;
#endif
{
mp_size size = src->size;
mp_srcptr src_ptr;
mp_ptr dst_ptr;
if (size >= 0)
{
/* As with infinite precision: one's complement, two's complement.
But this can be simplified using the identity -x = ~x + 1.
So we're going to compute (~~x) + 1 = x + 1! */
if (dst->alloc < size + 1)
_mpz_realloc (dst, size + 1);
src_ptr = src->d;
dst_ptr = dst->d;
if (size == 0)
{
/* Special case, as mpn_add wants the first arg's size >= the
second arg's size. */
dst_ptr[0] = 1;
dst->size = -1;
return;
}
{
mp_limb one = 1;
int cy;
cy = mpn_add (dst_ptr, src_ptr, size, &one, 1);
if (cy)
{
dst_ptr[size] = cy;
size++;
}
}
/* Store a negative size, to indicate ones-extension. */
dst->size = -size;
}
else
{
/* As with infinite precision: two's complement, then one's complement.
But that can be simplified using the identity -x = ~(x - 1).
So we're going to compute ~~(x - 1) = x - 1! */
size = -size;
if (dst->alloc < size)
_mpz_realloc (dst, size);
src_ptr = src->d;
dst_ptr = dst->d;
{
mp_limb one = 1;
size += mpn_sub (dst_ptr, src_ptr, size, &one, 1);
}
/* Store a positive size, to indicate zero-extension. */
dst->size = size;
}
}

117
gnu/lib/libgmp/mpz_div.c Normal file
View file

@ -0,0 +1,117 @@
/* mpz_div -- divide two integers and produce a quotient.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_div (MP_INT *quot, const MP_INT *num, const MP_INT *den)
#else
mpz_div (quot, num, den)
MP_INT *quot;
const MP_INT *num;
const MP_INT *den;
#endif
{
mp_srcptr np, dp;
mp_ptr qp, rp;
mp_size nsize = num->size;
mp_size dsize = den->size;
mp_size qsize, rsize;
mp_size sign_quotient = nsize ^ dsize;
unsigned normalization_steps;
nsize = ABS (nsize);
dsize = ABS (dsize);
/* Ensure space is enough for quotient. */
qsize = nsize - dsize + 1; /* qsize cannot be bigger than this. */
if (qsize <= 0)
{
quot->size = 0;
return;
}
if (quot->alloc < qsize)
_mpz_realloc (quot, qsize);
qp = quot->d;
np = num->d;
dp = den->d;
rp = (mp_ptr) alloca ((nsize + 1) * BYTES_PER_MP_LIMB);
count_leading_zeros (normalization_steps, dp[dsize - 1]);
/* Normalize the denominator and the numerator. */
if (normalization_steps != 0)
{
mp_ptr tp;
mp_limb ndigit;
/* Shift up the denominator setting the most significant bit of
the most significant word. Use temporary storage not to clobber
the original contents of the denominator. */
tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
(void) mpn_lshift (tp, dp, dsize, normalization_steps);
dp = tp;
/* Shift up the numerator, possibly introducing a new most
significant word. Move the shifted numerator in the remainder
meanwhile. */
ndigit = mpn_lshift (rp, np, nsize, normalization_steps);
if (ndigit != 0)
{
rp[nsize] = ndigit;
rsize = nsize + 1;
}
else
rsize = nsize;
}
else
{
/* The denominator is already normalized, as required.
Copy it to temporary space if it overlaps with the quotient. */
if (dp == qp)
{
dp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
MPN_COPY ((mp_ptr) dp, qp, dsize);
}
/* Move the numerator to the remainder. */
MPN_COPY (rp, np, nsize);
rsize = nsize;
}
qsize = rsize - dsize + mpn_div (qp, rp, rsize, dp, dsize);
/* Normalize the quotient. We may have at most one leading
zero-word, so no loop is needed. */
if (qsize > 0)
qsize -= (qp[qsize - 1] == 0);
if (sign_quotient < 0)
qsize = -qsize;
quot->size = qsize;
alloca (0);
}

53
gnu/lib/libgmp/mpz_div_2exp.c Normal file
View file

@ -0,0 +1,53 @@
/* mpz_div_2exp -- Divide a bignum by 2**CNT
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_div_2exp (MP_INT *w, const MP_INT *u, unsigned long int cnt)
#else
mpz_div_2exp (w, u, cnt)
MP_INT *w;
const MP_INT *u;
unsigned long int cnt;
#endif
{
mp_size usize = u->size;
mp_size wsize;
mp_size abs_usize = ABS (usize);
mp_size limb_cnt;
limb_cnt = cnt / BITS_PER_MP_LIMB;
wsize = abs_usize - limb_cnt;
if (wsize <= 0)
wsize = 0;
else
{
if (w->alloc < wsize)
_mpz_realloc (w, wsize);
wsize = mpn_rshift (w->d, u->d + limb_cnt, abs_usize - limb_cnt,
cnt % BITS_PER_MP_LIMB);
}
w->size = (usize >= 0) ? wsize : -wsize;
}

65
gnu/lib/libgmp/mpz_div_ui.c Normal file
View file

@ -0,0 +1,65 @@
/* mpz_div_ui(quot, dividend, divisor_limb)
-- Divide DIVIDEND by DIVISOR_LIMB and store the result in QUOT.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_div_ui (MP_INT *quot, const MP_INT *dividend, unsigned long int divisor_limb)
#else
mpz_div_ui (quot, dividend, divisor_limb)
MP_INT *quot;
const MP_INT *dividend;
unsigned long int divisor_limb;
#endif
{
mp_size sign_dividend;
mp_size dividend_size, quot_size;
mp_ptr dividend_ptr, quot_ptr;
sign_dividend = dividend->size;
dividend_size = ABS (dividend->size);
if (dividend_size == 0)
{
quot->size = 0;
return;
}
/* No need for temporary allocation and copying if QUOT == DIVIDEND as
the divisor is just one limb, and thus no intermediate remainders
need to be stored. */
if (quot->alloc < dividend_size)
_mpz_realloc (quot, dividend_size);
quot_ptr = quot->d;
dividend_ptr = dividend->d;
mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb);
/* The quotient is DIVIDEND_SIZE limbs, but the most significant
might be zero. Set QUOT_SIZE properly. */
quot_size = dividend_size - (quot_ptr[dividend_size - 1] == 0);
quot->size = sign_dividend >= 0 ? quot_size : -quot_size;
}

38
gnu/lib/libgmp/mpz_dm.c Normal file
View file

@ -0,0 +1,38 @@
/* mpz_divmod(quot,rem,dividend,divisor) -- Set QUOT to DIVIDEND/DIVISOR,
and REM to DIVIDEND mod DIVISOR.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_divmod (MP_INT *quot, MP_INT *rem, const MP_INT *num, const MP_INT *den)
#else
mpz_divmod (quot, rem, num, den)
MP_INT *quot;
MP_INT *rem;
const MP_INT *num;
const MP_INT *den;
#endif
#define COMPUTE_QUOTIENT
#include "mpz_dmincl.c"

81
gnu/lib/libgmp/mpz_dm_ui.c Normal file
View file

@ -0,0 +1,81 @@
/* mpz_divmod_ui(quot,rem,dividend,short_divisor) --
Set QUOT to DIVIDEND / SHORT_DIVISOR
and REM to DIVIDEND mod SHORT_DIVISOR.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_divmod_ui (MP_INT *quot, MP_INT *rem,
const MP_INT *dividend, unsigned long int divisor_limb)
#else
mpz_divmod_ui (quot, rem, dividend, divisor_limb)
MP_INT *quot;
MP_INT *rem;
const MP_INT *dividend;
unsigned long int divisor_limb;
#endif
{
mp_size sign_dividend;
mp_size dividend_size, quot_size;
mp_ptr dividend_ptr, quot_ptr;
mp_limb remainder_limb;
sign_dividend = dividend->size;
dividend_size = ABS (dividend->size);
if (dividend_size == 0)
{
quot->size = 0;
rem->size = 0;
return;
}
/* No need for temporary allocation and copying if QUOT == DIVIDEND as
the divisor is just one limb, and thus no intermediate remainders
need to be stored. */
if (quot->alloc < dividend_size)
_mpz_realloc (quot, dividend_size);
quot_ptr = quot->d;
dividend_ptr = dividend->d;
remainder_limb = mpn_divmod_1 (quot_ptr,
dividend_ptr, dividend_size, divisor_limb);
if (remainder_limb == 0)
rem->size = 0;
else
{
/* Store the single-limb remainder. We don't check if there's space
for just one limb, since no function ever makes zero space. */
rem->size = sign_dividend >= 0 ? 1 : -1;
rem->d[0] = remainder_limb;
}
/* The quotient is DIVIDEND_SIZE limbs, but the most significant
might be zero. Set QUOT_SIZE properly. */
quot_size = dividend_size - (quot_ptr[dividend_size - 1] == 0);
quot->size = sign_dividend >= 0 ? quot_size : -quot_size;
}

172
gnu/lib/libgmp/mpz_dmincl.c Normal file
View file

@ -0,0 +1,172 @@
/* mpz_dmincl.c -- include file for mpz_dm.c, mpz_mod.c, mdiv.c.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
/* THIS CODE IS OBSOLETE. IT WILL SOON BE REPLACED BY CLEANER CODE WITH
LESS MEMORY ALLOCATION OVERHEAD. */
/* If den == quot, den needs temporary storage.
If den == rem, den needs temporary storage.
If num == quot, num needs temporary storage.
If den has temporary storage, it can be normalized while being copied,
i.e no extra storage should be allocated. */
/* This is the function body of mdiv, mpz_divmod, and mpz_mod.
If COMPUTE_QUOTIENT is defined, the quotient is put in the MP_INT
object quot, otherwise that variable is not referenced at all.
The remainder is always computed, and the result is put in the MP_INT
object rem. */
{
mp_ptr np, dp;
mp_ptr qp, rp;
mp_size nsize = num->size;
mp_size dsize = den->size;
mp_size qsize, rsize;
mp_size sign_remainder = nsize;
#ifdef COMPUTE_QUOTIENT
mp_size sign_quotient = nsize ^ dsize;
#endif
unsigned normalization_steps;
nsize = ABS (nsize);
dsize = ABS (dsize);
/* Ensure space is enough for quotient and remainder. */
/* We need space for an extra limb in the remainder, because it's
up-shifted (normalized) below. */
rsize = nsize + 1;
if (rem->alloc < rsize)
_mpz_realloc (rem, rsize);
qsize = nsize - dsize + 1; /* qsize cannot be bigger than this. */
if (qsize <= 0)
{
#ifdef COMPUTE_QUOTIENT
quot->size = 0;
#endif
if (num != rem)
{
rem->size = num->size;
MPN_COPY (rem->d, num->d, nsize);
}
return;
}
#ifdef COMPUTE_QUOTIENT
if (quot->alloc < qsize)
_mpz_realloc (quot, qsize);
qp = quot->d;
#else
qp = (mp_ptr) alloca (qsize * BYTES_PER_MP_LIMB);
#endif
np = num->d;
dp = den->d;
rp = rem->d;
/* Make sure quot and num are different. Otherwise the numerator
would be successively overwritten by the quotient digits. */
if (qp == np)
{
np = (mp_ptr) alloca (nsize * BYTES_PER_MP_LIMB);
MPN_COPY (np, qp, nsize);
}
count_leading_zeros (normalization_steps, dp[dsize - 1]);
/* Normalize the denominator, i.e. make its most significant bit set by
shifting it NORMALIZATION_STEPS bits to the left. Also shift the
numerator the same number of steps (to keep the quotient the same!). */
if (normalization_steps != 0)
{
mp_ptr tp;
mp_limb ndigit;
/* Shift up the denominator setting the most significant bit of
the most significant word. Use temporary storage not to clobber
the original contents of the denominator. */
tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
(void) mpn_lshift (tp, dp, dsize, normalization_steps);
dp = tp;
/* Shift up the numerator, possibly introducing a new most
significant word. Move the shifted numerator in the remainder
meanwhile. */
ndigit = mpn_lshift (rp, np, nsize, normalization_steps);
if (ndigit != 0)
{
rp[nsize] = ndigit;
rsize = nsize + 1;
}
else
rsize = nsize;
}
else
{
#ifdef COMPUTE_QUOTIENT
if (rem == den || quot == den)
#else
if (rem == den)
#endif
{
mp_ptr tp;
tp = (mp_ptr) alloca (dsize * BYTES_PER_MP_LIMB);
MPN_COPY (tp, dp, dsize);
dp = tp;
}
/* Move the numerator to the remainder. */
if (rp != np)
MPN_COPY (rp, np, nsize);
rsize = nsize;
}
qsize = rsize - dsize + mpn_div (qp, rp, rsize, dp, dsize);
rsize = dsize;
/* Normalize the remainder. */
while (rsize > 0)
{
if (rp[rsize - 1] != 0)
break;
rsize--;
}
if (normalization_steps != 0)
rsize = mpn_rshift (rp, rp, rsize, normalization_steps);
rem->size = (sign_remainder >= 0) ? rsize : -rsize;
#ifdef COMPUTE_QUOTIENT
/* Normalize the quotient. We may have at most one leading
zero-word, so no loop is needed. */
if (qsize > 0)
qsize -= (qp[qsize - 1] == 0);
quot->size = (sign_quotient >= 0) ? qsize : -qsize;
#endif
alloca (0);
}

156
gnu/lib/libgmp/mpz_fac_ui.c Normal file
View file

@ -0,0 +1,156 @@
/* mpz_fac_ui(result, n) -- Set RESULT to N!.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#ifdef DBG
#include <stdio.h>
#endif
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_fac_ui (MP_INT *result, unsigned long int n)
#else
mpz_fac_ui (result, n)
MP_INT *result;
unsigned long int n;
#endif
{
#if SIMPLE_FAC
/* Be silly. Just multiply the numbers in ascending order. O(n**2). */
mp_limb k;
mpz_set_ui (result, (mp_limb) 1);
for (k = 2; k <= n; k++)
mpz_mul_ui (result, result, k);
#else
/* Be smarter. Multiply groups of numbers in ascending order until the
product doesn't fit in a limb. Multiply these partial products in a
balanced binary tree fashion, to make the operand have as equal sizes
as possible. (When the operands have about the same size, mpn_mul
becomes faster.) */
mp_limb k;
mp_limb p1, p0, p;
/* Stack of partial products, used to make the computation balanced
(i.e. make the sizes of the multiplication operands equal). The
topmost position of MP_STACK will contain a one-limb partial product,
the second topmost will contain a two-limb partial product, and so
on. MP_STACK[0] will contain a partial product with 2**t limbs.
To compute n! MP_STACK needs to be less than
log(n)**2/log(BITS_PER_MP_LIMB), so 30 is surely enough. */
#define MP_STACK_SIZE 30
MP_INT mp_stack[MP_STACK_SIZE];
/* TOP is an index into MP_STACK, giving the topmost element.
TOP_LIMIT_SO_FAR is the largets value it has taken so far. */
int top, top_limit_so_far;
/* Count of the total number of limbs put on MP_STACK so far. This
variable plays an essential role in making the compututation balanced.
See below. */
unsigned int tree_cnt;
top = top_limit_so_far = -1;
tree_cnt = 0;
p = 1;
for (k = 2; k <= n; k++)
{
/* Multiply the partial product in P with K. */
umul_ppmm (p1, p0, p, k);
/* Did we get overflow into the high limb, i.e. is the partial
product now more than one limb? */
if (p1 != 0)
{
tree_cnt++;
if (tree_cnt % 2 == 0)
{
mp_size i;
/* TREE_CNT is even (i.e. we have generated an even number of
one-limb partial products), which means that we have a
single-limb product on the top of MP_STACK. */
mpz_mul_ui (&mp_stack[top], &mp_stack[top], p);
/* If TREE_CNT is divisable by 4, 8,..., we have two
similar-sized partial products with 2, 4,... limbs at
the topmost two positions of MP_STACK. Multiply them
to form a new partial product with 4, 8,... limbs. */
for (i = 4; (tree_cnt & (i - 1)) == 0; i <<= 1)
{
mpz_mul (&mp_stack[top - 1],
&mp_stack[top], &mp_stack[top - 1]);
top--;
}
}
else
{
/* Put the single-limb partial product in P on the stack.
(The next time we get a single-limb product, we will
multiply the two together.) */
top++;
if (top > top_limit_so_far)
{
if (top > MP_STACK_SIZE)
abort();
/* The stack is now bigger than ever, initialize the top
element. */
mpz_init_set_ui (&mp_stack[top], p);
top_limit_so_far++;
}
else
mpz_set_ui (&mp_stack[top], p);
}
/* We ignored the last result from umul_ppmm. Put K in P as the
first component of the next single-limb partial product. */
p = k;
}
else
/* We didn't get overflow in umul_ppmm. Put p0 in P and try
with one more value of K. */
p = p0;
}
/* We have partial products in mp_stack[0..top], in descending order.
We also have a small partial product in p.
Their product is the final result. */
if (top < 0)
mpz_set_ui (result, p);
else
mpz_mul_ui (result, &mp_stack[top--], p);
while (top >= 0)
mpz_mul (result, result, &mp_stack[top--]);
/* Free the storage allocated for MP_STACK. */
for (top = top_limit_so_far; top >= 0; top--)
mpz_clear (&mp_stack[top]);
#endif
}

169
gnu/lib/libgmp/mpz_gcd.c Normal file
View file

@ -0,0 +1,169 @@
/* mpz_gcd -- Calculate the greatest common divisior of two integers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#ifndef BERKELEY_MP
void
#ifdef __STDC__
mpz_gcd (MP_INT *w, const MP_INT *u, const MP_INT *v)
#else
mpz_gcd (w, u, v)
MP_INT *w;
const MP_INT *u;
const MP_INT *v;
#endif
#else /* BERKELEY_MP */
void
#ifdef __STDC__
gcd (const MP_INT *u, const MP_INT *v, MP_INT *w)
#else
gcd (u, v, w)
const MP_INT *u;
const MP_INT *v;
MP_INT *w;
#endif
#endif /* BERKELEY_MP */
{
mp_size usize, vsize, wsize;
mp_ptr up_in, vp_in;
mp_ptr up, vp;
mp_ptr wp;
mp_size i;
mp_limb d;
int bcnt;
mp_size w_bcnt;
mp_limb cy_digit;
usize = ABS (u->size);
vsize = ABS (v->size);
/* GCD(0,v) == v. */
if (usize == 0)
{
if (w->alloc < vsize)
_mpz_realloc (w, vsize);
w->size = vsize;
MPN_COPY (w->d, v->d, vsize);
return;
}
/* GCD(0,u) == u. */
if (vsize == 0)
{
if (w->alloc < usize)
_mpz_realloc (w, usize);
w->size = usize;
MPN_COPY (w->d, u->d, usize);
return;
}
/* Make U odd by shifting it down as many bit positions as there
are zero bits. Put the result in temporary space. */
up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
up_in = u->d;
for (i = 0; (d = up_in[i]) == 0; i++)
;
count_leading_zeros (bcnt, d & -d);
bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
usize = mpn_rshift (up, up_in + i, usize - i, bcnt);
bcnt += i * BITS_PER_MP_LIMB;
w_bcnt = bcnt;
/* Make V odd by shifting it down as many bit positions as there
are zero bits. Put the result in temporary space. */
vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
vp_in = v->d;
for (i = 0; (d = vp_in[i]) == 0; i++)
;
count_leading_zeros (bcnt, d & -d);
bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
vsize = mpn_rshift (vp, vp_in + i, vsize - i, bcnt);
/* W_BCNT is set to the minimum of the number of zero bits in U and V.
Thus it represents the number of common 2 factors. */
bcnt += i * BITS_PER_MP_LIMB;
if (bcnt < w_bcnt)
w_bcnt = bcnt;
for (;;)
{
int cmp;
cmp = usize - vsize != 0 ? usize - vsize : mpn_cmp (up, vp, usize);
/* If U and V have become equal, we have found the GCD. */
if (cmp == 0)
break;
if (cmp > 0)
{
/* Replace U by (U - V) >> cnt, with cnt being the least value
making U odd again. */
usize += mpn_sub (up, up, usize, vp, vsize);
for (i = 0; (d = up[i]) == 0; i++)
;
count_leading_zeros (bcnt, d & -d);
bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
usize = mpn_rshift (up, up + i, usize - i, bcnt);
}
else
{
/* Replace V by (V - U) >> cnt, with cnt being the least value
making V odd again. */
vsize += mpn_sub (vp, vp, vsize, up, usize);
for (i = 0; (d = vp[i]) == 0; i++)
;
count_leading_zeros (bcnt, d & -d);
bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
vsize = mpn_rshift (vp, vp + i, vsize - i, bcnt);
}
}
/* GCD(U_IN, V_IN) now is U * 2**W_BCNT. */
wsize = usize + w_bcnt / BITS_PER_MP_LIMB + 1;
if (w->alloc < wsize)
_mpz_realloc (w, wsize);
wp = w->d;
MPN_ZERO (wp, w_bcnt / BITS_PER_MP_LIMB);
cy_digit = mpn_lshift (wp + w_bcnt / BITS_PER_MP_LIMB, up, usize,
w_bcnt % BITS_PER_MP_LIMB);
wsize = usize + w_bcnt / BITS_PER_MP_LIMB;
if (cy_digit != 0)
{
wp[wsize] = cy_digit;
wsize++;
}
w->size = wsize;
alloca (0);
}

80
gnu/lib/libgmp/mpz_gcdext.c Normal file
View file

@ -0,0 +1,80 @@
/* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that
g = as + bt.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
/* Botch: SLOW! */
void
#ifdef __STDC__
mpz_gcdext (MP_INT *g, MP_INT *s, MP_INT *t, const MP_INT *a, const MP_INT *b)
#else
mpz_gcdext (g, s, t, a, b)
MP_INT *g;
MP_INT *s;
MP_INT *t;
const MP_INT *a;
const MP_INT *b;
#endif
{
MP_INT s0, s1, q, r, x, d0, d1;
mpz_init_set_ui (&s0, 1);
mpz_init_set_ui (&s1, 0);
mpz_init (&q);
mpz_init (&r);
mpz_init (&x);
mpz_init_set (&d0, a);
mpz_init_set (&d1, b);
while (d1.size != 0)
{
mpz_divmod (&q, &r, &d0, &d1);
mpz_set (&d0, &d1);
mpz_set (&d1, &r);
mpz_mul (&x, &s1, &q);
mpz_sub (&x, &s0, &x);
mpz_set (&s0, &s1);
mpz_set (&s1, &x);
}
if (t != NULL)
{
mpz_mul (&x, &s0, a);
mpz_sub (&x, &d0, &x);
if (b->size == 0)
t->size = 0;
else
mpz_div (t, &x, b);
}
mpz_set (s, &s0);
mpz_set (g, &d0);
mpz_clear (&s0);
mpz_clear (&s1);
mpz_clear (&q);
mpz_clear (&r);
mpz_clear (&x);
mpz_clear (&d0);
mpz_clear (&d1);
}

40
gnu/lib/libgmp/mpz_get_si.c Normal file
View file

@ -0,0 +1,40 @@
/* mpz_get_si(integer) -- Return the least significant digit from INTEGER.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
signed long int
#ifdef __STDC__
mpz_get_si (const MP_INT *integer)
#else
mpz_get_si (integer)
const MP_INT *integer;
#endif
{
mp_size size = integer->size;
if (size > 0)
return integer->d[0] % ((mp_limb) 1 << (BITS_PER_MP_LIMB - 1));
else if (size < 0)
return -(integer->d[0] % ((mp_limb) 1 << (BITS_PER_MP_LIMB - 1)));
else
return 0;
}

39
gnu/lib/libgmp/mpz_get_str.c Normal file
View file

@ -0,0 +1,39 @@
/* mpz_get_str (string, base, mp_src) -- Convert the multiple precision
number MP_SRC to a string STRING of base BASE. If STRING is NULL
allocate space for the result. In any case, return a pointer to the
result. If STRING is not NULL, the caller must ensure enough space is
available to store the result.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
char *
#ifdef __STDC__
mpz_get_str (char *str, int base, const MP_INT *m)
#else
mpz_get_str (str, base, m)
char *str;
int base;
const MP_INT *m;
#endif
{
return _mpz_get_str (str, base, m);
}

36
gnu/lib/libgmp/mpz_get_ui.c Normal file
View file

@ -0,0 +1,36 @@
/* mpz_get_ui(integer) -- Return the least significant digit from INTEGER.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
unsigned long int
#ifdef __STDC__
mpz_get_ui (const MP_INT *integer)
#else
mpz_get_ui (integer)
const MP_INT *integer;
#endif
{
if (integer->size == 0)
return 0;
else
return integer->d[0];
}

35
gnu/lib/libgmp/mpz_init.c Normal file
View file

@ -0,0 +1,35 @@
/* mpz_init() -- Make a new multiple precision number with value 0.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_init (MP_INT *x)
#else
mpz_init (x)
MP_INT *x;
#endif
{
x->alloc = 1;
x->d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->alloc);
x->size = 0;
}

72
gnu/lib/libgmp/mpz_inp_raw.c Normal file
View file

@ -0,0 +1,72 @@
/* mpz_inp_raw -- Input a MP_INT in raw, but endianess, and wordsize
independent format (as output by mpz_out_raw).
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_inp_raw (MP_INT *x, FILE *file)
#else
mpz_inp_raw (x, file)
MP_INT *x;
FILE *file;
#endif
{
int i;
mp_size s;
mp_size xsize;
mp_ptr xp;
unsigned int c;
mp_limb x_digit;
mp_size x_index;
xsize = 0;
for (i = 4 - 1; i >= 0; i--)
{
c = fgetc (file);
xsize = (xsize << BITS_PER_CHAR) | c;
}
/* ??? Sign extend xsize for non-32 bit machines? */
x_index = (ABS (xsize) + BYTES_PER_MP_LIMB - 1) / BYTES_PER_MP_LIMB - 1;
if (x->alloc < x_index)
_mpz_realloc (x, x_index);
xp = x->d;
x->size = xsize / BYTES_PER_MP_LIMB;
x_digit = 0;
for (s = ABS (xsize) - 1; s >= 0; s--)
{
i = s % BYTES_PER_MP_LIMB;
c = fgetc (file);
x_digit = (x_digit << BITS_PER_CHAR) | c;
if (i == 0)
{
xp[x_index--] = x_digit;
x_digit = 0;
}
}
}

105
gnu/lib/libgmp/mpz_inp_str.c Normal file
View file

@ -0,0 +1,105 @@
/* mpz_inp_str(dest_integer, stream, base) -- Input a number in base
BASE from stdio stream STREAM and store the result in DEST_INTEGER.
Copyright (C) 1991, 1993 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include <ctype.h>
#include "gmp.h"
#include "gmp-impl.h"
static int
char_ok_for_base (c, base)
int c;
int base;
{
if (isdigit (c))
return (unsigned) c - '0' < base;
if (islower (c))
return (unsigned) c - 'a' + 10 < base;
if (isupper (c))
return (unsigned) c - 'A' + 10 < base;
return 0;
}
void
#ifdef __STDC__
mpz_inp_str (MP_INT *dest, FILE *stream, int base)
#else
mpz_inp_str (dest, stream, base)
MP_INT *dest;
FILE *stream;
int base;
#endif
{
char *str;
size_t str_size;
size_t i;
int c;
int negative = 0;
str_size = 100;
str = (char *) (*_mp_allocate_func) (str_size);
c = getc (stream);
if (c == '-')
{
negative = 1;
c = getc (stream);
}
/* If BASE is 0, try to find out the base by looking at the initial
characters. */
if (base == 0)
{
base = 10;
if (c == '0')
{
base = 8;
c = getc (stream);
if (c == 'x' || c == 'X')
{
base = 16;
c = getc (stream);
}
}
}
for (i = 0; char_ok_for_base (c, base); i++)
{
if (i >= str_size)
{
size_t old_str_size = str_size;
str_size = str_size * 3 / 2;
str = (char *) (*_mp_reallocate_func) (str, old_str_size, str_size);
}
str[i] = c;
c = getc (stream);
}
ungetc (c, stream);
str[i] = 0;
_mpz_set_str (dest, str, base);
if (negative)
dest->size = -dest->size;
(*_mp_free_func) (str, str_size);
}

242
gnu/lib/libgmp/mpz_ior.c Normal file
View file

@ -0,0 +1,242 @@
/* mpz_ior -- Logical inclusive or.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#define min(l,o) ((l) < (o) ? (l) : (o))
#define max(h,i) ((h) > (i) ? (h) : (i))
void
#ifdef __STDC__
mpz_ior (MP_INT *res, const MP_INT *op1, const MP_INT *op2)
#else
mpz_ior (res, op1, op2)
MP_INT *res;
const MP_INT *op1;
const MP_INT *op2;
#endif
{
mp_srcptr op1_ptr, op2_ptr;
mp_size op1_size, op2_size;
mp_ptr res_ptr;
mp_size res_size;
mp_size i;
op1_size = op1->size;
op2_size = op2->size;
op1_ptr = op1->d;
op2_ptr = op2->d;
res_ptr = res->d;
if (op1_size >= 0)
{
if (op2_size >= 0)
{
if (op1_size >= op2_size)
{
if (res->alloc < op1_size)
{
_mpz_realloc (res, op1_size);
op1_ptr = op1->d;
op2_ptr = op2->d;
res_ptr = res->d;
}
if (res_ptr != op1_ptr)
MPN_COPY (res_ptr + op2_size, op1_ptr + op2_size,
op1_size - op2_size);
for (i = op2_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] | op2_ptr[i];
res_size = op1_size;
}
else
{
if (res->alloc < op2_size)
{
_mpz_realloc (res, op2_size);
op1_ptr = op1->d;
op2_ptr = op2->d;
res_ptr = res->d;
}
if (res_ptr != op2_ptr)
MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size,
op2_size - op1_size);
for (i = op1_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] | op2_ptr[i];
res_size = op2_size;
}
res->size = res_size;
return;
}
else /* op2_size < 0 */
/* Fall through to the code at the end of the function. */
;
}
else
{
if (op2_size < 0)
{
mp_ptr opx;
mp_limb cy;
mp_limb one = 1;
/* Both operands are negative, so will be the result.
-((-OP1) | (-OP2)) = -(~(OP1 - 1) | ~(OP2 - 1)) =
= ~(~(OP1 - 1) | ~(OP2 - 1)) + 1 =
= ((OP1 - 1) & (OP2 - 1)) + 1 */
op1_size = -op1_size;
op2_size = -op2_size;
res_size = min (op1_size, op2_size);
/* Possible optimization: Decrease mpn_sub precision,
as we won't use the entire res of both. */
opx = (mp_ptr) alloca (op1_size * BYTES_PER_MP_LIMB);
op1_size += mpn_sub (opx, op1_ptr, op1_size, &one, 1);
op1_ptr = opx;
opx = (mp_ptr) alloca (op2_size * BYTES_PER_MP_LIMB);
op2_size += mpn_sub (opx, op2_ptr, op2_size, &one, 1);
op2_ptr = opx;
if (res->alloc < res_size)
{
_mpz_realloc (res, res_size);
res_ptr = res->d;
/* Don't re-read OP1_PTR and OP2_PTR. They point to
temporary space--never to the space RES->D used
to point to before reallocation. */
}
/* First loop finds the size of the result. */
for (i = res_size - 1; i >= 0; i--)
if ((op1_ptr[i] & op2_ptr[i]) != 0)
break;
res_size = i + 1;
/* Second loop computes the real result. */
for (i = res_size - 1; i >= 0; i--)
res_ptr[i] = op1_ptr[i] & op2_ptr[i];
if (res_size != 0)
{
cy = mpn_add (res_ptr, res_ptr, res_size, &one, 1);
if (cy)
{
res_ptr[res_size] = cy;
res_size++;
}
}
else
{
res_ptr[0] = 1;
res_size = 1;
}
res->size = -res_size;
return;
}
else
{
/* We should compute -OP1 | OP2. Swap OP1 and OP2 and fall
through to the code that handles OP1 | -OP2. */
{const MP_INT *t = op1; op1 = op2; op2 = t;}
{mp_srcptr t = op1_ptr; op1_ptr = op2_ptr; op2_ptr = t;}
{mp_size t = op1_size; op1_size = op2_size; op2_size = t;}
}
}
{
mp_ptr opx;
mp_limb cy;
mp_limb one = 1;
mp_size res_alloc;
/* Operand 2 negative, so will be the result.
-(OP1 | (-OP2)) = -(OP1 | ~(OP2 - 1)) =
= ~(OP1 | ~(OP2 - 1)) + 1 =
= (~OP1 & (OP2 - 1)) + 1 */
op2_size = -op2_size;
res_alloc = op2_size;
opx = (mp_ptr) alloca (op2_size * BYTES_PER_MP_LIMB);
op2_size += mpn_sub (opx, op2_ptr, op2_size, &one, 1);
op2_ptr = opx;
if (res->alloc < res_alloc)
{
_mpz_realloc (res, res_alloc);
op1_ptr = op1->d;
res_ptr = res->d;
/* Don't re-read OP2_PTR. It points to temporary space--never
to the space RES->D used to point to before reallocation. */
}
if (op1_size >= op2_size)
{
/* We can just ignore the part of OP1 that stretches above OP2,
because the result limbs are zero there. */
/* First loop finds the size of the result. */
for (i = op2_size - 1; i >= 0; i--)
if ((~op1_ptr[i] & op2_ptr[i]) != 0)
break;
res_size = i + 1;
}
else
{
res_size = op2_size;
/* Copy the part of OP2 that stretches above OP1, to RES. */
MPN_COPY (res_ptr + op1_size, op2_ptr + op1_size,
op2_size - op1_size);
}
/* Second loop computes the real result. */
for (i = res_size - 1; i >= 0; i--)
res_ptr[i] = ~op1_ptr[i] & op2_ptr[i];
if (res_size != 0)
{
cy = mpn_add (res_ptr, res_ptr, res_size, &one, 1);
if (cy)
{
res_ptr[res_size] = cy;
res_size++;
}
}
else
{
res_ptr[0] = 1;
res_size = 1;
}
res->size = -res_size;
alloca (0);
return;
}
}

45
gnu/lib/libgmp/mpz_iset.c Normal file
View file

@ -0,0 +1,45 @@
/* mpz_init_set (src_integer) -- Make a new multiple precision number with
a value copied from SRC_INTEGER.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_init_set (MP_INT *x, const MP_INT *src)
#else
mpz_init_set (x, src)
MP_INT *x;
const MP_INT *src;
#endif
{
mp_size size;
mp_size abs_size;
size = src->size;
abs_size = ABS (size);
x->alloc = abs_size == 0 ? 1 : abs_size;
x->d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->alloc);
MPN_COPY (x->d, src->d, abs_size);
x->size = size;
}

48
gnu/lib/libgmp/mpz_iset_si.c Normal file
View file

@ -0,0 +1,48 @@
/* mpz_init_set_si(val) -- Make a new multiple precision number with
value val.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_init_set_si (MP_INT *x, signed long int val)
#else
mpz_init_set_si (x, val)
MP_INT *x;
signed long int val;
#endif
{
x->alloc = 1;
x->d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->alloc);
if (val > 0)
{
x->d[0] = val;
x->size = 1;
}
else if (val < 0)
{
x->d[0] = -val;
x->size = -1;
}
else
x->size = 0;
}

42
gnu/lib/libgmp/mpz_iset_str.c Normal file
View file

@ -0,0 +1,42 @@
/* mpz_init_set_str(mpz, string, base) -- Initialize MPZ and set it to the
value in the \0-terminated ascii string STRING in base BASE. Return 0 if
the string was accepted, -1 if an error occured. If BASE == 0 determine
the base in the C standard way, i.e. 0xhh...h means base 16, 0oo...o
means base 8, otherwise assume base 10.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
int
#ifdef __STDC__
mpz_init_set_str (MP_INT *x, const char *str, int base)
#else
mpz_init_set_str (x, str, base)
MP_INT *x;
const char *str;
int base;
#endif
{
x->alloc = 1;
x->d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->alloc);
return _mpz_set_str (x, str, base);
}

43
gnu/lib/libgmp/mpz_iset_ui.c Normal file
View file

@ -0,0 +1,43 @@
/* mpz_init_set_ui(val) -- Make a new multiple precision number with
value val.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_init_set_ui (MP_INT *x, unsigned long int val)
#else
mpz_init_set_ui (x, val)
MP_INT *x;
unsigned long int val;
#endif
{
x->alloc = 1;
x->d = (mp_ptr) (*_mp_allocate_func) (BYTES_PER_MP_LIMB * x->alloc);
if (val > 0)
{
x->d[0] = val;
x->size = 1;
}
else
x->size = 0;
}

52
gnu/lib/libgmp/mpz_mdiv.c Normal file
View file

@ -0,0 +1,52 @@
/* mpz_mdiv -- Mathematical DIVision and MODulo, i.e. division that rounds
the quotient towards -infinity.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mdiv (MP_INT *quot,
const MP_INT *dividend, const MP_INT *divisor)
#else
mpz_mdiv (quot, dividend, divisor)
MP_INT *quot;
const MP_INT *dividend;
const MP_INT *divisor;
#endif
{
if ((dividend->size ^ divisor->size) >= 0)
{
/* When the dividend and the divisor has same sign, this function
gives same result as mpz_div. */
mpz_div (quot, dividend, divisor);
}
else
{
MP_INT rem;
MPZ_TMP_INIT (&rem, 1 + ABS (dividend->size));
mpz_divmod (quot, &rem, dividend, divisor);
if (rem.size != 0)
mpz_sub_ui (quot, quot, 1);
}
}

43
gnu/lib/libgmp/mpz_mdiv_ui.c Normal file
View file

@ -0,0 +1,43 @@
/* mpz_mdiv_ui -- Mathematical DIVision and MODulo, i.e. division that rounds
the quotient towards -infinity.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mdiv_ui (MP_INT *quot,
const MP_INT *dividend, unsigned long int divisor)
#else
mpz_mdiv_ui (quot, dividend, divisor)
MP_INT *quot;
const MP_INT *dividend;
unsigned long int divisor;
#endif
{
MP_INT rem;
MPZ_TMP_INIT (&rem, 1 + ABS (dividend->size));
mpz_divmod_ui (quot, &rem, dividend, divisor);
if (rem.size < 0)
mpz_sub_ui (quot, quot, 1);
}

64
gnu/lib/libgmp/mpz_mdm.c Normal file
View file

@ -0,0 +1,64 @@
/* mpz_mdivmod -- Mathematical DIVision and MODulo, i.e. division that rounds
the quotient towards -infinity, and with the remainder non-negative.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mdivmod (MP_INT *quot, MP_INT *rem,
const MP_INT *dividend, const MP_INT *divisor)
#else
mpz_mdivmod (quot, rem, dividend, divisor)
MP_INT *quot;
MP_INT *rem;
const MP_INT *dividend;
const MP_INT *divisor;
#endif
{
if ((dividend->size ^ divisor->size) >= 0)
{
/* When the dividend and the divisor has same sign, this function
gives same result as mpz_divmod. */
mpz_divmod (quot, rem, dividend, divisor);
}
else
{
MP_INT temp_divisor; /* N.B.: lives until function returns! */
/* We need the original value of the divisor after the quotient and
remainder have been preliminary calculated. We have to copy it to
temporary space if it's the same variable as either QUOT or REM. */
if (quot == divisor || rem == divisor)
{
MPZ_TMP_INIT (&temp_divisor, ABS (divisor->size));
mpz_set (&temp_divisor, divisor);
divisor = &temp_divisor;
}
mpz_divmod (quot, rem, dividend, divisor);
if (rem->size != 0)
{
mpz_sub_ui (quot, quot, 1);
mpz_add (rem, rem, divisor);
}
}
}

58
gnu/lib/libgmp/mpz_mdm_ui.c Normal file
View file

@ -0,0 +1,58 @@
/* mpz_mdivmod -- Mathematical DIVision and MODulo, i.e. division that rounds
the quotient towards -infinity, and with the remainder non-negative.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
unsigned long int
#ifdef __STDC__
mpz_mdivmod_ui (MP_INT *quot, MP_INT *rem,
const MP_INT *dividend, unsigned long int divisor)
#else
mpz_mdivmod_ui (quot, rem, dividend, divisor)
MP_INT *quot;
MP_INT *rem;
const MP_INT *dividend;
unsigned long int divisor;
#endif
{
MP_INT temp_rem; /* N.B.: lives until function returns! */
/* If the user doesn't want the remainder to be stored in an integer
object, allocate a scratch variable for it. */
if (rem == NULL)
{
MPZ_TMP_INIT (&temp_rem, 1 + ABS (dividend->size));
rem = &temp_rem;
}
mpz_divmod_ui (quot, rem, dividend, divisor);
if (rem->size < 0)
{
mpz_sub_ui (quot, quot, 1);
mpz_add_ui (rem, rem, divisor);
}
if (rem->size == 0)
return 0;
return rem->d[0];
}

60
gnu/lib/libgmp/mpz_mmod.c Normal file
View file

@ -0,0 +1,60 @@
/* mpz_mmod -- Mathematical MODulo, i.e. with the remainder
non-negative.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mmod (MP_INT *rem,
const MP_INT *dividend, const MP_INT *divisor)
#else
mpz_mmod (rem, dividend, divisor)
MP_INT *rem;
const MP_INT *dividend;
const MP_INT *divisor;
#endif
{
if ((dividend->size ^ divisor->size) >= 0)
{
/* When the dividend and the divisor has same sign, this function
gives same result as mpz_mod. */
mpz_mod (rem, dividend, divisor);
}
else
{
MP_INT temp_divisor; /* N.B.: lives until function returns! */
/* We need the original value of the divisor after the remainder has
been preliminary calculated. We have to copy it to temporary
space if it's the same variable as REM. */
if (rem == divisor)
{
MPZ_TMP_INIT (&temp_divisor, ABS (divisor->size));
mpz_set (&temp_divisor, divisor);
divisor = &temp_divisor;
}
mpz_mod (rem, dividend, divisor);
if (rem->size != 0)
mpz_add (rem, rem, divisor);
}
}

52
gnu/lib/libgmp/mpz_mmod_ui.c Normal file
View file

@ -0,0 +1,52 @@
/* mpz_mmod -- Mathematical MODulo, i.e. with the remainder
non-negative.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
unsigned long int
#ifdef __STDC__
mpz_mmod_ui (MP_INT *rem,
const MP_INT *dividend, unsigned long int divisor)
#else
mpz_mmod_ui (rem, dividend, divisor)
MP_INT *rem;
const MP_INT *dividend;
unsigned long int divisor;
#endif
{
MP_INT temp_rem; /* N.B.: lives until function returns! */
if (rem == NULL)
{
MPZ_TMP_INIT (&temp_rem, 1 + ABS (dividend->size));
rem = &temp_rem;
}
mpz_mod_ui (rem, dividend, divisor);
if (rem->size < 0)
mpz_add_ui (rem, rem, divisor);
if (rem->size == 0)
return 0;
return rem->d[0];
}

36
gnu/lib/libgmp/mpz_mod.c Normal file
View file

@ -0,0 +1,36 @@
/* mpz_mod(rem, dividend, divisor) -- Set REM to DIVIDEND mod DIVISOR.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_mod (MP_INT *rem, const MP_INT *num, const MP_INT *den)
#else
mpz_mod (rem, num, den)
MP_INT *rem;
const MP_INT *num;
const MP_INT *den;
#endif
#undef COMPUTE_QUOTIENT
#include "mpz_dmincl.c"

82
gnu/lib/libgmp/mpz_mod_2exp.c Normal file
View file

@ -0,0 +1,82 @@
/* mpz_mod_2exp -- divide a MP_INT by 2**n and produce a remainder.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mod_2exp (MP_INT *res, const MP_INT *in, unsigned long int cnt)
#else
mpz_mod_2exp (res, in, cnt)
MP_INT *res;
const MP_INT *in;
unsigned long int cnt;
#endif
{
mp_size in_size = ABS (in->size);
mp_size res_size;
mp_size limb_cnt = cnt / BITS_PER_MP_LIMB;
mp_srcptr in_ptr = in->d;
if (in_size > limb_cnt)
{
/* The input operand is (probably) greater than 2**CNT. */
mp_limb x;
x = in_ptr[limb_cnt] & (((mp_limb) 1 << cnt % BITS_PER_MP_LIMB) - 1);
if (x != 0)
{
res_size = limb_cnt + 1;
if (res->alloc < res_size)
_mpz_realloc (res, res_size);
res->d[limb_cnt] = x;
}
else
{
mp_size i;
for (i = limb_cnt - 1; i >= 0; i--)
if (in_ptr[i] != 0)
break;
res_size = i + 1;
if (res->alloc < res_size)
_mpz_realloc (res, res_size);
limb_cnt = res_size;
}
}
else
{
/* The input operand is smaller than 2**CNT. We perform a no-op,
apart from that we might need to copy IN to RES. */
res_size = in_size;
if (res->alloc < res_size)
_mpz_realloc (res, res_size);
limb_cnt = res_size;
}
if (res != in)
MPN_COPY (res->d, in->d, limb_cnt);
res->size = (in->size >= 0) ? res_size : -res_size;
}

65
gnu/lib/libgmp/mpz_mod_ui.c Normal file
View file

@ -0,0 +1,65 @@
/* mpz_mod_ui(rem, dividend, divisor_limb)
-- Set REM to DIVDEND mod DIVISOR_LIMB.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_mod_ui (MP_INT *rem, const MP_INT *dividend,
unsigned long int divisor_limb)
#else
mpz_mod_ui (rem, dividend, divisor_limb)
MP_INT *rem;
const MP_INT *dividend;
unsigned long int divisor_limb;
#endif
{
mp_size sign_dividend;
mp_size dividend_size;
mp_limb remainder_limb;
sign_dividend = dividend->size;
dividend_size = ABS (dividend->size);
if (dividend_size == 0)
{
rem->size = 0;
return;
}
/* No need for temporary allocation and copying if QUOT == DIVIDEND as
the divisor is just one limb, and thus no intermediate remainders
need to be stored. */
remainder_limb = mpn_mod_1 (dividend->d, dividend_size, divisor_limb);
if (remainder_limb == 0)
rem->size = 0;
else
{
/* Store the single-limb remainder. We don't check if there's space
for just one limb, since no function ever makes zero space. */
rem->size = sign_dividend >= 0 ? 1 : -1;
rem->d[0] = remainder_limb;
}
}

114
gnu/lib/libgmp/mpz_mul.c Normal file
View file

@ -0,0 +1,114 @@
/* mpz_mul -- Multiply two integers.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#ifndef BERKELEY_MP
void
#ifdef __STDC__
mpz_mul (MP_INT *w, const MP_INT *u, const MP_INT *v)
#else
mpz_mul (w, u, v)
MP_INT *w;
const MP_INT *u;
const MP_INT *v;
#endif
#else /* BERKELEY_MP */
void
#ifdef __STDC__
mult (const MP_INT *u, const MP_INT *v, MP_INT *w)
#else
mult (u, v, w)
const MP_INT *u;
const MP_INT *v;
MP_INT *w;
#endif
#endif /* BERKELEY_MP */
{
mp_size usize = u->size;
mp_size vsize = v->size;
mp_size wsize;
mp_size sign_product;
mp_ptr up, vp;
mp_ptr wp;
mp_ptr free_me = NULL;
size_t free_me_size;
sign_product = usize ^ vsize;
usize = ABS (usize);
vsize = ABS (vsize);
if (usize < vsize)
{
/* Swap U and V. */
{const MP_INT *t = u; u = v; v = t;}
{mp_size t = usize; usize = vsize; vsize = t;}
}
up = u->d;
vp = v->d;
wp = w->d;
/* Ensure W has space enough to store the result. */
wsize = usize + vsize;
if (w->alloc < wsize)
{
if (wp == up || wp == vp)
{
free_me = wp;
free_me_size = w->alloc;
}
else
(*_mp_free_func) (wp, w->alloc * BYTES_PER_MP_LIMB);
w->alloc = wsize;
wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
w->d = wp;
}
else
{
/* Make U and V not overlap with W. */
if (wp == up)
{
/* W and U are identical. Allocate temporary space for U. */
up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
/* Is V identical too? Keep it identical with U. */
if (wp == vp)
vp = up;
/* Copy to the temporary space. */
MPN_COPY (up, wp, usize);
}
else if (wp == vp)
{
/* W and V are identical. Allocate temporary space for V. */
vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
/* Copy to the temporary space. */
MPN_COPY (vp, wp, vsize);
}
}
wsize = mpn_mul (wp, up, usize, vp, vsize);
w->size = sign_product < 0 ? -wsize : wsize;
if (free_me != NULL)
(*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
alloca (0);
}

68
gnu/lib/libgmp/mpz_mul_2exp.c Normal file
View file

@ -0,0 +1,68 @@
/* mpz_mul_2exp -- Multiply a bignum by 2**CNT
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_mul_2exp (MP_INT *w, const MP_INT *u, unsigned long int cnt)
#else
mpz_mul_2exp (w, u, cnt)
MP_INT *w;
const MP_INT *u;
unsigned long int cnt;
#endif
{
mp_size usize = u->size;
mp_size abs_usize = ABS (usize);
mp_size wsize;
mp_size limb_cnt;
mp_ptr wp;
mp_limb wdigit;
if (usize == 0)
{
w->size = 0;
return;
}
limb_cnt = cnt / BITS_PER_MP_LIMB;
wsize = abs_usize + limb_cnt + 1;
if (w->alloc < wsize)
_mpz_realloc (w, wsize);
wp = w->d;
wdigit = mpn_lshift (wp + limb_cnt, u->d, abs_usize,
cnt % BITS_PER_MP_LIMB);
wsize = abs_usize + limb_cnt;
if (wdigit != 0)
{
wp[wsize] = wdigit;
wsize++;
}
/* Zero all whole digits at low end. Do it here and not before calling
mpn_lshift, not to loose for U == W. */
MPN_ZERO (wp, limb_cnt);
w->size = (usize >= 0) ? wsize : -wsize;
}

78
gnu/lib/libgmp/mpz_mul_ui.c Normal file
View file

@ -0,0 +1,78 @@
/* mpz_mul_ui(product, multiplier, small_multiplicand) -- Set
PRODUCT to MULTIPLICATOR times SMALL_MULTIPLICAND.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
void
#ifdef __STDC__
mpz_mul_ui (MP_INT *prod, const MP_INT *mult,
unsigned long int small_mult)
#else
mpz_mul_ui (prod, mult, small_mult)
MP_INT *prod;
const MP_INT *mult;
unsigned long int small_mult;
#endif
{
mp_size mult_size = mult->size;
mp_size sign_product = mult_size;
mp_size i;
mp_limb cy;
mp_size prod_size;
mp_srcptr mult_ptr;
mp_ptr prod_ptr;
mult_size = ABS (mult_size);
if (mult_size == 0 || small_mult == 0)
{
prod->size = 0;
return;
}
prod_size = mult_size + 1;
if (prod->alloc < prod_size)
_mpz_realloc (prod, prod_size);
mult_ptr = mult->d;
prod_ptr = prod->d;
cy = 0;
for (i = 0; i < mult_size; i++)
{
mp_limb p1, p0;
umul_ppmm (p1, p0, small_mult, mult_ptr[i]);
p0 += cy;
cy = p1 + (p0 < cy);
prod_ptr[i] = p0;
}
prod_size = mult_size;
if (cy != 0)
{
prod_ptr[mult_size] = cy;
prod_size++;
}
prod->size = sign_product > 0 ? prod_size : -prod_size;
}

46
gnu/lib/libgmp/mpz_neg.c Normal file
View file

@ -0,0 +1,46 @@
/* mpz_neg(MP_INT *dst, MP_INT *src) -- Assign the negated value of SRC to DST.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_neg (MP_INT *dst, const MP_INT *src)
#else
mpz_neg (dst, src)
MP_INT *dst;
const MP_INT *src;
#endif
{
mp_size src_size = src->size;
if (src != dst)
{
mp_size abs_src_size = ABS (src_size);
if (dst->alloc < abs_src_size)
_mpz_realloc (dst, abs_src_size);
MPN_COPY (dst->d, src->d, abs_src_size);
}
dst->size = -src_size;
}

55
gnu/lib/libgmp/mpz_out_raw.c Normal file
View file

@ -0,0 +1,55 @@
/* mpz_out_raw -- Output a MP_INT in raw, but endianess-independent format.
Copyright (C) 1991, 1992 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_out_raw (FILE *file, const MP_INT *x)
#else
mpz_out_raw (file, x)
FILE *file;
const MP_INT *x;
#endif
{
int i;
mp_size s;
mp_size xsize = x->size;
mp_srcptr xp = x->d;
mp_size out_size = xsize * BYTES_PER_MP_LIMB;
/* Make the size 4 bytes on all machines, to make the format portable. */
for (i = 4 - 1; i >= 0; i--)
fputc ((out_size >> (i * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR), file);
/* Output from the most significant digit to the least significant digit,
with each digit also output in decreasing significance order. */
for (s = ABS (xsize) - 1; s >= 0; s--)
{
mp_limb x_digit;
x_digit = xp[s];
for (i = BYTES_PER_MP_LIMB - 1; i >= 0; i--)
fputc ((x_digit >> (i * BITS_PER_CHAR)) % (1 << BITS_PER_CHAR), file);
}
}

45
gnu/lib/libgmp/mpz_out_str.c Normal file
View file

@ -0,0 +1,45 @@
/* mpz_out_str(stream, base, integer) -- Output to STREAM the multi prec.
integer INTEGER in base BASE.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <stdio.h>
#include "gmp.h"
#include "gmp-impl.h"
void
#ifdef __STDC__
mpz_out_str (FILE *stream, int base, const MP_INT *x)
#else
mpz_out_str (stream, base, x)
FILE *stream;
int base;
const MP_INT *x;
#endif
{
char *str;
size_t str_size;
str_size = ((size_t) (ABS (x->size) * BITS_PER_MP_LIMB
* __mp_bases[ABS (base)].chars_per_bit_exactly)) + 3;
str = (char *) alloca (str_size);
_mpz_get_str (str, base, x);
fputs (str, stream);
alloca (0);
}

118
gnu/lib/libgmp/mpz_perfsqr.c Normal file
View file

@ -0,0 +1,118 @@
/* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a pefect square,
zero otherwise.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#if BITS_PER_MP_LIMB == 32
static unsigned int primes[] = {3, 5, 7, 11, 13, 17, 19, 23, 29};
static unsigned long int residue_map[] =
{0x3, 0x13, 0x17, 0x23b, 0x161b, 0x1a317, 0x30af3, 0x5335f, 0x13d122f3};
#define PP 0xC0CFD797L /* 3 x 5 x 7 x 11 x 13 x ... x 29 */
#endif
/* sq_res_0x100[x mod 0x100] == 1 iff x mod 0x100 is a quadratic residue
modulo 0x100. */
static char sq_res_0x100[0x100] =
{
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
};
int
#ifdef __STDC__
mpz_perfect_square_p (const MP_INT *a)
#else
mpz_perfect_square_p (a)
const MP_INT *a;
#endif
{
mp_limb n1, n0;
mp_size i;
mp_size asize = a->size;
mp_srcptr aptr = a->d;
mp_limb rem;
mp_ptr root_ptr;
/* No negative numbers are perfect squares. */
if (asize < 0)
return 0;
/* The first test excludes 55/64 (85.9%) of the perfect square candidates
in O(1) time. */
if (sq_res_0x100[aptr[0] % 0x100] == 0)
return 0;
#if BITS_PER_MP_LIMB == 32
/* The second test excludes 30652543/30808063 (99.5%) of the remaining
perfect square candidates in O(n) time. */
/* Firstly, compute REM = A mod PP. */
n1 = aptr[asize - 1];
if (n1 >= PP)
{
n1 = 0;
i = asize - 1;
}
else
i = asize - 2;
for (; i >= 0; i--)
{
mp_limb dummy;
n0 = aptr[i];
udiv_qrnnd (dummy, n1, n1, n0, PP);
}
rem = n1;
/* We have A mod PP in REM. Now decide if REM is a quadratic residue
modulo the factors in PP. */
for (i = 0; i < (sizeof primes) / sizeof (int); i++)
{
unsigned int p;
p = primes[i];
rem %= p;
if ((residue_map[i] & (1L << rem)) == 0)
return 0;
}
#endif
/* For the third and last test, we finally compute the square root,
to make sure we've really got a perfect square. */
root_ptr = (mp_ptr) alloca ((asize + 1) / 2 * BYTES_PER_MP_LIMB);
/* Iff mpn_sqrt returns zero, the square is perfect. */
{
int retval = !mpn_sqrt (root_ptr, NULL, aptr, asize);
alloca (0);
return retval;
}
}

110
gnu/lib/libgmp/mpz_pow_ui.c Normal file
View file

@ -0,0 +1,110 @@
/* mpz_pow_ui(res, base, exp) -- Set RES to BASE**EXP.
Copyright (C) 1991 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.
The GNU MP 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with the GNU MP Library; see the file COPYING. If not, write to
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#ifndef BERKELEY_MP
void
#ifdef __STDC__
mpz_pow_ui (MP_INT *r, const MP_INT *b, unsigned long int e)
#else
mpz_pow_ui (r, b, e)
MP_INT *r;
const MP_INT *b;
unsigned long int e;
#endif
#else /* BERKELEY_MP */
void
#ifdef __STDC__
rpow (const MP_INT *b, signed short int e, MP_INT *r)
#else
rpow (b, e, r)
const MP_INT *b;
signed short int e;
MP_INT *r;
#endif
#endif /* BERKELEY_MP */
{
mp_ptr rp, bp, tp, xp;
mp_size rsize, bsize;
int cnt, i;
bsize = ABS (b->size);
/* Single out cases that give result == 0 or 1. These tests are here
to simplify the general code below, not to optimize. */
if (bsize == 0
#ifdef BERKELEY_MP
|| e < 0
#endif
)
{
r->size = 0;
return;
}
if (e == 0)
{
r->d[0] = 1;
r->size = 1;
return;
}
/* Count the number of leading zero bits of the base's most
significant limb. */
count_leading_zeros (cnt, b->d[bsize - 1]);
/* Over-estimate space requirements and allocate enough space for the
final result in two temporary areas. The two areas are used to
alternately hold the input and recieve the product for mpn_mul.
(This scheme is used to fulfill the requirements of mpn_mul; that
the product space may not be the same as any of the input operands.) */
rsize = bsize * e - cnt * e / BITS_PER_MP_LIMB;
rp = (mp_ptr) alloca (rsize * BYTES_PER_MP_LIMB);
tp = (mp_ptr) alloca (rsize * BYTES_PER_MP_LIMB);
bp = b->d;
MPN_COPY (rp, bp, bsize);
rsize = bsize;
count_leading_zeros (cnt, e);
for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
{
rsize = mpn_mul (tp, rp, rsize, rp, rsize);
xp = tp; tp = rp; rp = xp;
if ((e & ((mp_limb) 1 << i)) != 0)
{
rsize = mpn_mul (tp, rp, rsize, bp, bsize);
xp = tp; tp = rp; rp = xp;
}
}
/* Now then we know the exact space requirements, reallocate if
necessary. */
if (r->alloc < rsize)
_mpz_realloc (r, rsize);
MPN_COPY (r->d, rp, rsize);
r->size = (e & 1) == 0 || b->size >= 0 ? rsize : -rsize;
alloca (0);
}

Some files were not shown because too many files have changed in this diff Show more