/[dtapublic]/sf_code/esrgpcpj/shared/c_datd/gmp_ints.c
ViewVC logotype

Contents of /sf_code/esrgpcpj/shared/c_datd/gmp_ints.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (show annotations) (download)
Sat Oct 8 06:43:03 2016 UTC (7 years, 5 months ago) by dashley
File MIME type: text/plain
File size: 147987 byte(s)
Initial commit.
1 // $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 }
3889
3890
3891 //07/29/01: Visual inspection passed. Seems to work fine--not explicitly unit-tested
3892 //directly, but was tested from Tcl.
3893 void GMP_INTS_mpz_gcd(GMP_INTS_mpz_struct *result,
3894 const GMP_INTS_mpz_struct *arg1,
3895 const GMP_INTS_mpz_struct *arg2)
3896 {
3897 GMP_INTS_mpz_struct u, v, q, r;
3898 int loop_count;
3899
3900 //Eyeball the inputs carefully.
3901 assert(result != NULL);
3902 assert(result->n_allocd > 0);
3903 assert(result->limbs != NULL);
3904 assert(arg1 != NULL);
3905 assert(arg1->n_allocd > 0);
3906 assert(arg1->limbs != NULL);
3907 assert(arg2 != NULL);
3908 assert(arg2->n_allocd > 0);
3909 assert(arg2->limbs != NULL);
3910
3911 //Args are not allowed to be same object.
3912 assert(arg1 != arg2);
3913
3914 //If either input is error or taint, taint the output.
3915 {
3916 int taint;
3917
3918 taint = GMP_INTS_two_op_flags_map(arg1->flags, arg2->flags);
3919
3920 result->flags = 0;
3921 //"result" starts off with a clean slate. Must do this
3922 //after taint calculation in case locations of arg1 or arg2
3923 //are the same as result.
3924 if (taint)
3925 {
3926 result->flags = taint;
3927 return;
3928 }
3929 }
3930
3931 //If either input is zero, the result is 1.
3932 if (GMP_INTS_mpz_is_zero(arg1) || GMP_INTS_mpz_is_zero(arg2))
3933 {
3934 GMP_INTS_mpz_set_ui(result, 1);
3935 return;
3936 }
3937
3938 //Allocate space for locals.
3939 GMP_INTS_mpz_init(&u);
3940 GMP_INTS_mpz_init(&v);
3941 GMP_INTS_mpz_init(&q);
3942 GMP_INTS_mpz_init(&r);
3943
3944 //We are following Knuth Vol 2, p. 337, the modern Euclidian algorithm.
3945 //Note: There are faster algorithms for GCD, but because I hacked up the
3946 //GMP multiple-precision library so badly, those aren't included. This one
3947 //is logically correct but sub-optimal. Perhaps at a later time faster
3948 //algorithms will be re-included.
3949 //Copy inputs to u and v.
3950 GMP_INTS_mpz_copy(&u, arg1);
3951 GMP_INTS_mpz_copy(&v, arg2);
3952
3953 //Take the absolute value of each argument. We know that neither is zero,
3954 //but one or both might be negative.
3955 GMP_INTS_mpz_abs(&u);
3956 GMP_INTS_mpz_abs(&v);
3957
3958 //Begin Euclid's algorithm. There are really three possibilities:
3959 // a)We terminate normally.
3960 // b)Somehow we generate a math error and terminate based on flags.
3961 // c)Due to some unknown error in the math functions, we go on forever,
3962 // and the program locks up.
3963 GMP_INTS_mpz_set_ui(&r, 1);
3964 loop_count = 0;
3965 while (!GMP_INTS_mpz_is_zero(&r) && !q.flags && !r.flags && (loop_count < 100000))
3966 {
3967 loop_count++;
3968
3969 GMP_INTS_mpz_tdiv_qr(&q, &r, &u, &v);
3970 GMP_INTS_mpz_copy(&u, &v);
3971 GMP_INTS_mpz_copy(&v, &r);
3972 }
3973
3974 //Let's hope we didn't get out of the loop based on loop count.
3975 assert(loop_count != 100000);
3976
3977 //u now contains the answer.
3978 GMP_INTS_mpz_copy(result, &u);
3979
3980 //Deallocate space for locals.
3981 GMP_INTS_mpz_clear(&u);
3982 GMP_INTS_mpz_clear(&v);
3983 GMP_INTS_mpz_clear(&q);
3984 GMP_INTS_mpz_clear(&r);
3985 }
3986
3987
3988 /******************************************************************/
3989 /*** COMPARISON AND SIZING FUNCTIONS ***************************/
3990 /******************************************************************/
3991 //07/24/01: Visual inspection only, due to simplicity.
3992 int GMP_INTS_mpz_fits_uint_p (const GMP_INTS_mpz_struct *src)
3993 {
3994 GMP_INTS_size_t size;
3995 GMP_INTS_limb_t mpl;
3996
3997 //Eyeball the input parameter.
3998 assert(src != NULL);
3999 assert(src->n_allocd > 0);
4000 assert(src->limbs != NULL);
4001
4002 mpl = src->limbs[0];
4003 size = src->size;
4004 if (size < 0 || size > 1)
4005 return(0);
4006
4007 //The following line came from the original GNU code.
4008 //It isn't necessary in our case since limbs and ints are
4009 //both 32 bits, but it will do no harm.
4010 return (mpl <= (~(unsigned int) 0));
4011 }
4012
4013
4014 unsigned GMP_INTS_mpz_get_limb_zero(const GMP_INTS_mpz_struct *src)
4015 {
4016 //Eyeball the inputs.
4017 assert(src != NULL);
4018 assert(src->n_allocd > 0);
4019 assert(src->limbs != NULL);
4020
4021 if (!src->size)
4022 return(0);
4023 else
4024 return(src->limbs[0]);
4025 }
4026
4027
4028 //07/24/01: Visual inspection only. Understood the comparisons
4029 //and seems like they should work, but ... a little beyond my
4030 //comfort zone without testing. Trusting GNU on this one ...
4031 int GMP_INTS_mpz_fits_sint_p (const GMP_INTS_mpz_struct *src)
4032 {
4033 GMP_INTS_size_t size;
4034 GMP_INTS_limb_t mpl;
4035
4036 //Eyeball the input parameter.
4037 assert(src != NULL);
4038 assert(src->n_allocd > 0);
4039 assert(src->limbs != NULL);
4040
4041 mpl = src->limbs[0];
4042 size = src->size;
4043 if (size > 0)
4044 {
4045 if (size > 1)
4046 return 0;
4047 return mpl < ~((~(unsigned int) 0) >> 1);
4048 }
4049 else
4050 {
4051 if (size < -1)
4052 return 0;
4053 return mpl <= ~((~(unsigned int) 0) >> 1);
4054 }
4055 }
4056
4057
4058 //07/24/01: Visual inspection only. One issue that caught
4059 //my eye is that in one place function returned neg value if
4060 //< and pos value if >. Was within spec, but corrected because
4061 //it concerned me as I often test against -1 and 1. Seems
4062 //to invite accidents.
4063 int GMP_INTS_mpz_cmp (const GMP_INTS_mpz_struct *u,
4064 const GMP_INTS_mpz_struct *v)
4065 {
4066 GMP_INTS_size_t usize = u->size;
4067 GMP_INTS_size_t vsize = v->size;
4068 GMP_INTS_size_t size;
4069 GMP_INTS_limb_srcptr up, vp;
4070 int cmp;
4071
4072 //Eyeball the input parameters.
4073 assert(u != NULL);
4074 assert(u->n_allocd > 0);
4075 assert(u->limbs != NULL);
4076 assert(v != NULL);
4077 assert(v->n_allocd > 0);
4078 assert(v->limbs != NULL);
4079
4080 if (usize < vsize)
4081 return(-1);
4082 else if (usize > vsize)
4083 return(1);
4084
4085 if (usize == 0)
4086 return(0);
4087
4088 size = GMP_INTS_abs_of_size_t(usize);
4089
4090 up = u->limbs;
4091 vp = v->limbs;
4092
4093 cmp = GMP_INTS_mpn_cmp (up, vp, size);
4094
4095 if (cmp == 0)
4096 return(0);
4097
4098 if ((cmp < 0) == (usize < 0))
4099 return(1);
4100 else
4101 return(-1);
4102 }
4103
4104
4105 //07/24/01: Not visually inspected. Relying on
4106 //GNU ...
4107 int GMP_INTS_mpz_cmp_ui (const GMP_INTS_mpz_struct *u,
4108 unsigned long int v_digit)
4109 {
4110 GMP_INTS_size_t usize = u->size;
4111
4112 //Eyeball the input parameter.
4113 assert(u != NULL);
4114 assert(u->n_allocd > 0);
4115 assert(u->limbs != NULL);
4116
4117 if (usize == 0)
4118 return -(v_digit != 0);
4119
4120 if (usize == 1)
4121 {
4122 GMP_INTS_limb_t u_digit;
4123
4124 u_digit = u->limbs[0];
4125 if (u_digit > v_digit)
4126 return 1;
4127 if (u_digit < v_digit)
4128 return -1;
4129 return 0;
4130 }
4131
4132 return (usize > 0) ? 1 : -1;
4133 }
4134
4135
4136 //07/24/01: Not visually inspected. Relying on GNU.
4137 int GMP_INTS_mpz_cmp_si (const GMP_INTS_mpz_struct *u,
4138 signed long int v_digit)
4139 {
4140 GMP_INTS_size_t usize = u->size;
4141 GMP_INTS_size_t vsize;
4142 GMP_INTS_limb_t u_digit;
4143
4144 //Eyeball the input parameter.
4145 assert(u != NULL);
4146 assert(u->n_allocd > 0);
4147 assert(u->limbs != NULL);
4148
4149 vsize = 0;
4150 if (v_digit > 0)
4151 vsize = 1;
4152 else if (v_digit < 0)
4153 {
4154 vsize = -1;
4155 v_digit = -v_digit;
4156 }
4157
4158 if (usize != vsize)
4159 return usize - vsize;
4160
4161 if (usize == 0)
4162 return 0;
4163
4164 u_digit = u->limbs[0];
4165
4166 if (u_digit == (GMP_INTS_limb_t) (unsigned long) v_digit)
4167 return 0;
4168
4169 if (u_digit > (GMP_INTS_limb_t) (unsigned long) v_digit)
4170 return usize;
4171 else
4172 return -usize;
4173 }
4174
4175
4176 //07/24/01: Not visually inspected. Counting on GNU.
4177 int GMP_INTS_mpz_cmpabs (const GMP_INTS_mpz_struct *u,
4178 const GMP_INTS_mpz_struct *v)
4179 {
4180 GMP_INTS_size_t usize = u->size;
4181 GMP_INTS_size_t vsize = v->size;
4182 GMP_INTS_limb_srcptr up, vp;
4183 int cmp;
4184
4185 //Eyeball the input parameters.
4186 assert(u != NULL);
4187 assert(u->n_allocd > 0);
4188 assert(u->limbs != NULL);
4189 assert(v != NULL);
4190 assert(v->n_allocd > 0);
4191 assert(v->limbs != NULL);
4192
4193 usize = GMP_INTS_abs_of_size_t(usize);
4194 vsize = GMP_INTS_abs_of_size_t(vsize);
4195
4196 if (usize != vsize)
4197 return usize - vsize;
4198
4199 if (usize == 0)
4200 return 0;
4201
4202 up = u->limbs;
4203 vp = v->limbs;
4204
4205 cmp = GMP_INTS_mpn_cmp (up, vp, usize);
4206
4207 return cmp;
4208 }
4209
4210 //07/24/01: Not visually inspected. Counting on GNU.
4211 int GMP_INTS_mpz_cmpabs_ui(const GMP_INTS_mpz_struct *u,
4212 unsigned long int v_digit)
4213 {
4214 GMP_INTS_size_t usize = u->size;
4215
4216 //Eyeball the input parameter.
4217 assert(u != NULL);
4218 assert(u->n_allocd > 0);
4219 assert(u->limbs != NULL);
4220
4221 if (usize == 0)
4222 return -(v_digit != 0);
4223
4224 usize = GMP_INTS_abs_of_size_t(usize);
4225
4226 if (usize == 1)
4227 {
4228 GMP_INTS_limb_t u_digit;
4229
4230 u_digit = u->limbs[0];
4231 if (u_digit > v_digit)
4232 return 1;
4233 if (u_digit < v_digit)
4234 return -1;
4235 return 0;
4236 }
4237
4238 return 1;
4239 }
4240
4241
4242 /******************************************************************/
4243 /*** VERSION CONTROL IDENTITY FUNCTIONS ************************/
4244 /******************************************************************/
4245
4246 //07/18/01: Visual inspection only. Function deemed too
4247 //simple for unit testing.
4248 const char *GMP_INTS_cvcinfo(void)
4249 {
4250 return("$Header: /cvsroot/esrg/sfesrg/esrgpcpj/shared/c_datd/gmp_ints.c,v 1.22 2002/01/27 15:18:44 dtashley Exp $");
4251 }
4252
4253
4254 //07/18/01: Visual inspection only. Function deemed too
4255 //simple for unit testing.
4256 const char *GMP_INTS_hvcinfo(void)
4257 {
4258 return(GMP_INTS_H_VERSION);
4259 }
4260
4261
4262 //**************************************************************************
4263 // $Log: gmp_ints.c,v $
4264 // Revision 1.22 2002/01/27 15:18:44 dtashley
4265 // Minor typo corrected.
4266 //
4267 // Revision 1.21 2002/01/27 15:04:14 dtashley
4268 // HYREACH added to batch build (which required some changes in other programs).
4269 //
4270 // Revision 1.20 2001/08/16 19:49:40 dtashley
4271 // Beginning to prepare for v1.05 release.
4272 //
4273 // Revision 1.19 2001/08/15 06:56:05 dtashley
4274 // Substantial progress. Safety check-in.
4275 //
4276 // Revision 1.18 2001/08/07 10:42:48 dtashley
4277 // Completion of CFRATNUM extensions and DOS command-line utility.
4278 //
4279 // Revision 1.17 2001/07/30 02:51:18 dtashley
4280 // INTGCD extension and command-line utility finished up.
4281 //
4282 // Revision 1.16 2001/07/29 07:18:22 dtashley
4283 // Completion of ARBINT INTFAC extension.
4284 //
4285 // Revision 1.15 2001/07/25 23:40:02 dtashley
4286 // Completion of INTFAC program, many changes to handling of large
4287 // integers.
4288 //
4289 // Revision 1.14 2001/07/21 01:39:01 dtashley
4290 // Safety check-in. Major function to output an integer as rows of digits
4291 // has been completed. This was the last major function that needed to be
4292 // completed before useful command-line utilities can be constructed.
4293 //
4294 // Revision 1.13 2001/07/19 20:06:03 dtashley
4295 // Division finished. String formatting functions underway. Safety check-in.
4296 //
4297 // Revision 1.12 2001/07/18 21:53:09 dtashley
4298 // Division function finished and passes preliminary tests. Safety check-in.
4299 //
4300 // Revision 1.11 2001/07/17 22:30:14 dtashley
4301 // Safety check-in. Division function under construction.
4302 //
4303 // Revision 1.10 2001/07/16 17:46:46 dtashley
4304 // Multiplication finished, and only indirectly unit-tested. More detailed unit
4305 // test must follow, but expect no problems.
4306 //
4307 // Revision 1.9 2001/07/16 00:28:22 dtashley
4308 // Safety check-in. Addition and subtraction functions finished.
4309 //
4310 // Revision 1.8 2001/07/15 06:40:10 dtashley
4311 // Adaptation of GNU arbitrary-size integer package integrated into IjuScripter
4312 // and IjuConsole.
4313 //
4314 // Revision 1.7 2001/07/15 00:59:52 dtashley
4315 // Safety check-in. Commit before working on laptop.
4316 //
4317 // Revision 1.6 2001/07/14 07:03:37 dtashley
4318 // Safety check-in. Modifications and progress.
4319 //
4320 // Revision 1.5 2001/07/14 02:05:02 dtashley
4321 // Safety check-in. Almost ready for first unit-testing.
4322 //
4323 // Revision 1.4 2001/07/13 21:02:20 dtashley
4324 // Version control reporting changes.
4325 //
4326 // Revision 1.3 2001/07/13 06:54:57 dtashley
4327 // Safety check-in. Substantial progress and modifications.
4328 //
4329 // Revision 1.2 2001/07/13 00:57:08 dtashley
4330 // Safety check-in. Substantial progress on port.
4331 //
4332 // Revision 1.1 2001/07/12 05:07:02 dtashley
4333 // Initial checkin.
4334 //
4335 //**************************************************************************
4336 // End of CCMALLOC.C.

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25