/[dtapublic]/projs/trunk/shared_source/c_datd/gmp_ints.c
ViewVC logotype

Annotation of /projs/trunk/shared_source/c_datd/gmp_ints.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations) (download)
Sat Oct 8 06:43:03 2016 UTC (7 years, 9 months ago) by dashley
Original Path: sf_code/esrgpcpj/shared/c_datd/gmp_ints.c
File MIME type: text/plain
File size: 147987 byte(s)
Initial commit.
1 dashley 25 // $Header: /cvsroot/esrg/sfesrg/esrgpcpj/shared/c_datd/gmp_ints.c,v 1.22 2002/01/27 15:18:44 dtashley Exp $
2    
3     //--------------------------------------------------------------------------------
4     //Copyright 2001 David T. Ashley
5     //-------------------------------------------------------------------------------------------------
6     //This source code and any program in which it is compiled/used is provided under the GNU GENERAL
7     //PUBLIC LICENSE, Version 3, full license text below.
8     //-------------------------------------------------------------------------------------------------
9     // GNU GENERAL PUBLIC LICENSE
10     // Version 3, 29 June 2007
11     //
12     // Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
13     // Everyone is permitted to copy and distribute verbatim copies
14     // of this license document, but changing it is not allowed.
15     //
16     // Preamble
17     //
18     // The GNU General Public License is a free, copyleft license for
19     //software and other kinds of works.
20     //
21     // The licenses for most software and other practical works are designed
22     //to take away your freedom to share and change the works. By contrast,
23     //the GNU General Public License is intended to guarantee your freedom to
24     //share and change all versions of a program--to make sure it remains free
25     //software for all its users. We, the Free Software Foundation, use the
26     //GNU General Public License for most of our software; it applies also to
27     //any other work released this way by its authors. You can apply it to
28     //your programs, too.
29     //
30     // When we speak of free software, we are referring to freedom, not
31     //price. Our General Public Licenses are designed to make sure that you
32     //have the freedom to distribute copies of free software (and charge for
33     //them if you wish), that you receive source code or can get it if you
34     //want it, that you can change the software or use pieces of it in new
35     //free programs, and that you know you can do these things.
36     //
37     // To protect your rights, we need to prevent others from denying you
38     //these rights or asking you to surrender the rights. Therefore, you have
39     //certain responsibilities if you distribute copies of the software, or if
40     //you modify it: responsibilities to respect the freedom of others.
41     //
42     // For example, if you distribute copies of such a program, whether
43     //gratis or for a fee, you must pass on to the recipients the same
44     //freedoms that you received. You must make sure that they, too, receive
45     //or can get the source code. And you must show them these terms so they
46     //know their rights.
47     //
48     // Developers that use the GNU GPL protect your rights with two steps:
49     //(1) assert copyright on the software, and (2) offer you this License
50     //giving you legal permission to copy, distribute and/or modify it.
51     //
52     // For the developers' and authors' protection, the GPL clearly explains
53     //that there is no warranty for this free software. For both users' and
54     //authors' sake, the GPL requires that modified versions be marked as
55     //changed, so that their problems will not be attributed erroneously to
56     //authors of previous versions.
57     //
58     // Some devices are designed to deny users access to install or run
59     //modified versions of the software inside them, although the manufacturer
60     //can do so. This is fundamentally incompatible with the aim of
61     //protecting users' freedom to change the software. The systematic
62     //pattern of such abuse occurs in the area of products for individuals to
63     //use, which is precisely where it is most unacceptable. Therefore, we
64     //have designed this version of the GPL to prohibit the practice for those
65     //products. If such problems arise substantially in other domains, we
66     //stand ready to extend this provision to those domains in future versions
67     //of the GPL, as needed to protect the freedom of users.
68     //
69     // Finally, every program is threatened constantly by software patents.
70     //States should not allow patents to restrict development and use of
71     //software on general-purpose computers, but in those that do, we wish to
72     //avoid the special danger that patents applied to a free program could
73     //make it effectively proprietary. To prevent this, the GPL assures that
74     //patents cannot be used to render the program non-free.
75     //
76     // The precise terms and conditions for copying, distribution and
77     //modification follow.
78     //
79     // TERMS AND CONDITIONS
80     //
81     // 0. Definitions.
82     //
83     // "This License" refers to version 3 of the GNU General Public License.
84     //
85     // "Copyright" also means copyright-like laws that apply to other kinds of
86     //works, such as semiconductor masks.
87     //
88     // "The Program" refers to any copyrightable work licensed under this
89     //License. Each licensee is addressed as "you". "Licensees" and
90     //"recipients" may be individuals or organizations.
91     //
92     // To "modify" a work means to copy from or adapt all or part of the work
93     //in a fashion requiring copyright permission, other than the making of an
94     //exact copy. The resulting work is called a "modified version" of the
95     //earlier work or a work "based on" the earlier work.
96     //
97     // A "covered work" means either the unmodified Program or a work based
98     //on the Program.
99     //
100     // To "propagate" a work means to do anything with it that, without
101     //permission, would make you directly or secondarily liable for
102     //infringement under applicable copyright law, except executing it on a
103     //computer or modifying a private copy. Propagation includes copying,
104     //distribution (with or without modification), making available to the
105     //public, and in some countries other activities as well.
106     //
107     // To "convey" a work means any kind of propagation that enables other
108     //parties to make or receive copies. Mere interaction with a user through
109     //a computer network, with no transfer of a copy, is not conveying.
110     //
111     // An interactive user interface displays "Appropriate Legal Notices"
112     //to the extent that it includes a convenient and prominently visible
113     //feature that (1) displays an appropriate copyright notice, and (2)
114     //tells the user that there is no warranty for the work (except to the
115     //extent that warranties are provided), that licensees may convey the
116     //work under this License, and how to view a copy of this License. If
117     //the interface presents a list of user commands or options, such as a
118     //menu, a prominent item in the list meets this criterion.
119     //
120     // 1. Source Code.
121     //
122     // The "source code" for a work means the preferred form of the work
123     //for making modifications to it. "Object code" means any non-source
124     //form of a work.
125     //
126     // A "Standard Interface" means an interface that either is an official
127     //standard defined by a recognized standards body, or, in the case of
128     //interfaces specified for a particular programming language, one that
129     //is widely used among developers working in that language.
130     //
131     // The "System Libraries" of an executable work include anything, other
132     //than the work as a whole, that (a) is included in the normal form of
133     //packaging a Major Component, but which is not part of that Major
134     //Component, and (b) serves only to enable use of the work with that
135     //Major Component, or to implement a Standard Interface for which an
136     //implementation is available to the public in source code form. A
137     //"Major Component", in this context, means a major essential component
138     //(kernel, window system, and so on) of the specific operating system
139     //(if any) on which the executable work runs, or a compiler used to
140     //produce the work, or an object code interpreter used to run it.
141     //
142     // The "Corresponding Source" for a work in object code form means all
143     //the source code needed to generate, install, and (for an executable
144     //work) run the object code and to modify the work, including scripts to
145     //control those activities. However, it does not include the work's
146     //System Libraries, or general-purpose tools or generally available free
147     //programs which are used unmodified in performing those activities but
148     //which are not part of the work. For example, Corresponding Source
149     //includes interface definition files associated with source files for
150     //the work, and the source code for shared libraries and dynamically
151     //linked subprograms that the work is specifically designed to require,
152     //such as by intimate data communication or control flow between those
153     //subprograms and other parts of the work.
154     //
155     // The Corresponding Source need not include anything that users
156     //can regenerate automatically from other parts of the Corresponding
157     //Source.
158     //
159     // The Corresponding Source for a work in source code form is that
160     //same work.
161     //
162     // 2. Basic Permissions.
163     //
164     // All rights granted under this License are granted for the term of
165     //copyright on the Program, and are irrevocable provided the stated
166     //conditions are met. This License explicitly affirms your unlimited
167     //permission to run the unmodified Program. The output from running a
168     //covered work is covered by this License only if the output, given its
169     //content, constitutes a covered work. This License acknowledges your
170     //rights of fair use or other equivalent, as provided by copyright law.
171     //
172     // You may make, run and propagate covered works that you do not
173     //convey, without conditions so long as your license otherwise remains
174     //in force. You may convey covered works to others for the sole purpose
175     //of having them make modifications exclusively for you, or provide you
176     //with facilities for running those works, provided that you comply with
177     //the terms of this License in conveying all material for which you do
178     //not control copyright. Those thus making or running the covered works
179     //for you must do so exclusively on your behalf, under your direction
180     //and control, on terms that prohibit them from making any copies of
181     //your copyrighted material outside their relationship with you.
182     //
183     // Conveying under any other circumstances is permitted solely under
184     //the conditions stated below. Sublicensing is not allowed; section 10
185     //makes it unnecessary.
186     //
187     // 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
188     //
189     // No covered work shall be deemed part of an effective technological
190     //measure under any applicable law fulfilling obligations under article
191     //11 of the WIPO copyright treaty adopted on 20 December 1996, or
192     //similar laws prohibiting or restricting circumvention of such
193     //measures.
194     //
195     // When you convey a covered work, you waive any legal power to forbid
196     //circumvention of technological measures to the extent such circumvention
197     //is effected by exercising rights under this License with respect to
198     //the covered work, and you disclaim any intention to limit operation or
199     //modification of the work as a means of enforcing, against the work's
200     //users, your or third parties' legal rights to forbid circumvention of
201     //technological measures.
202     //
203     // 4. Conveying Verbatim Copies.
204     //
205     // You may convey verbatim copies of the Program's source code as you
206     //receive it, in any medium, provided that you conspicuously and
207     //appropriately publish on each copy an appropriate copyright notice;
208     //keep intact all notices stating that this License and any
209     //non-permissive terms added in accord with section 7 apply to the code;
210     //keep intact all notices of the absence of any warranty; and give all
211     //recipients a copy of this License along with the Program.
212     //
213     // You may charge any price or no price for each copy that you convey,
214     //and you may offer support or warranty protection for a fee.
215     //
216     // 5. Conveying Modified Source Versions.
217     //
218     // You may convey a work based on the Program, or the modifications to
219     //produce it from the Program, in the form of source code under the
220     //terms of section 4, provided that you also meet all of these conditions:
221     //
222     // a) The work must carry prominent notices stating that you modified
223     // it, and giving a relevant date.
224     //
225     // b) The work must carry prominent notices stating that it is
226     // released under this License and any conditions added under section
227     // 7. This requirement modifies the requirement in section 4 to
228     // "keep intact all notices".
229     //
230     // c) You must license the entire work, as a whole, under this
231     // License to anyone who comes into possession of a copy. This
232     // License will therefore apply, along with any applicable section 7
233     // additional terms, to the whole of the work, and all its parts,
234     // regardless of how they are packaged. This License gives no
235     // permission to license the work in any other way, but it does not
236     // invalidate such permission if you have separately received it.
237     //
238     // d) If the work has interactive user interfaces, each must display
239     // Appropriate Legal Notices; however, if the Program has interactive
240     // interfaces that do not display Appropriate Legal Notices, your
241     // work need not make them do so.
242     //
243     // A compilation of a covered work with other separate and independent
244     //works, which are not by their nature extensions of the covered work,
245     //and which are not combined with it such as to form a larger program,
246     //in or on a volume of a storage or distribution medium, is called an
247     //"aggregate" if the compilation and its resulting copyright are not
248     //used to limit the access or legal rights of the compilation's users
249     //beyond what the individual works permit. Inclusion of a covered work
250     //in an aggregate does not cause this License to apply to the other
251     //parts of the aggregate.
252     //
253     // 6. Conveying Non-Source Forms.
254     //
255     // You may convey a covered work in object code form under the terms
256     //of sections 4 and 5, provided that you also convey the
257     //machine-readable Corresponding Source under the terms of this License,
258     //in one of these ways:
259     //
260     // a) Convey the object code in, or embodied in, a physical product
261     // (including a physical distribution medium), accompanied by the
262     // Corresponding Source fixed on a durable physical medium
263     // customarily used for software interchange.
264     //
265     // b) Convey the object code in, or embodied in, a physical product
266     // (including a physical distribution medium), accompanied by a
267     // written offer, valid for at least three years and valid for as
268     // long as you offer spare parts or customer support for that product
269     // model, to give anyone who possesses the object code either (1) a
270     // copy of the Corresponding Source for all the software in the
271     // product that is covered by this License, on a durable physical
272     // medium customarily used for software interchange, for a price no
273     // more than your reasonable cost of physically performing this
274     // conveying of source, or (2) access to copy the
275     // Corresponding Source from a network server at no charge.
276     //
277     // c) Convey individual copies of the object code with a copy of the
278     // written offer to provide the Corresponding Source. This
279     // alternative is allowed only occasionally and noncommercially, and
280     // only if you received the object code with such an offer, in accord
281     // with subsection 6b.
282     //
283     // d) Convey the object code by offering access from a designated
284     // place (gratis or for a charge), and offer equivalent access to the
285     // Corresponding Source in the same way through the same place at no
286     // further charge. You need not require recipients to copy the
287     // Corresponding Source along with the object code. If the place to
288     // copy the object code is a network server, the Corresponding Source
289     // may be on a different server (operated by you or a third party)
290     // that supports equivalent copying facilities, provided you maintain
291     // clear directions next to the object code saying where to find the
292     // Corresponding Source. Regardless of what server hosts the
293     // Corresponding Source, you remain obligated to ensure that it is
294     // available for as long as needed to satisfy these requirements.
295     //
296     // e) Convey the object code using peer-to-peer transmission, provided
297     // you inform other peers where the object code and Corresponding
298     // Source of the work are being offered to the general public at no
299     // charge under subsection 6d.
300     //
301     // A separable portion of the object code, whose source code is excluded
302     //from the Corresponding Source as a System Library, need not be
303     //included in conveying the object code work.
304     //
305     // A "User Product" is either (1) a "consumer product", which means any
306     //tangible personal property which is normally used for personal, family,
307     //or household purposes, or (2) anything designed or sold for incorporation
308     //into a dwelling. In determining whether a product is a consumer product,
309     //doubtful cases shall be resolved in favor of coverage. For a particular
310     //product received by a particular user, "normally used" refers to a
311     //typical or common use of that class of product, regardless of the status
312     //of the particular user or of the way in which the particular user
313     //actually uses, or expects or is expected to use, the product. A product
314     //is a consumer product regardless of whether the product has substantial
315     //commercial, industrial or non-consumer uses, unless such uses represent
316     //the only significant mode of use of the product.
317     //
318     // "Installation Information" for a User Product means any methods,
319     //procedures, authorization keys, or other information required to install
320     //and execute modified versions of a covered work in that User Product from
321     //a modified version of its Corresponding Source. The information must
322     //suffice to ensure that the continued functioning of the modified object
323     //code is in no case prevented or interfered with solely because
324     //modification has been made.
325     //
326     // If you convey an object code work under this section in, or with, or
327     //specifically for use in, a User Product, and the conveying occurs as
328     //part of a transaction in which the right of possession and use of the
329     //User Product is transferred to the recipient in perpetuity or for a
330     //fixed term (regardless of how the transaction is characterized), the
331     //Corresponding Source conveyed under this section must be accompanied
332     //by the Installation Information. But this requirement does not apply
333     //if neither you nor any third party retains the ability to install
334     //modified object code on the User Product (for example, the work has
335     //been installed in ROM).
336     //
337     // The requirement to provide Installation Information does not include a
338     //requirement to continue to provide support service, warranty, or updates
339     //for a work that has been modified or installed by the recipient, or for
340     //the User Product in which it has been modified or installed. Access to a
341     //network may be denied when the modification itself materially and
342     //adversely affects the operation of the network or violates the rules and
343     //protocols for communication across the network.
344     //
345     // Corresponding Source conveyed, and Installation Information provided,
346     //in accord with this section must be in a format that is publicly
347     //documented (and with an implementation available to the public in
348     //source code form), and must require no special password or key for
349     //unpacking, reading or copying.
350     //
351     // 7. Additional Terms.
352     //
353     // "Additional permissions" are terms that supplement the terms of this
354     //License by making exceptions from one or more of its conditions.
355     //Additional permissions that are applicable to the entire Program shall
356     //be treated as though they were included in this License, to the extent
357     //that they are valid under applicable law. If additional permissions
358     //apply only to part of the Program, that part may be used separately
359     //under those permissions, but the entire Program remains governed by
360     //this License without regard to the additional permissions.
361     //
362     // When you convey a copy of a covered work, you may at your option
363     //remove any additional permissions from that copy, or from any part of
364     //it. (Additional permissions may be written to require their own
365     //removal in certain cases when you modify the work.) You may place
366     //additional permissions on material, added by you to a covered work,
367     //for which you have or can give appropriate copyright permission.
368     //
369     // Notwithstanding any other provision of this License, for material you
370     //add to a covered work, you may (if authorized by the copyright holders of
371     //that material) supplement the terms of this License with terms:
372     //
373     // a) Disclaiming warranty or limiting liability differently from the
374     // terms of sections 15 and 16 of this License; or
375     //
376     // b) Requiring preservation of specified reasonable legal notices or
377     // author attributions in that material or in the Appropriate Legal
378     // Notices displayed by works containing it; or
379     //
380     // c) Prohibiting misrepresentation of the origin of that material, or
381     // requiring that modified versions of such material be marked in
382     // reasonable ways as different from the original version; or
383     //
384     // d) Limiting the use for publicity purposes of names of licensors or
385     // authors of the material; or
386     //
387     // e) Declining to grant rights under trademark law for use of some
388     // trade names, trademarks, or service marks; or
389     //
390     // f) Requiring indemnification of licensors and authors of that
391     // material by anyone who conveys the material (or modified versions of
392     // it) with contractual assumptions of liability to the recipient, for
393     // any liability that these contractual assumptions directly impose on
394     // those licensors and authors.
395     //
396     // All other non-permissive additional terms are considered "further
397     //restrictions" within the meaning of section 10. If the Program as you
398     //received it, or any part of it, contains a notice stating that it is
399     //governed by this License along with a term that is a further
400     //restriction, you may remove that term. If a license document contains
401     //a further restriction but permits relicensing or conveying under this
402     //License, you may add to a covered work material governed by the terms
403     //of that license document, provided that the further restriction does
404     //not survive such relicensing or conveying.
405     //
406     // If you add terms to a covered work in accord with this section, you
407     //must place, in the relevant source files, a statement of the
408     //additional terms that apply to those files, or a notice indicating
409     //where to find the applicable terms.
410     //
411     // Additional terms, permissive or non-permissive, may be stated in the
412     //form of a separately written license, or stated as exceptions;
413     //the above requirements apply either way.
414     //
415     // 8. Termination.
416     //
417     // You may not propagate or modify a covered work except as expressly
418     //provided under this License. Any attempt otherwise to propagate or
419     //modify it is void, and will automatically terminate your rights under
420     //this License (including any patent licenses granted under the third
421     //paragraph of section 11).
422     //
423     // However, if you cease all violation of this License, then your
424     //license from a particular copyright holder is reinstated (a)
425     //provisionally, unless and until the copyright holder explicitly and
426     //finally terminates your license, and (b) permanently, if the copyright
427     //holder fails to notify you of the violation by some reasonable means
428     //prior to 60 days after the cessation.
429     //
430     // Moreover, your license from a particular copyright holder is
431     //reinstated permanently if the copyright holder notifies you of the
432     //violation by some reasonable means, this is the first time you have
433     //received notice of violation of this License (for any work) from that
434     //copyright holder, and you cure the violation prior to 30 days after
435     //your receipt of the notice.
436     //
437     // Termination of your rights under this section does not terminate the
438     //licenses of parties who have received copies or rights from you under
439     //this License. If your rights have been terminated and not permanently
440     //reinstated, you do not qualify to receive new licenses for the same
441     //material under section 10.
442     //
443     // 9. Acceptance Not Required for Having Copies.
444     //
445     // You are not required to accept this License in order to receive or
446     //run a copy of the Program. Ancillary propagation of a covered work
447     //occurring solely as a consequence of using peer-to-peer transmission
448     //to receive a copy likewise does not require acceptance. However,
449     //nothing other than this License grants you permission to propagate or
450     //modify any covered work. These actions infringe copyright if you do
451     //not accept this License. Therefore, by modifying or propagating a
452     //covered work, you indicate your acceptance of this License to do so.
453     //
454     // 10. Automatic Licensing of Downstream Recipients.
455     //
456     // Each time you convey a covered work, the recipient automatically
457     //receives a license from the original licensors, to run, modify and
458     //propagate that work, subject to this License. You are not responsible
459     //for enforcing compliance by third parties with this License.
460     //
461     // An "entity transaction" is a transaction transferring control of an
462     //organization, or substantially all assets of one, or subdividing an
463     //organization, or merging organizations. If propagation of a covered
464     //work results from an entity transaction, each party to that
465     //transaction who receives a copy of the work also receives whatever
466     //licenses to the work the party's predecessor in interest had or could
467     //give under the previous paragraph, plus a right to possession of the
468     //Corresponding Source of the work from the predecessor in interest, if
469     //the predecessor has it or can get it with reasonable efforts.
470     //
471     // You may not impose any further restrictions on the exercise of the
472     //rights granted or affirmed under this License. For example, you may
473     //not impose a license fee, royalty, or other charge for exercise of
474     //rights granted under this License, and you may not initiate litigation
475     //(including a cross-claim or counterclaim in a lawsuit) alleging that
476     //any patent claim is infringed by making, using, selling, offering for
477     //sale, or importing the Program or any portion of it.
478     //
479     // 11. Patents.
480     //
481     // A "contributor" is a copyright holder who authorizes use under this
482     //License of the Program or a work on which the Program is based. The
483     //work thus licensed is called the contributor's "contributor version".
484     //
485     // A contributor's "essential patent claims" are all patent claims
486     //owned or controlled by the contributor, whether already acquired or
487     //hereafter acquired, that would be infringed by some manner, permitted
488     //by this License, of making, using, or selling its contributor version,
489     //but do not include claims that would be infringed only as a
490     //consequence of further modification of the contributor version. For
491     //purposes of this definition, "control" includes the right to grant
492     //patent sublicenses in a manner consistent with the requirements of
493     //this License.
494     //
495     // Each contributor grants you a non-exclusive, worldwide, royalty-free
496     //patent license under the contributor's essential patent claims, to
497     //make, use, sell, offer for sale, import and otherwise run, modify and
498     //propagate the contents of its contributor version.
499     //
500     // In the following three paragraphs, a "patent license" is any express
501     //agreement or commitment, however denominated, not to enforce a patent
502     //(such as an express permission to practice a patent or covenant not to
503     //sue for patent infringement). To "grant" such a patent license to a
504     //party means to make such an agreement or commitment not to enforce a
505     //patent against the party.
506     //
507     // If you convey a covered work, knowingly relying on a patent license,
508     //and the Corresponding Source of the work is not available for anyone
509     //to copy, free of charge and under the terms of this License, through a
510     //publicly available network server or other readily accessible means,
511     //then you must either (1) cause the Corresponding Source to be so
512     //available, or (2) arrange to deprive yourself of the benefit of the
513     //patent license for this particular work, or (3) arrange, in a manner
514     //consistent with the requirements of this License, to extend the patent
515     //license to downstream recipients. "Knowingly relying" means you have
516     //actual knowledge that, but for the patent license, your conveying the
517     //covered work in a country, or your recipient's use of the covered work
518     //in a country, would infringe one or more identifiable patents in that
519     //country that you have reason to believe are valid.
520     //
521     // If, pursuant to or in connection with a single transaction or
522     //arrangement, you convey, or propagate by procuring conveyance of, a
523     //covered work, and grant a patent license to some of the parties
524     //receiving the covered work authorizing them to use, propagate, modify
525     //or convey a specific copy of the covered work, then the patent license
526     //you grant is automatically extended to all recipients of the covered
527     //work and works based on it.
528     //
529     // A patent license is "discriminatory" if it does not include within
530     //the scope of its coverage, prohibits the exercise of, or is
531     //conditioned on the non-exercise of one or more of the rights that are
532     //specifically granted under this License. You may not convey a covered
533     //work if you are a party to an arrangement with a third party that is
534     //in the business of distributing software, under which you make payment
535     //to the third party based on the extent of your activity of conveying
536     //the work, and under which the third party grants, to any of the
537     //parties who would receive the covered work from you, a discriminatory
538     //patent license (a) in connection with copies of the covered work
539     //conveyed by you (or copies made from those copies), or (b) primarily
540     //for and in connection with specific products or compilations that
541     //contain the covered work, unless you entered into that arrangement,
542     //or that patent license was granted, prior to 28 March 2007.
543     //
544     // Nothing in this License shall be construed as excluding or limiting
545     //any implied license or other defenses to infringement that may
546     //otherwise be available to you under applicable patent law.
547     //
548     // 12. No Surrender of Others' Freedom.
549     //
550     // If conditions are imposed on you (whether by court order, agreement or
551     //otherwise) that contradict the conditions of this License, they do not
552     //excuse you from the conditions of this License. If you cannot convey a
553     //covered work so as to satisfy simultaneously your obligations under this
554     //License and any other pertinent obligations, then as a consequence you may
555     //not convey it at all. For example, if you agree to terms that obligate you
556     //to collect a royalty for further conveying from those to whom you convey
557     //the Program, the only way you could satisfy both those terms and this
558     //License would be to refrain entirely from conveying the Program.
559     //
560     // 13. Use with the GNU Affero General Public License.
561     //
562     // Notwithstanding any other provision of this License, you have
563     //permission to link or combine any covered work with a work licensed
564     //under version 3 of the GNU Affero General Public License into a single
565     //combined work, and to convey the resulting work. The terms of this
566     //License will continue to apply to the part which is the covered work,
567     //but the special requirements of the GNU Affero General Public License,
568     //section 13, concerning interaction through a network will apply to the
569     //combination as such.
570     //
571     // 14. Revised Versions of this License.
572     //
573     // The Free Software Foundation may publish revised and/or new versions of
574     //the GNU General Public License from time to time. Such new versions will
575     //be similar in spirit to the present version, but may differ in detail to
576     //address new problems or concerns.
577     //
578     // Each version is given a distinguishing version number. If the
579     //Program specifies that a certain numbered version of the GNU General
580     //Public License "or any later version" applies to it, you have the
581     //option of following the terms and conditions either of that numbered
582     //version or of any later version published by the Free Software
583     //Foundation. If the Program does not specify a version number of the
584     //GNU General Public License, you may choose any version ever published
585     //by the Free Software Foundation.
586     //
587     // If the Program specifies that a proxy can decide which future
588     //versions of the GNU General Public License can be used, that proxy's
589     //public statement of acceptance of a version permanently authorizes you
590     //to choose that version for the Program.
591     //
592     // Later license versions may give you additional or different
593     //permissions. However, no additional obligations are imposed on any
594     //author or copyright holder as a result of your choosing to follow a
595     //later version.
596     //
597     // 15. Disclaimer of Warranty.
598     //
599     // THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
600     //APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
601     //HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
602     //OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
603     //THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
604     //PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
605     //IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
606     //ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
607     //
608     // 16. Limitation of Liability.
609     //
610     // IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
611     //WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
612     //THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
613     //GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
614     //USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
615     //DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
616     //PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
617     //EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
618     //SUCH DAMAGES.
619     //
620     // 17. Interpretation of Sections 15 and 16.
621     //
622     // If the disclaimer of warranty and limitation of liability provided
623     //above cannot be given local legal effect according to their terms,
624     //reviewing courts shall apply local law that most closely approximates
625     //an absolute waiver of all civil liability in connection with the
626     //Program, unless a warranty or assumption of liability accompanies a
627     //copy of the Program in return for a fee.
628     //
629     // END OF TERMS AND CONDITIONS
630     //
631     // How to Apply These Terms to Your New Programs
632     //
633     // If you develop a new program, and you want it to be of the greatest
634     //possible use to the public, the best way to achieve this is to make it
635     //free software which everyone can redistribute and change under these terms.
636     //
637     // To do so, attach the following notices to the program. It is safest
638     //to attach them to the start of each source file to most effectively
639     //state the exclusion of warranty; and each file should have at least
640     //the "copyright" line and a pointer to where the full notice is found.
641     //
642     // <one line to give the program's name and a brief idea of what it does.>
643     // Copyright (C) <year> <name of author>
644     //
645     // This program is free software: you can redistribute it and/or modify
646     // it under the terms of the GNU General Public License as published by
647     // the Free Software Foundation, either version 3 of the License, or
648     // (at your option) any later version.
649     //
650     // This program is distributed in the hope that it will be useful,
651     // but WITHOUT ANY WARRANTY; without even the implied warranty of
652     // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
653     // GNU General Public License for more details.
654     //
655     // You should have received a copy of the GNU General Public License
656     // along with this program. If not, see <http://www.gnu.org/licenses/>.
657     //
658     //Also add information on how to contact you by electronic and paper mail.
659     //
660     // If the program does terminal interaction, make it output a short
661     //notice like this when it starts in an interactive mode:
662     //
663     // <program> Copyright (C) <year> <name of author>
664     // This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
665     // This is free software, and you are welcome to redistribute it
666     // under certain conditions; type `show c' for details.
667     //
668     //The hypothetical commands `show w' and `show c' should show the appropriate
669     //parts of the General Public License. Of course, your program's commands
670     //might be different; for a GUI interface, you would use an "about box".
671     //
672     // You should also get your employer (if you work as a programmer) or school,
673     //if any, to sign a "copyright disclaimer" for the program, if necessary.
674     //For more information on this, and how to apply and follow the GNU GPL, see
675     //<http://www.gnu.org/licenses/>.
676     //
677     // The GNU General Public License does not permit incorporating your program
678     //into proprietary programs. If your program is a subroutine library, you
679     //may consider it more useful to permit linking proprietary applications with
680     //the library. If this is what you want to do, use the GNU Lesser General
681     //Public License instead of this License. But first, please read
682     //<http://www.gnu.org/philosophy/why-not-lgpl.html>.
683     //-------------------------------------------------------------------------------------------------
684     //--------------------------------------------------------------------------------
685     #define MODULE_GMP_INTS
686    
687     #include <assert.h>
688     #include <ctype.h>
689     #include <process.h>
690     #include <stddef.h>
691     #include <stdio.h>
692     #include <string.h>
693     #include <time.h>
694    
695    
696     /* Only included the guarded allocation header if we are compiling
697     ** a DOS console type application. Other types of applications have
698     ** other ways of protecting for out of memory. Including the
699     ** header would do no harm in these cases, but do no good, either.
700     */
701     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
702     #include "ccmalloc.h"
703     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
704     #include "tclalloc.h"
705     #else
706     /* Do nothing. */
707     #endif
708    
709     #include "bstrfunc.h"
710     #include "charfunc.h"
711     #include "fcmiof.h"
712     #include "gmp_ints.h"
713     #include "intfunc.h"
714    
715    
716     /******************************************************************/
717     /*** CUSTOM ALLOCATION FUNCTIONS *******************************/
718     /******************************************************************/
719     /* We need to decide here on how memory not on the stack will be
720     ** allocated (i.e. what will become of the standard functions
721     ** like malloc, free, etc.). This is dependent on the type of
722     ** application we're making. The possible types are so far are:
723     ** APP_TYPE_SIMPLE_DOS_CONSOLE :
724     ** Simple DOS console application.
725     ** APP_TYPE_IJUSCRIPTER_IJUCONSOLE:
726     ** The Tcl tool build.
727     **
728     ** The custom allocation functions here are a "portal" or "wrapper"
729     ** for how the integer and rational number functions should
730     ** get memory.
731     **
732     ** The functions below are standard, except that the GNU MP team
733     ** built more generality into what allocation schemes could be
734     ** used by including size information in some calls that don't
735     ** normally get it. That is why there are some extra calls below
736     ** where the information is discarded. Other than that, these are
737     ** standard allocation calls.
738     */
739     //07/15/01: Visual inspection only. Function deemed too
740     //simple for unit testing.
741     void *GMP_INTS_malloc( size_t size )
742     {
743     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
744     return(CCMALLOC_malloc(size));
745     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
746     return(TclpAlloc(size));
747     #else
748     return(malloc(size));
749     #endif
750     }
751    
752    
753     //07/15/01: Visual inspection only. Function deemed too
754     //simple for unit testing.
755     void *GMP_INTS_calloc( size_t num, size_t size )
756     {
757     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
758     return(CCMALLOC_calloc(num, size));
759     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
760     return(TclpCalloc(num, size));
761     #else
762     return(calloc(num, size));
763     #endif
764     }
765    
766    
767     //07/15/01: Visual inspection only. Function deemed too
768     //simple for unit testing.
769     void *GMP_INTS_realloc( void *memblock, size_t size )
770     {
771     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
772     return(CCMALLOC_realloc(memblock, size));
773     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
774     return(TclpRealloc(memblock, size));
775     #else
776     return(realloc(memblock, size));
777     #endif
778     }
779    
780    
781     //07/15/01: Visual inspection only. Function deemed too
782     //simple for unit testing.
783     void *GMP_INTS_realloc_w_size( void *memblock,
784     size_t old_size,
785     size_t size )
786     {
787     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
788     return(CCMALLOC_realloc(memblock, size));
789     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
790     return(TclpRealloc(memblock, size));
791     #else
792     return(realloc(memblock, size));
793     #endif
794     }
795    
796    
797     //07/15/01: Visual inspection only. Function deemed too
798     //simple for unit testing.
799     void GMP_INTS_free( void *memblock )
800     {
801     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
802     CCMALLOC_free(memblock);
803     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
804     TclpFree(memblock);
805     #else
806     free(memblock);
807     #endif
808     }
809    
810    
811     //07/15/01: Visual inspection only. Function deemed too
812     //simple for unit testing.
813     void GMP_INTS_free_w_size( void *memblock, size_t size )
814     {
815     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
816     CCMALLOC_free(memblock);
817     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
818     TclpFree(memblock);
819     #else
820     free(memblock);
821     #endif
822     }
823    
824    
825     /******************************************************************/
826     /*** PORTABILITY CHECK FUNCTIONS *******************************/
827     /******************************************************************/
828     //Because there is the risk that Microsoft Visual C++ might
829     //change in the future, the following function can be called
830     //to see if the assumptions about data sizes are valid. This
831     //function returns TRUE if there is a problem, or FALSE
832     //otherwise.
833    
834     //07/15/01: Unit testing complete.
835     int GMP_INTS_data_sizes_are_wrong(void)
836     {
837     int i;
838     GMP_INTS_limb_t tv;
839     _int64 tv64;
840    
841     //Check the number of bit rolls required to get the limb
842     //to go to zero again. This had better be 32.
843     tv = 1;
844     i = 0;
845     while (tv)
846     {
847     tv <<= 1;
848     i++;
849     }
850     if (i != 32)
851     return(1);
852    
853     //Check that an _int64 is really and truly 64 bits.
854     tv64 = 1;
855     i = 0;
856     while (tv64)
857     {
858     tv64 <<= 1;
859     i++;
860     }
861     if (i != 64)
862     return(1);
863    
864     //Room for additional tests here if needed later.
865    
866     return(0);
867     }
868    
869    
870     /******************************************************************/
871     /*** ERROR STRING IDENTIFICATION AND PROCESSING FUNCTIONS *******/
872     /******************************************************************/
873    
874     int GMP_INTS_identify_nan_string(const char *s)
875     {
876     assert(s != NULL);
877    
878     if (!strcmp(s, GMP_INTS_EF_INTOVF_POS_STRING))
879     return(0);
880     else if (!strcmp(s, GMP_INTS_EF_INTOVF_NEG_STRING))
881     return(1);
882     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_POS_STRING))
883     return(2);
884     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING))
885     return(3);
886     else
887     return(-1);
888     }
889    
890    
891     const char *GMP_INTS_supply_nan_string(int idx)
892     {
893     assert((idx >= 0) && (idx <= 3));
894    
895     if (idx==0)
896     return(GMP_INTS_EF_INTOVF_POS_STRING);
897     else if (idx==1)
898     return(GMP_INTS_EF_INTOVF_NEG_STRING);
899     else if (idx==2)
900     return(GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
901     else
902     return(GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
903     }
904    
905    
906     /******************************************************************/
907     /*** DEBUG PRINTING FUNCTIONS **********************************/
908     /******************************************************************/
909     //These functions are for printing out integers and limbs
910     //and groups of limbs for unit testing and debugging.
911    
912     //07/15/01: Exempt from testing because debug/development
913     //function.
914     void GMP_INTS_print_limb_group(FILE *stream,
915     GMP_INTS_limb_srcptr lg,
916     GMP_INTS_size_t n,
917     char *desc)
918     {
919     int i;
920    
921     assert(stream != NULL);
922     assert(n >= 0);
923     assert(desc != NULL);
924    
925     if (!lg)
926     {
927     fprintf(stream, " %s: NULL\n", desc);
928     }
929     else
930     {
931     for (i=n-1; i>=0; i--)
932     {
933     fprintf(stream, " %s[%2d]: 0x%8X\n", desc, i, lg[i]);
934     }
935     }
936     }
937    
938    
939     void GMP_INTS_mpz_print_int(FILE *stream,
940     const GMP_INTS_mpz_struct *arg,
941     char *desc)
942     {
943     int i;
944    
945     assert(stream != NULL);
946     assert(arg != NULL);
947     assert(desc != NULL);
948    
949     fprintf(stream, "Printing integer:\n %s\n", desc);
950    
951     fprintf(stream, " flags: %d\n", arg->flags);
952     fprintf(stream, " ptr value to body: %p\n", arg);
953     fprintf(stream, " n_allocd: %d\n", arg->n_allocd);
954     fprintf(stream, " size: %d\n", arg->size);
955     fprintf(stream, " limbs (ptr val): %p\n", arg->limbs);
956    
957     for (i=arg->n_allocd-1; i>=0; i--)
958     {
959     fprintf(stream, " limb[%3d]: %8X\n", i, arg->limbs[i]);
960     }
961     }
962    
963    
964     /******************************************************************/
965     /*** LOW-LEVEL MACRO REPLACEMENTS ******************************/
966     /******************************************************************/
967     //The functions in this category are replacements for macros.
968     //Clarity was gained at the expense of speed.
969    
970     int GMP_INTS_mpz_get_flags (const GMP_INTS_mpz_struct *arg)
971     {
972     assert(arg != NULL);
973     assert(arg->n_allocd > 0);
974    
975     return(arg->flags);
976     }
977    
978    
979     //07/15/01: Visual inspection only. Function deemed too
980     //simple for unit testing.
981     GMP_INTS_size_t GMP_INTS_abs_of_size_t(GMP_INTS_size_t arg)
982     {
983     //Be sure that the bit pattern does not represent the maximum
984     //negative argument. Negating this would give the result of
985     //zero, which is not what is intended.
986     assert(arg != 0x80000000);
987    
988     if (arg < 0)
989     return(-arg);
990     else
991     return(arg);
992     }
993    
994    
995     //07/15/01: Visual inspection only. Function deemed too
996     //simple for unit testing.
997     int GMP_INTS_mpz_sgn(const GMP_INTS_mpz_struct *arg)
998     {
999     assert(arg != NULL);
1000     assert(arg->n_allocd > 0);
1001    
1002     if (arg->size > 0)
1003     return(1);
1004     else if (arg->size == 0)
1005     return(0);
1006     else
1007     return(-1);
1008     }
1009    
1010    
1011     //07/15/01: Visual inspection only. Function deemed too
1012     //simple for unit testing.
1013     int GMP_INTS_mpz_is_neg(const GMP_INTS_mpz_struct *arg)
1014     {
1015     assert(arg != NULL);
1016     assert(arg->n_allocd > 0);
1017    
1018     if (GMP_INTS_mpz_sgn(arg) == -1)
1019     return(1);
1020     else
1021     return(0);
1022     }
1023    
1024    
1025     //07/15/01: Visual inspection only. Function deemed too
1026     //simple for unit testing.
1027     int GMP_INTS_mpz_is_zero(const GMP_INTS_mpz_struct *arg)
1028     {
1029     assert(arg != NULL);
1030     assert(arg->n_allocd > 0);
1031    
1032     if (GMP_INTS_mpz_sgn(arg) == 0)
1033     return(1);
1034     else
1035     return(0);
1036     }
1037    
1038    
1039     //07/15/01: Visual inspection only. Function deemed too
1040     //simple for unit testing.
1041     int GMP_INTS_mpz_is_pos(const GMP_INTS_mpz_struct *arg)
1042     {
1043     assert(arg != NULL);
1044     assert(arg->n_allocd > 0);
1045    
1046     if (GMP_INTS_mpz_sgn(arg) == 1)
1047     return(1);
1048     else
1049     return(0);
1050     }
1051    
1052    
1053     //07/15/01: Visual inspection only. Function deemed too
1054     //simple for unit testing.
1055     int GMP_INTS_mpz_is_odd(const GMP_INTS_mpz_struct *arg)
1056     {
1057     assert(arg != NULL);
1058     assert(arg->n_allocd > 0);
1059    
1060     if (arg->size == 0)
1061     return 0;
1062     else if ((arg->limbs[0] & 0x1) != 0)
1063     return 1;
1064     else
1065     return 0;
1066     }
1067    
1068    
1069     //07/15/01: Visual inspection only. Function deemed too
1070     //simple for unit testing.
1071     int GMP_INTS_mpz_is_even(const GMP_INTS_mpz_struct *arg)
1072     {
1073     assert(arg != NULL);
1074     assert(arg->n_allocd > 0);
1075    
1076     if (GMP_INTS_mpz_is_odd(arg))
1077     return 0;
1078     else
1079     return 1;
1080     }
1081    
1082    
1083     void GMP_INTS_mpz_negate(GMP_INTS_mpz_struct *arg)
1084     {
1085     //Eyeball the input parameters.
1086     assert(arg != NULL);
1087     assert(arg->n_allocd > 0);
1088     assert(arg->limbs != NULL);
1089    
1090     arg->size = -(arg->size);
1091     }
1092    
1093    
1094     //07/15/01: Visual inspection only. Function deemed too
1095     //simple for unit testing.
1096     void GMP_INTS_mpn_normalize(GMP_INTS_limb_ptr limb_array,
1097     GMP_INTS_size_t *idx)
1098     {
1099     assert(limb_array != NULL);
1100     assert(idx != NULL);
1101     assert(idx >= 0);
1102    
1103     while (*idx > 0)
1104     {
1105     if (limb_array[*idx - 1] != 0)
1106     break;
1107     (*idx)--;
1108     }
1109     }
1110    
1111    
1112     //07/15/01: Visual inspection only. Function deemed too
1113     //simple for unit testing.
1114     void GMP_INTS_mpn_copy_limbs(GMP_INTS_limb_ptr dest,
1115     GMP_INTS_limb_srcptr src,
1116     GMP_INTS_size_t n)
1117     {
1118     GMP_INTS_size_t i;
1119    
1120     assert(dest != NULL);
1121     assert(src != NULL);
1122     assert(n >= 0);
1123    
1124     for (i=0; i<n; i++)
1125     dest[i] = src[i];
1126     }
1127    
1128    
1129     /******************************************************************/
1130     /*** LOW-LEVEL ARITHMETIC FUNCTIONS ****************************/
1131     /******************************************************************/
1132    
1133     int GMP_INTS_two_op_flags_map(int flags1, int flags2)
1134     {
1135     int cf;
1136    
1137     if (!flags1 && !flags2)
1138     {
1139     return(0);
1140     }
1141     else
1142     {
1143     cf = flags1 | flags2;
1144    
1145     if (cf & (GMP_INTS_EF_INTOVF_POS | GMP_INTS_EF_INTOVF_TAINT_POS))
1146     {
1147     //In either case here, the result will be tainted
1148     //positive.
1149     return(GMP_INTS_EF_INTOVF_TAINT_POS);
1150     }
1151     else if (cf & (GMP_INTS_EF_INTOVF_NEG | GMP_INTS_EF_INTOVF_TAINT_NEG))
1152     {
1153     //In either case here, the result will be tainted
1154     //negative.
1155     return(GMP_INTS_EF_INTOVF_TAINT_NEG);
1156     }
1157     else
1158     {
1159     //This case is where the flags ints are non-zero, but
1160     //no known bits are set. This is surely some kind of
1161     //an internal software error. In debug mode, want to
1162     //alert to error. In actual operation, want to just
1163     //pretend an ordinary positive taint.
1164     assert(0);
1165     return(GMP_INTS_EF_INTOVF_TAINT_POS);
1166     }
1167     }
1168     }
1169    
1170    
1171     GMP_INTS_limb_t GMP_INTS_mpn_add_1 (GMP_INTS_limb_ptr res_ptr,
1172     GMP_INTS_limb_srcptr s1_ptr,
1173     GMP_INTS_size_t s1_size,
1174     GMP_INTS_limb_t s2_limb)
1175     {
1176     GMP_INTS_limb_t x;
1177    
1178     assert(res_ptr != NULL);
1179     assert(s1_ptr != NULL);
1180     assert(s1_size > 0);
1181    
1182     x = *s1_ptr++;
1183     s2_limb = x + s2_limb;
1184     *res_ptr++ = s2_limb;
1185     //Since limbs are unsigned, the test below tests if there
1186     //was a carry, i.e. a positive rollover.
1187     if (s2_limb < x)
1188     {
1189     while (--s1_size != 0)
1190     {
1191     x = *s1_ptr++ + 1;
1192     *res_ptr++ = x;
1193     if (x != 0)
1194     goto fin;
1195     }
1196    
1197     return 1;
1198     }
1199    
1200     fin:
1201     if (res_ptr != s1_ptr)
1202     {
1203     GMP_INTS_size_t i;
1204     for (i = 0; i < s1_size - 1; i++)
1205     {
1206     res_ptr[i] = s1_ptr[i];
1207     }
1208     }
1209    
1210     return 0;
1211     }
1212    
1213    
1214     GMP_INTS_limb_t GMP_INTS_mpn_sub_1(GMP_INTS_limb_ptr res_ptr,
1215     GMP_INTS_limb_srcptr s1_ptr,
1216     GMP_INTS_size_t s1_size,
1217     GMP_INTS_limb_t s2_limb)
1218     {
1219     GMP_INTS_limb_t x;
1220    
1221     assert(res_ptr != NULL);
1222     assert(s1_ptr != NULL);
1223     assert(s1_size > 0);
1224    
1225     x = *s1_ptr++;
1226     s2_limb = x - s2_limb;
1227     *res_ptr++ = s2_limb;
1228     //The test below detects a borrow.
1229     if (s2_limb > x)
1230     {
1231     while (--s1_size != 0)
1232     {
1233     x = *s1_ptr++;
1234     *res_ptr++ = x - 1;
1235     if (x != 0)
1236     goto fin;
1237     }
1238    
1239     return 1;
1240     }
1241    
1242     fin:
1243     if (res_ptr != s1_ptr)
1244     {
1245     GMP_INTS_size_t i;
1246     for (i = 0; i < s1_size - 1; i++)
1247     {
1248     res_ptr[i] = s1_ptr[i];
1249     }
1250     }
1251     return 0;
1252     }
1253    
1254    
1255     //07/15/01: Am willing to skip unit-testing on this.
1256     //Understand the logic (i.e. passes visual inspection),
1257     //and comes from GNU-MP. Hope any defects here will be
1258     //caught in testing of GMP_INTS_mpz_mul() and other
1259     //higher-level functions.
1260     GMP_INTS_limb_t GMP_INTS_mpn_mul_1 (GMP_INTS_limb_ptr res_ptr,
1261     GMP_INTS_limb_srcptr s1_ptr,
1262     GMP_INTS_size_t s1_size,
1263     GMP_INTS_limb_t s2_limb)
1264     {
1265     GMP_INTS_limb_t cy_limb;
1266     GMP_INTS_size_t j;
1267     GMP_INTS_limb_t prod_high, prod_low;
1268     unsigned _int64 temp;
1269    
1270     assert(res_ptr != NULL);
1271     assert(s1_ptr != NULL);
1272     assert(s1_size > 0);
1273    
1274     /* The loop counter and index J goes from -S1_SIZE to -1. This way
1275     the loop becomes faster. */
1276     j = -s1_size;
1277    
1278     /* Offset the base pointers to compensate for the negative indices. */
1279     s1_ptr -= j;
1280     res_ptr -= j;
1281    
1282     cy_limb = 0;
1283     do
1284     {
1285     //The original code here was the following macro:
1286     //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
1287     //Will use the 64-bit data type of MSVC++ to achieve
1288     //the same effect.
1289     //
1290     //NOTE AS OF 07/13/01: I have looked at the assembly-
1291     //language, and the lines below are a real sore spot.
1292     //The multiply is fairly direct (although there is a
1293     //function call), but the shift does not behave as
1294     //expected--there is a function call and a loop to
1295     //go through the 32 iterations. After logical testing,
1296     //may want to clean this out--this would surely
1297     //result in a speed increase. This is a sore spot.
1298     temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
1299     prod_low = (GMP_INTS_limb_t)temp;
1300     prod_high = (GMP_INTS_limb_t)(temp >> 32);
1301    
1302     prod_low += cy_limb;
1303     cy_limb = (prod_low < cy_limb) + prod_high;
1304    
1305     res_ptr[j] = prod_low;
1306     }
1307     while (++j != 0);
1308    
1309     return cy_limb;
1310     }
1311    
1312    
1313     //07/15/01: Am willing to skip unit-testing on this.
1314     //Understand the logic (i.e. passes visual inspection),
1315     //and comes from GNU-MP. Hope any defects here will be
1316     //caught in testing of GMP_INTS_mpz_add() and other
1317     //higher-level functions.
1318     GMP_INTS_limb_t GMP_INTS_mpn_add_n(GMP_INTS_limb_ptr res_ptr,
1319     GMP_INTS_limb_srcptr s1_ptr,
1320     GMP_INTS_limb_srcptr s2_ptr,
1321     GMP_INTS_size_t size)
1322     {
1323     GMP_INTS_limb_t x, y, cy;
1324     GMP_INTS_size_t j;
1325    
1326     assert(res_ptr != NULL);
1327     assert(s1_ptr != NULL);
1328     assert(s2_ptr != NULL);
1329    
1330     /* The loop counter and index J goes from -SIZE to -1. This way
1331     the loop becomes faster. */
1332     j = -size;
1333    
1334     /* Offset the base pointers to compensate for the negative indices. */
1335     s1_ptr -= j;
1336     s2_ptr -= j;
1337     res_ptr -= j;
1338    
1339     cy = 0;
1340     do
1341     {
1342     y = s2_ptr[j];
1343     x = s1_ptr[j];
1344     y += cy; /* add previous carry to one addend */
1345     cy = (y < cy); /* get out carry from that addition */
1346     y = x + y; /* add other addend */
1347     cy = (y < x) + cy; /* get out carry from that add, combine */
1348     res_ptr[j] = y;
1349     }
1350     while (++j != 0);
1351    
1352     return cy;
1353     }
1354    
1355    
1356     //07/15/01: Am willing to skip unit-testing on this.
1357     //Understand the logic (i.e. passes visual inspection),
1358     //and comes from GNU-MP. Hope any defects here will be
1359     //caught in testing of GMP_INTS_mpz_mul() and other
1360     //higher-level functions.
1361     GMP_INTS_limb_t GMP_INTS_mpn_addmul_1 (GMP_INTS_limb_ptr res_ptr,
1362     GMP_INTS_limb_srcptr s1_ptr,
1363     GMP_INTS_size_t s1_size,
1364     GMP_INTS_limb_t s2_limb)
1365     {
1366     GMP_INTS_limb_t cy_limb;
1367     GMP_INTS_size_t j;
1368     GMP_INTS_limb_t prod_high, prod_low;
1369     GMP_INTS_limb_t x;
1370     unsigned _int64 temp;
1371    
1372     //Eyeball the inputs carefully.
1373     assert(res_ptr != NULL);
1374     assert(s1_ptr != NULL);
1375     assert(s1_size > 0);
1376    
1377     /* The loop counter and index J goes from -SIZE to -1. This way
1378     the loop becomes faster. */
1379     j = -s1_size;
1380    
1381     /* Offset the base pointers to compensate for the negative indices. */
1382     res_ptr -= j;
1383     s1_ptr -= j;
1384    
1385     cy_limb = 0;
1386     do
1387     {
1388     //The original code here was the following macro:
1389     //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
1390     //Will use the 64-bit data type of MSVC++ to achieve
1391     //the same effect.
1392     //
1393     //NOTE AS OF 07/14/01: I have not looked at the assembly-
1394     //language, but the assembly-language generated by what
1395     //is below is suspected to have performance problems.
1396     //May want to come back to this.
1397     temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
1398     prod_low = (GMP_INTS_limb_t)temp;
1399     prod_high = (GMP_INTS_limb_t)(temp >> 32);
1400    
1401     prod_low += cy_limb;
1402     cy_limb = (prod_low < cy_limb) + prod_high;
1403    
1404     x = res_ptr[j];
1405     prod_low = x + prod_low;
1406     cy_limb += (prod_low < x);
1407     res_ptr[j] = prod_low;
1408     }
1409     while (++j != 0);
1410    
1411     return cy_limb;
1412     }
1413    
1414    
1415     //07/15/01: Am willing to skip unit-testing on this.
1416     //Understand the logic (i.e. passes visual inspection),
1417     //and comes from GNU-MP.
1418     GMP_INTS_limb_t GMP_INTS_mpn_add (GMP_INTS_limb_ptr res_ptr,
1419     GMP_INTS_limb_srcptr s1_ptr,
1420     GMP_INTS_size_t s1_size,
1421     GMP_INTS_limb_srcptr s2_ptr,
1422     GMP_INTS_size_t s2_size)
1423     {
1424     GMP_INTS_limb_t cy_limb = 0;
1425    
1426     assert(res_ptr != NULL);
1427     assert(s1_ptr != NULL);
1428     assert(s2_ptr != NULL);
1429    
1430     //Numbers apparently must be arranged with sizes so that
1431     //LIMBS(s1) >= LIMBS(s2).
1432     //Add the parts up to the most significant limb of S2.
1433     if (s2_size != 0)
1434     cy_limb = GMP_INTS_mpn_add_n (res_ptr,
1435     s1_ptr,
1436     s2_ptr,
1437     s2_size);
1438    
1439     //Process the carry result, and propagate the carries up through
1440     //the parts of S1 that don't exist in S2, i.e. propagate the
1441     //carries upward in S1.
1442     if (s1_size - s2_size != 0)
1443     cy_limb = GMP_INTS_mpn_add_1 (res_ptr + s2_size,
1444     s1_ptr + s2_size,
1445     s1_size - s2_size,
1446     cy_limb);
1447     return cy_limb;
1448     }
1449    
1450    
1451     //07/15/01: Am willing to skip unit-testing on this.
1452     //Understand the logic (i.e. passes visual inspection),
1453     //and comes from GNU-MP.
1454     GMP_INTS_limb_t GMP_INTS_mpn_sub_n(GMP_INTS_limb_ptr res_ptr,
1455     GMP_INTS_limb_srcptr s1_ptr,
1456     GMP_INTS_limb_srcptr s2_ptr,
1457     GMP_INTS_size_t size)
1458     {
1459     GMP_INTS_limb_t x, y, cy;
1460     GMP_INTS_size_t j;
1461    
1462     assert(res_ptr != NULL);
1463     assert(s1_ptr != NULL);
1464     assert(s2_ptr != NULL);
1465    
1466     /* The loop counter and index J goes from -SIZE to -1. This way
1467     the loop becomes faster. */
1468     j = -size;
1469    
1470     /* Offset the base pointers to compensate for the negative indices. */
1471     s1_ptr -= j;
1472     s2_ptr -= j;
1473     res_ptr -= j;
1474    
1475     cy = 0;
1476     do
1477     {
1478     y = s2_ptr[j];
1479     x = s1_ptr[j];
1480     y += cy; /* add previous carry to subtrahend */
1481     cy = (y < cy); /* get out carry from that addition */
1482     y = x - y; /* main subtract */
1483     cy = (y > x) + cy; /* get out carry from the subtract, combine */
1484     res_ptr[j] = y;
1485     }
1486     while (++j != 0);
1487    
1488     return cy;
1489     }
1490    
1491    
1492     //07/17/01: Am willing to skip unit-testing on this.
1493     //Understand the logic (i.e. passes visual inspection),
1494     //and comes from GNU-MP.
1495     GMP_INTS_limb_t GMP_INTS_mpn_sub (GMP_INTS_limb_ptr res_ptr,
1496     GMP_INTS_limb_srcptr s1_ptr,
1497     GMP_INTS_size_t s1_size,
1498     GMP_INTS_limb_srcptr s2_ptr,
1499     GMP_INTS_size_t s2_size)
1500     {
1501     GMP_INTS_limb_t cy_limb = 0;
1502    
1503     assert(res_ptr != NULL);
1504     assert(s1_ptr != NULL);
1505     assert(s2_ptr != NULL);
1506    
1507     if (s2_size != 0)
1508     cy_limb = GMP_INTS_mpn_sub_n(res_ptr,
1509     s1_ptr,
1510     s2_ptr,
1511     s2_size);
1512    
1513     if (s1_size - s2_size != 0)
1514     cy_limb = GMP_INTS_mpn_sub_1(res_ptr + s2_size,
1515     s1_ptr + s2_size,
1516     s1_size - s2_size,
1517     cy_limb);
1518     return cy_limb;
1519     }
1520    
1521    
1522     //07/17/01: Am willing to skip unit-testing on this.
1523     //Understand the logic (i.e. passes visual inspection),
1524     //and comes from GNU-MP.
1525     GMP_INTS_limb_t GMP_INTS_mpn_lshift(GMP_INTS_limb_ptr wp,
1526     GMP_INTS_limb_srcptr up,
1527     GMP_INTS_size_t usize,
1528     unsigned int cnt)
1529     {
1530     GMP_INTS_limb_t high_limb, low_limb;
1531     unsigned sh_1, sh_2;
1532     GMP_INTS_size_t i;
1533     GMP_INTS_limb_t retval;
1534    
1535     assert(wp != NULL);
1536     assert(up != NULL);
1537     assert(usize > 0);
1538     assert(cnt > 0);
1539    
1540     sh_1 = cnt;
1541    
1542     wp += 1;
1543     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
1544     //This automatically implies that can't call this function to shift more
1545     //than 31 places.
1546     i = usize - 1;
1547     low_limb = up[i];
1548     retval = low_limb >> sh_2; //Return value is the amount shifted
1549     //off the top.
1550     high_limb = low_limb;
1551     while (--i >= 0)
1552     {
1553     low_limb = up[i];
1554     wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
1555     high_limb = low_limb;
1556     }
1557     wp[i] = high_limb << sh_1;
1558    
1559     return retval;
1560     }
1561    
1562    
1563     //07/17/01: Am willing to skip unit-testing on this.
1564     //Understand the logic more or less (i.e. passes visual inspection),
1565     //and comes from GNU-MP.
1566     /* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
1567     and store the USIZE least significant limbs of the result at WP.
1568     The bits shifted out to the right are returned.
1569    
1570     Argument constraints:
1571     1. 0 < CNT < BITS_PER_MP_LIMB
1572     2. If the result is to be written over the input, WP must be <= UP.
1573     */
1574     GMP_INTS_limb_t GMP_INTS_mpn_rshift (GMP_INTS_limb_ptr wp,
1575     GMP_INTS_limb_srcptr up,
1576     GMP_INTS_size_t usize,
1577     unsigned int cnt)
1578     {
1579     GMP_INTS_limb_t high_limb, low_limb;
1580     unsigned sh_1, sh_2;
1581     GMP_INTS_size_t i;
1582     GMP_INTS_limb_t retval;
1583    
1584     assert(wp != NULL);
1585     assert(up != NULL);
1586     assert(usize > 0);
1587     assert(cnt > 0);
1588    
1589     sh_1 = cnt;
1590    
1591     wp -= 1;
1592     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
1593     high_limb = up[0];
1594     retval = high_limb << sh_2;
1595     low_limb = high_limb;
1596    
1597     for (i = 1; i < usize; i++)
1598     {
1599     high_limb = up[i];
1600     wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
1601     low_limb = high_limb;
1602     }
1603     wp[i] = low_limb >> sh_1;
1604    
1605     return retval;
1606     }
1607    
1608    
1609     //07/17/01: Am willing to skip unit-testing on this.
1610     //Understand the logic (i.e. passes visual inspection),
1611     //and comes from GNU-MP.
1612     int GMP_INTS_mpn_cmp (GMP_INTS_limb_srcptr op1_ptr,
1613     GMP_INTS_limb_srcptr op2_ptr,
1614     GMP_INTS_size_t size)
1615     {
1616     GMP_INTS_size_t i;
1617     GMP_INTS_limb_t op1_word, op2_word;
1618    
1619     assert(op1_ptr != NULL);
1620     assert(op2_ptr != NULL);
1621    
1622     for (i = size - 1; i >= 0; i--)
1623     {
1624     op1_word = op1_ptr[i];
1625     op2_word = op2_ptr[i];
1626     if (op1_word != op2_word)
1627     goto diff;
1628     }
1629     return 0;
1630    
1631     diff:
1632     //This can *not* be simplified to
1633     // op2_word - op2_word
1634     //since that expression might give signed overflow.
1635     return (op1_word > op2_word) ? 1 : -1;
1636     }
1637    
1638    
1639     //07/15/01: Am willing to skip unit-testing on this.
1640     //Understand the logic (i.e. passes visual inspection),
1641     //and comes from GNU-MP. Hope any defects here will be
1642     //caught in testing of GMP_INTS_mpz_mul() and other
1643     //higher-level functions.
1644     void GMP_INTS_mpn_mul_basecase (GMP_INTS_limb_ptr prodp,
1645     GMP_INTS_limb_srcptr up,
1646     GMP_INTS_size_t usize,
1647     GMP_INTS_limb_srcptr vp,
1648     GMP_INTS_size_t vsize)
1649     {
1650     assert(prodp != NULL);
1651     assert(up != NULL);
1652     assert(usize > 0);
1653     assert(vp != NULL);
1654     assert(vsize > 0);
1655    
1656     /* We first multiply by the low order one or two limbs, as the result can
1657     be stored, not added, to PROD. We also avoid a loop for zeroing this
1658     way. */
1659     prodp[usize] = GMP_INTS_mpn_mul_1 (prodp, up, usize, vp[0]);
1660     prodp++;
1661     vp++;
1662     vsize--;
1663    
1664     /* For each iteration in the loop, multiply U with one limb from V, and
1665     add the result to PROD. */
1666     while (vsize != 0)
1667     {
1668     prodp[usize] = GMP_INTS_mpn_addmul_1 (prodp, up, usize, vp[0]);
1669     prodp++,
1670     vp++,
1671     vsize--;
1672     }
1673     }
1674    
1675    
1676     //07/15/01: No unit testing possible--this is a passthrough.
1677     //In the original GNU MP code, there were several multiplication
1678     //algorithms, and this function would select one based on the
1679     //size of the operands and other considerations. The code has been
1680     //pared so that only simple multiplication is used, which is why
1681     //this function contains only a single pass-thru function call.
1682     void GMP_INTS_mpn_mul_n (GMP_INTS_limb_ptr p,
1683     GMP_INTS_limb_srcptr a,
1684     GMP_INTS_limb_srcptr b,
1685     GMP_INTS_size_t n)
1686     {
1687     GMP_INTS_mpn_mul_basecase (p, a, n, b, n);
1688     }
1689    
1690    
1691     //07/16/01: Visual inspection OK. Will not perform unit testing.
1692     GMP_INTS_limb_t GMP_INTS_mpn_mul(GMP_INTS_limb_ptr prodp,
1693     GMP_INTS_limb_srcptr up,
1694     GMP_INTS_size_t un,
1695     GMP_INTS_limb_srcptr vp,
1696     GMP_INTS_size_t vn)
1697     {
1698     //This is a gutted version of the GNU MP function. The GNU
1699     //MP function considered the case of a square, and also
1700     //better algorithms that pay off with large operands.
1701     //This gutted version uses only basic multiplication
1702     //(O(N**2)).
1703    
1704     //Eyeball the input parameters.
1705     assert(prodp != NULL);
1706     assert(up != NULL);
1707     assert(un >= 0);
1708     assert(vp != NULL);
1709     assert(vn >= 0);
1710    
1711     /* Basic long multiplication. */
1712     GMP_INTS_mpn_mul_basecase (prodp, up, un, vp, vn);
1713    
1714     //Return the most significant limb (which might be zero).
1715     //This is different than
1716     //most other functions, which return the spillover.
1717     return prodp[un + vn - 1];
1718     }
1719    
1720    
1721     /******************************************************************/
1722     /*** LIMB SPACE REALLOCATION FUNCTIONS *************************/
1723     /******************************************************************/
1724    
1725     void *GMP_INTS_mpz_realloc (GMP_INTS_mpz_struct *m,
1726     GMP_INTS_size_t new_size)
1727     {
1728     /* Never allocate zero space. */
1729     if (new_size <= 0)
1730     new_size = 1;
1731    
1732     m->limbs = (GMP_INTS_limb_ptr)
1733     GMP_INTS_realloc_w_size (m->limbs,
1734     m->n_allocd * sizeof(GMP_INTS_limb_t),
1735     new_size * sizeof(GMP_INTS_limb_t));
1736     m->n_allocd = new_size;
1737    
1738     return (void *) m->limbs;
1739     }
1740    
1741    
1742     /******************************************************************/
1743     /*** PUBLIC INITIALIZATION AND MEMORY MANAGEMENT FUNCTIONS *****/
1744     /******************************************************************/
1745    
1746     void GMP_INTS_mpz_init (GMP_INTS_mpz_struct *x)
1747     {
1748     assert(x != NULL);
1749    
1750     //The structure (the header block) exists in the
1751     //caller's area. Most likely it is a local variable.
1752     //This is OK, because it doesn't take up much space.
1753    
1754     //Start off with no errors.
1755     x->flags = 0;
1756    
1757     //Allocate space for one limb, which is the most
1758     //basic amount. This will grow, almost certainly.
1759     x->limbs = GMP_INTS_malloc(sizeof(GMP_INTS_limb_t));
1760    
1761     //Indicate that one limb was allocated.
1762     x->n_allocd = 1;
1763    
1764     //Set the size to 0. This signals a value of zero.
1765     x->size = 0;
1766     }
1767    
1768    
1769     void GMP_INTS_mpz_clear (GMP_INTS_mpz_struct *x)
1770     {
1771     //Be sure the passed pointer is not NULL.
1772     assert(x != NULL);
1773    
1774     //Be sure that the amount allocated is also above zero.
1775     //Anything else represents a logical error.
1776     assert(x->n_allocd > 0);
1777    
1778     //Be sure that the pointer to the allocated limbs
1779     //is not NULL. Anything else would be a logical
1780     //error.
1781     assert(x->limbs != NULL);
1782    
1783     //Deallocate the space for the limbs. The pointer is
1784     //set NULL and the allocated amount set to zero
1785     // so in case clear is called again it will be
1786     //a detectable error.
1787     GMP_INTS_free_w_size(x->limbs,
1788     x->n_allocd * sizeof(GMP_INTS_limb_t));
1789     x->limbs = NULL;
1790     x->n_allocd = 0;
1791     }
1792    
1793    
1794     /******************************************************************/
1795     /*** PUBLIC ASSIGNMENT FUNCTIONS *******************************/
1796     /******************************************************************/
1797    
1798     void GMP_INTS_mpz_copy( GMP_INTS_mpz_struct *dst,
1799     const GMP_INTS_mpz_struct *src)
1800     {
1801     GMP_INTS_size_t i, n;
1802    
1803     //Eyeball the input parameters.
1804     assert(dst != NULL);
1805     assert(dst->n_allocd > 0);
1806     assert(dst->limbs != NULL);
1807     assert(src != NULL);
1808     assert(src->n_allocd > 0);
1809     assert(src->limbs != NULL);
1810    
1811     //Source and destination may not be the same.
1812     assert(src != dst);
1813    
1814     //Figure out the real size of the source. We need to take the absolute
1815     //value.
1816     n = GMP_INTS_abs_of_size_t(src->size);
1817    
1818     //Reallocate the destination to be bigger if necessary.
1819     if (dst->n_allocd < n)
1820     {
1821     GMP_INTS_mpz_realloc (dst, n);
1822     }
1823    
1824     //Copy the non-dynamic fields in the header.
1825     dst->flags = src->flags;
1826     dst->size = src->size;
1827    
1828     //Copy the limbs.
1829     for (i=0; i<n; i++)
1830     dst->limbs[i] = src->limbs[i];
1831     }
1832    
1833    
1834     void GMP_INTS_mpz_set_ui (GMP_INTS_mpz_struct *dest,
1835     unsigned long int val)
1836     {
1837     assert(dest != NULL);
1838    
1839     /* We don't check if the allocation is enough, since the rest of the
1840     package ensures it's at least 1, which is what we need here. */
1841    
1842     dest->flags = 0;
1843     //A set operation resets any errors.
1844    
1845     if (val > 0)
1846     {
1847     dest->limbs[0] = val;
1848     dest->size = 1;
1849     }
1850     else
1851     {
1852     dest->size = 0;
1853     }
1854     }
1855    
1856    
1857     void GMP_INTS_mpz_set_si (GMP_INTS_mpz_struct *dest,
1858     signed long int val)
1859     {
1860     assert(dest != NULL);
1861    
1862     /* We don't check if the allocation is enough, since the rest of the
1863     package ensures it's at least 1, which is what we need here. */
1864    
1865     dest->flags = 0;
1866     //A set operation resets any errors.
1867    
1868     if (val > 0)
1869     {
1870     dest->limbs[0] = val;
1871     dest->size = 1;
1872     }
1873     else if (val < 0)
1874     {
1875     dest->limbs[0] = (unsigned long) -val;
1876     dest->size = -1;
1877     }
1878     else
1879     {
1880     dest->size = 0;
1881     }
1882     }
1883    
1884    
1885     void GMP_INTS_mpz_set_simple_char_str(GMP_INTS_mpz_struct *z,
1886     const char *s)
1887     {
1888     int sign=1;
1889     int digval;
1890     GMP_INTS_mpz_struct digvalz, k10;
1891    
1892     //Eyeball the arguments.
1893     assert(z != NULL);
1894     assert(z->n_allocd > 0);
1895     assert(z->limbs != NULL);
1896     assert(s != NULL);
1897    
1898     //Set the arbitrary integer to zero. This will also kill
1899     //any error flags.
1900     GMP_INTS_mpz_set_ui(z, 0);
1901    
1902     //Allocate an integer for our private use to hold each digit
1903     //value.
1904     GMP_INTS_mpz_init(&digvalz);
1905    
1906     //Allocate the constant 10, which we will use often.
1907     GMP_INTS_mpz_init(&k10);
1908     GMP_INTS_mpz_set_ui(&k10, 10);
1909    
1910     //As long as there are are digits and no flags set, keep
1911     //multiplying and adding the value of the digit. Non-
1912     //digits are simply ignored.
1913     while (!(z->flags) && (*s))
1914     {
1915     if (*s == '-')
1916     {
1917     sign = -sign;
1918     }
1919     else
1920     {
1921     digval = CHARFUNC_digit_to_val(*s);
1922     if (digval >= 0)
1923     {
1924     GMP_INTS_mpz_set_ui(&digvalz, digval);
1925     GMP_INTS_mpz_mul(z, z, &k10);
1926     GMP_INTS_mpz_add(z, z, &digvalz);
1927     }
1928     }
1929     s++;
1930     }
1931    
1932     //Adjust the final sign of the result.
1933     if (sign < 0)
1934     z->size = -(z->size);
1935    
1936     //Deallocate our temporary integers.
1937     GMP_INTS_mpz_clear(&digvalz);
1938     GMP_INTS_mpz_clear(&k10);
1939     }
1940    
1941    
1942     void GMP_INTS_mpz_set_sci_not_num(GMP_INTS_mpz_struct *z,
1943     int *failure,
1944     const char *s)
1945     {
1946     int parse_failure;
1947     //Return code from the floating point parsing
1948     //function.
1949     char mant_sign;
1950     //Sign character, if any, from the mantissa,
1951     //or N otherwise.
1952     size_t mant_bdp;
1953     //The index to the start of the mantissa before
1954     //the decimal point.
1955     size_t mant_bdp_len;
1956     //The length of the mantissa before the decimal
1957     //point. Zero means not defined, i.e. that
1958     //no characters were parsed and interpreted as
1959     //that part of a floating point number.
1960     size_t mant_adp;
1961     size_t mant_adp_len;
1962     //Similar fields for after the decimal point.
1963     char exp_sign;
1964     //Sign of the exponent, if any, or N otherwise.
1965     size_t exp;
1966     size_t exp_len;
1967     //Similar fields as to the mantissa, but for the
1968     //exponent.
1969     size_t si;
1970     //Iteration variable.
1971     int exponent_val;
1972     //The value of the exponent. We can't accept
1973     //an exponent outside the range of a 24-bit
1974     //signed integer. The 24-bit limit is arbitrary.
1975     //For one thing, it gives room to detect overflow
1976     //as are adding and multiplying by 10.
1977    
1978     //Eyeball the input parameters.
1979     assert(z != NULL);
1980     assert(z->n_allocd > 0);
1981     assert(z->limbs != NULL);
1982     assert(failure != NULL);
1983     assert(s != NULL);
1984    
1985     //Start off believing no failure.
1986     *failure = 0;
1987    
1988     //Set the output to zero. This is the default case for some
1989     //steps below.
1990     GMP_INTS_mpz_set_ui(z, 0);
1991    
1992     //Attempt to parse the number as a general number
1993     //in scientific notation.
1994     BSTRFUNC_parse_gen_sci_not_num(s,
1995     &parse_failure,
1996     &mant_sign,
1997     &mant_bdp,
1998     &mant_bdp_len,
1999     &mant_adp,
2000     &mant_adp_len,
2001     &exp_sign,
2002     &exp,
2003     &exp_len);
2004    
2005     //If it wouldn't parse as a general number, can't go further.
2006     if (parse_failure)
2007     {
2008     *failure = 1;
2009     return;
2010     }
2011     else if (!exp_len && !mant_adp_len)
2012     {
2013     //There was no exponent, and no portion after
2014     //the decimal point. Can just parse as an integer.
2015     char *temp_buf;
2016    
2017     //Allocate the temporary buffer to be one character longer
2018     //than the length specified for the parsed mantissa.
2019     temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
2020    
2021     //Copy from the parsed area into the temporary buffer.
2022     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
2023     temp_buf[si-mant_bdp] = s[si];
2024     temp_buf[mant_bdp_len] = 0;
2025    
2026     //Set the arbitrary integer to the value of the character
2027     //string.
2028     GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
2029    
2030     //If the number parsed as negative, invert.
2031     if (mant_sign == '-')
2032     z->size = -z->size;
2033    
2034     //Deallocate the temporary buffer.
2035     GMP_INTS_free(temp_buf);
2036     }
2037     else if (!exp_len && mant_adp_len)
2038     {
2039     char *temp_buf;
2040    
2041     //In this case, there are digits after the decimal point,
2042     //but no exponent specified. The only way this makes
2043     //sense is if all of the digits are zero--otherwise it
2044     //cannot be an integer.
2045     for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
2046     {
2047     if (s[si] != '0')
2048     {
2049     *failure = 1;
2050     return;
2051     }
2052     }
2053    
2054     //We're clean. They are only zeros. Execute as per
2055     //integer code.
2056    
2057     //Allocate the temporary buffer to be one character longer
2058     //than the length specified for the parsed mantissa.
2059     temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
2060    
2061     //Copy from the parsed area into the temporary buffer.
2062     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
2063     temp_buf[si-mant_bdp] = s[si];
2064     temp_buf[mant_bdp_len] = 0;
2065    
2066     //Set the arbitrary integer to the value of the character
2067     //string.
2068     GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
2069    
2070     //If the number parsed as negative, invert.
2071     if (mant_sign == '-')
2072     z->size = -z->size;
2073    
2074     //Deallocate the temporary buffer.
2075     GMP_INTS_free(temp_buf);
2076     }
2077     else if (exp_len)
2078     {
2079     //This is the most difficult case, where an exponent
2080     //is specified. There are several complex subcases,
2081     //such as:
2082     // a)If the exponent is too positive or too negative,
2083     // we can't use it. In general, we won't tackle
2084     // an exponent that won't fit in a signed 24-bit
2085     // integer. This provides a range of from
2086     // -8,388,608 to +8,388,607. This dwarfs the
2087     // 100,000 or so digit preprocessor limit,
2088     // and should be adequate for any practical
2089     // application.
2090     // b)If the exponent is zero, we ignore it.
2091     // c)If the exponent is positive, it has to
2092     // be large enough to overcome any
2093     // digits past the decimal point, otherwise
2094     // we don't end up with an integer.
2095     // d)If the exponent is negative, there have to
2096     // be enough digits so that an integer remains
2097     // after the exponent is applied. This
2098     // generally requires trailing zeros on the
2099     // part before the decimal point.
2100    
2101     //First, tackle the exponent. Process the
2102     //exponent into a signed integer. We have to
2103     //balk at anything outside of 24 bits. The
2104     //procedure used automatically handles
2105     //leading zeros correctly.
2106     exponent_val = 0;
2107     for (si=exp; si<(exp+exp_len); si++)
2108     {
2109     int val;
2110    
2111     val = CHARFUNC_digit_to_val(s[si]);
2112    
2113     assert(val >= 0 && val <= 9);
2114    
2115     exponent_val *= 10;
2116     exponent_val += val;
2117    
2118     if (((exp_sign=='-') && (exponent_val>8388608))
2119     ||
2120     ((exp_sign != '-') && (exponent_val>8388607)))
2121     {
2122     *failure = 1;
2123     return;
2124     }
2125     }
2126    
2127     //If we're here, the exponent has been computed and
2128     //is within 24 bits. However, we need to adjust for
2129     //the sign.
2130     if (exp_sign == '-')
2131     exponent_val = -exponent_val;
2132    
2133     //We need to make accurate assertions about the
2134     //portion of the number, if any, after the decimal point.
2135     //This means that we need to effectively discard
2136     //trailing zeros. To do this, we do not need to
2137     //relocate the string, we can just back off the index
2138     //to bypass any trailing zeros.
2139     while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))
2140     mant_adp_len--;
2141    
2142     //We also need to make accurate assertions about the
2143     //portion of the number, if any, before the decimal
2144     //point. It is known that the parsing function
2145     //isn't tolerant of multiple zeros, but zero is a
2146     //special case. Let's advance the pointer to the
2147     //part before the decimal point so that zero will
2148     //have zero length.
2149     while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))
2150     {
2151     mant_bdp++;
2152     mant_bdp_len--;
2153     }
2154    
2155     //If we are dealing with zero, who cares about the
2156     //exponent? Just return the value of zero.
2157     if (!mant_bdp_len && !mant_adp_len)
2158     {
2159     *failure = 0;
2160     GMP_INTS_mpz_set_ui(z, 0);
2161     return;
2162     }
2163    
2164     //Beyond this point, we have something non-zero.
2165     //If the exponent is positive, it must be at least
2166     //as large as the number of digits beyond the
2167     //decimal point in order to form an integer. If the
2168     //exponent is zero, there must be no digits after the
2169     //decimal point. If the exponent is negative, there
2170     //must be no digits after the decimal point, and the
2171     //trailing zeros on the part before the decimal point
2172     //must be adequate to handle the right decimal shift.
2173     if (exponent_val == 0)
2174     {
2175     if (mant_adp_len)
2176     {
2177     *failure = 1;
2178     return;
2179     }
2180     }
2181     else if (exponent_val > 0)
2182     {
2183     if ((int)mant_adp_len > exponent_val)
2184     {
2185     *failure = 1;
2186     return;
2187     }
2188     }
2189     else //exponent_val < 0
2190     {
2191     if (mant_adp_len)
2192     {
2193     *failure = 1;
2194     return;
2195     }
2196     else
2197     {
2198     //Count the number of trailing zeros on the part
2199     //before the decimal point.
2200     size_t trailing_zero_count;
2201     int idx;
2202    
2203     trailing_zero_count = 0;
2204    
2205     for(idx = mant_bdp + mant_bdp_len - 1;
2206     (mant_bdp_len != 0) && (idx >= (int)mant_bdp);
2207     idx--)
2208     {
2209     if (s[idx] == '0')
2210     trailing_zero_count++;
2211     else
2212     break;
2213     }
2214    
2215     //Check on the assertion about trailing zeros.
2216     if ((int)trailing_zero_count < -exponent_val)
2217     {
2218     *failure = 1;
2219     return;
2220     }
2221     }
2222     }
2223    
2224     {
2225     //Create a string long enough to hold the digits
2226     //before the decimal point plus the ones after and
2227     //convert that to an arbitrary integer.
2228     //Form a power of 10 which is 10 exponentiated to
2229     //the absolute value of the exponent. If the
2230     //exponent was positive, multiply by it. If the
2231     //exponent was negative, divide by it.
2232     char *conv_str;
2233     size_t sidx;
2234     GMP_INTS_mpz_struct power_of_ten, k10, trash;
2235    
2236     GMP_INTS_mpz_init(&power_of_ten);
2237     GMP_INTS_mpz_init(&k10);
2238     GMP_INTS_mpz_init(&trash);
2239    
2240     conv_str = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + mant_adp_len + 1));
2241    
2242     sidx=0;
2243    
2244     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
2245     {
2246     conv_str[sidx] = s[si];
2247     sidx++;
2248     }
2249     for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
2250     {
2251     conv_str[sidx] = s[si];
2252     sidx++;
2253     }
2254     conv_str[sidx] = 0;
2255    
2256     assert(sidx == (mant_bdp_len + mant_adp_len));
2257    
2258     GMP_INTS_mpz_set_simple_char_str(z, conv_str);
2259    
2260     GMP_INTS_mpz_set_ui(&k10, 10);
2261    
2262     if (exponent_val > 0)
2263     GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, exponent_val-mant_adp_len);
2264     else
2265     GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, -exponent_val);
2266    
2267     if (exponent_val >= 0)
2268     {
2269     GMP_INTS_mpz_mul(z, z, &power_of_ten);
2270     }
2271     else
2272     {
2273     GMP_INTS_mpz_tdiv_qr (&k10,
2274     &trash,
2275     z,
2276     &power_of_ten);
2277     GMP_INTS_mpz_copy(z, &k10);
2278     }
2279    
2280     //If the argument had a minus sign, invert.
2281     if (mant_sign == '-')
2282     z->size = -z->size;
2283    
2284     GMP_INTS_free(conv_str);
2285    
2286     GMP_INTS_mpz_clear(&trash);
2287     GMP_INTS_mpz_clear(&k10);
2288     GMP_INTS_mpz_clear(&power_of_ten);
2289    
2290     //Finally, if the arbitrary integer has overflowed, this is
2291     //a parse failure. Must declare as such.
2292     if (z->flags)
2293     *failure = 1;
2294     }
2295     }
2296     else
2297     {
2298     *failure = 1;
2299     return;
2300     }
2301     }
2302    
2303    
2304     void GMP_INTS_mpz_set_general_int(GMP_INTS_mpz_struct *z,
2305     int *failure,
2306     const char *s)
2307     {
2308     //Eyeball the input parameters.
2309     assert(z != NULL);
2310     assert(z->n_allocd > 0);
2311     assert(z->limbs != NULL);
2312     assert(failure != NULL);
2313     assert(s != NULL);
2314    
2315     //Try to parse it as a simple integer.
2316     if (BSTRFUNC_is_sint_wo_commas(s))
2317     {
2318     GMP_INTS_mpz_set_simple_char_str(z, s);
2319     *failure = 0;
2320     return;
2321     }
2322     //If that didn't work, try to parse it as a simple
2323     //integer with commas.
2324     else if (BSTRFUNC_is_sint_w_commas(s))
2325     {
2326     GMP_INTS_mpz_set_simple_char_str(z, s);
2327     *failure = 0;
2328     return;
2329     }
2330     //If neither of those worked, try to parse it as
2331     //something containing scientific notation.
2332     else
2333     {
2334     GMP_INTS_mpz_set_sci_not_num(z, failure, s);
2335    
2336     if (!*failure)
2337     {
2338     //We were able to parse it that way.
2339     //Everything is set up, just return.
2340     return;
2341     }
2342     else
2343     {
2344     //We're out of options. All parsing failed.
2345     GMP_INTS_mpz_set_ui(z, 0);
2346     *failure = 1;
2347     return;
2348     }
2349     }
2350     }
2351    
2352    
2353     void GMP_INTS_mpz_parse_into_uint32(unsigned *result,
2354     int *failure,
2355     char *s)
2356     {
2357     GMP_INTS_mpz_struct arb_int;
2358    
2359     //Eyeball the input parameters.
2360     assert(result != NULL);
2361     assert(failure != NULL);
2362     assert(s != NULL);
2363    
2364     //Allocate space for the one arbitrary integer we need.
2365     GMP_INTS_mpz_init(&arb_int);
2366    
2367     //Try to parse the string into an arbitrary length integer
2368     //using all methods known to man.
2369     GMP_INTS_mpz_set_general_int(&arb_int, failure, s);
2370    
2371     //If the parse failed, we must declare failure and return
2372     //0.
2373     if (*failure)
2374     {
2375     *result = 0;
2376     *failure = 1;
2377     }
2378     else
2379     {
2380     //We might have success, but it might be negative or
2381     //too big.
2382     if (arb_int.size == 1)
2383     {
2384     *result = arb_int.limbs[0];
2385     *failure = 0;
2386     }
2387     else if (arb_int.size == 0)
2388     {
2389     *result = 0;
2390     *failure = 0;
2391     }
2392     else
2393     {
2394     *result = 0;
2395     *failure = 1;
2396     }
2397     }
2398    
2399     //Deallocate the arbitrary integer.
2400     GMP_INTS_mpz_clear(&arb_int);
2401     }
2402    
2403    
2404     void GMP_INTS_mpz_swap(GMP_INTS_mpz_struct *a,
2405     GMP_INTS_mpz_struct *b)
2406     {
2407     GMP_INTS_mpz_struct temp;
2408    
2409     //Eyeball the input parameters.
2410     assert(a != NULL);
2411     assert(a->n_allocd > 0);
2412     assert(a->limbs != NULL);
2413     assert(b != NULL);
2414     assert(b->n_allocd > 0);
2415     assert(b->limbs != NULL);
2416    
2417     //Make the swap via memory copy.
2418     memcpy(&temp, a, sizeof(GMP_INTS_mpz_struct));
2419     memcpy(a, b, sizeof(GMP_INTS_mpz_struct));
2420     memcpy(b, &temp, sizeof(GMP_INTS_mpz_struct));
2421     }
2422    
2423    
2424     /******************************************************************/
2425     /*** PUBLIC ARITHMETIC FUNCTIONS *******************************/
2426     /******************************************************************/
2427    
2428     //07/15/01: Unit test and visual inspection passed.
2429     void GMP_INTS_mpz_add ( GMP_INTS_mpz_struct *w,
2430     const GMP_INTS_mpz_struct *u,
2431     const GMP_INTS_mpz_struct *v)
2432     {
2433     GMP_INTS_limb_srcptr up, vp;
2434     GMP_INTS_limb_ptr wp;
2435     GMP_INTS_size_t usize, vsize, wsize;
2436     GMP_INTS_size_t abs_usize;
2437     GMP_INTS_size_t abs_vsize;
2438    
2439     //Look at the input parameters carefully.
2440     assert(w != NULL);
2441     assert(u != NULL);
2442     assert(v != NULL);
2443     assert(w->n_allocd > 0);
2444     assert(u->n_allocd > 0);
2445     assert(v->n_allocd > 0);
2446     assert(w->limbs != NULL);
2447     assert(u->limbs != NULL);
2448     assert(v->limbs != NULL);
2449    
2450     //Handle the case of a tainted result. If either of the
2451     //two inputs are either direct overflows or tainted by
2452     //an overflow, mark the result tainted and do not perform
2453     //any arithmetic operation.
2454     {
2455     int taint;
2456    
2457     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
2458    
2459     w->flags = 0;
2460     //"w" starts off with a clean slate. Must do this
2461     //after taint calculation in case locations of u or v
2462     //are the same as w.
2463     if (taint)
2464     {
2465     w->flags = taint;
2466     return;
2467     }
2468     }
2469    
2470     usize = u->size;
2471     vsize = v->size;
2472     abs_usize = GMP_INTS_abs_of_size_t(usize);
2473     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
2474    
2475     //Arrange things so that U has at least as many
2476     //limbs as V, i.e. limbs(U) >= limbs(V);
2477     if (abs_usize < abs_vsize)
2478     {
2479     const GMP_INTS_mpz_struct *tmp_ptr;
2480     GMP_INTS_size_t tmp_size;
2481    
2482     //Swap U and V. This does no harm, because we are
2483     //manipulating only local variables. This does not
2484     //affect the caller.
2485     tmp_ptr = u;
2486     u = v;
2487     v = tmp_ptr;
2488     tmp_size = usize;
2489     usize = vsize;
2490     vsize = tmp_size;
2491     tmp_size = abs_usize;
2492     abs_usize = abs_vsize;
2493     abs_vsize = tmp_size;
2494     }
2495    
2496     /* True: ABS_USIZE >= ABS_VSIZE. */
2497    
2498     /* If not space for w (and possible carry), increase space. */
2499     wsize = abs_usize + 1;
2500     if (w->n_allocd < wsize)
2501     GMP_INTS_mpz_realloc(w, wsize);
2502    
2503     //These pointers must be obtained after realloc. At this point,
2504     //u or v may be the same as w.
2505     up = u->limbs;
2506     vp = v->limbs;
2507     wp = w->limbs;
2508    
2509     if ((usize ^ vsize) < 0)
2510     {
2511     //U and V have different sign. Need to compare them to determine
2512     //which operand to subtract from which.
2513    
2514     //This test is right since ABS_USIZE >= ABS_VSIZE.
2515     //If the equality case is ruled out, then U has more limbs
2516     //than V, which means that it is bigger in magnitude.
2517     if (abs_usize != abs_vsize)
2518     {
2519     GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
2520     wsize = abs_usize;
2521    
2522     //Normalize the result. This was formerly a macro.
2523     //To normalize in this context means to trim the size
2524     //down to eliminate any leading zero limbs that came
2525     //about because the size of the result of an operation
2526     //was overestimated.
2527     GMP_INTS_mpn_normalize(wp, &wsize);
2528    
2529     if (usize < 0)
2530     wsize = -wsize;
2531     }
2532     else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
2533     {
2534     GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
2535     wsize = abs_usize;
2536     GMP_INTS_mpn_normalize(wp, &wsize);
2537     if (usize >= 0)
2538     wsize = -wsize;
2539     }
2540     else
2541     {
2542     GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
2543     wsize = abs_usize;
2544     GMP_INTS_mpn_normalize(wp, &wsize);
2545     if (usize < 0)
2546     wsize = -wsize;
2547     }
2548     }
2549     else
2550     {
2551     /* U and V have same sign. Add them. */
2552     GMP_INTS_limb_t cy_limb
2553     = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
2554     wp[abs_usize] = cy_limb;
2555     wsize = abs_usize + cy_limb;
2556     if (usize < 0)
2557     wsize = -wsize;
2558     }
2559    
2560     w->size = wsize;
2561    
2562     //Handle the case of an overflowed result. If the result
2563     //of the addition is too big or too small, mark it as
2564     //overflowed.
2565     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2566     {
2567     w->flags = GMP_INTS_EF_INTOVF_POS;
2568     }
2569     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2570     {
2571     w->flags = GMP_INTS_EF_INTOVF_NEG;
2572     }
2573     }
2574    
2575    
2576     //07/15/01: Unit testing skipped because of recursive
2577     //nature. Visual inspection OK.
2578     void GMP_INTS_mpz_add_ui ( GMP_INTS_mpz_struct *w,
2579     const GMP_INTS_mpz_struct *u,
2580     unsigned long int v)
2581     {
2582     //The GNU MP version of this is quite efficient, and this
2583     //makes sense since it is a common operation. However,
2584     //for simplicity just define this recursively in terms
2585     //of the ADD function. This can always be made quicker
2586     //later (by changing back to the GNU MP version).
2587     GMP_INTS_mpz_struct temp;
2588    
2589     //Eyeball the inputs carefully.
2590     assert(w != NULL);
2591     assert(w->n_allocd > 0);
2592     assert(w->limbs != NULL);
2593     assert(u != NULL);
2594     assert(u->n_allocd > 0);
2595     assert(u->limbs != NULL);
2596    
2597     //Create a temporary integer.
2598     GMP_INTS_mpz_init(&temp);
2599    
2600     //Set the temporary integer to the value of the input
2601     //argument.
2602     GMP_INTS_mpz_set_ui(&temp, v);
2603    
2604     //Do the actual addition. This recursive definition
2605     //is inherently wasteful, but I'm after clarity, not
2606     //warp speed.
2607     GMP_INTS_mpz_add(w, u, &temp);
2608    
2609     //Destroy the temporary integer (this will reclaim the
2610     //memory).
2611     GMP_INTS_mpz_clear(&temp);
2612     }
2613    
2614    
2615     //07/15/01: Visual inspection passed. Not unit tested
2616     //because of symmetry with GMP_INTS_mpz_add().
2617     void GMP_INTS_mpz_sub ( GMP_INTS_mpz_struct *w,
2618     const GMP_INTS_mpz_struct *u,
2619     const GMP_INTS_mpz_struct *v)
2620     {
2621     GMP_INTS_limb_srcptr up, vp;
2622     GMP_INTS_limb_ptr wp;
2623     GMP_INTS_size_t usize, vsize, wsize;
2624     GMP_INTS_size_t abs_usize;
2625     GMP_INTS_size_t abs_vsize;
2626    
2627     //Look at the input parameters carefully.
2628     assert(w != NULL);
2629     assert(u != NULL);
2630     assert(v != NULL);
2631     assert(w->n_allocd > 0);
2632     assert(u->n_allocd > 0);
2633     assert(v->n_allocd > 0);
2634     assert(w->limbs != NULL);
2635     assert(u->limbs != NULL);
2636     assert(v->limbs != NULL);
2637    
2638     //Handle the case of a tainted result. If either of the
2639     //two inputs are either direct overflows or tainted by
2640     //an overflow, mark the result tainted and do not perform
2641     //any arithmetic operation.
2642     {
2643     int taint;
2644    
2645     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
2646    
2647     w->flags = 0;
2648     //"w" starts off with a clean slate. Must do this
2649     //after taint calculation in case locations of u or v
2650     //are the same as w.
2651     if (taint)
2652     {
2653     w->flags = taint;
2654     return;
2655     }
2656     }
2657    
2658     usize = u->size;
2659     vsize = -(v->size); /* The "-" makes the difference from mpz_add */
2660     abs_usize = GMP_INTS_abs_of_size_t(usize);
2661     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
2662    
2663     if (abs_usize < abs_vsize)
2664     {
2665     const GMP_INTS_mpz_struct *tmp_ptr;
2666     GMP_INTS_size_t tmp_size;
2667    
2668     //Swap U and V. This does no harm, because we are
2669     //manipulating only local variables. This does not
2670     //affect the caller.
2671     tmp_ptr = u;
2672     u = v;
2673     v = tmp_ptr;
2674     tmp_size = usize;
2675     usize = vsize;
2676     vsize = tmp_size;
2677     tmp_size = abs_usize;
2678     abs_usize = abs_vsize;
2679     abs_vsize = tmp_size;
2680     }
2681    
2682     /* True: ABS_USIZE >= ABS_VSIZE. */
2683    
2684     /* If not space for w (and possible carry), increase space. */
2685     wsize = abs_usize + 1;
2686     if (w->n_allocd < wsize)
2687     GMP_INTS_mpz_realloc (w, wsize);
2688    
2689     /* These must be after realloc (u or v may be the same as w). */
2690     up = u->limbs;
2691     vp = v->limbs;
2692     wp = w->limbs;
2693    
2694     if ((usize ^ vsize) < 0)
2695     {
2696     //U and V have different sign. Need to compare them to determine
2697     //which operand to subtract from which.
2698    
2699     //This test is right since ABS_USIZE >= ABS_VSIZE.
2700     if (abs_usize != abs_vsize)
2701     {
2702     GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
2703     wsize = abs_usize;
2704     GMP_INTS_mpn_normalize(wp, &wsize);
2705     if (usize < 0)
2706     wsize = -wsize;
2707     }
2708     else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
2709     {
2710     GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
2711     wsize = abs_usize;
2712     GMP_INTS_mpn_normalize(wp, &wsize);
2713     if (usize >= 0)
2714     wsize = -wsize;
2715     }
2716     else
2717     {
2718     GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
2719     wsize = abs_usize;
2720     GMP_INTS_mpn_normalize (wp, &wsize);
2721     if (usize < 0)
2722     wsize = -wsize;
2723     }
2724     }
2725     else
2726     {
2727     /* U and V have same sign. Add them. */
2728     GMP_INTS_limb_t cy_limb
2729     = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
2730     wp[abs_usize] = cy_limb;
2731     wsize = abs_usize + cy_limb;
2732     if (usize < 0)
2733     wsize = -wsize;
2734     }
2735    
2736     w->size = wsize;
2737    
2738     //Handle the case of an overflowed result. If the result
2739     //of the addition is too big or too small, mark it as
2740     //overflowed.
2741     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2742     {
2743     w->flags = GMP_INTS_EF_INTOVF_POS;
2744     }
2745     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2746     {
2747     w->flags = GMP_INTS_EF_INTOVF_NEG;
2748     }
2749     }
2750    
2751    
2752     //07/15/01: Unit testing skipped because of recursive
2753     //nature. Visual inspection OK.
2754     void GMP_INTS_mpz_sub_ui ( GMP_INTS_mpz_struct *w,
2755     const GMP_INTS_mpz_struct *u,
2756     unsigned long int v)
2757     {
2758     //The GNU MP version of this is quite efficient, and this
2759     //makes sense since it is a common operation. However,
2760     //for simplicity just define this recursively in terms
2761     //of the SUB function. This can always be made quicker
2762     //later (by changing back to the GNU MP version).
2763     GMP_INTS_mpz_struct temp;
2764    
2765     //Eyeball the inputs carefully.
2766     assert(w != NULL);
2767     assert(w->n_allocd > 0);
2768     assert(w->limbs != NULL);
2769     assert(u != NULL);
2770     assert(u->n_allocd > 0);
2771     assert(u->limbs != NULL);
2772    
2773     //Create a temporary integer.
2774     GMP_INTS_mpz_init(&temp);
2775    
2776     //Set the temporary integer to the value of the input
2777     //argument.
2778     GMP_INTS_mpz_set_ui(&temp, v);
2779    
2780     //Do the actual subtraction. This recursive definition
2781     //is inherently wasteful, but I'm after clarity, not
2782     //warp speed.
2783     GMP_INTS_mpz_sub(w, u, &temp);
2784    
2785     //Destroy the temporary integer (this will reclaim the
2786     //memory).
2787     GMP_INTS_mpz_clear(&temp);
2788     }
2789    
2790    
2791     void GMP_INTS_mpz_mul ( GMP_INTS_mpz_struct *w,
2792     const GMP_INTS_mpz_struct *u,
2793     const GMP_INTS_mpz_struct *v)
2794     {
2795     GMP_INTS_size_t usize = u->size;
2796     GMP_INTS_size_t vsize = v->size;
2797     GMP_INTS_size_t wsize;
2798     GMP_INTS_size_t sign_product;
2799     GMP_INTS_limb_ptr up, vp;
2800     GMP_INTS_limb_ptr wp;
2801     GMP_INTS_limb_ptr free_me = NULL;
2802     GMP_INTS_size_t free_me_size;
2803     GMP_INTS_limb_t cy_limb;
2804    
2805     //Eyeball the inputs.
2806     assert(w != NULL);
2807     assert(w->n_allocd > 0);
2808     assert(w->limbs != NULL);
2809     assert(u != NULL);
2810     assert(u->n_allocd > 0);
2811     assert(u->limbs != NULL);
2812     assert(v != NULL);
2813     assert(v->n_allocd > 0);
2814     assert(v->limbs != NULL);
2815    
2816     //Handle the case of a tainted result. If either of the
2817     //two inputs are either direct overflows or tainted by
2818     //an overflow, mark the result tainted and do not perform
2819     //any arithmetic operation.
2820     {
2821     int taint;
2822    
2823     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
2824    
2825     w->flags = 0;
2826     //"w" starts off with a clean slate. Must do this
2827     //after taint calculation in case locations of u or v
2828     //are the same as w.
2829     if (taint)
2830     {
2831     w->flags = taint;
2832     return;
2833     }
2834     }
2835    
2836     sign_product = usize ^ vsize;
2837     usize = GMP_INTS_abs_of_size_t(usize);
2838     vsize = GMP_INTS_abs_of_size_t(vsize);
2839    
2840     //Handle the case of a certain result overflow (why do the math when
2841     //the result is certain?). In general, when multiplying two inputs
2842     //whose sizes are M limbs and N limbs, the size of the result will be
2843     //either M+N or M+N-1 limbs. If M+N-1 > MAX_ALLOWED, then can declare
2844     //an early overflow.
2845     if ((usize + vsize - 1) > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2846     {
2847     if (sign_product < 0)
2848     w->flags = GMP_INTS_EF_INTOVF_NEG;
2849     else
2850     w->flags = GMP_INTS_EF_INTOVF_POS;
2851    
2852     return;
2853     }
2854    
2855    
2856     if (usize < vsize)
2857     {
2858     //Temporary variables just for the swap.
2859     const GMP_INTS_mpz_struct *tmp_ptr;
2860     GMP_INTS_size_t tmp_size;
2861    
2862     //Swap U and V.
2863     tmp_ptr = u;
2864     u = v;
2865     v = tmp_ptr;
2866     tmp_size = usize;
2867     usize = vsize;
2868     vsize = tmp_size;
2869     }
2870    
2871     //Grab pointers to the arrays of limbs.
2872     up = u->limbs;
2873     vp = v->limbs;
2874     wp = w->limbs;
2875    
2876     /* Ensure W has space enough to store the result. */
2877     wsize = usize + vsize;
2878     if (w->n_allocd < wsize)
2879     {
2880     if (wp == up || wp == vp)
2881     {
2882     free_me = wp;
2883     free_me_size = w->n_allocd;
2884     }
2885     else
2886     {
2887     GMP_INTS_free_w_size (wp, w->n_allocd * sizeof(GMP_INTS_limb_t));
2888     }
2889    
2890     w->n_allocd = wsize;
2891     wp = (GMP_INTS_limb_ptr)
2892     GMP_INTS_malloc (wsize * sizeof(GMP_INTS_limb_t));
2893     w->limbs = wp;
2894     }
2895     else
2896     {
2897     /* Make U and V not overlap with W. */
2898     if (wp == up)
2899     {
2900     /* W and U are identical. Allocate temporary space for U. */
2901     up = (GMP_INTS_limb_ptr)
2902     _alloca(usize * sizeof(GMP_INTS_limb_t));
2903     /* Is V identical too? Keep it identical with U. */
2904     if (wp == vp)
2905     vp = up;
2906     /* Copy to the temporary space. */
2907     GMP_INTS_mpn_copy_limbs(up, wp, usize);
2908     }
2909     else if (wp == vp)
2910     {
2911     /* W and V are identical. Allocate temporary space for V. */
2912     vp = (GMP_INTS_limb_ptr)
2913     _alloca(vsize * sizeof(GMP_INTS_limb_t));
2914     /* Copy to the temporary space. */
2915     GMP_INTS_mpn_copy_limbs(vp, wp, vsize);
2916     }
2917     }
2918    
2919     if (vsize == 0)
2920     {
2921     wsize = 0;
2922     }
2923     else
2924     {
2925     cy_limb = GMP_INTS_mpn_mul (wp, up, usize, vp, vsize);
2926     wsize = usize + vsize;
2927     wsize -= cy_limb == 0;
2928     }
2929    
2930     w->size = sign_product < 0 ? -wsize : wsize;
2931    
2932     if (free_me != NULL)
2933     GMP_INTS_free_w_size (free_me, free_me_size * sizeof(GMP_INTS_limb_t));
2934    
2935     //Final check for overflow.
2936     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2937     w->flags = GMP_INTS_EF_INTOVF_POS;
2938     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2939     w->flags = GMP_INTS_EF_INTOVF_NEG;
2940     }
2941    
2942    
2943     //07/15/01: Unit testing skipped because of recursive
2944     //nature. Visual inspection OK.
2945     void GMP_INTS_mpz_mul_si ( GMP_INTS_mpz_struct *w,
2946     const GMP_INTS_mpz_struct *u,
2947     long int v)
2948     {
2949     GMP_INTS_mpz_struct temp;
2950    
2951     //Eyeball the inputs carefully.
2952     assert(w != NULL);
2953     assert(w->n_allocd > 0);
2954     assert(w->limbs != NULL);
2955     assert(u != NULL);
2956     assert(u->n_allocd > 0);
2957     assert(u->limbs != NULL);
2958    
2959     //Create a temporary integer.
2960     GMP_INTS_mpz_init(&temp);
2961    
2962     //Set the temporary integer to the value of the input
2963     //argument.
2964     GMP_INTS_mpz_set_si(&temp, v);
2965    
2966     //Do the actual multiplication. This recursive definition
2967     //is inherently wasteful, but I'm after clarity, not
2968     //warp speed.
2969     GMP_INTS_mpz_mul(w, u, &temp);
2970    
2971     //Destroy the temporary integer (this will reclaim the
2972     //memory).
2973     GMP_INTS_mpz_clear(&temp);
2974     }
2975    
2976    
2977     //07/15/01: Unit testing skipped because of recursive
2978     //nature. Visual inspection OK.
2979     void GMP_INTS_mpz_mul_ui ( GMP_INTS_mpz_struct *w,
2980     const GMP_INTS_mpz_struct *u,
2981     unsigned long int v)
2982     {
2983     GMP_INTS_mpz_struct temp;
2984    
2985     //Eyeball the inputs carefully.
2986     assert(w != NULL);
2987     assert(w->size >= 0);
2988     assert(w->limbs != NULL);
2989     assert(u != NULL);
2990     assert(u->size >= 0);
2991     assert(u->limbs != NULL);
2992    
2993     //Create a temporary integer.
2994     GMP_INTS_mpz_init(&temp);
2995    
2996     //Set the temporary integer to the value of the input
2997     //argument.
2998     GMP_INTS_mpz_set_ui(&temp, v);
2999    
3000     //Do the actual multiplication. This recursive definition
3001     //is inherently wasteful, but I'm after clarity, not
3002     //warp speed.
3003     GMP_INTS_mpz_mul(w, u, &temp);
3004    
3005     //Destroy the temporary integer (this will reclaim the
3006     //memory).
3007     GMP_INTS_mpz_clear(&temp);
3008     }
3009    
3010    
3011     void GMP_INTS_mpz_tdiv_qr ( GMP_INTS_mpz_struct *quot,
3012     GMP_INTS_mpz_struct *rem,
3013     const GMP_INTS_mpz_struct *num,
3014     const GMP_INTS_mpz_struct *den)
3015     {
3016     GMP_INTS_size_t abs_num_size,
3017     abs_den_size,
3018     quotient_sign,
3019     remainder_sign,
3020     numerator_bitsize,
3021     denominator_bitsize,
3022     division_loop_count,
3023     division_loop_count_mod_32,
3024     division_loop_count_div_32,
3025     division_counter,
3026     i;
3027     GMP_INTS_limb_t temp_limb;
3028     GMP_INTS_limb_ptr trial_divisor;
3029    
3030     //Eyeball the input parameters.
3031     assert(quot != NULL);
3032     assert(quot->n_allocd > 0);
3033     assert(quot->limbs != NULL);
3034     assert(rem != NULL);
3035     assert(rem->n_allocd > 0);
3036     assert(rem->limbs != NULL);
3037     assert(num != NULL);
3038     assert(num->n_allocd > 0);
3039     assert(num->limbs != NULL);
3040     assert(den != NULL);
3041     assert(den->n_allocd > 0);
3042     assert(den->limbs != NULL);
3043    
3044     //We require for this function that the numerator, denominator, quotient, and
3045     //remainder all be distinct.
3046     assert(quot != rem);
3047     assert(quot != num);
3048     assert(quot != den);
3049     assert(rem != num);
3050     assert(rem != den);
3051     assert(num != den);
3052    
3053     //The GNU code was probably very efficient, but exceeded
3054     //my abilities to analyze. This is the classic
3055     //division algorithm.
3056    
3057     //First, start off with the quotient and remainder having
3058     //no error flags set. These will be set if appropriate.
3059     quot->flags = 0;
3060     rem->flags = 0;
3061    
3062     //First, handle tainted inputs. If the numerator or denominator
3063     //are bad or tainted, the quotient and remainder get tainted
3064     //automatically.
3065     {
3066     int taint;
3067    
3068     taint = GMP_INTS_two_op_flags_map(num->flags, den->flags);
3069    
3070     if (taint)
3071     {
3072     quot->flags = taint;
3073     rem->flags = taint;
3074     return;
3075     }
3076     }
3077    
3078     //The second possible cause for taint is if the divisor is
3079     //zero. This will get both the value of positive overflow.
3080     if (den->size == 0)
3081     {
3082     quot->flags = GMP_INTS_EF_INTOVF_POS;
3083     rem->flags = GMP_INTS_EF_INTOVF_POS;
3084     return;
3085     }
3086    
3087     //Handle the special case of a numerator of zero. If the numerator
3088     //is zero, the quotient and remainder are zero automatically.
3089     if (num->size == 0)
3090     {
3091     GMP_INTS_mpz_set_ui(quot, 0);
3092     GMP_INTS_mpz_set_ui(rem, 0);
3093     return;
3094     }
3095    
3096     //Generally, nothing else can go wrong as far as taint. The
3097     //value of the quotient is confined to be no larger than the
3098     //numerator, and the value of the remainder is confined to
3099     //be no larger than denominator-1. So, generally, if the
3100     //inputs are in size bounds, the outputs will be also.
3101    
3102     //Figure out how large in limbs the numerator and denominator actually
3103     //are.
3104     abs_num_size = GMP_INTS_abs_of_size_t(num->size);
3105     abs_den_size = GMP_INTS_abs_of_size_t(den->size);
3106    
3107     //Figure out the sign of things. We want the following relationship
3108     //to be true:
3109     // num/den = quot + rem/den.
3110     //The way to achieve this is to assign the sign of the quotient in the traditional
3111     //way, then to assign the remainder to have the same sign as the numerator.
3112     quotient_sign = num->size ^ den->size;
3113     remainder_sign = num->size;
3114    
3115     //The remainder starts off with the absolute value of the numerator, and then
3116     //we subtract from it as part of the division loop.
3117     GMP_INTS_mpz_copy(rem, num);
3118     //We know after the copy that the amount of space allocated in the remainder
3119     //MUST be at least as large as the absolute value of the numerator. So from
3120     //this point forward we use the space.
3121     assert(rem->n_allocd >= abs_num_size);
3122    
3123     //Figure out the number of significant bits in the numerator and denominator.
3124     //This determines the loop count over which we do the shift division loop.
3125     numerator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_num_size;
3126    
3127     i = abs_num_size - 1;
3128    
3129     //We need to be extra careful here. One failure mode is that an integer
3130     //data structure is corrupted and the "size" field reflects limbs
3131     //that are zero. Need to watch that this kind of failure doesn't
3132     //cause memory access errors.
3133     assert(num->limbs[i] != 0);
3134     if (num->limbs[i] == 0)
3135     {
3136     quot->flags = GMP_INTS_EF_INTOVF_POS;
3137     rem->flags = GMP_INTS_EF_INTOVF_POS;
3138     return;
3139     }
3140    
3141     temp_limb = 0x80000000;
3142    
3143     while (((num->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
3144     {
3145     numerator_bitsize--;
3146     temp_limb >>= 1;
3147     }
3148    
3149     denominator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_den_size;
3150    
3151     i = abs_den_size - 1;
3152    
3153     //We need to be extra careful here. One failure mode is that an integer
3154     //data structure is corrupted and the "size" field reflects limbs
3155     //that are zero. Need to watch that this kind of failure doesn't
3156     //cause memory access errors.
3157     assert(den->limbs[i] != 0);
3158     if (den->limbs[i] == 0)
3159     {
3160     quot->flags = GMP_INTS_EF_INTOVF_POS;
3161     rem->flags = GMP_INTS_EF_INTOVF_POS;
3162     return;
3163     }
3164    
3165     temp_limb = 0x80000000;
3166    
3167     while (((den->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
3168     {
3169     denominator_bitsize--;
3170     temp_limb >>= 1;
3171     }
3172    
3173     //The quotient starts off with the value of zero, but we consistently may
3174     //mask 1 into it and shift left. We need to be sure that we have as much
3175     //shift space there as is in the numerator. For this purpose we need to
3176     //prepare a block of clear memory as large as the numerator's.
3177     if (quot->n_allocd < abs_num_size)
3178     {
3179     GMP_INTS_mpz_realloc(quot, abs_num_size); //Make it big enough.
3180     }
3181     //Now, zero the memory.
3182     for (i=0; i<abs_num_size; i++)
3183     quot->limbs[i] = 0;
3184    
3185     //Determine the division loop count. This is the difference
3186     //in bit sizes between the numerator and denominator. It is
3187     //possible for this number to be negative, which means that the
3188     //main division loop will be executed zero times. This gives the
3189     //right results.
3190     division_loop_count = numerator_bitsize - denominator_bitsize;
3191    
3192     //We need to calculate some important numbers from the division loop
3193     //count. We need to know this number MOD 32 (which tells how far to
3194     //shift the divisor bitwise to line up with the numerator), and we
3195     //also need this number DIV 32 for the limb-wise shift.
3196     division_loop_count_mod_32 = division_loop_count % 32;
3197     division_loop_count_div_32 = division_loop_count / 32;
3198    
3199     //We now need a shift register in which we shift the denominator up
3200     //for repeated comparisons. We should dynamically allocate this to
3201     //be the same size as the numerator. Using _alloca() is OK, as one
3202     //of the unit tests is to be sure that _alloca() will handle integer
3203     //of the maximum allowed size.
3204     trial_divisor = _alloca(abs_num_size * sizeof(GMP_INTS_limb_t));
3205    
3206     //Our trial divisor needs to start off with the divisor shifted up
3207     //so that the most significant bit is aligned with the numerator.
3208     for (i = 0; i < abs_num_size; i++)
3209     {
3210     if ((division_loop_count < 0) || (i < division_loop_count_div_32))
3211     {
3212     trial_divisor[i] = 0;
3213     }
3214     else
3215     {
3216     if ((i-division_loop_count_div_32) < abs_den_size)
3217     trial_divisor[i] = den->limbs[i - division_loop_count_div_32];
3218     else
3219     trial_divisor[i] = 0;
3220     }
3221     }
3222    
3223     //The code above planted the limbs in the right place. Now need to shift bits
3224     //upward by the remaining number.
3225     if ((division_loop_count > 0) && (division_loop_count_mod_32 > 0))
3226     {
3227     //There is an existing function we can call to do the left shift.
3228     GMP_INTS_mpn_lshift(trial_divisor,
3229     trial_divisor,
3230     abs_num_size,
3231     division_loop_count_mod_32);
3232     }
3233    
3234    
3235     //Everything is ready to go. Now begin the division loop itself. It is possible
3236     //for the loop to execute zero times, which will happen if the denominator is longer
3237     //in bits than the numerator. In such cases, we can't execute this loop even once
3238     //because the math assumes that the numerator is at least as long as the denominator.
3239     for (division_counter = 0; division_counter < division_loop_count+1; division_counter++)
3240     {
3241     //Shift the quotient left one bit.
3242     GMP_INTS_mpn_lshift(quot->limbs,
3243     quot->limbs,
3244     abs_num_size,
3245     1);
3246    
3247     //If the remainder is at least as large as the trial divisor, subtract the trial
3248     //divisor from the remainder and mask in the quotient.
3249     if (GMP_INTS_mpn_cmp(rem->limbs,
3250     trial_divisor,
3251     abs_num_size) >= 0)
3252     {
3253     GMP_INTS_mpn_sub(rem->limbs,
3254     rem->limbs,
3255     abs_num_size,
3256     trial_divisor,
3257     abs_num_size);
3258     quot->limbs[0] |= 1;
3259     }
3260    
3261     //Shift the trial divisor right one bit.
3262     GMP_INTS_mpn_rshift(trial_divisor,
3263     trial_divisor,
3264     abs_num_size,
3265     1);
3266     } //End for each iteration of the division loop.
3267    
3268     //Normalize the quotient and the remainder. The normalization
3269     //process is to bring the sizes down if we have leading
3270     //zeros.
3271     quot->size = abs_num_size;
3272     GMP_INTS_mpn_normalize(quot->limbs, &(quot->size));
3273     rem->size = abs_num_size;
3274     GMP_INTS_mpn_normalize(rem->limbs, &(rem->size));
3275    
3276     //Adjust the signs as required.
3277     if (quotient_sign < 0)
3278     quot->size = -(quot->size);
3279     if (remainder_sign < 0)
3280     rem->size = -(rem->size);
3281     }
3282    
3283    
3284     void GMP_INTS_mpz_fac_ui(GMP_INTS_mpz_struct *result,
3285     unsigned long int n)
3286     {
3287     //Just multiply the numbers in ascending order. The original
3288     //GNU library contained a much more elegant algorithm, but
3289     //this is more direct.
3290    
3291     unsigned long int k;
3292    
3293     GMP_INTS_mpz_set_ui (result, 1L);
3294    
3295     for (k = 2; (k <= n) && !(result->flags); k++)
3296     GMP_INTS_mpz_mul_ui (result, result, k);
3297     }
3298    
3299    
3300     /******************************************************************/
3301     /*** PUBLIC CONVERSION AND OUTPUT FUNCTIONS ********************/
3302     /******************************************************************/
3303     //07/18/01: Visual inspection OK. Function returns
3304     //reasonable values even out to 100,000 digits--seems OK.
3305     int GMP_INTS_mpz_size_in_base_10(const GMP_INTS_mpz_struct *arg)
3306     {
3307     _int64 n;
3308    
3309     //Eyeball the input parameter.
3310     assert(arg != NULL);
3311     assert(arg->n_allocd > 0);
3312     assert(arg->limbs != NULL);
3313    
3314     //Get the number of limbs occupied by the integer.
3315     //Because even the digit zero takes some space,
3316     //don't accept zero for an answer.
3317     n = GMP_INTS_abs_of_size_t(arg->size);
3318     if (n==0)
3319     n = 1;
3320    
3321     //Convert this to the number of bits. Generously
3322     //ignore any unused leading bits.
3323     n *= 32;
3324    
3325     //Used a slightly high best rational approximation in F_{65535}
3326     //to go from the number of bits to the number of
3327     //digits. The division discards, so bump the result
3328     //up by 1 to compensate for possible truncation. The number
3329     //we are aproximating is ln(2)/ln(10).
3330     n *= 12655;
3331     n /= 42039;
3332     n++;
3333    
3334     //Compensate for possible commas in the result. Again,
3335     //consider truncation.
3336     n *= 4;
3337     n /= 3;
3338     n++;
3339    
3340     //Compensate for the minus sign, the trailing zero,
3341     //cosmic rays striking the computer from the martian
3342     //listening post camoflaged on the moon, and the
3343     //possibility that we might need to put text in the
3344     //string if any flag is set.
3345     n += 100;
3346    
3347     //And that should be a good return value.
3348     return((int) n);
3349     }
3350    
3351    
3352     //07/19/01: Visual inspection and unit test is OK.
3353     void GMP_INTS_mpz_to_string(char *out,
3354     const GMP_INTS_mpz_struct *in)
3355     {
3356     //Eyeball the input parameters.
3357     assert(out != NULL);
3358     assert(in != NULL);
3359     assert(in->n_allocd > 0);
3360     assert(in->limbs != NULL);
3361    
3362     //If any of the flags are set, stuff in the text.
3363     if (in->flags)
3364     {
3365     if (in->flags & GMP_INTS_EF_INTOVF_POS)
3366     {
3367     strcpy(out, GMP_INTS_EF_INTOVF_POS_STRING);
3368     }
3369     else if (in->flags & GMP_INTS_EF_INTOVF_NEG)
3370     {
3371     strcpy(out, GMP_INTS_EF_INTOVF_NEG_STRING);
3372     }
3373     else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_POS)
3374     {
3375     strcpy(out, GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
3376     }
3377     else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_NEG)
3378     {
3379     strcpy(out, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
3380     }
3381     else
3382     {
3383     strcpy(out, "INTERNAL_ERROR");
3384     }
3385     }
3386     else
3387     {
3388     //Ordinary integer conversion.
3389     GMP_INTS_mpz_struct num, den, quot, rem, k10;
3390    
3391     //Allocate space for the temporary integers.
3392     GMP_INTS_mpz_init(&num);
3393     GMP_INTS_mpz_init(&den);
3394     GMP_INTS_mpz_init(&quot);
3395     GMP_INTS_mpz_init(&rem);
3396     GMP_INTS_mpz_init(&k10);
3397    
3398     //Assign the constant 10.
3399     GMP_INTS_mpz_set_ui(&k10, 10);
3400    
3401     //If the integer is zero, assign that.
3402     if (in->size == 0)
3403     {
3404     strcpy(out, "0");
3405     }
3406     else
3407     {
3408     //We have to do a full conversion. The algorithm
3409     //is division by 10, each time obtaining the least
3410     //significant digit, until finally the quotient is
3411     //zero.
3412     char *ptr;
3413    
3414     ptr = out;
3415    
3416     GMP_INTS_mpz_copy(&num, in);
3417     GMP_INTS_mpz_copy(&den, &k10);
3418     do
3419     {
3420     #if 0
3421     printf("Values before division:\n");
3422     FCMIOF_hline();
3423     GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
3424     FCMIOF_hline();
3425     GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
3426     FCMIOF_hline();
3427     GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
3428     FCMIOF_hline();
3429     GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
3430     FCMIOF_hline();
3431    
3432     if (num.size > 1)
3433     FCMIOF_hline();
3434     #endif
3435    
3436     GMP_INTS_mpz_tdiv_qr(&quot, &rem, &num, &den);
3437     #if 0
3438     printf("Values after division:\n");
3439     FCMIOF_hline();
3440     GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
3441     FCMIOF_hline();
3442     GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
3443     FCMIOF_hline();
3444     GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
3445     FCMIOF_hline();
3446     GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
3447     FCMIOF_hline();
3448     #endif
3449    
3450     if (rem.size != 0)
3451     {
3452     *ptr = '0' + (char)(rem.limbs[0]);
3453     }
3454     else
3455     {
3456     *ptr = '0';
3457     }
3458     ptr++;
3459     GMP_INTS_mpz_copy(&num, &quot);
3460     //printf("digit\n");
3461     }
3462     while (!GMP_INTS_mpz_is_zero(&quot));
3463    
3464     //Finally, if the input was negative, tack on the
3465     //minus sign.
3466     if (GMP_INTS_mpz_is_neg(in))
3467     {
3468     *ptr = '-';
3469     ptr++;
3470     }
3471    
3472     //Finally, tack on the trailing zero terminator.
3473     *ptr = 0;
3474     ptr++;
3475    
3476     //Reverse the string.
3477     BSTRFUNC_str_reverse(out);
3478     }
3479    
3480     //Deallocate the integers.
3481     GMP_INTS_mpz_clear(&num);
3482     GMP_INTS_mpz_clear(&den);
3483     GMP_INTS_mpz_clear(&quot);
3484     GMP_INTS_mpz_clear(&rem);
3485     GMP_INTS_mpz_clear(&k10);
3486     }
3487     }
3488    
3489    
3490     void GMP_INTS_mpz_long_int_format_to_stream(FILE *s,
3491     const GMP_INTS_mpz_struct *i,
3492     const char *desc)
3493     {
3494     int line_len;
3495     int digits_per_line;
3496     char *digits;
3497     int num_digits;
3498     int nlines;
3499     int cur_line;
3500     int number_desc_width;
3501    
3502     //Eyeball the inputs, make sure the caller isn't doing
3503     //something stupid.
3504     assert(s != NULL);
3505     assert(i != NULL);
3506     assert(i->n_allocd > 0);
3507     assert(i->limbs != NULL);
3508     assert(desc != NULL);
3509    
3510     //Obtain the line length assumed for formatted output.
3511     line_len = FCMIOF_get_line_len();
3512    
3513     //The description width allowed is 20.
3514     number_desc_width = 20;
3515    
3516     /* The number of digits per line that we assume must be a multiple of
3517     ** three. The formula below was not examined very carefully, but it
3518     ** works fine for a line length of 78. If line length is changed,
3519     ** this formula may need to be examined very carefully and rewritten.
3520     */
3521     digits_per_line = INTFUNC_max(3, ((((line_len-42)*3)/4)/3)*3);
3522     assert(digits_per_line >= 3);
3523    
3524     /* We now need to get a digit string corresponding to this
3525     ** number. First, need to figure out how much and
3526     ** allocate the space.
3527     */
3528     digits = GMP_INTS_malloc(GMP_INTS_mpz_size_in_base_10(i) * sizeof(char));
3529     GMP_INTS_mpz_to_string(digits, i);
3530    
3531     //If the number is negative, delete the leading minus sign.
3532     //The rest of the display algorithm needs an unsigned
3533     //series of digits.
3534     if (*digits == '-')
3535     {
3536     int i = 0;
3537    
3538     do
3539     {
3540     digits[i] = digits[i+1];
3541     i++;
3542     }
3543     while(digits[i-1]);
3544     }
3545    
3546     //Figure out how many digits in the string representation.
3547     num_digits = strlen(digits);
3548    
3549     /* As the first order of business, figure out how many lines the beast
3550     ** will require.
3551     */
3552     if (i->flags)
3553     {
3554     nlines = 1; /* Only one line required for NAN verbeage. */
3555     }
3556     else if (GMP_INTS_mpz_is_zero(i))
3557     {
3558     nlines = 1; /* The zero value requires one line. */
3559     }
3560     else
3561     {
3562     /* In any other case, have a formula.
3563     */
3564     nlines = 1 + (num_digits - 1) / digits_per_line;
3565     }
3566    
3567     /* Iterate through each line, spitting out whatever is appropriate. */
3568     for (cur_line = 0; cur_line < nlines; cur_line++)
3569     {
3570     int cur_digit_on_line;
3571    
3572     /* If this is the first line, spit out the description, right-aligned.
3573     ** Otherwise, spit spaces.
3574     */
3575     if (!cur_line)
3576     {
3577     /* First line. */
3578     int len;
3579    
3580     len = strlen(desc);
3581    
3582     if (len <= number_desc_width)
3583     {
3584     /* Description is shorter or equal, pad on left. */
3585     FCMIOF_stream_repchar(s, ' ', number_desc_width - len);
3586     fprintf(s, "%s", desc);
3587     }
3588     else
3589     {
3590     /* Description is too long, truncate. */
3591     int i;
3592    
3593     for (i=0; i<number_desc_width; i++)
3594     fprintf(s, "%c", desc[i]);
3595     }
3596    
3597     fprintf(s, ": ");
3598    
3599     /* If the number is negative, throw in a minus sign. */
3600     if (GMP_INTS_mpz_is_neg(i) && !(i->flags))
3601     {
3602     fprintf(s, "- ");
3603     }
3604     else
3605     {
3606     fprintf(s, " ");
3607     }
3608     }
3609     else
3610     {
3611     /* Every line but first line. */
3612     FCMIOF_stream_repchar(s, ' ', number_desc_width+4);
3613     }
3614    
3615     for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++)
3616     {
3617     int idx_into_string;
3618     /* Index into the string which is our digit of interest.
3619     */
3620    
3621     /* Compute the index. The equation is based on the ordering
3622     ** of presentation, for example,
3623     **
3624     ** 7 6
3625     ** 5 4 3
3626     ** 2 1 0.
3627     **
3628     ** With a little thought, the equation should make sense.
3629     ** The index won't always be used to index into the string.
3630     */
3631     idx_into_string =
3632     ((((nlines-1) - cur_line) * digits_per_line)
3633     +
3634     (digits_per_line - 1 - cur_digit_on_line));
3635    
3636     /* Print the appropriate digit or a space. The NAN case and the
3637     ** zero case need to be treated specially.
3638     */
3639     if (i->flags)
3640     {
3641     /* Not a number. Everything is blank, except spell out
3642     ** description of condition at the end of the string of
3643     ** digits.
3644     */
3645     int index_from_right;
3646     int virtual_index;
3647    
3648     index_from_right = digits_per_line - 1 - cur_digit_on_line;
3649     //The index calculated above is calculated so that the
3650     //final position on the line has index [0].
3651     assert(index_from_right >= 0 && index_from_right < digits_per_line);
3652    
3653     //Now, calculate the "virtual index". The virtual index
3654     //is the actual number of characters from the right, taking
3655     //into account commas.
3656     virtual_index = index_from_right + index_from_right/3;
3657    
3658     if (((index_from_right % 3) == 2) && cur_digit_on_line)
3659     {
3660     //We are one position past a comma. This means
3661     //that we might need a "fill" character to go
3662     //where the comma should have gone.
3663    
3664     if (virtual_index + 1 < num_digits)
3665     {
3666     //The character we should print exists.
3667     fprintf(s, "%c", digits[num_digits - 2 - virtual_index]);
3668     }
3669     else
3670     {
3671     //The character doesn't exist, because the error
3672     //string is apparently too short. Must print a
3673     //space, instead.
3674     fprintf(s, " ");
3675     }
3676     }
3677    
3678     //We've done the fill character, if the position we're in
3679     //is one past a comma. Now, do the ordinary character
3680     //corresponding to a digit position.
3681     if (virtual_index < num_digits)
3682     {
3683     //The character we should print exists.
3684     fprintf(s, "%c", digits[num_digits - 1 - virtual_index]);
3685     }
3686     else
3687     {
3688     //The character doesn't exist, because the error
3689     //string is apparently too short. Must print a
3690     //space, instead.
3691     fprintf(s, " ");
3692     }
3693     }
3694     else if (GMP_INTS_mpz_is_zero(i))
3695     {
3696     /* This is the zero case. For zero, there is only one line,
3697     ** and every character except the last one is a blank.
3698     */
3699     if (cur_digit_on_line == (digits_per_line - 1))
3700     {
3701     fprintf(s, "0");
3702     }
3703     else
3704     {
3705     fprintf(s, " ");
3706     }
3707     }
3708     else
3709     {
3710     /* This is a valid number which is not zero. Need to print
3711     ** the digits.
3712     */
3713    
3714     if (idx_into_string < num_digits)
3715     {
3716     int actual_index;
3717    
3718     actual_index = num_digits - 1 - idx_into_string;
3719     //This is a string reversal mapping. The original
3720     //code stored strings least significant digit first,
3721     //but this code uses most significant digit first.
3722     assert((actual_index >= 0) && (actual_index < num_digits));
3723     fprintf(s, "%c", digits[actual_index]);
3724     }
3725     else
3726     {
3727     fprintf(s, " ");
3728     }
3729     } /* End of digit case.
3730    
3731     /* Now handle the commas. The rules for commas are straightforward.
3732     ** a)NAN never has a comma.
3733     ** b)Zeros never have a comma.
3734     ** c)The final line, last digit never has a comma.
3735     ** d)Everything else in multiples of three ...
3736     */
3737     if (!(idx_into_string % 3) && (idx_into_string))
3738     {
3739     if (i->flags)
3740     {
3741     //fprintf(s, " ");
3742     }
3743     else if (!num_digits)
3744     {
3745     fprintf(s, " ");
3746     }
3747     else
3748     {
3749     if (idx_into_string < num_digits)
3750     {
3751     fprintf(s, ",");
3752     }
3753     else
3754     {
3755     fprintf(s, " ");
3756     }
3757     }
3758     }
3759     } /* End for each digit on the current line. */
3760    
3761     /* For the first line, print out an informative message
3762     ** advising of the number of digits. For all other lines
3763     ** print nothing.
3764     */
3765     if (!cur_line && !(i->flags))
3766     {
3767     if (nlines == 1)
3768     fprintf(s, " ");
3769    
3770     if (num_digits <= 1)
3771     {
3772     fprintf(s, " ( 1 digit )\n");
3773     }
3774     else if (num_digits < 1000)
3775     {
3776     fprintf(s, " (%7d digits)\n", num_digits);
3777     }
3778     else
3779     {
3780     fprintf(s, " (%3d,%03d digits)\n", num_digits / 1000, num_digits % 1000);
3781     }
3782     }
3783     else
3784     {
3785     fprintf(s, "\n");
3786     }
3787     } /* End for each line. */
3788    
3789     //Deallocate the string space.
3790     GMP_INTS_free(digits);
3791     }
3792    
3793    
3794     void GMP_INTS_mpz_arb_int_raw_to_stream(FILE *s,
3795     const GMP_INTS_mpz_struct *i)
3796     {
3797     int size_reqd;
3798     char *digits;
3799    
3800     //Eyeball the input parameters.
3801     assert(s != NULL);
3802     assert(i != NULL);
3803     assert(i->n_allocd > 0);
3804     assert(i->limbs != NULL);
3805    
3806     size_reqd = GMP_INTS_mpz_size_in_base_10(i);
3807     digits = GMP_INTS_malloc(size_reqd * sizeof(char));
3808     GMP_INTS_mpz_to_string(digits, i);
3809     fprintf(s, "%s", digits);
3810     GMP_INTS_free(digits);
3811     }
3812    
3813    
3814     //07/24/01: Passed visual inspection and unit tests.
3815     void GMP_INTS_mpz_pow_ui( GMP_INTS_mpz_struct *result,
3816     const GMP_INTS_mpz_struct *base,
3817     unsigned exponent)
3818     {
3819     GMP_INTS_mpz_struct temp;
3820     //Temporary location to hold the base raised to
3821     //a binary power (repeated squaring).
3822    
3823     //Eyeball the input parameters.
3824     assert(result != NULL);
3825     assert(result->n_allocd > 0);
3826     assert(result->limbs != NULL);
3827     assert(base != NULL);
3828     assert(base->n_allocd > 0);
3829     assert(base->limbs != NULL);
3830    
3831     //For this function, the base and the result may not
3832     //be the same object.
3833     assert(result != base);
3834    
3835     //If the base is tained, the output is tainted by association.
3836     {
3837     int taint;
3838    
3839     taint = GMP_INTS_two_op_flags_map(base->flags, 0);
3840    
3841     if (taint)
3842     {
3843     result->flags = taint;
3844     return;
3845     }
3846     }
3847    
3848     //Allocate our temporary variable and set it to the base.
3849     GMP_INTS_mpz_init(&temp);
3850     GMP_INTS_mpz_copy(&temp, base);
3851    
3852     //The result begins with the value of 1.
3853     GMP_INTS_mpz_set_ui(result, 1);
3854    
3855     //Loop through, processing each bit of the exponent. This is a fairly effective
3856     //algorithm, but not the optimal one (Knuth points this out).
3857     while (exponent && !result->flags)
3858     {
3859     if (exponent & 0x1)
3860     {
3861     GMP_INTS_mpz_mul(result, result, &temp);
3862     }
3863    
3864     //Square the temporary variable. Because squaring of arb integer
3865     //may be very expensive, the test against 1 (i.e. last iteration)
3866     //certainly pays for itself.
3867     if (exponent != 1)
3868     GMP_INTS_mpz_mul(&temp, &temp, &temp);
3869    
3870     exponent >>= 1;
3871     }
3872    
3873     //Deallocate our temporary variable.
3874     GMP_INTS_mpz_clear(&temp);
3875     }
3876    
3877    
3878     void GMP_INTS_mpz_abs(GMP_INTS_mpz_struct *arg)
3879     {
3880     //Eyeball the input parameter.
3881     assert(arg != NULL);
3882     assert(arg->n_allocd > 0);
3883     assert(arg->limbs != NULL);
3884    
3885     //Take the absolute value.
3886     if (arg->size < 0)
3887     arg->size = -arg->size;
3888