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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations) (download)
Sat Oct 8 06:43:03 2016 UTC (6 years, 5 months ago) by dashley
Original Path: sf_code/esrgpcpj/shared/c_datd/gmp_ralg.c
File MIME type: text/plain
File size: 132737 byte(s)
Initial commit.
1 dashley 25 // $Header: /cvsroot/esrg/sfesrg/esrgpcpj/shared/c_datd/gmp_ralg.c,v 1.10 2002/01/27 17:58:15 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_RALG
686    
687     #include <assert.h>
688     #include <stdio.h>
689     #include <string.h>
690     #include <time.h>
691    
692    
693     #include "fcmiof.h"
694     #include "gmp_ints.h"
695     #include "gmp_rats.h"
696     #include "gmp_ralg.h"
697     #include "intfunc.h"
698    
699     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
700     #include "ccmalloc.h"
701     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
702     #include "tclalloc.h"
703     #else
704     /* Do nothing. */
705     #endif
706    
707    
708     /******************************************************************/
709     /*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/
710     /******************************************************************/
711     //08/16/01: Visual inspection OK.
712     void GMP_RALG_cfdecomp_init(
713     GMP_RALG_cf_app_struct *decomp,
714     int *failure,
715     GMP_INTS_mpz_struct *num,
716     GMP_INTS_mpz_struct *den)
717     {
718     int loop_counter, i;
719     GMP_INTS_mpz_struct arb_temp1, arb_temp2;
720    
721     //Eyeball the input parameters. The rest of the eyeballing
722     //will occur as functions are called to manipulate the
723     //numerator and denominator passed in.
724     assert(decomp != NULL);
725     assert(failure != NULL);
726     assert(num != NULL);
727     assert(den != NULL);
728    
729     //Allocate space for temporary integers.
730     GMP_INTS_mpz_init(&arb_temp1);
731     GMP_INTS_mpz_init(&arb_temp2);
732    
733     //Begin believing no failure.
734     *failure = 0;
735    
736     //Initialize the copy of the numerator and denominator and
737     //copy these into the structure.
738     GMP_INTS_mpz_init(&(decomp->num));
739     GMP_INTS_mpz_copy(&(decomp->num), num);
740     GMP_INTS_mpz_init(&(decomp->den));
741     GMP_INTS_mpz_copy(&(decomp->den), den);
742    
743     //Allocate the space for the first increment of the
744     //growable areas. We need to use different memory
745     //allocation functions depending on whether we're
746     //in a Tcl build or a DOS command-line utility
747     //build.
748     //Space for partial quotients.
749     decomp->a =
750     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
751     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
752     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
753     (GMP_INTS_mpz_struct *)
754     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
755     #else
756     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
757     #endif
758    
759     //Dividends.
760     decomp->dividend =
761     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
762     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
763     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
764     (GMP_INTS_mpz_struct *)
765     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
766     #else
767     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
768     #endif
769    
770     //Divisors.
771     decomp->divisor =
772     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
773     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
774     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
775     (GMP_INTS_mpz_struct *)
776     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
777     #else
778     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
779     #endif
780    
781     //Remainders.
782     decomp->remainder =
783     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
784     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
785     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
786     (GMP_INTS_mpz_struct *)
787     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
788     #else
789     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
790     #endif
791    
792     //Convergent numerators.
793     decomp->p =
794     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
795     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
796     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
797     (GMP_INTS_mpz_struct *)
798     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
799     #else
800     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
801     #endif
802    
803     //Convergent denominators.
804     decomp->q =
805     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
806     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
807     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
808     (GMP_INTS_mpz_struct *)
809     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
810     #else
811     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
812     #endif
813    
814     //Now the number of allocated slots is what we just allocated.
815     decomp->nallocd = GMP_RALG_CF_ALLOC_INCREMENT;
816    
817     //The number of slots actually used is zero, to start with.
818     decomp->n = 0;
819    
820     //There are a number of conditions that will lead to an error
821     //where we can't successfully form the continued fraction
822     //decomposition. These errors are:
823     // a)Either component is NAN.
824     // b)Zero denominator.
825     // c)Either component negative.
826     //In these cases, we'll pretend we got 0/1 for the number
827     //and set things accordingly, and we'll set the failure
828     //flag for the caller.
829     //
830     if (GMP_INTS_mpz_get_flags(num)
831     ||
832     GMP_INTS_mpz_get_flags(den)
833     ||
834     GMP_INTS_mpz_is_zero(den)
835     ||
836     GMP_INTS_mpz_is_neg(num)
837     ||
838     GMP_INTS_mpz_is_neg(den))
839     {
840     *failure = 1;
841     decomp->n = 1;
842    
843     GMP_INTS_mpz_set_ui(&(decomp->num), 0);
844     GMP_INTS_mpz_set_ui(&(decomp->den), 1);
845    
846     GMP_INTS_mpz_init(decomp->dividend);
847     GMP_INTS_mpz_set_ui(decomp->dividend, 0);
848    
849     GMP_INTS_mpz_init(decomp->divisor);
850     GMP_INTS_mpz_set_ui(decomp->divisor, 1);
851    
852     GMP_INTS_mpz_init(decomp->a);
853     GMP_INTS_mpz_set_ui(decomp->a, 0);
854    
855     GMP_INTS_mpz_init(decomp->remainder);
856     GMP_INTS_mpz_set_ui(decomp->remainder, 0);
857    
858     GMP_INTS_mpz_init(decomp->p);
859     GMP_INTS_mpz_set_ui(decomp->p, 0);
860    
861     GMP_INTS_mpz_init(decomp->q);
862     GMP_INTS_mpz_set_ui(decomp->q, 1);
863    
864     goto return_point;
865     }
866    
867     //If we're here there are not any errors that we
868     //are willing to detect. We should be clear
869     //for a continued fraction decomposition.
870     loop_counter = 0;
871     do
872     {
873     //Increment the count of "rows", because we're
874     //about to add one.
875     decomp->n++;
876    
877     //If we have used up all the space available
878     //for integers, we have to allocate more.
879     if (decomp->n > decomp->nallocd)
880     {
881     //We now have more allocated space.
882     decomp->nallocd += GMP_RALG_CF_ALLOC_INCREMENT;
883    
884     //Be absolutely sure we have not made a greivous
885     //error.
886     assert(decomp->n <= decomp->nallocd);
887    
888     //Space for dividends.
889     decomp->dividend =
890     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
891     CCMALLOC_realloc(
892     decomp->dividend,
893     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
894     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
895     (GMP_INTS_mpz_struct *)
896     TclpRealloc((char *)decomp->dividend,
897     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
898     #else
899     realloc(decomp->dividend, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
900     #endif
901    
902     //Space for divisors.
903     decomp->divisor =
904     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
905     CCMALLOC_realloc(
906     decomp->divisor,
907     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
908     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
909     (GMP_INTS_mpz_struct *)
910     TclpRealloc((char *)decomp->divisor,
911     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
912     #else
913     realloc(decomp->divisor, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
914     #endif
915    
916     //Space for partial quotients.
917     decomp->a =
918     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
919     CCMALLOC_realloc(
920     decomp->a,
921     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
922     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
923     (GMP_INTS_mpz_struct *)
924     TclpRealloc((char *)decomp->a,
925     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
926     #else
927     realloc(decomp->a, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
928     #endif
929    
930     //Space for remainders.
931     decomp->remainder =
932     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
933     CCMALLOC_realloc(
934     decomp->remainder,
935     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
936     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
937     (GMP_INTS_mpz_struct *)
938     TclpRealloc((char *)decomp->remainder,
939     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
940     #else
941     realloc(decomp->remainder, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
942     #endif
943    
944     //Space for partial quotient numerators.
945     decomp->p =
946     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
947     CCMALLOC_realloc(
948     decomp->p,
949     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
950     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
951     (GMP_INTS_mpz_struct *)
952     TclpRealloc((char *)decomp->p,
953     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
954     #else
955     realloc(decomp->p, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
956     #endif
957    
958     //Space for partial quotient denominators.
959     decomp->q =
960     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
961     CCMALLOC_realloc(
962     decomp->q,
963     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
964     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
965     (GMP_INTS_mpz_struct *)
966     TclpRealloc((char *)decomp->q,
967     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
968     #else
969     realloc(decomp->q, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
970     #endif
971     }
972    
973     //At this point, we have enough space to do the next round of operations.
974     //Set up an index variable.
975     i = decomp->n - 1;
976    
977     //Initialize all of the integers at this round.
978     GMP_INTS_mpz_init(decomp->dividend + i);
979     GMP_INTS_mpz_init(decomp->divisor + i);
980     GMP_INTS_mpz_init(decomp->a + i);
981     GMP_INTS_mpz_init(decomp->remainder + i);
982     GMP_INTS_mpz_init(decomp->p + i);
983     GMP_INTS_mpz_init(decomp->q + i);
984    
985     //Perform the next round of continued fraction decomposition. This
986     //is standard stuff.
987     if (i==0)
988     {
989     //In the 0th round, we essentially perform initial assignments.
990     GMP_INTS_mpz_copy(decomp->dividend, &(decomp->num));
991     GMP_INTS_mpz_copy(decomp->divisor, &(decomp->den));
992    
993     //Make the division to get quotient and remainder.
994     GMP_INTS_mpz_tdiv_qr(decomp->a, decomp->remainder, decomp->dividend, decomp->divisor);
995    
996     //The convergents in the first round are always the quotient over 1.
997     GMP_INTS_mpz_copy(decomp->p, decomp->a);
998     GMP_INTS_mpz_set_ui(decomp->q, 1);
999     }
1000     else if (i==1)
1001     {
1002     //In the 1st round, the only point of caution is that we have to
1003     //consider p(k-2)/q(k-2) carefully.
1004     GMP_INTS_mpz_copy(decomp->dividend + 1, decomp->divisor + 0);
1005     GMP_INTS_mpz_copy(decomp->divisor + 1, decomp->remainder + 0);
1006    
1007     //Make the division to get quotient and remainder.
1008     GMP_INTS_mpz_tdiv_qr(decomp->a + 1,
1009     decomp->remainder + 1,
1010     decomp->dividend + 1,
1011     decomp->divisor + 1);
1012    
1013     //Need to compute the numerator of the convergent. This will be:
1014     // a(1) p(0) + p(-1) = a(1)p(0) + 1.
1015     GMP_INTS_mpz_mul(decomp->p + 1, decomp->a + 1, decomp->p + 0);
1016     GMP_INTS_mpz_set_ui(&arb_temp1, 1);
1017     GMP_INTS_mpz_add(decomp->p + 1, decomp->p + 1, &arb_temp1);
1018    
1019     //Need to compute the denominator of the convergent. This will
1020     //be a(1)q(0) + q(-1) = a(1) q(0) = a(1).
1021     GMP_INTS_mpz_copy(decomp->q + 1, decomp->a + 1);
1022     }
1023     else
1024     {
1025     //In the general case, it is a simple formula.
1026     //Rotate in the dividend and divisor from the previous round.
1027     GMP_INTS_mpz_copy(decomp->dividend + i, decomp->divisor + i - 1);
1028     GMP_INTS_mpz_copy(decomp->divisor + i, decomp->remainder + i - 1);
1029    
1030     //Make the division to get quotient and remainder.
1031     GMP_INTS_mpz_tdiv_qr(decomp->a + i,
1032     decomp->remainder + i,
1033     decomp->dividend + i,
1034     decomp->divisor + i);
1035    
1036     //Need to compute the numerator of the convergent. This will be:
1037     // p(i) = a(i)p(i-1) + p(i-2)
1038     GMP_INTS_mpz_mul(decomp->p + i, decomp->a + i, decomp->p + i - 1);
1039     GMP_INTS_mpz_add(decomp->p + i, decomp->p + i, decomp->p + i - 2);
1040    
1041     //Need to compute the numerator of the convergent. This will be:
1042     // q(i) = q(i)q(i-1) + q(i-2)
1043     GMP_INTS_mpz_mul(decomp->q + i, decomp->a + i, decomp->q + i - 1);
1044     GMP_INTS_mpz_add(decomp->q + i, decomp->q + i, decomp->q + i - 2);
1045     }
1046    
1047     loop_counter++;
1048     } while(!GMP_INTS_mpz_is_zero(decomp->remainder + decomp->n - 1) && loop_counter < 100000);
1049    
1050     //In debug builds, be sure we did not terminate based on the loop counter.
1051     assert(loop_counter != 100000);
1052    
1053     return_point:
1054    
1055     //Deallocate space for temporary integers.
1056     GMP_INTS_mpz_clear(&arb_temp1);
1057     GMP_INTS_mpz_clear(&arb_temp2);
1058     }
1059    
1060    
1061     //08/16/01: Visual inspection OK.
1062     void GMP_RALG_cfdecomp_destroy(
1063     GMP_RALG_cf_app_struct *decomp
1064     )
1065     {
1066     int i;
1067    
1068     //Eyeball the input parameters. Other eyeballing
1069     //will be done as integers are destroyed.
1070     assert(decomp != NULL);
1071    
1072     //First, destroy the things that are bound directly
1073     //to the record.
1074     GMP_INTS_mpz_clear(&(decomp->num));
1075     GMP_INTS_mpz_clear(&(decomp->den));
1076    
1077     //Now, destroy every integer which is allocated.
1078     for (i=0; i < decomp->n; i++)
1079     {
1080     GMP_INTS_mpz_clear(decomp->dividend + i);
1081     GMP_INTS_mpz_clear(decomp->divisor + i);
1082     GMP_INTS_mpz_clear(decomp->a + i);
1083     GMP_INTS_mpz_clear(decomp->remainder + i);
1084     GMP_INTS_mpz_clear(decomp->p + i);
1085     GMP_INTS_mpz_clear(decomp->q + i);
1086     }
1087    
1088     //Now, destroy the arrays of integers.
1089     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1090     CCMALLOC_free(decomp->dividend);
1091     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1092     TclpFree((char *)decomp->dividend);
1093     #else
1094     free(decomp->dividend);
1095     #endif
1096     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1097     CCMALLOC_free(decomp->divisor);
1098     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1099     TclpFree((char *)decomp->divisor);
1100     #else
1101     free(decomp->divisor);
1102     #endif
1103     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1104     CCMALLOC_free(decomp->a);
1105     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1106     TclpFree((char *)decomp->a);
1107     #else
1108     free(decomp->a);
1109     #endif
1110     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1111     CCMALLOC_free(decomp->remainder);
1112     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1113     TclpFree((char *)decomp->remainder);
1114     #else
1115     free(decomp->remainder);
1116     #endif
1117     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1118     CCMALLOC_free(decomp->p);
1119     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1120     TclpFree((char *)decomp->p);
1121     #else
1122     free(decomp->p);
1123     #endif
1124     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1125     CCMALLOC_free(decomp->q);
1126     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1127     TclpFree((char *)decomp->q);
1128     #else
1129     free(decomp->q);
1130     #endif
1131     }
1132    
1133    
1134     /******************************************************************/
1135     /*** FORMATTED OUTPUT FUNCTIONS *********************************/
1136     /******************************************************************/
1137     //08/16/01: Visual inspection OK.
1138     void GMP_RALG_cfdecomp_emit(
1139     FILE *s,
1140     char *banner,
1141     GMP_RALG_cf_app_struct *decomp,
1142     int nf,
1143     int dap,
1144     const GMP_INTS_mpz_struct *dap_denominator)
1145     {
1146     int i;
1147    
1148     GMP_INTS_mpz_struct arb_temp, arb_quotient, arb_remainder;
1149    
1150     //Eyeball the input parameters. The banner is allowed to
1151     //be null, so can't check that.
1152     assert(s != NULL);
1153     assert(decomp != NULL);
1154    
1155     //Allocate our temporary integers.
1156     GMP_INTS_mpz_init(&arb_temp);
1157     GMP_INTS_mpz_init(&arb_quotient);
1158     GMP_INTS_mpz_init(&arb_remainder);
1159    
1160     //If banner requested and noformat option not used.
1161     if (banner && !nf)
1162     {
1163     FCMIOF_stream_bannerheading(s, banner, 1);
1164     }
1165    
1166     //Dump the input numerator.
1167     if (!nf)
1168     {
1169     GMP_INTS_mpz_long_int_format_to_stream(s,
1170     &(decomp->num),
1171     "Input Numerator");
1172     }
1173     else
1174     {
1175     GMP_INTS_mpz_arb_int_raw_to_stream(s, &(decomp->num));
1176     fprintf(s, "\n");
1177     }
1178    
1179     //Separator if not in unformatted mode.
1180     if (!nf)
1181     FCMIOF_stream_hline(s);
1182    
1183     //Dump the input denominator.
1184     if (!nf)
1185     {
1186     GMP_INTS_mpz_long_int_format_to_stream(s,
1187     &(decomp->den),
1188     "Input Denominator");
1189     }
1190     else
1191     {
1192     GMP_INTS_mpz_arb_int_raw_to_stream(s, &(decomp->den));
1193     fprintf(s, "\n");
1194     }
1195    
1196     //Separator if not in unformatted mode.
1197     if (!nf)
1198     FCMIOF_stream_hline(s);
1199    
1200     for (i=0; i<decomp->n; i++)
1201     {
1202     char strbuf[100];
1203     //Buffer to prepare description.
1204    
1205     //Print out the dividend at each round.
1206     if (!nf)
1207     {
1208     sprintf(strbuf, "dividend(%d)", i);
1209     GMP_INTS_mpz_long_int_format_to_stream(s,
1210     decomp->dividend + i,
1211     strbuf);
1212     }
1213     else
1214     {
1215     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->dividend+i);
1216     fprintf(s, "\n");
1217     }
1218    
1219     //Separator if not in unformatted mode.
1220     if (!nf)
1221     FCMIOF_stream_hline(s);
1222    
1223     //Print out the divisor at each round.
1224     if (!nf)
1225     {
1226     sprintf(strbuf, "divisor(%d)", i);
1227     GMP_INTS_mpz_long_int_format_to_stream(s,
1228     decomp->divisor + i,
1229     strbuf);
1230     }
1231     else
1232     {
1233     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->divisor+i);
1234     fprintf(s, "\n");
1235     }
1236    
1237     //Separator if not in unformatted mode.
1238     if (!nf)
1239     FCMIOF_stream_hline(s);
1240    
1241     //Print out partial quotient at each round.
1242     if (!nf)
1243     {
1244     sprintf(strbuf, "a(%d)", i);
1245     GMP_INTS_mpz_long_int_format_to_stream(s,
1246     decomp->a + i,
1247     strbuf);
1248     }
1249     else
1250     {
1251     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->a+i);
1252     fprintf(s, "\n");
1253     }
1254    
1255     //Separator if not in unformatted mode.
1256     if (!nf)
1257     FCMIOF_stream_hline(s);
1258    
1259     //It doesn't make any sense to print out the
1260     //remainder, because this becomes the divisor
1261     //for the next round. It is just wasted output
1262     //lines.
1263    
1264     //Print out the convergent numerator at
1265     //each round.
1266     if (!nf)
1267     {
1268     sprintf(strbuf, "p(%d)", i);
1269     GMP_INTS_mpz_long_int_format_to_stream(s,
1270     decomp->p + i,
1271     strbuf);
1272     }
1273     else
1274     {
1275     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->p+i);
1276     fprintf(s, "\n");
1277     }
1278    
1279     //Separator if not in unformatted mode.
1280     if (!nf)
1281     FCMIOF_stream_hline(s);
1282    
1283     //Print out the convergent denominator at
1284     //each round.
1285     if (!nf)
1286     {
1287     sprintf(strbuf, "q(%d)", i);
1288     GMP_INTS_mpz_long_int_format_to_stream(s,
1289     decomp->q + i,
1290     strbuf);
1291     }
1292     else
1293     {
1294     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->q+i);
1295     fprintf(s, "\n");
1296     }
1297    
1298     //Separator if not in unformatted mode.
1299     if (!nf)
1300     FCMIOF_stream_hline(s);
1301    
1302     if (dap)
1303     {
1304     //Calculate the DAP numerator
1305     GMP_INTS_mpz_mul(&arb_temp, dap_denominator, decomp->p + i);
1306     GMP_INTS_mpz_tdiv_qr(&arb_quotient, &arb_remainder,
1307     &arb_temp, decomp->q + i);
1308    
1309     //Print DAP numerator.
1310     if (!nf)
1311     {
1312     sprintf(strbuf, "dap_h(%d)", i);
1313     GMP_INTS_mpz_long_int_format_to_stream(s,
1314     &arb_quotient,
1315     strbuf);
1316     }
1317     else
1318     {
1319     GMP_INTS_mpz_arb_int_raw_to_stream(s, &arb_quotient);
1320     fprintf(s, "\n");
1321     }
1322    
1323     //Separator if not in unformatted mode.
1324     if (!nf)
1325     FCMIOF_stream_hline(s);
1326    
1327     //Print DAP denominator.
1328     if (!nf)
1329     {
1330     sprintf(strbuf, "dap_k(%d)", i);
1331     GMP_INTS_mpz_long_int_format_to_stream(s,
1332     dap_denominator,
1333     strbuf);
1334     }
1335     else
1336     {
1337     GMP_INTS_mpz_arb_int_raw_to_stream(s, dap_denominator);
1338     fprintf(s, "\n");
1339     }
1340    
1341     //Separator if not in unformatted mode.
1342     if (!nf)
1343     FCMIOF_stream_hline(s);
1344     }
1345     }
1346    
1347     //Deallocate our temporary integers.
1348     GMP_INTS_mpz_clear(&arb_temp);
1349     GMP_INTS_mpz_clear(&arb_quotient);
1350     GMP_INTS_mpz_clear(&arb_remainder);
1351     }
1352    
1353    
1354     /******************************************************************/
1355     /*** FAREY SERIES PREDECESSOR AND SUCCESSOR FUNCTIONS ***********/
1356     /******************************************************************/
1357     //08/16/01: Visual inspection OK.
1358     void GMP_RALG_farey_predecessor(
1359     GMP_RATS_mpq_struct *result,
1360     const GMP_RATS_mpq_struct *plus_two,
1361     const GMP_RATS_mpq_struct *plus_one,
1362     const GMP_INTS_mpz_struct *N)
1363     {
1364     GMP_RATS_mpq_struct result_copy;
1365     //Used to hold return value in case the result
1366     //is the same as either of the input arguments.
1367     GMP_INTS_mpz_struct temp1, temp2, floor_func;
1368     //Temporary integers.
1369    
1370     assert(result != NULL);
1371     assert(plus_two != NULL);
1372     assert(plus_one != NULL);
1373     assert(N != NULL);
1374    
1375     //Initialize the variables used.
1376     GMP_RATS_mpq_init(&result_copy);
1377     GMP_INTS_mpz_init(&temp1);
1378     GMP_INTS_mpz_init(&temp2);
1379     GMP_INTS_mpz_init(&floor_func);
1380    
1381     //Numerator of the term in the floor function.
1382     GMP_INTS_mpz_add(&temp1, &(plus_two->den), N);
1383    
1384     //Term in the floor function. This is used to
1385     //calculate both numerator and denominator, so we save it.
1386     GMP_INTS_mpz_tdiv_qr(&floor_func, &temp2, &temp1, &(plus_one->den));
1387    
1388     //Product of result of floor function and numerator--now
1389     //forming the numerator of the output.
1390     GMP_INTS_mpz_mul(&temp2, &floor_func, &(plus_one->num));
1391    
1392     //Final result assigned to numerator.
1393     GMP_INTS_mpz_sub(&(result_copy.num), &temp2, &(plus_two->num));
1394    
1395     //Product of result of floor function and denominator--now
1396     //forming the denominator of the output.
1397     GMP_INTS_mpz_mul(&temp2, &floor_func, &(plus_one->den));
1398    
1399     //Final result assigned to denominator.
1400     GMP_INTS_mpz_sub(&(result_copy.den), &temp2, &(plus_two->den));
1401    
1402     //Copy the result to the object owned by the caller.
1403     GMP_RATS_mpq_copy(result, &result_copy);
1404    
1405     //Deallocate dynamic memory.
1406     GMP_RATS_mpq_clear(&result_copy);
1407     GMP_INTS_mpz_clear(&temp1);
1408     GMP_INTS_mpz_clear(&temp2);
1409     GMP_INTS_mpz_clear(&floor_func);
1410     }
1411    
1412    
1413     //08/16/01: Visual inspection OK.
1414     void GMP_RALG_farey_successor(
1415     GMP_RATS_mpq_struct *result,
1416     const GMP_RATS_mpq_struct *minus_two,
1417     const GMP_RATS_mpq_struct *minus_one,
1418     const GMP_INTS_mpz_struct *N)
1419     {
1420     GMP_RATS_mpq_struct result_copy;
1421     //Used to hold return value in case the result
1422     //is the same as either of the input arguments.
1423     GMP_INTS_mpz_struct temp1, temp2, floor_func;
1424     //Temporary integers.
1425    
1426     assert(result != NULL);
1427     assert(minus_two != NULL);
1428     assert(minus_one != NULL);
1429     assert(N != NULL);
1430    
1431     //Initialize the variables used.
1432     GMP_RATS_mpq_init(&result_copy);
1433     GMP_INTS_mpz_init(&temp1);
1434     GMP_INTS_mpz_init(&temp2);
1435     GMP_INTS_mpz_init(&floor_func);
1436    
1437     //Numerator of the term in the floor function.
1438     GMP_INTS_mpz_add(&temp1, &(minus_two->den), N);
1439    
1440     //Term in the floor function. This is used to
1441     //calculate both numerator and denominator, so we save it.
1442     GMP_INTS_mpz_tdiv_qr(&floor_func, &temp2, &temp1, &(minus_one->den));
1443    
1444     //Product of result of floor function and numerator--now
1445     //forming the numerator of the output.
1446     GMP_INTS_mpz_mul(&temp2, &floor_func, &(minus_one->num));
1447    
1448     //Final result assigned to numerator.
1449     GMP_INTS_mpz_sub(&(result_copy.num), &temp2, &(minus_two->num));
1450    
1451     //Product of result of floor function and denominator--now
1452     //forming the denominator of the output.
1453     GMP_INTS_mpz_mul(&temp2, &floor_func, &(minus_one->den));
1454    
1455     //Final result assigned to denominator.
1456     GMP_INTS_mpz_sub(&(result_copy.den), &temp2, &(minus_two->den));
1457    
1458     //Copy the result to the object owned by the caller.
1459     GMP_RATS_mpq_copy(result, &result_copy);
1460    
1461     //Deallocate dynamic memory.
1462     GMP_RATS_mpq_clear(&result_copy);
1463     GMP_INTS_mpz_clear(&temp1);
1464     GMP_INTS_mpz_clear(&temp2);
1465     GMP_INTS_mpz_clear(&floor_func);
1466     }
1467    
1468    
1469     //08/16/01: Visual inspection OK.
1470     void GMP_RALG_enclosing_farey_neighbors(
1471     const GMP_RATS_mpq_struct *rn_in,
1472     const GMP_INTS_mpz_struct *N,
1473     const GMP_RALG_cf_app_struct *cf_rep,
1474     int *equality,
1475     GMP_RATS_mpq_struct *left,
1476     GMP_RATS_mpq_struct *right)
1477     {
1478     GMP_RATS_mpq_struct rn_abs;
1479     //Absolute value of rational number supplied.
1480     GMP_RATS_mpq_struct previous_convergent;
1481     //Convergent before the one that has the denominator
1482     //not exceeding the order of the series. Need to fudge
1483     //a little bit because don't have -1-th order convergents
1484     //tabulated.
1485     GMP_RATS_mpq_struct other_neighbor;
1486     //The other neighbor besides the highest-order convergent
1487     //without denominator too large.
1488     GMP_INTS_mpz_struct temp1, temp2, temp3, temp4;
1489     //Temporary integers.
1490     int ho_conv;
1491     //Index of highest-ordered convergent that does not have
1492     //denominator too large.
1493    
1494     //Eyeball the parameters.
1495     assert(rn_in != NULL);
1496     assert(N != NULL);
1497     assert(cf_rep != NULL);
1498     assert(equality != NULL);
1499     assert(left != NULL);
1500     assert(right != NULL);
1501    
1502     //Allocate dynamic variables.
1503     GMP_RATS_mpq_init(&rn_abs);
1504     GMP_RATS_mpq_init(&previous_convergent);
1505     GMP_RATS_mpq_init(&other_neighbor);
1506     GMP_INTS_mpz_init(&temp1);
1507     GMP_INTS_mpz_init(&temp2);
1508     GMP_INTS_mpz_init(&temp3);
1509     GMP_INTS_mpz_init(&temp4);
1510    
1511     //Zero is a troublesome case, because it requires us to
1512     //cross signs. Split this case out explicitly.
1513     if (GMP_INTS_mpz_is_zero(&(rn_in->num)))
1514     {
1515     *equality = 1; //0/1 a member of Farey series of any order.
1516     GMP_INTS_mpz_set_si(&(left->num), -1);
1517     GMP_INTS_mpz_copy(&(left->den), N);
1518     GMP_INTS_mpz_set_si(&(right->num), 1);
1519     GMP_INTS_mpz_copy(&(right->den), N);
1520     }
1521     else
1522     {
1523     //Make a copy of the rational number in. As a condition of
1524     //using this function, it must be normalized, but it still
1525     //may be negative. Our strategy is to treat the number as
1526     //positive, find the neighbors, then if it was negative
1527     //complement and re-order the neighbors. In other words,
1528     //find neighbors to a negative number by symmetry, not
1529     //by forming the CF representation of a negative number.
1530     //Also, we can't touch the input parameter.
1531     GMP_RATS_mpq_copy(&rn_abs, rn_in);
1532     GMP_INTS_mpz_abs(&(rn_abs.num));
1533    
1534     //Find the index of the highest-ordered convergent
1535     //with a denominator not exceeding the denominator of
1536     //the rational number supplied. The zero'th order
1537     //convergent has a denominator of 1, so that one
1538     //at least is safe.
1539    
1540     //Assign either the "left" or right
1541     //neighbor to be the highest-ordered
1542     //convergent with a denominator not exceeding the
1543     //denominator of the rational number input. I say
1544     //"either" because the properties of convergents let
1545     //us know based on the oddness or evenness of the order
1546     //which side it is on.
1547     ho_conv = 0;
1548     while (((ho_conv + 1) < cf_rep->n) && (GMP_INTS_mpz_cmp(cf_rep->q + ho_conv + 1, N) <= 0))
1549     {
1550     #if 0
1551     //Some questions about this loop--debugging output.
1552     printf("ho_conv : %d\n", ho_conv);
1553     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1554     cf_rep->q + ho_conv + 1,
1555     "decomp_den");
1556     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1557     &(rn_abs.den),
1558     "rn_in_den");
1559     printf("Compare result: %d\n\n", GMP_INTS_mpz_cmp(cf_rep->q + ho_conv + 1, &(rn_abs.den)));
1560     #endif
1561    
1562     ho_conv++;
1563     }
1564    
1565     if (INTFUNC_is_even(ho_conv))
1566     {
1567     GMP_INTS_mpz_copy(&(left->num), cf_rep->p + ho_conv);
1568     GMP_INTS_mpz_copy(&(left->den), cf_rep->q + ho_conv);
1569     }
1570     else
1571     {
1572     GMP_INTS_mpz_copy(&(right->num), cf_rep->p + ho_conv);
1573     GMP_INTS_mpz_copy(&(right->den), cf_rep->q + ho_conv);
1574     }
1575    
1576     //Now, we need to calculate the other neighbor based
1577     //on the standard formula. This is a little tricky
1578     //because we don't have the -1-th order convergents
1579     //tabulated so we have to fudge a little bit.
1580     if (ho_conv == 0)
1581     {
1582     GMP_RATS_mpq_set_si(&previous_convergent, 1, 0);
1583     }
1584     else
1585     {
1586     GMP_INTS_mpz_copy(&(previous_convergent.num), cf_rep->p + ho_conv - 1);
1587     GMP_INTS_mpz_copy(&(previous_convergent.den), cf_rep->q + ho_conv - 1);
1588     }
1589    
1590     //Calculate the other neighbor according to the standard
1591     //formula.
1592     GMP_INTS_mpz_sub(&temp1, N, &(previous_convergent.den));
1593     GMP_INTS_mpz_tdiv_qr(&temp2, &temp3, &temp1, cf_rep->q + ho_conv);
1594     //temp2 now contains term from floor() function in formula.
1595     GMP_INTS_mpz_mul(&temp1, &temp2, cf_rep->p + ho_conv);
1596     GMP_INTS_mpz_add(&(other_neighbor.num), &temp1, &(previous_convergent.num));
1597     GMP_INTS_mpz_mul(&temp1, &temp2, cf_rep->q + ho_conv);
1598     GMP_INTS_mpz_add(&(other_neighbor.den), &temp1, &(previous_convergent.den));
1599    
1600     //Copy the other neighbor into the right slot.
1601     if (INTFUNC_is_even(ho_conv))
1602     {
1603     GMP_RATS_mpq_copy(right, &other_neighbor);
1604     }
1605     else
1606     {
1607     GMP_RATS_mpq_copy(left, &other_neighbor);
1608     }
1609    
1610     //Set the equality flag. We have equality if and only
1611     //if the denominator of the rational number is less than
1612     //or equal to N.
1613     if (GMP_INTS_mpz_cmp(&(rn_abs.den), N) <= 0)
1614     {
1615     *equality = 1;
1616     }
1617     else
1618     {
1619     *equality = 0;
1620     }
1621    
1622     //In the event of equality, we don't really have the true
1623     //neighbors. If the final convergent is even-ordered,
1624     //the left needs to be replaced. If the final convergent
1625     //is odd-ordered, the right needs to be replaced.
1626     if (*equality)
1627     {
1628     if (INTFUNC_is_even(ho_conv))
1629     {
1630     //Left needs to be replaced.
1631     GMP_RALG_farey_predecessor(
1632     left,
1633     right,
1634     &rn_abs,
1635     N);
1636     }
1637     else
1638     {
1639     //Right needs to be replaced.
1640     GMP_RALG_farey_successor(
1641     right,
1642     left,
1643     &rn_abs,
1644     N);
1645     }
1646     }
1647    
1648     //OK, we should be all done. The final catch is that if
1649     //the number supplied was negative, we need to invert
1650     //and re-order the neighbors.
1651     if (GMP_INTS_mpz_is_neg(&(rn_in->num)))
1652     {
1653     GMP_RATS_mpq_swap(left, right);
1654     GMP_INTS_mpz_negate(&(left->num));
1655     GMP_INTS_mpz_negate(&(right->num));
1656     }
1657     } //End if (rn==0) else clause
1658    
1659     //Deallocate dynamic variables.
1660     GMP_RATS_mpq_clear(&rn_abs);
1661     GMP_RATS_mpq_clear(&previous_convergent);
1662     GMP_RATS_mpq_clear(&other_neighbor);
1663     GMP_INTS_mpz_clear(&temp1);
1664     GMP_INTS_mpz_clear(&temp2);
1665     GMP_INTS_mpz_clear(&temp3);
1666     GMP_INTS_mpz_clear(&temp4);
1667     }
1668    
1669    
1670    
1671     //08/16/01: Visual inspection OK. Did not fully inspect the
1672     //iterative part of this function. Unit testing will be
1673     //careful, expect that to catch any anomalies.
1674     void GMP_RALG_consecutive_fab_terms(
1675     const GMP_RATS_mpq_struct *rn_in,
1676     const GMP_INTS_mpz_struct *kmax,
1677     const GMP_INTS_mpz_struct *hmax,
1678     int n_left_in,
1679     int n_right_in,
1680     GMP_RALG_fab_neighbor_collection_struct *result
1681     )
1682     {
1683     int error_flag, equality_flag;
1684     char *estring_kmax_neg = "KMAX is zero, negative, or NAN.";
1685     char *estring_hmax_neg = "HMAX is negative or NAN.";
1686     char *estring_general = "Unspecified general error in GMP_RALG_consecutive_fab_terms().";
1687    
1688     GMP_RATS_mpq_struct q_temp1, q_temp2, q_temp3, q_temp4,
1689     left_neighbor, right_neighbor,
1690     left_neighbor_abs, right_neighbor_abs,
1691     hmax_over_one_neg, corner_point_neg,
1692     abs_norm_recip_rn;
1693    
1694     //Eyeball input parameters.
1695     assert(rn_in != NULL);
1696     assert(kmax != NULL);
1697     assert(n_left_in >= 0);
1698     assert(n_left_in <= 0x00FFFFFF);
1699     assert(n_right_in >= 0);
1700     assert(n_right_in <= 0x00FFFFFF);
1701     assert(result != NULL);
1702    
1703     //Allocate all of the dynamic memory we'll need for this function. It will be
1704     //released at the end.
1705     GMP_RATS_mpq_init(&q_temp1);
1706     GMP_RATS_mpq_init(&q_temp2);
1707     GMP_RATS_mpq_init(&q_temp3);
1708     GMP_RATS_mpq_init(&q_temp4);
1709     GMP_RATS_mpq_init(&left_neighbor);
1710     GMP_RATS_mpq_init(&right_neighbor);
1711     GMP_RATS_mpq_init(&left_neighbor_abs);
1712     GMP_RATS_mpq_init(&right_neighbor_abs);
1713     GMP_RATS_mpq_init(&hmax_over_one_neg);
1714     GMP_RATS_mpq_init(&corner_point_neg);
1715     GMP_RATS_mpq_init(&abs_norm_recip_rn);
1716    
1717     //Zero out the result block. This is the easiest way to give many variables
1718     //default values of 0, FALSE, and NULL.
1719     memset(result, 0, sizeof(GMP_RALG_fab_neighbor_collection_struct));
1720    
1721     //Allocate all integer and rational number structures in the result block.
1722     GMP_RATS_mpq_init(&(result->rn_in));
1723     GMP_INTS_mpz_init(&(result->kmax_in));
1724     GMP_INTS_mpz_init(&(result->hmax_in));
1725     GMP_RATS_mpq_init(&(result->hmax_over_one));
1726     GMP_RATS_mpq_init(&(result->corner_point));
1727     GMP_RATS_mpq_init(&(result->corner_point_minus_one));
1728     GMP_RATS_mpq_init(&(result->corner_point_plus_one));
1729     GMP_RATS_mpq_init(&(result->norm_rn));
1730     GMP_RATS_mpq_init(&(result->abs_norm_rn));
1731    
1732     //Fill in the rational number, exactly as passed.
1733     GMP_RATS_mpq_copy(&(result->rn_in), rn_in);
1734    
1735     //Fill in the number of left and right neighbors that the caller wants.
1736     //However, let's of course say nothing less than zero and nothing more
1737     //than 10000 neighbors on either side.
1738     result->n_left_in = INTFUNC_min(INTFUNC_max(0, n_left_in), 10000);
1739     result->n_right_in = INTFUNC_min(INTFUNC_max(0, n_right_in), 10000);
1740    
1741     //Fill in the value of KMAX, exactly as passed. If it is not at least
1742     //the value of 1 or if error flags, croak.
1743     GMP_INTS_mpz_copy(&(result->kmax_in), kmax);
1744     if (GMP_INTS_mpz_get_flags(kmax) || GMP_INTS_mpz_is_zero(kmax) || GMP_INTS_mpz_is_neg(kmax))
1745     {
1746     result->error =
1747     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1748     CCMALLOC_malloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1749     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1750     TclpAlloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1751     #else
1752     malloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1753     #endif
1754     strcpy(result->error, estring_kmax_neg);
1755     goto return_point;
1756     }
1757    
1758     //Decide whether the caller intends to use HMAX. Neg is error, but zero
1759     //is a signal that don't intend to use.
1760     if (hmax)
1761     {
1762     GMP_INTS_mpz_copy(&(result->hmax_in), hmax);
1763     if (GMP_INTS_mpz_get_flags(hmax) || GMP_INTS_mpz_is_neg(hmax))
1764     {
1765     result->error =
1766     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1767     CCMALLOC_malloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1768     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1769     TclpAlloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1770     #else
1771     malloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1772     #endif
1773     strcpy(result->error, estring_hmax_neg);
1774     goto return_point;
1775     }
1776     else if (GMP_INTS_mpz_is_pos(hmax))
1777     {
1778     result->hmax_supplied = 1;
1779     }
1780     }
1781    
1782     //If HMAX has been supplied, assign and normalize the
1783     //corner point.
1784     if (result->hmax_supplied)
1785     {
1786     GMP_INTS_mpz_copy(&(result->corner_point.num), &(result->hmax_in));
1787     GMP_INTS_mpz_copy(&(result->corner_point.den), &(result->kmax_in));
1788     GMP_RATS_mpq_normalize(&(result->corner_point));
1789     }
1790    
1791     //If HMAX has been supplied, we want to get the continued
1792     //fraction representation of both the corner point and its
1793     //reciprocal. This is because we're going to need to
1794     //find its adjacent points so we can easily crawl
1795     //around a rectangular region of the integer lattice.
1796     if (result->hmax_supplied)
1797     {
1798     //CF representation of corner point.
1799     GMP_RALG_cfdecomp_init(&(result->corner_point_cf_rep),
1800     &error_flag,
1801     &(result->corner_point.num),
1802     &(result->corner_point.den));
1803     result->corner_point_cf_rep_formed = 1;
1804    
1805     //If there was an error forming the CF representation
1806     //of the corner point, bail out.
1807     if (error_flag)
1808     {
1809     result->error =
1810     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1811     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1812     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1813     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1814     #else
1815     malloc(sizeof(char) * (strlen(estring_general) + 1));
1816     #endif
1817     strcpy(result->error, estring_general);
1818     goto return_point;
1819     }
1820    
1821     //CF representation of reciprocal of corner point.
1822     GMP_RALG_cfdecomp_init(&(result->corner_point_recip_cf_rep),
1823     &error_flag,
1824     &(result->corner_point.den),
1825     &(result->corner_point.num));
1826     result->corner_point_recip_cf_rep_formed = 1;
1827    
1828     //If there was an error forming the CF representation
1829     //of the reciprocal of the corner point, bail out.
1830     if (error_flag)
1831     {
1832     result->error =
1833     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1834     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1835     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1836     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1837     #else
1838     malloc(sizeof(char) * (strlen(estring_general) + 1));
1839     #endif
1840     strcpy(result->error, estring_general);
1841     goto return_point;
1842     }
1843     }
1844    
1845     //Normalize the rational number supplied.
1846     GMP_RATS_mpq_copy(&(result->norm_rn), rn_in);
1847     GMP_RATS_mpq_normalize(&(result->norm_rn));
1848    
1849     //Form the absolute value of the normalized
1850     //version, and set the neg flag.
1851     GMP_RATS_mpq_copy(&(result->abs_norm_rn),&(result->norm_rn));
1852     if (GMP_INTS_mpz_is_neg(&(result->abs_norm_rn.num)))
1853     {
1854     GMP_INTS_mpz_negate(&(result->abs_norm_rn.num));
1855     result->rn_is_neg = 1;
1856     }
1857    
1858     //Form the continued fraction representation of the
1859     //absolute value of the rational number supplied.
1860     //This is always required, because we cannot get any
1861     //neighbors without it.
1862     GMP_RALG_cfdecomp_init(&(result->rn_abs_cf_rep),
1863     &error_flag,
1864     &(result->abs_norm_rn.num),
1865     &(result->abs_norm_rn.den));
1866     result->rn_abs_cf_rep_formed = 1;
1867    
1868     //If there was an error forming the CF representation
1869     //of the absolute value of rational number supplied, bail out.
1870     if (error_flag)
1871     {
1872     result->error =
1873     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1874     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1875     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1876     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1877     #else
1878     malloc(sizeof(char) * (strlen(estring_general) + 1));
1879     #endif
1880     strcpy(result->error, estring_general);
1881     goto return_point;
1882     }
1883    
1884     //Set the equality flag. The rational number supplied is
1885     //in the series of interest if and only if, in its normalized
1886     //form, K <= KMAX, and if HMAX was supplied, H <= HMAX.
1887     if (GMP_INTS_mpz_cmp(&(result->abs_norm_rn.den), kmax) <= 0)
1888     {
1889     if (result->hmax_supplied)
1890     {
1891     if (GMP_INTS_mpz_cmp(&(result->abs_norm_rn.num), hmax) <= 0)
1892     {
1893     result->equality = 1;
1894     }
1895     else
1896     {
1897     result->equality = 0;
1898     }
1899     }
1900     else
1901     {
1902     result->equality = 1;
1903     }
1904     }
1905     else
1906     {
1907     result->equality = 0;
1908     }
1909    
1910     //The final cause of error is if the rational number
1911     //supplied is outside the interval [-HMAX/1, HMAX/1].
1912     //In such cases, simply refuse to calculate
1913     //any approximations. This error can only occur
1914     //if HMAX is specified. If only KMAX is specified,
1915     //this error cannot occur.
1916     if (result->hmax_supplied)
1917     {
1918     //Form the rational number HMAX/1. We will use it for
1919     //a comparison.
1920     GMP_INTS_mpz_copy(&(result->hmax_over_one.num), hmax);
1921     GMP_INTS_mpz_set_ui(&(result->hmax_over_one.den), 1);
1922    
1923     //If the comparison shows that the absolute value of
1924     //the rational number in is larger than HMAX over 1,
1925     //then declare an error.
1926     if (GMP_RATS_mpq_cmp(&(result->abs_norm_rn),&(result->hmax_over_one),NULL) > 0)
1927     {
1928     result->error =
1929     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1930     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1931     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1932     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1933     #else
1934     malloc(sizeof(char) * (strlen(estring_general) + 1));
1935     #endif
1936     strcpy(result->error, estring_general);
1937     goto return_point;
1938     }
1939     }
1940    
1941     //If we're here, we're very much clean. The only thing
1942     //that could go wrong is an overflow.
1943    
1944     //Allocate space for the left and right arrays of
1945     //neighbors.
1946     if (result->n_left_in)
1947     {
1948     result->lefts =
1949     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1950     CCMALLOC_malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1951     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1952     (GMP_RALG_fab_neighbor_struct *)
1953     TclpAlloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1954     #else
1955     malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1956     #endif
1957     }
1958    
1959     if (result->n_right_in)
1960     {
1961     result->rights =
1962     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1963     CCMALLOC_malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1964     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1965     (GMP_RALG_fab_neighbor_struct *)
1966     TclpAlloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1967     #else
1968     malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1969     #endif
1970     }
1971    
1972     //If the rational number supplied is above the corner
1973     //point, we want to form the continued fraction representation
1974     //of its reciprocal.
1975     if (result->hmax_supplied)
1976     {
1977     if (GMP_RATS_mpq_cmp(&(result->abs_norm_rn),&(result->corner_point),NULL) > 0)
1978     {
1979     GMP_RALG_cfdecomp_init(&(result->rn_abs_recip_cf_rep),
1980     &error_flag,
1981     &(result->abs_norm_rn.den),
1982     &(result->abs_norm_rn.num));
1983     result->rn_abs_recip_cf_rep_formed = 1;
1984     }
1985     }
1986    
1987     //If HMAX has been supplied, we want to calculate the points just below and above
1988     //the corner point. The reason we want to do this is because we need to gracefully
1989     //"round the corner" in either direction.
1990     //
1991     //Calculate the point just to the left of the corner point.
1992     if (result->hmax_supplied)
1993     {
1994     GMP_RALG_enclosing_farey_neighbors(
1995     &(result->corner_point),
1996     &(result->kmax_in),
1997     &(result->corner_point_cf_rep),
1998     &equality_flag,
1999     &(result->corner_point_minus_one),
2000     &(q_temp1)
2001     );
2002     }
2003    
2004     //Calculate the point just to the right of the corner point. This is
2005     //where HMAX is the dominant constraint. We need to find the left
2006     //Farey neighbor to the reciprocal of the corner point in the Farey
2007     //series of order HMAX, then take its reciprocal. There is the possibility
2008     //if KMAX=1 that this point will have a denominator of zero, but we
2009     //will accept that as a number here, since we should never hit it,
2010     //as it will be beyond HMAX/1.
2011     if (result->hmax_supplied)
2012     {
2013     GMP_RATS_mpq_copy(&q_temp1, &(result->corner_point));
2014     GMP_INTS_mpz_swap(&(q_temp1.num), &(q_temp1.den));
2015     GMP_RALG_enclosing_farey_neighbors(
2016     &q_temp1,
2017     &(result->hmax_in),
2018     &(result->corner_point_recip_cf_rep),
2019     &equality_flag,
2020     &(result->corner_point_plus_one),
2021     &(q_temp2)
2022     );
2023     GMP_INTS_mpz_swap(&(result->corner_point_plus_one.num), &(result->corner_point_plus_one.den));
2024     }
2025    
2026     //Calculate the complement of HMAX/1. Nothing that we generate can go beyond
2027     //this to the south.
2028     if (result->hmax_supplied)
2029     {
2030     GMP_RATS_mpq_copy(&(hmax_over_one_neg), &(result->hmax_over_one));
2031     GMP_INTS_mpz_negate(&(hmax_over_one_neg.num));
2032     }
2033    
2034     //Also calculate the complement of HMAX/KMAX.
2035     if (result->hmax_supplied)
2036     {
2037     GMP_RATS_mpq_copy(&(corner_point_neg), &(result->corner_point));
2038     GMP_INTS_mpz_negate(&(corner_point_neg.num));
2039     }
2040    
2041     //Form the reciprocal of the absolute value of the normalized value of
2042     //the rational number in.
2043     GMP_RATS_mpq_copy(&abs_norm_recip_rn, &(result->abs_norm_rn));
2044     GMP_RATS_mpq_swap_components(&abs_norm_recip_rn);
2045    
2046     //OK, now we get down to brass tacks. Iterate first to get the
2047     //left neighbors. The ordinary complexity of this is also compounded
2048     //by the fact that we must handle negative numbers as well--everything
2049     //from -HMAX/1 to HMAX/1.
2050     //
2051     //PSEUDO-CODE:
2052     // a)If the rational number to approximate is <= -HMAX/1 or there are no
2053     // left neighbors requested, terminate with no neighbors on the left.
2054     //
2055     // b)Find the right neighbor of the rational number supplied.
2056     //
2057     // c)Find the left neighbor of the rational number supplied.
2058     //
2059     // d)While (queued_count < count)
2060     //
2061     // e)Queue the left neighbor, queued_count++
2062     //
2063     // f)If (queued_count >= count), break.
2064     //
2065     // g)If (left_neighbor <= -HMAX/1), break
2066     //
2067     // h)Advance the frame one to the left.
2068     //
2069     //**************************************************************************
2070     // a)If the rational number to approximate is <= -HMAX/1 or there are no
2071     // left neighbors requested, terminate with no neighbors on the left.
2072     //**************************************************************************
2073     if ((result->hmax_supplied && GMP_RATS_mpq_cmp(&(result->norm_rn), &hmax_over_one_neg, NULL) <= 0)
2074     || (n_left_in <= 0))
2075     goto done_with_left_neighbors;
2076    
2077     //**************************************************************************
2078     // b)Find the right neighbor of the rational number supplied.
2079     //**************************************************************************
2080     // c)Find the left neighbor of the rational number supplied.
2081     //**************************************************************************
2082     if (!result->hmax_supplied)
2083     {
2084     //In this case, the notion of corner point is meaningless, because
2085     //there is no constraint on H. We can just go on our merry way. Get
2086     //the two neighbors.
2087     GMP_RALG_enclosing_farey_neighbors(
2088     &(result->norm_rn),
2089     &(result->kmax_in),
2090     &(result->rn_abs_cf_rep),
2091     &equality_flag,
2092     &left_neighbor,
2093     &right_neighbor
2094     );
2095     //The enclosing Farey neighbor function is prohibited from identifying the
2096     //rational number itself as a Farey term. If the number is in the Farey
2097     //series, we must replace the right neighbor, otherwise we cannot apply
2098     //the standard recursive formulas.
2099     if (equality_flag)
2100     {
2101     GMP_RATS_mpq_copy(&right_neighbor, &(result->norm_rn));
2102     }
2103     }
2104     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) < 0)
2105     {
2106     //The rational number specified is negative and below the negative corner point.
2107     //This means that HMAX is the dominant constraint. We need to find the
2108     //neighbors in the Farey series of order HMAX, then reorder and invert, etc.
2109     GMP_RALG_enclosing_farey_neighbors(
2110     &abs_norm_recip_rn,
2111     &(result->hmax_in),
2112     &(result->rn_abs_recip_cf_rep),
2113     &equality_flag,
2114     &left_neighbor,
2115     &right_neighbor
2116     );
2117    
2118     //Again, if the number specified was already in the series of interest,
2119     //we need to swap in the right neighbor.
2120     if (equality_flag)
2121     {
2122     GMP_RATS_mpq_copy(&right_neighbor, &abs_norm_recip_rn);
2123     }
2124    
2125     //Take the reciprocal of both neighbors, which will put them out of order,
2126     //then negate them, which will put them back in order.
2127     GMP_RATS_mpq_swap_components(&left_neighbor);
2128     GMP_INTS_mpz_negate(&(left_neighbor.num));
2129     GMP_RATS_mpq_swap_components(&right_neighbor);
2130     GMP_INTS_mpz_negate(&(right_neighbor.num));
2131     }
2132     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) == 0)
2133     {
2134     //The rational number specified is the negative corner point. In this case
2135     //Because we can never return the corner point itself as a left neighbor,
2136     //we need to set the left value to be the negative of one past, and the right
2137     //to be the negative of the corner point.
2138     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_plus_one));
2139     GMP_INTS_mpz_negate(&(left_neighbor.num));
2140     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point));
2141     GMP_INTS_mpz_negate(&(right_neighbor.num));
2142     }
2143     else if ((GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) > 0)
2144     &&
2145     (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) < 0))
2146     {
2147     //The rational number specified is in the area dominated by the KMAX constraint
2148     //between -HMAX/KMAX and HMAX/KMAX. The ordinary Farey neighbor function will
2149     //handle this correctly. Again, we need to adjust the output if the number
2150     //is already formable, because the Farey neighbor function is prohibited from
2151     //returning the number itself as a neighbor.
2152     GMP_RALG_enclosing_farey_neighbors(
2153     &(result->norm_rn),
2154     &(result->kmax_in),
2155     &(result->rn_abs_cf_rep),
2156     &equality_flag,
2157     &left_neighbor,
2158     &right_neighbor
2159     );
2160     //The enclosing Farey neighbor function is prohibited from identifying the
2161     //rational number itself as a Farey term. If the number is in the Farey
2162     //series, we must replace the right neighbor, otherwise we cannot apply
2163     //the standard recursive formulas.
2164     if (equality_flag)
2165     {
2166     GMP_RATS_mpq_copy(&right_neighbor, &(result->norm_rn));
2167     }
2168     }
2169     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) == 0)
2170     {
2171     //The rational number specified is the corner point. In this case
2172     //because we can never return the corner point itself as a left neighbor,
2173     //we need to set the left value to be one before, and the right
2174     //to be the corner point.
2175     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_minus_one));
2176     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point));
2177     }
2178     else
2179     {
2180     //The only possibility left is that the number is positive and above the
2181     //corner point where HMAX is the dominant constraint.
2182     GMP_RALG_enclosing_farey_neighbors(
2183     &abs_norm_recip_rn,
2184     &(result->hmax_in),
2185     &(result->rn_abs_recip_cf_rep),
2186     &equality_flag,
2187     &left_neighbor,
2188     &right_neighbor
2189     );
2190    
2191     //Again, if the number specified was already in the series of interest,
2192     //we need to swap in the neighbor. This time, however, it must be the
2193     //left neighbor because taking the reciprocals will reverse the order.
2194     if (equality_flag)
2195     {
2196     GMP_RATS_mpq_copy(&left_neighbor, &abs_norm_recip_rn);
2197     }
2198    
2199     //Take the reciprocal of both neighbors, which will put them out of order,
2200     //then swap them, which will put them back in order.
2201     GMP_RATS_mpq_swap_components(&left_neighbor);
2202     GMP_RATS_mpq_swap_components(&right_neighbor);
2203     GMP_RATS_mpq_swap(&left_neighbor, &right_neighbor);
2204     }
2205    
2206     #if 0
2207     //Print out the left neighbor and right neighbor determined, for debugging.
2208     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2209     &(left_neighbor.num),
2210     "left_neigh_num");
2211     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2212     &(left_neighbor.den),
2213     "left_neigh_den");
2214     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2215     &(right_neighbor.num),
2216     "right_neigh_num");
2217     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2218     &(right_neighbor.den),
2219     "right_neigh_den");
2220     #endif
2221    
2222    
2223     //**************************************************************************
2224     // d)While (queued_count < count)
2225     //**************************************************************************
2226     while (result->n_left_out < result->n_left_in)
2227     {
2228     //**************************************************************************
2229     // e)Queue the left neighbor, queued_count++
2230     //**************************************************************************
2231     (result->lefts + result->n_left_out)->index = -(result->n_left_out + 1);
2232     GMP_RATS_mpq_init(&((result->lefts + result->n_left_out)->neighbor));
2233     GMP_RATS_mpq_copy(&((result->lefts + result->n_left_out)->neighbor), &left_neighbor);
2234     (result->n_left_out)++;
2235    
2236     //**************************************************************************
2237     // f)If (queued_count >= count), break.
2238     //**************************************************************************
2239     //By the way, this step is to save unnecessary contortions once we've met
2240     //the quota.
2241     if (result->n_left_out >= result->n_left_in)
2242     break;
2243    
2244     //**************************************************************************
2245     // g)If (left_neighbor <= -HMAX/1), break
2246     //**************************************************************************
2247     //This breaks us when we've queued the most negative number we can--can't go
2248     //further. This only applies for cases where KMAX is also specified.
2249     if (result->hmax_supplied
2250     &&
2251     GMP_RATS_mpq_cmp(&left_neighbor, &hmax_over_one_neg, NULL) <= 0)
2252     break;
2253    
2254     //**************************************************************************
2255     // h)Advance the frame one to the left.
2256     //**************************************************************************
2257     //Advancing one frame to the left is a complicated affair, requiring several
2258     //subcases. We break it up into regions which are best visualized using
2259     //a graph of the integer lattice with dots for each rational number.
2260     if (!(result->hmax_supplied))
2261     {
2262     //This is the case where we're are looking only in the
2263     //Farey series.
2264     if (GMP_INTS_mpz_is_pos(&left_neighbor.num))
2265     {
2266     //In this case, the left neighbor and right neighbor
2267     //are both positive, and we can apply the standard
2268     //formulas.
2269     GMP_RALG_farey_predecessor(&q_temp1,
2270     &right_neighbor,
2271     &left_neighbor,
2272     &(result->kmax_in));
2273     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2274     GMP_RATS_mpq_copy(&left_neighbor, &q_temp1);
2275     }
2276     else if (GMP_INTS_mpz_is_zero(&left_neighbor.num))
2277     {
2278     //In this case, we are making the transition from positive
2279     //to negative.
2280     GMP_INTS_mpz_set_si(&(left_neighbor.num), -1);
2281     GMP_INTS_mpz_copy(&(left_neighbor.den), &(result->kmax_in));
2282     GMP_RATS_mpq_set_si(&right_neighbor, 0, 1);
2283     }
2284     else
2285     {
2286     //Here, we are purely negative and decreasing. Need to negate
2287     //the numbers, find successor, then negate.
2288     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2289     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2290     GMP_INTS_mpz_abs(&(q_temp1.num));
2291     GMP_INTS_mpz_abs(&(q_temp2.num));
2292     GMP_RALG_farey_successor(&q_temp3,
2293     &q_temp2,
2294     &q_temp1,
2295     &(result->kmax_in));
2296     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2297     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
2298     GMP_INTS_mpz_negate(&(left_neighbor.num));
2299     }
2300     }
2301     else if (GMP_RATS_mpq_cmp(&left_neighbor, &(result->corner_point), NULL) > 0)
2302     {
2303     //We are above the top corner point. In this case HMAX is the dominant
2304     //constraint, and we find our food by taking reciprocals and applying
2305     //the standard relationships in the Farey series of order HMAX.
2306     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2307     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2308     GMP_RATS_mpq_swap_components(&q_temp1);
2309     GMP_RATS_mpq_swap_components(&q_temp2);
2310     GMP_RALG_farey_successor(&q_temp3,
2311     &q_temp2,
2312     &q_temp1,
2313     &(result->hmax_in));
2314     GMP_RATS_mpq_swap_components(&q_temp3);
2315     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2316     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
2317     }
2318     else if (GMP_RATS_mpq_cmp(&left_neighbor, &(result->corner_point), NULL) == 0)
2319     {
2320     //We are precisely at the corner point. This is where we round the corner.
2321     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2322     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_minus_one));
2323     }
2324     else if (GMP_INTS_mpz_is_pos(&left_neighbor.num))
2325     {
2326     //In this region we are going straight down the Farey series.
2327     GMP_RALG_farey_predecessor(&q_temp1,
2328     &right_neighbor,
2329     &left_neighbor,
2330     &(result->kmax_in));
2331     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2332     GMP_RATS_mpq_copy(&left_neighbor, &q_temp1);
2333     }
2334     else if (GMP_INTS_mpz_is_zero(&left_neighbor.num))
2335     {
2336     //In this case, we are making the transition from positive
2337     //to negative.
2338     GMP_INTS_mpz_set_si(&(left_neighbor.num), -1);
2339     GMP_INTS_mpz_copy(&(left_neighbor.den), &(result->kmax_in));
2340     GMP_RATS_mpq_set_si(&right_neighbor, 0, 1);
2341     }
2342     else if (GMP_RATS_mpq_cmp(&left_neighbor, &corner_point_neg, NULL) > 0)
2343     {
2344     //Here, we are purely negative and decreasing. Need to negate
2345     //the numbers, find successor, then negate.
2346     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2347     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2348     GMP_INTS_mpz_abs(&(q_temp1.num));
2349     GMP_INTS_mpz_abs(&(q_temp2.num));
2350     GMP_RALG_farey_successor(&q_temp3,
2351     &q_temp2,
2352     &q_temp1,
2353     &(result->kmax_in));
2354     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2355     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
2356     GMP_INTS_mpz_negate(&(left_neighbor.num));
2357     }
2358     else if (GMP_RATS_mpq_cmp(&left_neighbor, &corner_point_neg, NULL) == 0)
2359     {
2360     //This is where we are rounding the negative corner.
2361     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2362     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_plus_one));
2363     GMP_INTS_mpz_negate(&(left_neighbor.num));
2364     }
2365     else
2366     {
2367     //Here we're going in the negative direction along the "bottom" of the
2368     //rectangle.
2369     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2370     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2371     GMP_INTS_mpz_abs(&(q_temp1.num));
2372     GMP_INTS_mpz_abs(&(q_temp2.num));
2373     GMP_RATS_mpq_swap_components(&q_temp1);
2374     GMP_RATS_mpq_swap_components(&q_temp2);
2375     GMP_RALG_farey_predecessor(&q_temp3,
2376     &q_temp2,
2377     &q_temp1,
2378     &(result->hmax_in));
2379     GMP_RATS_mpq_swap_components(&q_temp3);
2380     GMP_INTS_mpz_negate(&(q_temp3.num));
2381     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
2382     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
2383     }
2384     }
2385    
2386    
2387     done_with_left_neighbors: ;
2388    
2389     //**************************************************************************
2390     // a)If the rational number to approximate is >= HMAX/1 or there are no
2391     // right neighbors requested, terminate with no neighbors on the right.
2392     //**************************************************************************
2393     if ((result->hmax_supplied && GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->hmax_over_one), NULL) >= 0)
2394     || (n_right_in <= 0))
2395     goto done_with_right_neighbors;
2396    
2397     //**************************************************************************
2398     // b)Find the right neighbor of the rational number supplied.
2399     //**************************************************************************
2400     // c)Find the left neighbor of the rational number supplied.
2401     //**************************************************************************
2402     if (!result->hmax_supplied)
2403     {
2404     //In this case, the notion of corner point is meaningless, because
2405     //there is no constraint on H. We can just go on our merry way. Get
2406     //the two neighbors.
2407     GMP_RALG_enclosing_farey_neighbors(
2408     &(result->norm_rn),
2409     &(result->kmax_in),
2410     &(result->rn_abs_cf_rep),
2411     &equality_flag,
2412     &left_neighbor,
2413     &right_neighbor
2414     );
2415     //The enclosing Farey neighbor function is prohibited from identifying the
2416     //rational number itself as a Farey term. If the number is in the Farey
2417     //series, we must replace the left neighbor, otherwise we cannot apply
2418     //the standard recursive formulas.
2419     if (equality_flag)
2420     {
2421     GMP_RATS_mpq_copy(&left_neighbor, &(result->norm_rn));
2422     }
2423     }
2424     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) < 0)
2425     {
2426     //The rational number specified is negative and below the negative corner point.
2427     //This means that HMAX is the dominant constraint. We need to find the
2428     //neighbors in the Farey series of order HMAX, then reorder and invert, etc.
2429     GMP_RALG_enclosing_farey_neighbors(
2430     &abs_norm_recip_rn,
2431     &(result->hmax_in),
2432     &(result->rn_abs_recip_cf_rep),
2433     &equality_flag,
2434     &left_neighbor,
2435     &right_neighbor
2436     );
2437    
2438     //Again, if the number specified was already in the series of interest,
2439     //we need to swap in the left neighbor.
2440     if (equality_flag)
2441     {
2442     GMP_RATS_mpq_copy(&left_neighbor, &abs_norm_recip_rn);
2443     }
2444    
2445     //Take the reciprocal of both neighbors, which will put them out of order,
2446     //then negate them, which will put them back in order.
2447     GMP_RATS_mpq_swap_components(&left_neighbor);
2448     GMP_INTS_mpz_negate(&(left_neighbor.num));
2449     GMP_RATS_mpq_swap_components(&right_neighbor);
2450     GMP_INTS_mpz_negate(&(right_neighbor.num));
2451     }
2452     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) == 0)
2453     {
2454     //The rational number specified is the negative corner point. In this case
2455     //Because we can never return the corner point itself as a left neighbor,
2456     //we need to set the right value to be the negative of one before, and the left
2457     //to be the negative of the corner point.
2458     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point));
2459     GMP_INTS_mpz_negate(&(left_neighbor.num));
2460     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_minus_one));
2461     GMP_INTS_mpz_negate(&(right_neighbor.num));
2462     }
2463     else if ((GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) > 0)
2464     &&
2465     (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) < 0))
2466     {
2467     //The rational number specified is in the area dominated by the KMAX constraint
2468     //between -HMAX/KMAX and HMAX/KMAX. The ordinary Farey neighbor function will
2469     //handle this correctly. Again, we need to adjust the output if the number
2470     //is already formable, because the Farey neighbor function is prohibited from
2471     //returning the number itself as a neighbor.
2472     GMP_RALG_enclosing_farey_neighbors(
2473     &(result->norm_rn),
2474     &(result->kmax_in),
2475     &(result->rn_abs_cf_rep),
2476     &equality_flag,
2477     &left_neighbor,
2478     &right_neighbor
2479     );
2480     //The enclosing Farey neighbor function is prohibited from identifying the
2481     //rational number itself as a Farey term. If the number is in the Farey
2482     //series, we must replace the left neighbor, otherwise we cannot apply
2483     //the standard recursive formulas.
2484     if (equality_flag)
2485     {
2486     GMP_RATS_mpq_copy(&left_neighbor, &(result->norm_rn));
2487     }
2488     }
2489     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) == 0)
2490     {
2491     //The rational number specified is the positive corner point. In this case.
2492     //because we can never return the corner point itself as a right neighbor,
2493     //we need to set the right value to be one after, and the left
2494     //to be the corner point.
2495     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point));
2496     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_plus_one));
2497     }
2498     else
2499     {
2500     //The only possibility left is that the number is positive and at or above the
2501     //corner point where HMAX is the dominant constraint.
2502     GMP_RALG_enclosing_farey_neighbors(
2503     &abs_norm_recip_rn,
2504     &(result->hmax_in),
2505     &(result->rn_abs_recip_cf_rep),
2506     &equality_flag,
2507     &left_neighbor,
2508     &right_neighbor
2509     );
2510    
2511     //Again, if the number specified was already in the series of interest,
2512     //we need to swap in the neighbor. This time, however, it must be the
2513     //right neighbor because taking the reciprocals will reverse the order.
2514     if (equality_flag)
2515     {
2516     GMP_RATS_mpq_copy(&right_neighbor, &abs_norm_recip_rn);
2517     }
2518    
2519     //Take the reciprocal of both neighbors, which will put them out of order,
2520     //then swap them, which will put them back in order.
2521     GMP_RATS_mpq_swap_components(&left_neighbor);
2522     GMP_RATS_mpq_swap_components(&right_neighbor);
2523     GMP_RATS_mpq_swap(&left_neighbor, &right_neighbor);
2524     }
2525    
2526     #if 0
2527     //Print out the left neighbor and right neighbor determined, for debugging.
2528     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2529     &(left_neighbor.num),
2530     "left_neigh_num");
2531     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2532     &(left_neighbor.den),
2533     "left_neigh_den");
2534     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2535     &(right_neighbor.num),
2536     "right_neigh_num");
2537     GMP_INTS_mpz_long_int_format_to_stream(stdout,
2538     &(right_neighbor.den),
2539     "right_neigh_den");
2540     #endif
2541    
2542    
2543     //**************************************************************************
2544     // d)While (queued_count < count)
2545     //**************************************************************************
2546     while (result->n_right_out < result->n_right_in)
2547     {
2548     //**************************************************************************
2549     // e)Queue the right neighbor, queued_count++
2550     //**************************************************************************
2551     (result->rights + result->n_right_out)->index = (result->n_right_out + 1);
2552     GMP_RATS_mpq_init(&((result->rights + result->n_right_out)->neighbor));
2553     GMP_RATS_mpq_copy(&((result->rights + result->n_right_out)->neighbor), &right_neighbor);
2554     (result->n_right_out)++;
2555    
2556     //**************************************************************************
2557     // f)If (queued_count >= count), break.
2558     //**************************************************************************
2559     //By the way, this step is to save unnecessary contortions once we've met
2560     //the quota.
2561     if (result->n_right_out >= result->n_right_in)
2562     break;
2563    
2564     //**************************************************************************
2565     // g)If (right_neighbor >= HMAX/1), break
2566     //**************************************************************************
2567     //This breaks us when we've queued the most positive number we can--can't go
2568     //further. This only applies for cases where KMAX is also specified.
2569     if (result->hmax_supplied
2570     &&
2571     GMP_RATS_mpq_cmp(&right_neighbor, &(result->hmax_over_one), NULL) >= 0)
2572     break;
2573    
2574     //**************************************************************************
2575     // h)Advance the frame one to the right.
2576     //**************************************************************************
2577     //Advancing one frame to the right is a complicated affair, requiring several
2578     //subcases. We break it up into regions which are best visualized using
2579     //a graph of the integer lattice with dots for each rational number.
2580     if (!(result->hmax_supplied))
2581     {
2582     //This is the case where we're are looking only in the
2583     //Farey series.
2584     if (GMP_INTS_mpz_is_neg(&right_neighbor.num))
2585     {
2586     //Neg and increasing towards zero. Can handle by symmetry.
2587     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2588     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2589     GMP_INTS_mpz_abs(&(q_temp1.num));
2590     GMP_INTS_mpz_abs(&(q_temp2.num));
2591     GMP_RALG_farey_predecessor(&q_temp3,
2592     &q_temp1,
2593     &q_temp2,
2594     &(result->kmax_in));
2595     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2596     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2597     GMP_INTS_mpz_negate(&(right_neighbor.num));
2598     }
2599     else if (GMP_INTS_mpz_is_zero(&right_neighbor.num))
2600     {
2601     //Right term just hit zero and are increasing.
2602     //Left will become 0/1, right becomes 1/KMAX.
2603     GMP_RATS_mpq_set_si(&left_neighbor, 0, 1);
2604     GMP_INTS_mpz_set_si(&(right_neighbor.num), 1);
2605     GMP_INTS_mpz_copy(&(right_neighbor.den), &(result->kmax_in));
2606     }
2607     else
2608     {
2609     //Are above zero and increasing. Can use standard Farey
2610     //successor formula.
2611     GMP_RALG_farey_successor(&q_temp1,
2612     &left_neighbor,
2613     &right_neighbor,
2614     &(result->kmax_in));
2615     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2616     GMP_RATS_mpq_copy(&right_neighbor, &q_temp1);
2617     }
2618     }
2619     else if (GMP_RATS_mpq_cmp(&right_neighbor, &corner_point_neg, NULL) < 0)
2620     {
2621     //We are below the negative corner point and moving towards
2622     //zero, with HMAX the dominant constraint. We can proceed by
2623     //symmetry, producing a Farey successor and negating and inverting.
2624     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2625     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2626     GMP_INTS_mpz_abs(&(q_temp1.num));
2627     GMP_INTS_mpz_abs(&(q_temp2.num));
2628     GMP_RATS_mpq_swap_components(&q_temp1);
2629     GMP_RATS_mpq_swap_components(&q_temp2);
2630     GMP_RALG_farey_successor(&q_temp3,
2631     &q_temp1,
2632     &q_temp2,
2633     &(result->hmax_in));
2634     GMP_RATS_mpq_swap_components(&q_temp3);
2635     GMP_INTS_mpz_negate(&(q_temp3.num));
2636     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2637     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2638     }
2639     else if (GMP_RATS_mpq_cmp(&right_neighbor, &corner_point_neg, NULL) == 0)
2640     {
2641     //We are at the bottom negative corner point and need to "round the corner".
2642     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2643     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_minus_one));
2644     GMP_INTS_mpz_negate(&(right_neighbor.num));
2645     }
2646     else if (GMP_INTS_mpz_is_neg(&right_neighbor.num))
2647     {
2648     //In this region we are going straight up the Farey series approaching
2649     //zero. Need to negate to use standard relationships.
2650     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2651     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2652     GMP_INTS_mpz_abs(&(q_temp1.num));
2653     GMP_INTS_mpz_abs(&(q_temp2.num));
2654     GMP_RALG_farey_predecessor(&q_temp3,
2655     &q_temp1,
2656     &q_temp2,
2657     &(result->kmax_in));
2658     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2659     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2660     GMP_INTS_mpz_negate(&(right_neighbor.num));
2661     }
2662     else if (GMP_INTS_mpz_is_zero(&right_neighbor.num))
2663     {
2664     //Zero crossover.
2665     GMP_RATS_mpq_set_si(&left_neighbor, 0, 1);
2666     GMP_INTS_mpz_set_si(&(right_neighbor.num), 1);
2667     GMP_INTS_mpz_copy(&(right_neighbor.den), &(result->kmax_in));
2668     }
2669     else if (GMP_RATS_mpq_cmp(&right_neighbor, &(result->corner_point), NULL) < 0)
2670     {
2671     //Below corner point. Standard relationship applies.
2672     GMP_RALG_farey_successor(&q_temp1,
2673     &left_neighbor,
2674     &right_neighbor,
2675     &(result->kmax_in));
2676     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2677     GMP_RATS_mpq_copy(&right_neighbor, &q_temp1);
2678     }
2679     else if (GMP_RATS_mpq_cmp(&right_neighbor, &(result->corner_point), NULL) == 0)
2680     {
2681     //At the positive corner point.
2682     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2683     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_plus_one));
2684     }
2685     else
2686     {
2687     //Above the positive corner point and heading for HMAX/1.
2688     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2689     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2690     GMP_RATS_mpq_swap_components(&q_temp1);
2691     GMP_RATS_mpq_swap_components(&q_temp2);
2692     GMP_RALG_farey_predecessor(&q_temp3,
2693     &q_temp1,
2694     &q_temp2,
2695     &(result->hmax_in));
2696     GMP_RATS_mpq_swap_components(&q_temp3);
2697     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2698     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2699     }
2700     }
2701    
2702     done_with_right_neighbors: ;
2703    
2704     //This is a single return point so we catch all the dynamic memory
2705     //deallocation.
2706     return_point:
2707     GMP_RATS_mpq_clear(&q_temp1);
2708     GMP_RATS_mpq_clear(&q_temp2);
2709     GMP_RATS_mpq_clear(&q_temp3);
2710     GMP_RATS_mpq_clear(&q_temp4);
2711     GMP_RATS_mpq_clear(&left_neighbor);
2712     GMP_RATS_mpq_clear(&right_neighbor);
2713     GMP_RATS_mpq_clear(&left_neighbor_abs);
2714     GMP_RATS_mpq_clear(&right_neighbor_abs);
2715     GMP_RATS_mpq_clear(&hmax_over_one_neg);
2716     GMP_RATS_mpq_clear(&corner_point_neg);
2717     GMP_RATS_mpq_clear(&abs_norm_recip_rn);
2718     }
2719    
2720    
2721     //08/16/01: Visual inspection OK.
2722     void GMP_RALG_consecutive_fab_terms_result_free(
2723     GMP_RALG_fab_neighbor_collection_struct *arg
2724     )
2725     {
2726     int i;
2727    
2728     //Eyeball the input.
2729     assert(arg != NULL);
2730    
2731     //Deallocate all rational numbers and integers that MUST be allocated, i.e. they are
2732     //never conditional.
2733     GMP_RATS_mpq_clear(&(arg->rn_in));
2734     GMP_INTS_mpz_clear(&(arg->kmax_in));
2735     GMP_INTS_mpz_clear(&(arg->hmax_in));
2736     GMP_RATS_mpq_clear(&(arg->hmax_over_one));
2737     GMP_RATS_mpq_clear(&(arg->corner_point));
2738     GMP_RATS_mpq_clear(&(arg->corner_point_minus_one));
2739     GMP_RATS_mpq_clear(&(arg->corner_point_plus_one));
2740     GMP_RATS_mpq_clear(&(arg->norm_rn));
2741     GMP_RATS_mpq_clear(&(arg->abs_norm_rn));
2742    
2743     //Destroy any continued fraction representations that were
2744     //formed.
2745     if (arg->rn_abs_cf_rep_formed)
2746     {
2747     GMP_RALG_cfdecomp_destroy(&(arg->rn_abs_cf_rep));
2748     }
2749     if (arg->rn_abs_recip_cf_rep_formed)
2750     {
2751     GMP_RALG_cfdecomp_destroy(&(arg->rn_abs_recip_cf_rep));
2752     }
2753     if(arg->corner_point_cf_rep_formed)
2754     {
2755     GMP_RALG_cfdecomp_destroy(&(arg->corner_point_cf_rep));
2756     }
2757     if(arg->corner_point_recip_cf_rep_formed)
2758     {
2759     GMP_RALG_cfdecomp_destroy(&(arg->corner_point_recip_cf_rep));
2760     }
2761    
2762     //Walk through the lefts, freeing up any allocated rational numbers.
2763     for (i=0; i < arg->n_left_out; i++)
2764     {
2765     GMP_RATS_mpq_clear(&(arg->lefts[i].neighbor));
2766     }
2767    
2768     //Walk through the rights, freeing up any allocated rational numbers.
2769     for (i=0; i < arg->n_right_out; i++)
2770     {
2771     GMP_RATS_mpq_clear(&(arg->rights[i].neighbor));
2772     }
2773    
2774     //Deallocate any area assigned for lefts.
2775     if (arg->lefts)
2776     {
2777     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
2778     CCMALLOC_free(arg->lefts);
2779     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
2780     TclpFree((char *)arg->lefts);
2781     #else
2782     free(arg->lefts);
2783     #endif
2784    
2785     arg->lefts = NULL;
2786     }
2787    
2788     //Deallocate any area assigned for rights.
2789     if (arg->rights)
2790     {
2791     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
2792     CCMALLOC_free(arg->rights);
2793     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
2794     TclpFree((char *)arg->rights);
2795     #else
2796     free(arg->rights);
2797     #endif
2798    
2799     arg->rights = NULL;
2800     }
2801     }
2802    
2803    
2804     //08/16/01: Visual inspection OK.
2805     void GMP_RALG_consecutive_fab_terms_result_dump(
2806     FILE *s,
2807     GMP_RALG_fab_neighbor_collection_struct *arg
2808     )
2809     {
2810     int i;
2811     char buf[250];
2812    
2813     //Eyeball the input parameters.
2814     assert(s != NULL);
2815     assert(arg != NULL);
2816    
2817     //Announce intent.
2818     FCMIOF_stream_bannerheading(s,
2819     "Emitting Neighbor Decomp For Analysis",
2820     1);
2821    
2822     //Dump the fields, in order.
2823     GMP_INTS_mpz_long_int_format_to_stream(s,
2824     &(arg->rn_in.num),
2825     "rn_in_num");
2826     GMP_INTS_mpz_long_int_format_to_stream(s,
2827     &(arg->rn_in.den),
2828     "rn_in_den");
2829     GMP_INTS_mpz_long_int_format_to_stream(s,
2830     &(arg->kmax_in),
2831     "kmax_in");
2832     fprintf(s, " hmax_supplied: %12d\n", arg->hmax_supplied);
2833     GMP_INTS_mpz_long_int_format_to_stream(s,
2834     &(arg->hmax_in),
2835     "hmax_in");
2836     if (arg->error)
2837     {
2838     fprintf(s, " error: %s\n", arg->error);
2839     }
2840     else
2841     {
2842     fprintf(s, " error: NULL\n");
2843     }
2844    
2845     GMP_INTS_mpz_long_int_format_to_stream(s,
2846     &(arg->corner_point.num),
2847     "corner_point_num");
2848     GMP_INTS_mpz_long_int_format_to_stream(s,
2849     &(arg->corner_point.den),
2850     "corner_point_den");
2851    
2852     GMP_INTS_mpz_long_int_format_to_stream(s,
2853     &(arg->corner_point_minus_one.num),
2854     "cor_pt_minus_one_num");
2855     GMP_INTS_mpz_long_int_format_to_stream(s,
2856     &(arg->corner_point_minus_one.den),
2857     "cor_pt_minus_one_den");
2858    
2859     GMP_INTS_mpz_long_int_format_to_stream(s,
2860     &(arg->corner_point_plus_one.num),
2861     "cor_pt_plus_one_num");
2862     GMP_INTS_mpz_long_int_format_to_stream(s,
2863     &(arg->corner_point_plus_one.den),
2864     "cor_pt_plus_one_den");
2865    
2866     GMP_INTS_mpz_long_int_format_to_stream(s,
2867     &(arg->hmax_over_one.num),
2868     "hmax/1_num");
2869     GMP_INTS_mpz_long_int_format_to_stream(s,
2870     &(arg->hmax_over_one.den),
2871     "hmax/1_den");
2872    
2873     GMP_INTS_mpz_long_int_format_to_stream(s,
2874     &(arg->norm_rn.num),
2875     "norm_rn_num");
2876     GMP_INTS_mpz_long_int_format_to_stream(s,
2877     &(arg->norm_rn.den),
2878     "norm_rn_den");
2879    
2880     GMP_INTS_mpz_long_int_format_to_stream(s,
2881     &(arg->abs_norm_rn.num),
2882     "abs_norm_rn_num");
2883     GMP_INTS_mpz_long_int_format_to_stream(s,
2884     &(arg->abs_norm_rn.den),
2885     "abs_norm_rn_den");
2886    
2887     fprintf(s, " rn_is_neg: %12d\n", arg->rn_is_neg);
2888    
2889     fprintf(s, " n_left_in: %12d\n", arg->n_left_in);
2890    
2891     fprintf(s, " n_right_in: %12d\n", arg->n_right_in);
2892    
2893     fprintf(s, "rn_abs_cf_rep_formed: %12d\n", arg->rn_abs_cf_rep_formed);
2894    
2895     if (arg->rn_abs_cf_rep_formed)
2896     {
2897     GMP_RALG_cfdecomp_emit(s, "Abs Of RN In CF Decomp", &(arg->rn_abs_cf_rep), 0, 0, NULL);
2898     }
2899    
2900     fprintf(s, "rn_abs_recip_cf_rep_formed: %12d\n", arg->rn_abs_recip_cf_rep_formed);
2901    
2902     if (arg->rn_abs_recip_cf_rep_formed)
2903     {
2904     GMP_RALG_cfdecomp_emit(s, "Abs Of Recip Of RN In CF Decomp", &(arg->rn_abs_recip_cf_rep), 0, 0, NULL);
2905     }
2906    
2907     fprintf(s, "corner_point_cf_rep_formed: %12d\n", arg->corner_point_cf_rep_formed);
2908    
2909     if (arg->corner_point_cf_rep_formed)
2910     {
2911     GMP_RALG_cfdecomp_emit(s, "Corner Point CF Decomp", &(arg->corner_point_cf_rep), 0, 0, NULL);
2912     }
2913    
2914     fprintf(s, "cor_pt_recip_cf_rep_formed: %12d\n", arg->corner_point_recip_cf_rep_formed);
2915    
2916     if (arg->corner_point_recip_cf_rep_formed)
2917     {
2918     GMP_RALG_cfdecomp_emit(s, "Corner Point Recip CF Decomp", &(arg->corner_point_recip_cf_rep), 0, 0, NULL);
2919     }
2920    
2921     fprintf(s, " equality: %12d\n", arg->equality);
2922    
2923     fprintf(s, " n_left_out: %12d\n", arg->n_left_out);
2924    
2925     fprintf(s, " n_right_out: %12d\n", arg->n_right_out);
2926    
2927     for (i=0; i < arg->n_left_out; i++)
2928     {
2929     sprintf(buf, "Contents Of Left Neighbor Array [%d]", i);
2930     FCMIOF_stream_bannerheading(s,
2931     buf,
2932     0);
2933     fprintf(s, " index: %12d\n", arg->lefts[i].index);
2934     GMP_INTS_mpz_long_int_format_to_stream(s,
2935     &(arg->lefts[i].neighbor.num),
2936     "neighbor_num");
2937     GMP_INTS_mpz_long_int_format_to_stream(s,
2938     &(arg->lefts[i].neighbor.den),
2939     "neighbor_den");
2940     }
2941    
2942     for (i=0; i < arg->n_right_out; i++)
2943     {
2944     sprintf(buf, "Contents Of Right Neighbor Array [%d]", i);
2945     FCMIOF_stream_bannerheading(s,
2946     buf,
2947     0);
2948     fprintf(s, " index: %12d\n", arg->rights[i].index);
2949     GMP_INTS_mpz_long_int_format_to_stream(s,
2950     &(arg->rights[i].neighbor.num),
2951     "neighbor_num");
2952     GMP_INTS_mpz_long_int_format_to_stream(s,
2953     &(arg->rights[i].neighbor.den),
2954     "neighbor_den");
2955     }
2956    
2957     FCMIOF_stream_hline(s);
2958     }
2959    
2960    
2961     /******************************************************************/
2962     /*** VERSION CONTROL REPORTING FUNCTIONS ************************/
2963     /******************************************************************/
2964    
2965     //08/16/01: Visual inspection OK.
2966     const char *GMP_RALG_cvcinfo(void)
2967     {
2968     return("$Header: /cvsroot/esrg/sfesrg/esrgpcpj/shared/c_datd/gmp_ralg.c,v 1.10 2002/01/27 17:58:15 dtashley Exp $");
2969     }
2970    
2971    
2972     //08/16/01: Visual inspection OK.
2973     const char *GMP_RALG_hvcinfo(void)
2974     {
2975     return(GMP_RALG_H_VERSION);
2976     }
2977    
2978    
2979     //**************************************************************************
2980     // $Log: gmp_ralg.c,v $
2981     // Revision 1.10 2002/01/27 17:58:15 dtashley
2982     // CRC32, other programs modified to work under new directory structure.
2983     //
2984     // Revision 1.9 2001/08/18 18:33:13 dtashley
2985     // Preparing for release of v1.05.
2986     //
2987     // Revision 1.8 2001/08/16 19:49:40 dtashley
2988     // Beginning to prepare for v1.05 release.
2989     //
2990     // Revision 1.7 2001/08/15 06:56:05 dtashley
2991     // Substantial progress. Safety check-in.
2992     //
2993     // Revision 1.6 2001/08/12 10:20:58 dtashley
2994     // Safety check-in. Substantial progress.
2995     //
2996     // Revision 1.5 2001/08/07 10:42:48 dtashley
2997     // Completion of CFRATNUM extensions and DOS command-line utility.
2998     //
2999     // Revision 1.4 2001/07/13 21:02:20 dtashley
3000     // Version control reporting changes.
3001     //
3002     // Revision 1.3 2001/07/13 20:44:42 dtashley
3003     // Changes, CVS keyword expansion test.
3004     //
3005     // Revision 1.2 2001/07/13 00:57:08 dtashley
3006     // Safety check-in. Substantial progress on port.
3007     //
3008     // Revision 1.1 2001/07/12 05:42:06 dtashley
3009     // Initial checkin.
3010     //
3011     //**************************************************************************
3012     // End of GMP_RALG.C.

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25