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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 29 - (show annotations) (download)
Sat Oct 8 07:08:47 2016 UTC (6 years, 5 months ago) by dashley
Original Path: to_be_filed/sf_code/esrgpcpj/shared/c_datd/gmp_ralg.c
File MIME type: text/plain
File size: 132737 byte(s)
Directories relocated.
1 // $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