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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 56 - (hide annotations) (download)
Sat Oct 29 01:53:01 2016 UTC (7 years, 8 months ago) by dashley
File MIME type: text/plain
File size: 95292 byte(s)
License and property (keyword) changes.
1 dashley 56 //$Header$
2 dashley 25 //-------------------------------------------------------------------------------------------------
3 dashley 56 //This file is part of "David T. Ashley's Shared Source Code", a set of shared components
4     //integrated into many of David T. Ashley's projects.
5 dashley 25 //-------------------------------------------------------------------------------------------------
6 dashley 56 //This source code and any program in which it is compiled/used is provided under the MIT License,
7     //reproduced below.
8     //-------------------------------------------------------------------------------------------------
9     //Permission is hereby granted, free of charge, to any person obtaining a copy of
10     //this software and associated documentation files(the "Software"), to deal in the
11     //Software without restriction, including without limitation the rights to use,
12     //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
13     //Software, and to permit persons to whom the Software is furnished to do so,
14     //subject to the following conditions :
15 dashley 25 //
16 dashley 56 //The above copyright notice and this permission notice shall be included in all
17     //copies or substantial portions of the Software.
18 dashley 25 //
19 dashley 56 //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20     //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21     //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
22     //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23     //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24     //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25     //SOFTWARE.
26 dashley 25 //-------------------------------------------------------------------------------------------------
27     #define MODULE_GMP_RALG
28    
29     #include <assert.h>
30     #include <stdio.h>
31     #include <string.h>
32     #include <time.h>
33    
34    
35     #include "fcmiof.h"
36     #include "gmp_ints.h"
37     #include "gmp_rats.h"
38     #include "gmp_ralg.h"
39     #include "intfunc.h"
40    
41     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
42     #include "ccmalloc.h"
43     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
44     #include "tclalloc.h"
45     #else
46     /* Do nothing. */
47     #endif
48    
49    
50     /******************************************************************/
51     /*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/
52     /******************************************************************/
53     //08/16/01: Visual inspection OK.
54     void GMP_RALG_cfdecomp_init(
55     GMP_RALG_cf_app_struct *decomp,
56     int *failure,
57     GMP_INTS_mpz_struct *num,
58     GMP_INTS_mpz_struct *den)
59     {
60     int loop_counter, i;
61     GMP_INTS_mpz_struct arb_temp1, arb_temp2;
62    
63     //Eyeball the input parameters. The rest of the eyeballing
64     //will occur as functions are called to manipulate the
65     //numerator and denominator passed in.
66     assert(decomp != NULL);
67     assert(failure != NULL);
68     assert(num != NULL);
69     assert(den != NULL);
70    
71     //Allocate space for temporary integers.
72     GMP_INTS_mpz_init(&arb_temp1);
73     GMP_INTS_mpz_init(&arb_temp2);
74    
75     //Begin believing no failure.
76     *failure = 0;
77    
78     //Initialize the copy of the numerator and denominator and
79     //copy these into the structure.
80     GMP_INTS_mpz_init(&(decomp->num));
81     GMP_INTS_mpz_copy(&(decomp->num), num);
82     GMP_INTS_mpz_init(&(decomp->den));
83     GMP_INTS_mpz_copy(&(decomp->den), den);
84    
85     //Allocate the space for the first increment of the
86     //growable areas. We need to use different memory
87     //allocation functions depending on whether we're
88     //in a Tcl build or a DOS command-line utility
89     //build.
90     //Space for partial quotients.
91     decomp->a =
92     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
93     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
94     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
95     (GMP_INTS_mpz_struct *)
96     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
97     #else
98     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
99     #endif
100    
101     //Dividends.
102     decomp->dividend =
103     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
104     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
105     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
106     (GMP_INTS_mpz_struct *)
107     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
108     #else
109     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
110     #endif
111    
112     //Divisors.
113     decomp->divisor =
114     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
115     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
116     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
117     (GMP_INTS_mpz_struct *)
118     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
119     #else
120     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
121     #endif
122    
123     //Remainders.
124     decomp->remainder =
125     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
126     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
127     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
128     (GMP_INTS_mpz_struct *)
129     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
130     #else
131     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
132     #endif
133    
134     //Convergent numerators.
135     decomp->p =
136     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
137     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
138     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
139     (GMP_INTS_mpz_struct *)
140     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
141     #else
142     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
143     #endif
144    
145     //Convergent denominators.
146     decomp->q =
147     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
148     CCMALLOC_malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
149     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
150     (GMP_INTS_mpz_struct *)
151     TclpAlloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
152     #else
153     malloc(sizeof(GMP_INTS_mpz_struct) * GMP_RALG_CF_ALLOC_INCREMENT);
154     #endif
155    
156     //Now the number of allocated slots is what we just allocated.
157     decomp->nallocd = GMP_RALG_CF_ALLOC_INCREMENT;
158    
159     //The number of slots actually used is zero, to start with.
160     decomp->n = 0;
161    
162     //There are a number of conditions that will lead to an error
163     //where we can't successfully form the continued fraction
164     //decomposition. These errors are:
165     // a)Either component is NAN.
166     // b)Zero denominator.
167     // c)Either component negative.
168     //In these cases, we'll pretend we got 0/1 for the number
169     //and set things accordingly, and we'll set the failure
170     //flag for the caller.
171     //
172     if (GMP_INTS_mpz_get_flags(num)
173     ||
174     GMP_INTS_mpz_get_flags(den)
175     ||
176     GMP_INTS_mpz_is_zero(den)
177     ||
178     GMP_INTS_mpz_is_neg(num)
179     ||
180     GMP_INTS_mpz_is_neg(den))
181     {
182     *failure = 1;
183     decomp->n = 1;
184    
185     GMP_INTS_mpz_set_ui(&(decomp->num), 0);
186     GMP_INTS_mpz_set_ui(&(decomp->den), 1);
187    
188     GMP_INTS_mpz_init(decomp->dividend);
189     GMP_INTS_mpz_set_ui(decomp->dividend, 0);
190    
191     GMP_INTS_mpz_init(decomp->divisor);
192     GMP_INTS_mpz_set_ui(decomp->divisor, 1);
193    
194     GMP_INTS_mpz_init(decomp->a);
195     GMP_INTS_mpz_set_ui(decomp->a, 0);
196    
197     GMP_INTS_mpz_init(decomp->remainder);
198     GMP_INTS_mpz_set_ui(decomp->remainder, 0);
199    
200     GMP_INTS_mpz_init(decomp->p);
201     GMP_INTS_mpz_set_ui(decomp->p, 0);
202    
203     GMP_INTS_mpz_init(decomp->q);
204     GMP_INTS_mpz_set_ui(decomp->q, 1);
205    
206     goto return_point;
207     }
208    
209     //If we're here there are not any errors that we
210     //are willing to detect. We should be clear
211     //for a continued fraction decomposition.
212     loop_counter = 0;
213     do
214     {
215     //Increment the count of "rows", because we're
216     //about to add one.
217     decomp->n++;
218    
219     //If we have used up all the space available
220     //for integers, we have to allocate more.
221     if (decomp->n > decomp->nallocd)
222     {
223     //We now have more allocated space.
224     decomp->nallocd += GMP_RALG_CF_ALLOC_INCREMENT;
225    
226     //Be absolutely sure we have not made a greivous
227     //error.
228     assert(decomp->n <= decomp->nallocd);
229    
230     //Space for dividends.
231     decomp->dividend =
232     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
233     CCMALLOC_realloc(
234     decomp->dividend,
235     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
236     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
237     (GMP_INTS_mpz_struct *)
238     TclpRealloc((char *)decomp->dividend,
239     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
240     #else
241     realloc(decomp->dividend, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
242     #endif
243    
244     //Space for divisors.
245     decomp->divisor =
246     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
247     CCMALLOC_realloc(
248     decomp->divisor,
249     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
250     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
251     (GMP_INTS_mpz_struct *)
252     TclpRealloc((char *)decomp->divisor,
253     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
254     #else
255     realloc(decomp->divisor, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
256     #endif
257    
258     //Space for partial quotients.
259     decomp->a =
260     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
261     CCMALLOC_realloc(
262     decomp->a,
263     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
264     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
265     (GMP_INTS_mpz_struct *)
266     TclpRealloc((char *)decomp->a,
267     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
268     #else
269     realloc(decomp->a, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
270     #endif
271    
272     //Space for remainders.
273     decomp->remainder =
274     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
275     CCMALLOC_realloc(
276     decomp->remainder,
277     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
278     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
279     (GMP_INTS_mpz_struct *)
280     TclpRealloc((char *)decomp->remainder,
281     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
282     #else
283     realloc(decomp->remainder, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
284     #endif
285    
286     //Space for partial quotient numerators.
287     decomp->p =
288     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
289     CCMALLOC_realloc(
290     decomp->p,
291     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
292     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
293     (GMP_INTS_mpz_struct *)
294     TclpRealloc((char *)decomp->p,
295     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
296     #else
297     realloc(decomp->p, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
298     #endif
299    
300     //Space for partial quotient denominators.
301     decomp->q =
302     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
303     CCMALLOC_realloc(
304     decomp->q,
305     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
306     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
307     (GMP_INTS_mpz_struct *)
308     TclpRealloc((char *)decomp->q,
309     sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
310     #else
311     realloc(decomp->q, sizeof(GMP_INTS_mpz_struct) * decomp->nallocd);
312     #endif
313     }
314    
315     //At this point, we have enough space to do the next round of operations.
316     //Set up an index variable.
317     i = decomp->n - 1;
318    
319     //Initialize all of the integers at this round.
320     GMP_INTS_mpz_init(decomp->dividend + i);
321     GMP_INTS_mpz_init(decomp->divisor + i);
322     GMP_INTS_mpz_init(decomp->a + i);
323     GMP_INTS_mpz_init(decomp->remainder + i);
324     GMP_INTS_mpz_init(decomp->p + i);
325     GMP_INTS_mpz_init(decomp->q + i);
326    
327     //Perform the next round of continued fraction decomposition. This
328     //is standard stuff.
329     if (i==0)
330     {
331     //In the 0th round, we essentially perform initial assignments.
332     GMP_INTS_mpz_copy(decomp->dividend, &(decomp->num));
333     GMP_INTS_mpz_copy(decomp->divisor, &(decomp->den));
334    
335     //Make the division to get quotient and remainder.
336     GMP_INTS_mpz_tdiv_qr(decomp->a, decomp->remainder, decomp->dividend, decomp->divisor);
337    
338     //The convergents in the first round are always the quotient over 1.
339     GMP_INTS_mpz_copy(decomp->p, decomp->a);
340     GMP_INTS_mpz_set_ui(decomp->q, 1);
341     }
342     else if (i==1)
343     {
344     //In the 1st round, the only point of caution is that we have to
345     //consider p(k-2)/q(k-2) carefully.
346     GMP_INTS_mpz_copy(decomp->dividend + 1, decomp->divisor + 0);
347     GMP_INTS_mpz_copy(decomp->divisor + 1, decomp->remainder + 0);
348    
349     //Make the division to get quotient and remainder.
350     GMP_INTS_mpz_tdiv_qr(decomp->a + 1,
351     decomp->remainder + 1,
352     decomp->dividend + 1,
353     decomp->divisor + 1);
354    
355     //Need to compute the numerator of the convergent. This will be:
356     // a(1) p(0) + p(-1) = a(1)p(0) + 1.
357     GMP_INTS_mpz_mul(decomp->p + 1, decomp->a + 1, decomp->p + 0);
358     GMP_INTS_mpz_set_ui(&arb_temp1, 1);
359     GMP_INTS_mpz_add(decomp->p + 1, decomp->p + 1, &arb_temp1);
360    
361     //Need to compute the denominator of the convergent. This will
362     //be a(1)q(0) + q(-1) = a(1) q(0) = a(1).
363     GMP_INTS_mpz_copy(decomp->q + 1, decomp->a + 1);
364     }
365     else
366     {
367     //In the general case, it is a simple formula.
368     //Rotate in the dividend and divisor from the previous round.
369     GMP_INTS_mpz_copy(decomp->dividend + i, decomp->divisor + i - 1);
370     GMP_INTS_mpz_copy(decomp->divisor + i, decomp->remainder + i - 1);
371    
372     //Make the division to get quotient and remainder.
373     GMP_INTS_mpz_tdiv_qr(decomp->a + i,
374     decomp->remainder + i,
375     decomp->dividend + i,
376     decomp->divisor + i);
377    
378     //Need to compute the numerator of the convergent. This will be:
379     // p(i) = a(i)p(i-1) + p(i-2)
380     GMP_INTS_mpz_mul(decomp->p + i, decomp->a + i, decomp->p + i - 1);
381     GMP_INTS_mpz_add(decomp->p + i, decomp->p + i, decomp->p + i - 2);
382    
383     //Need to compute the numerator of the convergent. This will be:
384     // q(i) = q(i)q(i-1) + q(i-2)
385     GMP_INTS_mpz_mul(decomp->q + i, decomp->a + i, decomp->q + i - 1);
386     GMP_INTS_mpz_add(decomp->q + i, decomp->q + i, decomp->q + i - 2);
387     }
388    
389     loop_counter++;
390     } while(!GMP_INTS_mpz_is_zero(decomp->remainder + decomp->n - 1) && loop_counter < 100000);
391    
392     //In debug builds, be sure we did not terminate based on the loop counter.
393     assert(loop_counter != 100000);
394    
395     return_point:
396    
397     //Deallocate space for temporary integers.
398     GMP_INTS_mpz_clear(&arb_temp1);
399     GMP_INTS_mpz_clear(&arb_temp2);
400     }
401    
402    
403     //08/16/01: Visual inspection OK.
404     void GMP_RALG_cfdecomp_destroy(
405     GMP_RALG_cf_app_struct *decomp
406     )
407     {
408     int i;
409    
410     //Eyeball the input parameters. Other eyeballing
411     //will be done as integers are destroyed.
412     assert(decomp != NULL);
413    
414     //First, destroy the things that are bound directly
415     //to the record.
416     GMP_INTS_mpz_clear(&(decomp->num));
417     GMP_INTS_mpz_clear(&(decomp->den));
418    
419     //Now, destroy every integer which is allocated.
420     for (i=0; i < decomp->n; i++)
421     {
422     GMP_INTS_mpz_clear(decomp->dividend + i);
423     GMP_INTS_mpz_clear(decomp->divisor + i);
424     GMP_INTS_mpz_clear(decomp->a + i);
425     GMP_INTS_mpz_clear(decomp->remainder + i);
426     GMP_INTS_mpz_clear(decomp->p + i);
427     GMP_INTS_mpz_clear(decomp->q + i);
428     }
429    
430     //Now, destroy the arrays of integers.
431     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
432     CCMALLOC_free(decomp->dividend);
433     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
434     TclpFree((char *)decomp->dividend);
435     #else
436     free(decomp->dividend);
437     #endif
438     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
439     CCMALLOC_free(decomp->divisor);
440     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
441     TclpFree((char *)decomp->divisor);
442     #else
443     free(decomp->divisor);
444     #endif
445     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
446     CCMALLOC_free(decomp->a);
447     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
448     TclpFree((char *)decomp->a);
449     #else
450     free(decomp->a);
451     #endif
452     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
453     CCMALLOC_free(decomp->remainder);
454     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
455     TclpFree((char *)decomp->remainder);
456     #else
457     free(decomp->remainder);
458     #endif
459     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
460     CCMALLOC_free(decomp->p);
461     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
462     TclpFree((char *)decomp->p);
463     #else
464     free(decomp->p);
465     #endif
466     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
467     CCMALLOC_free(decomp->q);
468     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
469     TclpFree((char *)decomp->q);
470     #else
471     free(decomp->q);
472     #endif
473     }
474    
475    
476     /******************************************************************/
477     /*** FORMATTED OUTPUT FUNCTIONS *********************************/
478     /******************************************************************/
479     //08/16/01: Visual inspection OK.
480     void GMP_RALG_cfdecomp_emit(
481     FILE *s,
482     char *banner,
483     GMP_RALG_cf_app_struct *decomp,
484     int nf,
485     int dap,
486     const GMP_INTS_mpz_struct *dap_denominator)
487     {
488     int i;
489    
490     GMP_INTS_mpz_struct arb_temp, arb_quotient, arb_remainder;
491    
492     //Eyeball the input parameters. The banner is allowed to
493     //be null, so can't check that.
494     assert(s != NULL);
495     assert(decomp != NULL);
496    
497     //Allocate our temporary integers.
498     GMP_INTS_mpz_init(&arb_temp);
499     GMP_INTS_mpz_init(&arb_quotient);
500     GMP_INTS_mpz_init(&arb_remainder);
501    
502     //If banner requested and noformat option not used.
503     if (banner && !nf)
504     {
505     FCMIOF_stream_bannerheading(s, banner, 1);
506     }
507    
508     //Dump the input numerator.
509     if (!nf)
510     {
511     GMP_INTS_mpz_long_int_format_to_stream(s,
512     &(decomp->num),
513     "Input Numerator");
514     }
515     else
516     {
517     GMP_INTS_mpz_arb_int_raw_to_stream(s, &(decomp->num));
518     fprintf(s, "\n");
519     }
520    
521     //Separator if not in unformatted mode.
522     if (!nf)
523     FCMIOF_stream_hline(s);
524    
525     //Dump the input denominator.
526     if (!nf)
527     {
528     GMP_INTS_mpz_long_int_format_to_stream(s,
529     &(decomp->den),
530     "Input Denominator");
531     }
532     else
533     {
534     GMP_INTS_mpz_arb_int_raw_to_stream(s, &(decomp->den));
535     fprintf(s, "\n");
536     }
537    
538     //Separator if not in unformatted mode.
539     if (!nf)
540     FCMIOF_stream_hline(s);
541    
542     for (i=0; i<decomp->n; i++)
543     {
544     char strbuf[100];
545     //Buffer to prepare description.
546    
547     //Print out the dividend at each round.
548     if (!nf)
549     {
550     sprintf(strbuf, "dividend(%d)", i);
551     GMP_INTS_mpz_long_int_format_to_stream(s,
552     decomp->dividend + i,
553     strbuf);
554     }
555     else
556     {
557     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->dividend+i);
558     fprintf(s, "\n");
559     }
560    
561     //Separator if not in unformatted mode.
562     if (!nf)
563     FCMIOF_stream_hline(s);
564    
565     //Print out the divisor at each round.
566     if (!nf)
567     {
568     sprintf(strbuf, "divisor(%d)", i);
569     GMP_INTS_mpz_long_int_format_to_stream(s,
570     decomp->divisor + i,
571     strbuf);
572     }
573     else
574     {
575     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->divisor+i);
576     fprintf(s, "\n");
577     }
578    
579     //Separator if not in unformatted mode.
580     if (!nf)
581     FCMIOF_stream_hline(s);
582    
583     //Print out partial quotient at each round.
584     if (!nf)
585     {
586     sprintf(strbuf, "a(%d)", i);
587     GMP_INTS_mpz_long_int_format_to_stream(s,
588     decomp->a + i,
589     strbuf);
590     }
591     else
592     {
593     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->a+i);
594     fprintf(s, "\n");
595     }
596    
597     //Separator if not in unformatted mode.
598     if (!nf)
599     FCMIOF_stream_hline(s);
600    
601     //It doesn't make any sense to print out the
602     //remainder, because this becomes the divisor
603     //for the next round. It is just wasted output
604     //lines.
605    
606     //Print out the convergent numerator at
607     //each round.
608     if (!nf)
609     {
610     sprintf(strbuf, "p(%d)", i);
611     GMP_INTS_mpz_long_int_format_to_stream(s,
612     decomp->p + i,
613     strbuf);
614     }
615     else
616     {
617     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->p+i);
618     fprintf(s, "\n");
619     }
620    
621     //Separator if not in unformatted mode.
622     if (!nf)
623     FCMIOF_stream_hline(s);
624    
625     //Print out the convergent denominator at
626     //each round.
627     if (!nf)
628     {
629     sprintf(strbuf, "q(%d)", i);
630     GMP_INTS_mpz_long_int_format_to_stream(s,
631     decomp->q + i,
632     strbuf);
633     }
634     else
635     {
636     GMP_INTS_mpz_arb_int_raw_to_stream(s, decomp->q+i);
637     fprintf(s, "\n");
638     }
639    
640     //Separator if not in unformatted mode.
641     if (!nf)
642     FCMIOF_stream_hline(s);
643    
644     if (dap)
645     {
646     //Calculate the DAP numerator
647     GMP_INTS_mpz_mul(&arb_temp, dap_denominator, decomp->p + i);
648     GMP_INTS_mpz_tdiv_qr(&arb_quotient, &arb_remainder,
649     &arb_temp, decomp->q + i);
650    
651     //Print DAP numerator.
652     if (!nf)
653     {
654     sprintf(strbuf, "dap_h(%d)", i);
655     GMP_INTS_mpz_long_int_format_to_stream(s,
656     &arb_quotient,
657     strbuf);
658     }
659     else
660     {
661     GMP_INTS_mpz_arb_int_raw_to_stream(s, &arb_quotient);
662     fprintf(s, "\n");
663     }
664    
665     //Separator if not in unformatted mode.
666     if (!nf)
667     FCMIOF_stream_hline(s);
668    
669     //Print DAP denominator.
670     if (!nf)
671     {
672     sprintf(strbuf, "dap_k(%d)", i);
673     GMP_INTS_mpz_long_int_format_to_stream(s,
674     dap_denominator,
675     strbuf);
676     }
677     else
678     {
679     GMP_INTS_mpz_arb_int_raw_to_stream(s, dap_denominator);
680     fprintf(s, "\n");
681     }
682    
683     //Separator if not in unformatted mode.
684     if (!nf)
685     FCMIOF_stream_hline(s);
686     }
687     }
688    
689     //Deallocate our temporary integers.
690     GMP_INTS_mpz_clear(&arb_temp);
691     GMP_INTS_mpz_clear(&arb_quotient);
692     GMP_INTS_mpz_clear(&arb_remainder);
693     }
694    
695    
696     /******************************************************************/
697     /*** FAREY SERIES PREDECESSOR AND SUCCESSOR FUNCTIONS ***********/
698     /******************************************************************/
699     //08/16/01: Visual inspection OK.
700     void GMP_RALG_farey_predecessor(
701     GMP_RATS_mpq_struct *result,
702     const GMP_RATS_mpq_struct *plus_two,
703     const GMP_RATS_mpq_struct *plus_one,
704     const GMP_INTS_mpz_struct *N)
705     {
706     GMP_RATS_mpq_struct result_copy;
707     //Used to hold return value in case the result
708     //is the same as either of the input arguments.
709     GMP_INTS_mpz_struct temp1, temp2, floor_func;
710     //Temporary integers.
711    
712     assert(result != NULL);
713     assert(plus_two != NULL);
714     assert(plus_one != NULL);
715     assert(N != NULL);
716    
717     //Initialize the variables used.
718     GMP_RATS_mpq_init(&result_copy);
719     GMP_INTS_mpz_init(&temp1);
720     GMP_INTS_mpz_init(&temp2);
721     GMP_INTS_mpz_init(&floor_func);
722    
723     //Numerator of the term in the floor function.
724     GMP_INTS_mpz_add(&temp1, &(plus_two->den), N);
725    
726     //Term in the floor function. This is used to
727     //calculate both numerator and denominator, so we save it.
728     GMP_INTS_mpz_tdiv_qr(&floor_func, &temp2, &temp1, &(plus_one->den));
729    
730     //Product of result of floor function and numerator--now
731     //forming the numerator of the output.
732     GMP_INTS_mpz_mul(&temp2, &floor_func, &(plus_one->num));
733    
734     //Final result assigned to numerator.
735     GMP_INTS_mpz_sub(&(result_copy.num), &temp2, &(plus_two->num));
736    
737     //Product of result of floor function and denominator--now
738     //forming the denominator of the output.
739     GMP_INTS_mpz_mul(&temp2, &floor_func, &(plus_one->den));
740    
741     //Final result assigned to denominator.
742     GMP_INTS_mpz_sub(&(result_copy.den), &temp2, &(plus_two->den));
743    
744     //Copy the result to the object owned by the caller.
745     GMP_RATS_mpq_copy(result, &result_copy);
746    
747     //Deallocate dynamic memory.
748     GMP_RATS_mpq_clear(&result_copy);
749     GMP_INTS_mpz_clear(&temp1);
750     GMP_INTS_mpz_clear(&temp2);
751     GMP_INTS_mpz_clear(&floor_func);
752     }
753    
754    
755     //08/16/01: Visual inspection OK.
756     void GMP_RALG_farey_successor(
757     GMP_RATS_mpq_struct *result,
758     const GMP_RATS_mpq_struct *minus_two,
759     const GMP_RATS_mpq_struct *minus_one,
760     const GMP_INTS_mpz_struct *N)
761     {
762     GMP_RATS_mpq_struct result_copy;
763     //Used to hold return value in case the result
764     //is the same as either of the input arguments.
765     GMP_INTS_mpz_struct temp1, temp2, floor_func;
766     //Temporary integers.
767    
768     assert(result != NULL);
769     assert(minus_two != NULL);
770     assert(minus_one != NULL);
771     assert(N != NULL);
772    
773     //Initialize the variables used.
774     GMP_RATS_mpq_init(&result_copy);
775     GMP_INTS_mpz_init(&temp1);
776     GMP_INTS_mpz_init(&temp2);
777     GMP_INTS_mpz_init(&floor_func);
778    
779     //Numerator of the term in the floor function.
780     GMP_INTS_mpz_add(&temp1, &(minus_two->den), N);
781    
782     //Term in the floor function. This is used to
783     //calculate both numerator and denominator, so we save it.
784     GMP_INTS_mpz_tdiv_qr(&floor_func, &temp2, &temp1, &(minus_one->den));
785    
786     //Product of result of floor function and numerator--now
787     //forming the numerator of the output.
788     GMP_INTS_mpz_mul(&temp2, &floor_func, &(minus_one->num));
789    
790     //Final result assigned to numerator.
791     GMP_INTS_mpz_sub(&(result_copy.num), &temp2, &(minus_two->num));
792    
793     //Product of result of floor function and denominator--now
794     //forming the denominator of the output.
795     GMP_INTS_mpz_mul(&temp2, &floor_func, &(minus_one->den));
796    
797     //Final result assigned to denominator.
798     GMP_INTS_mpz_sub(&(result_copy.den), &temp2, &(minus_two->den));
799    
800     //Copy the result to the object owned by the caller.
801     GMP_RATS_mpq_copy(result, &result_copy);
802    
803     //Deallocate dynamic memory.
804     GMP_RATS_mpq_clear(&result_copy);
805     GMP_INTS_mpz_clear(&temp1);
806     GMP_INTS_mpz_clear(&temp2);
807     GMP_INTS_mpz_clear(&floor_func);
808     }
809    
810    
811     //08/16/01: Visual inspection OK.
812     void GMP_RALG_enclosing_farey_neighbors(
813     const GMP_RATS_mpq_struct *rn_in,
814     const GMP_INTS_mpz_struct *N,
815     const GMP_RALG_cf_app_struct *cf_rep,
816     int *equality,
817     GMP_RATS_mpq_struct *left,
818     GMP_RATS_mpq_struct *right)
819     {
820     GMP_RATS_mpq_struct rn_abs;
821     //Absolute value of rational number supplied.
822     GMP_RATS_mpq_struct previous_convergent;
823     //Convergent before the one that has the denominator
824     //not exceeding the order of the series. Need to fudge
825     //a little bit because don't have -1-th order convergents
826     //tabulated.
827     GMP_RATS_mpq_struct other_neighbor;
828     //The other neighbor besides the highest-order convergent
829     //without denominator too large.
830     GMP_INTS_mpz_struct temp1, temp2, temp3, temp4;
831     //Temporary integers.
832     int ho_conv;
833     //Index of highest-ordered convergent that does not have
834     //denominator too large.
835    
836     //Eyeball the parameters.
837     assert(rn_in != NULL);
838     assert(N != NULL);
839     assert(cf_rep != NULL);
840     assert(equality != NULL);
841     assert(left != NULL);
842     assert(right != NULL);
843    
844     //Allocate dynamic variables.
845     GMP_RATS_mpq_init(&rn_abs);
846     GMP_RATS_mpq_init(&previous_convergent);
847     GMP_RATS_mpq_init(&other_neighbor);
848     GMP_INTS_mpz_init(&temp1);
849     GMP_INTS_mpz_init(&temp2);
850     GMP_INTS_mpz_init(&temp3);
851     GMP_INTS_mpz_init(&temp4);
852    
853     //Zero is a troublesome case, because it requires us to
854     //cross signs. Split this case out explicitly.
855     if (GMP_INTS_mpz_is_zero(&(rn_in->num)))
856     {
857     *equality = 1; //0/1 a member of Farey series of any order.
858     GMP_INTS_mpz_set_si(&(left->num), -1);
859     GMP_INTS_mpz_copy(&(left->den), N);
860     GMP_INTS_mpz_set_si(&(right->num), 1);
861     GMP_INTS_mpz_copy(&(right->den), N);
862     }
863     else
864     {
865     //Make a copy of the rational number in. As a condition of
866     //using this function, it must be normalized, but it still
867     //may be negative. Our strategy is to treat the number as
868     //positive, find the neighbors, then if it was negative
869     //complement and re-order the neighbors. In other words,
870     //find neighbors to a negative number by symmetry, not
871     //by forming the CF representation of a negative number.
872     //Also, we can't touch the input parameter.
873     GMP_RATS_mpq_copy(&rn_abs, rn_in);
874     GMP_INTS_mpz_abs(&(rn_abs.num));
875    
876     //Find the index of the highest-ordered convergent
877     //with a denominator not exceeding the denominator of
878     //the rational number supplied. The zero'th order
879     //convergent has a denominator of 1, so that one
880     //at least is safe.
881    
882     //Assign either the "left" or right
883     //neighbor to be the highest-ordered
884     //convergent with a denominator not exceeding the
885     //denominator of the rational number input. I say
886     //"either" because the properties of convergents let
887     //us know based on the oddness or evenness of the order
888     //which side it is on.
889     ho_conv = 0;
890     while (((ho_conv + 1) < cf_rep->n) && (GMP_INTS_mpz_cmp(cf_rep->q + ho_conv + 1, N) <= 0))
891     {
892     #if 0
893     //Some questions about this loop--debugging output.
894     printf("ho_conv : %d\n", ho_conv);
895     GMP_INTS_mpz_long_int_format_to_stream(stdout,
896     cf_rep->q + ho_conv + 1,
897     "decomp_den");
898     GMP_INTS_mpz_long_int_format_to_stream(stdout,
899     &(rn_abs.den),
900     "rn_in_den");
901     printf("Compare result: %d\n\n", GMP_INTS_mpz_cmp(cf_rep->q + ho_conv + 1, &(rn_abs.den)));
902     #endif
903    
904     ho_conv++;
905     }
906    
907     if (INTFUNC_is_even(ho_conv))
908     {
909     GMP_INTS_mpz_copy(&(left->num), cf_rep->p + ho_conv);
910     GMP_INTS_mpz_copy(&(left->den), cf_rep->q + ho_conv);
911     }
912     else
913     {
914     GMP_INTS_mpz_copy(&(right->num), cf_rep->p + ho_conv);
915     GMP_INTS_mpz_copy(&(right->den), cf_rep->q + ho_conv);
916     }
917    
918     //Now, we need to calculate the other neighbor based
919     //on the standard formula. This is a little tricky
920     //because we don't have the -1-th order convergents
921     //tabulated so we have to fudge a little bit.
922     if (ho_conv == 0)
923     {
924     GMP_RATS_mpq_set_si(&previous_convergent, 1, 0);
925     }
926     else
927     {
928     GMP_INTS_mpz_copy(&(previous_convergent.num), cf_rep->p + ho_conv - 1);
929     GMP_INTS_mpz_copy(&(previous_convergent.den), cf_rep->q + ho_conv - 1);
930     }
931    
932     //Calculate the other neighbor according to the standard
933     //formula.
934     GMP_INTS_mpz_sub(&temp1, N, &(previous_convergent.den));
935     GMP_INTS_mpz_tdiv_qr(&temp2, &temp3, &temp1, cf_rep->q + ho_conv);
936     //temp2 now contains term from floor() function in formula.
937     GMP_INTS_mpz_mul(&temp1, &temp2, cf_rep->p + ho_conv);
938     GMP_INTS_mpz_add(&(other_neighbor.num), &temp1, &(previous_convergent.num));
939     GMP_INTS_mpz_mul(&temp1, &temp2, cf_rep->q + ho_conv);
940     GMP_INTS_mpz_add(&(other_neighbor.den), &temp1, &(previous_convergent.den));
941    
942     //Copy the other neighbor into the right slot.
943     if (INTFUNC_is_even(ho_conv))
944     {
945     GMP_RATS_mpq_copy(right, &other_neighbor);
946     }
947     else
948     {
949     GMP_RATS_mpq_copy(left, &other_neighbor);
950     }
951    
952     //Set the equality flag. We have equality if and only
953     //if the denominator of the rational number is less than
954     //or equal to N.
955     if (GMP_INTS_mpz_cmp(&(rn_abs.den), N) <= 0)
956     {
957     *equality = 1;
958     }
959     else
960     {
961     *equality = 0;
962     }
963    
964     //In the event of equality, we don't really have the true
965     //neighbors. If the final convergent is even-ordered,
966     //the left needs to be replaced. If the final convergent
967     //is odd-ordered, the right needs to be replaced.
968     if (*equality)
969     {
970     if (INTFUNC_is_even(ho_conv))
971     {
972     //Left needs to be replaced.
973     GMP_RALG_farey_predecessor(
974     left,
975     right,
976     &rn_abs,
977     N);
978     }
979     else
980     {
981     //Right needs to be replaced.
982     GMP_RALG_farey_successor(
983     right,
984     left,
985     &rn_abs,
986     N);
987     }
988     }
989    
990     //OK, we should be all done. The final catch is that if
991     //the number supplied was negative, we need to invert
992     //and re-order the neighbors.
993     if (GMP_INTS_mpz_is_neg(&(rn_in->num)))
994     {
995     GMP_RATS_mpq_swap(left, right);
996     GMP_INTS_mpz_negate(&(left->num));
997     GMP_INTS_mpz_negate(&(right->num));
998     }
999     } //End if (rn==0) else clause
1000    
1001     //Deallocate dynamic variables.
1002     GMP_RATS_mpq_clear(&rn_abs);
1003     GMP_RATS_mpq_clear(&previous_convergent);
1004     GMP_RATS_mpq_clear(&other_neighbor);
1005     GMP_INTS_mpz_clear(&temp1);
1006     GMP_INTS_mpz_clear(&temp2);
1007     GMP_INTS_mpz_clear(&temp3);
1008     GMP_INTS_mpz_clear(&temp4);
1009     }
1010    
1011    
1012    
1013     //08/16/01: Visual inspection OK. Did not fully inspect the
1014     //iterative part of this function. Unit testing will be
1015     //careful, expect that to catch any anomalies.
1016     void GMP_RALG_consecutive_fab_terms(
1017     const GMP_RATS_mpq_struct *rn_in,
1018     const GMP_INTS_mpz_struct *kmax,
1019     const GMP_INTS_mpz_struct *hmax,
1020     int n_left_in,
1021     int n_right_in,
1022     GMP_RALG_fab_neighbor_collection_struct *result
1023     )
1024     {
1025     int error_flag, equality_flag;
1026     char *estring_kmax_neg = "KMAX is zero, negative, or NAN.";
1027     char *estring_hmax_neg = "HMAX is negative or NAN.";
1028     char *estring_general = "Unspecified general error in GMP_RALG_consecutive_fab_terms().";
1029    
1030     GMP_RATS_mpq_struct q_temp1, q_temp2, q_temp3, q_temp4,
1031     left_neighbor, right_neighbor,
1032     left_neighbor_abs, right_neighbor_abs,
1033     hmax_over_one_neg, corner_point_neg,
1034     abs_norm_recip_rn;
1035    
1036     //Eyeball input parameters.
1037     assert(rn_in != NULL);
1038     assert(kmax != NULL);
1039     assert(n_left_in >= 0);
1040     assert(n_left_in <= 0x00FFFFFF);
1041     assert(n_right_in >= 0);
1042     assert(n_right_in <= 0x00FFFFFF);
1043     assert(result != NULL);
1044    
1045     //Allocate all of the dynamic memory we'll need for this function. It will be
1046     //released at the end.
1047     GMP_RATS_mpq_init(&q_temp1);
1048     GMP_RATS_mpq_init(&q_temp2);
1049     GMP_RATS_mpq_init(&q_temp3);
1050     GMP_RATS_mpq_init(&q_temp4);
1051     GMP_RATS_mpq_init(&left_neighbor);
1052     GMP_RATS_mpq_init(&right_neighbor);
1053     GMP_RATS_mpq_init(&left_neighbor_abs);
1054     GMP_RATS_mpq_init(&right_neighbor_abs);
1055     GMP_RATS_mpq_init(&hmax_over_one_neg);
1056     GMP_RATS_mpq_init(&corner_point_neg);
1057     GMP_RATS_mpq_init(&abs_norm_recip_rn);
1058    
1059     //Zero out the result block. This is the easiest way to give many variables
1060     //default values of 0, FALSE, and NULL.
1061     memset(result, 0, sizeof(GMP_RALG_fab_neighbor_collection_struct));
1062    
1063     //Allocate all integer and rational number structures in the result block.
1064     GMP_RATS_mpq_init(&(result->rn_in));
1065     GMP_INTS_mpz_init(&(result->kmax_in));
1066     GMP_INTS_mpz_init(&(result->hmax_in));
1067     GMP_RATS_mpq_init(&(result->hmax_over_one));
1068     GMP_RATS_mpq_init(&(result->corner_point));
1069     GMP_RATS_mpq_init(&(result->corner_point_minus_one));
1070     GMP_RATS_mpq_init(&(result->corner_point_plus_one));
1071     GMP_RATS_mpq_init(&(result->norm_rn));
1072     GMP_RATS_mpq_init(&(result->abs_norm_rn));
1073    
1074     //Fill in the rational number, exactly as passed.
1075     GMP_RATS_mpq_copy(&(result->rn_in), rn_in);
1076    
1077     //Fill in the number of left and right neighbors that the caller wants.
1078     //However, let's of course say nothing less than zero and nothing more
1079     //than 10000 neighbors on either side.
1080     result->n_left_in = INTFUNC_min(INTFUNC_max(0, n_left_in), 10000);
1081     result->n_right_in = INTFUNC_min(INTFUNC_max(0, n_right_in), 10000);
1082    
1083     //Fill in the value of KMAX, exactly as passed. If it is not at least
1084     //the value of 1 or if error flags, croak.
1085     GMP_INTS_mpz_copy(&(result->kmax_in), kmax);
1086     if (GMP_INTS_mpz_get_flags(kmax) || GMP_INTS_mpz_is_zero(kmax) || GMP_INTS_mpz_is_neg(kmax))
1087     {
1088     result->error =
1089     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1090     CCMALLOC_malloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1091     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1092     TclpAlloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1093     #else
1094     malloc(sizeof(char) * (strlen(estring_kmax_neg) + 1));
1095     #endif
1096     strcpy(result->error, estring_kmax_neg);
1097     goto return_point;
1098     }
1099    
1100     //Decide whether the caller intends to use HMAX. Neg is error, but zero
1101     //is a signal that don't intend to use.
1102     if (hmax)
1103     {
1104     GMP_INTS_mpz_copy(&(result->hmax_in), hmax);
1105     if (GMP_INTS_mpz_get_flags(hmax) || GMP_INTS_mpz_is_neg(hmax))
1106     {
1107     result->error =
1108     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1109     CCMALLOC_malloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1110     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1111     TclpAlloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1112     #else
1113     malloc(sizeof(char) * (strlen(estring_hmax_neg) + 1));
1114     #endif
1115     strcpy(result->error, estring_hmax_neg);
1116     goto return_point;
1117     }
1118     else if (GMP_INTS_mpz_is_pos(hmax))
1119     {
1120     result->hmax_supplied = 1;
1121     }
1122     }
1123    
1124     //If HMAX has been supplied, assign and normalize the
1125     //corner point.
1126     if (result->hmax_supplied)
1127     {
1128     GMP_INTS_mpz_copy(&(result->corner_point.num), &(result->hmax_in));
1129     GMP_INTS_mpz_copy(&(result->corner_point.den), &(result->kmax_in));
1130     GMP_RATS_mpq_normalize(&(result->corner_point));
1131     }
1132    
1133     //If HMAX has been supplied, we want to get the continued
1134     //fraction representation of both the corner point and its
1135     //reciprocal. This is because we're going to need to
1136     //find its adjacent points so we can easily crawl
1137     //around a rectangular region of the integer lattice.
1138     if (result->hmax_supplied)
1139     {
1140     //CF representation of corner point.
1141     GMP_RALG_cfdecomp_init(&(result->corner_point_cf_rep),
1142     &error_flag,
1143     &(result->corner_point.num),
1144     &(result->corner_point.den));
1145     result->corner_point_cf_rep_formed = 1;
1146    
1147     //If there was an error forming the CF representation
1148     //of the corner point, bail out.
1149     if (error_flag)
1150     {
1151     result->error =
1152     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1153     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1154     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1155     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1156     #else
1157     malloc(sizeof(char) * (strlen(estring_general) + 1));
1158     #endif
1159     strcpy(result->error, estring_general);
1160     goto return_point;
1161     }
1162    
1163     //CF representation of reciprocal of corner point.
1164     GMP_RALG_cfdecomp_init(&(result->corner_point_recip_cf_rep),
1165     &error_flag,
1166     &(result->corner_point.den),
1167     &(result->corner_point.num));
1168     result->corner_point_recip_cf_rep_formed = 1;
1169    
1170     //If there was an error forming the CF representation
1171     //of the reciprocal of the corner point, bail out.
1172     if (error_flag)
1173     {
1174     result->error =
1175     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1176     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1177     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1178     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1179     #else
1180     malloc(sizeof(char) * (strlen(estring_general) + 1));
1181     #endif
1182     strcpy(result->error, estring_general);
1183     goto return_point;
1184     }
1185     }
1186    
1187     //Normalize the rational number supplied.
1188     GMP_RATS_mpq_copy(&(result->norm_rn), rn_in);
1189     GMP_RATS_mpq_normalize(&(result->norm_rn));
1190    
1191     //Form the absolute value of the normalized
1192     //version, and set the neg flag.
1193     GMP_RATS_mpq_copy(&(result->abs_norm_rn),&(result->norm_rn));
1194     if (GMP_INTS_mpz_is_neg(&(result->abs_norm_rn.num)))
1195     {
1196     GMP_INTS_mpz_negate(&(result->abs_norm_rn.num));
1197     result->rn_is_neg = 1;
1198     }
1199    
1200     //Form the continued fraction representation of the
1201     //absolute value of the rational number supplied.
1202     //This is always required, because we cannot get any
1203     //neighbors without it.
1204     GMP_RALG_cfdecomp_init(&(result->rn_abs_cf_rep),
1205     &error_flag,
1206     &(result->abs_norm_rn.num),
1207     &(result->abs_norm_rn.den));
1208     result->rn_abs_cf_rep_formed = 1;
1209    
1210     //If there was an error forming the CF representation
1211     //of the absolute value of rational number supplied, bail out.
1212     if (error_flag)
1213     {
1214     result->error =
1215     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1216     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1217     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1218     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1219     #else
1220     malloc(sizeof(char) * (strlen(estring_general) + 1));
1221     #endif
1222     strcpy(result->error, estring_general);
1223     goto return_point;
1224     }
1225    
1226     //Set the equality flag. The rational number supplied is
1227     //in the series of interest if and only if, in its normalized
1228     //form, K <= KMAX, and if HMAX was supplied, H <= HMAX.
1229     if (GMP_INTS_mpz_cmp(&(result->abs_norm_rn.den), kmax) <= 0)
1230     {
1231     if (result->hmax_supplied)
1232     {
1233     if (GMP_INTS_mpz_cmp(&(result->abs_norm_rn.num), hmax) <= 0)
1234     {
1235     result->equality = 1;
1236     }
1237     else
1238     {
1239     result->equality = 0;
1240     }
1241     }
1242     else
1243     {
1244     result->equality = 1;
1245     }
1246     }
1247     else
1248     {
1249     result->equality = 0;
1250     }
1251    
1252     //The final cause of error is if the rational number
1253     //supplied is outside the interval [-HMAX/1, HMAX/1].
1254     //In such cases, simply refuse to calculate
1255     //any approximations. This error can only occur
1256     //if HMAX is specified. If only KMAX is specified,
1257     //this error cannot occur.
1258     if (result->hmax_supplied)
1259     {
1260     //Form the rational number HMAX/1. We will use it for
1261     //a comparison.
1262     GMP_INTS_mpz_copy(&(result->hmax_over_one.num), hmax);
1263     GMP_INTS_mpz_set_ui(&(result->hmax_over_one.den), 1);
1264    
1265     //If the comparison shows that the absolute value of
1266     //the rational number in is larger than HMAX over 1,
1267     //then declare an error.
1268     if (GMP_RATS_mpq_cmp(&(result->abs_norm_rn),&(result->hmax_over_one),NULL) > 0)
1269     {
1270     result->error =
1271     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1272     CCMALLOC_malloc(sizeof(char) * (strlen(estring_general) + 1));
1273     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1274     TclpAlloc(sizeof(char) * (strlen(estring_general) + 1));
1275     #else
1276     malloc(sizeof(char) * (strlen(estring_general) + 1));
1277     #endif
1278     strcpy(result->error, estring_general);
1279     goto return_point;
1280     }
1281     }
1282    
1283     //If we're here, we're very much clean. The only thing
1284     //that could go wrong is an overflow.
1285    
1286     //Allocate space for the left and right arrays of
1287     //neighbors.
1288     if (result->n_left_in)
1289     {
1290     result->lefts =
1291     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1292     CCMALLOC_malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1293     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1294     (GMP_RALG_fab_neighbor_struct *)
1295     TclpAlloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1296     #else
1297     malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_left_in);
1298     #endif
1299     }
1300    
1301     if (result->n_right_in)
1302     {
1303     result->rights =
1304     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
1305     CCMALLOC_malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1306     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
1307     (GMP_RALG_fab_neighbor_struct *)
1308     TclpAlloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1309     #else
1310     malloc(sizeof(GMP_RALG_fab_neighbor_struct) * result->n_right_in);
1311     #endif
1312     }
1313    
1314     //If the rational number supplied is above the corner
1315     //point, we want to form the continued fraction representation
1316     //of its reciprocal.
1317     if (result->hmax_supplied)
1318     {
1319     if (GMP_RATS_mpq_cmp(&(result->abs_norm_rn),&(result->corner_point),NULL) > 0)
1320     {
1321     GMP_RALG_cfdecomp_init(&(result->rn_abs_recip_cf_rep),
1322     &error_flag,
1323     &(result->abs_norm_rn.den),
1324     &(result->abs_norm_rn.num));
1325     result->rn_abs_recip_cf_rep_formed = 1;
1326     }
1327     }
1328    
1329     //If HMAX has been supplied, we want to calculate the points just below and above
1330     //the corner point. The reason we want to do this is because we need to gracefully
1331     //"round the corner" in either direction.
1332     //
1333     //Calculate the point just to the left of the corner point.
1334     if (result->hmax_supplied)
1335     {
1336     GMP_RALG_enclosing_farey_neighbors(
1337     &(result->corner_point),
1338     &(result->kmax_in),
1339     &(result->corner_point_cf_rep),
1340     &equality_flag,
1341     &(result->corner_point_minus_one),
1342     &(q_temp1)
1343     );
1344     }
1345    
1346     //Calculate the point just to the right of the corner point. This is
1347     //where HMAX is the dominant constraint. We need to find the left
1348     //Farey neighbor to the reciprocal of the corner point in the Farey
1349     //series of order HMAX, then take its reciprocal. There is the possibility
1350     //if KMAX=1 that this point will have a denominator of zero, but we
1351     //will accept that as a number here, since we should never hit it,
1352     //as it will be beyond HMAX/1.
1353     if (result->hmax_supplied)
1354     {
1355     GMP_RATS_mpq_copy(&q_temp1, &(result->corner_point));
1356     GMP_INTS_mpz_swap(&(q_temp1.num), &(q_temp1.den));
1357     GMP_RALG_enclosing_farey_neighbors(
1358     &q_temp1,
1359     &(result->hmax_in),
1360     &(result->corner_point_recip_cf_rep),
1361     &equality_flag,
1362     &(result->corner_point_plus_one),
1363     &(q_temp2)
1364     );
1365     GMP_INTS_mpz_swap(&(result->corner_point_plus_one.num), &(result->corner_point_plus_one.den));
1366     }
1367    
1368     //Calculate the complement of HMAX/1. Nothing that we generate can go beyond
1369     //this to the south.
1370     if (result->hmax_supplied)
1371     {
1372     GMP_RATS_mpq_copy(&(hmax_over_one_neg), &(result->hmax_over_one));
1373     GMP_INTS_mpz_negate(&(hmax_over_one_neg.num));
1374     }
1375    
1376     //Also calculate the complement of HMAX/KMAX.
1377     if (result->hmax_supplied)
1378     {
1379     GMP_RATS_mpq_copy(&(corner_point_neg), &(result->corner_point));
1380     GMP_INTS_mpz_negate(&(corner_point_neg.num));
1381     }
1382    
1383     //Form the reciprocal of the absolute value of the normalized value of
1384     //the rational number in.
1385     GMP_RATS_mpq_copy(&abs_norm_recip_rn, &(result->abs_norm_rn));
1386     GMP_RATS_mpq_swap_components(&abs_norm_recip_rn);
1387    
1388     //OK, now we get down to brass tacks. Iterate first to get the
1389     //left neighbors. The ordinary complexity of this is also compounded
1390     //by the fact that we must handle negative numbers as well--everything
1391     //from -HMAX/1 to HMAX/1.
1392     //
1393     //PSEUDO-CODE:
1394     // a)If the rational number to approximate is <= -HMAX/1 or there are no
1395     // left neighbors requested, terminate with no neighbors on the left.
1396     //
1397     // b)Find the right neighbor of the rational number supplied.
1398     //
1399     // c)Find the left neighbor of the rational number supplied.
1400     //
1401     // d)While (queued_count < count)
1402     //
1403     // e)Queue the left neighbor, queued_count++
1404     //
1405     // f)If (queued_count >= count), break.
1406     //
1407     // g)If (left_neighbor <= -HMAX/1), break
1408     //
1409     // h)Advance the frame one to the left.
1410     //
1411     //**************************************************************************
1412     // a)If the rational number to approximate is <= -HMAX/1 or there are no
1413     // left neighbors requested, terminate with no neighbors on the left.
1414     //**************************************************************************
1415     if ((result->hmax_supplied && GMP_RATS_mpq_cmp(&(result->norm_rn), &hmax_over_one_neg, NULL) <= 0)
1416     || (n_left_in <= 0))
1417     goto done_with_left_neighbors;
1418    
1419     //**************************************************************************
1420     // b)Find the right neighbor of the rational number supplied.
1421     //**************************************************************************
1422     // c)Find the left neighbor of the rational number supplied.
1423     //**************************************************************************
1424     if (!result->hmax_supplied)
1425     {
1426     //In this case, the notion of corner point is meaningless, because
1427     //there is no constraint on H. We can just go on our merry way. Get
1428     //the two neighbors.
1429     GMP_RALG_enclosing_farey_neighbors(
1430     &(result->norm_rn),
1431     &(result->kmax_in),
1432     &(result->rn_abs_cf_rep),
1433     &equality_flag,
1434     &left_neighbor,
1435     &right_neighbor
1436     );
1437     //The enclosing Farey neighbor function is prohibited from identifying the
1438     //rational number itself as a Farey term. If the number is in the Farey
1439     //series, we must replace the right neighbor, otherwise we cannot apply
1440     //the standard recursive formulas.
1441     if (equality_flag)
1442     {
1443     GMP_RATS_mpq_copy(&right_neighbor, &(result->norm_rn));
1444     }
1445     }
1446     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) < 0)
1447     {
1448     //The rational number specified is negative and below the negative corner point.
1449     //This means that HMAX is the dominant constraint. We need to find the
1450     //neighbors in the Farey series of order HMAX, then reorder and invert, etc.
1451     GMP_RALG_enclosing_farey_neighbors(
1452     &abs_norm_recip_rn,
1453     &(result->hmax_in),
1454     &(result->rn_abs_recip_cf_rep),
1455     &equality_flag,
1456     &left_neighbor,
1457     &right_neighbor
1458     );
1459    
1460     //Again, if the number specified was already in the series of interest,
1461     //we need to swap in the right neighbor.
1462     if (equality_flag)
1463     {
1464     GMP_RATS_mpq_copy(&right_neighbor, &abs_norm_recip_rn);
1465     }
1466    
1467     //Take the reciprocal of both neighbors, which will put them out of order,
1468     //then negate them, which will put them back in order.
1469     GMP_RATS_mpq_swap_components(&left_neighbor);
1470     GMP_INTS_mpz_negate(&(left_neighbor.num));
1471     GMP_RATS_mpq_swap_components(&right_neighbor);
1472     GMP_INTS_mpz_negate(&(right_neighbor.num));
1473     }
1474     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) == 0)
1475     {
1476     //The rational number specified is the negative corner point. In this case
1477     //Because we can never return the corner point itself as a left neighbor,
1478     //we need to set the left value to be the negative of one past, and the right
1479     //to be the negative of the corner point.
1480     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_plus_one));
1481     GMP_INTS_mpz_negate(&(left_neighbor.num));
1482     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point));
1483     GMP_INTS_mpz_negate(&(right_neighbor.num));
1484     }
1485     else if ((GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) > 0)
1486     &&
1487     (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) < 0))
1488     {
1489     //The rational number specified is in the area dominated by the KMAX constraint
1490     //between -HMAX/KMAX and HMAX/KMAX. The ordinary Farey neighbor function will
1491     //handle this correctly. Again, we need to adjust the output if the number
1492     //is already formable, because the Farey neighbor function is prohibited from
1493     //returning the number itself as a neighbor.
1494     GMP_RALG_enclosing_farey_neighbors(
1495     &(result->norm_rn),
1496     &(result->kmax_in),
1497     &(result->rn_abs_cf_rep),
1498     &equality_flag,
1499     &left_neighbor,
1500     &right_neighbor
1501     );
1502     //The enclosing Farey neighbor function is prohibited from identifying the
1503     //rational number itself as a Farey term. If the number is in the Farey
1504     //series, we must replace the right neighbor, otherwise we cannot apply
1505     //the standard recursive formulas.
1506     if (equality_flag)
1507     {
1508     GMP_RATS_mpq_copy(&right_neighbor, &(result->norm_rn));
1509     }
1510     }
1511     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) == 0)
1512     {
1513     //The rational number specified is the corner point. In this case
1514     //because we can never return the corner point itself as a left neighbor,
1515     //we need to set the left value to be one before, and the right
1516     //to be the corner point.
1517     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_minus_one));
1518     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point));
1519     }
1520     else
1521     {
1522     //The only possibility left is that the number is positive and above the
1523     //corner point where HMAX is the dominant constraint.
1524     GMP_RALG_enclosing_farey_neighbors(
1525     &abs_norm_recip_rn,
1526     &(result->hmax_in),
1527     &(result->rn_abs_recip_cf_rep),
1528     &equality_flag,
1529     &left_neighbor,
1530     &right_neighbor
1531     );
1532    
1533     //Again, if the number specified was already in the series of interest,
1534     //we need to swap in the neighbor. This time, however, it must be the
1535     //left neighbor because taking the reciprocals will reverse the order.
1536     if (equality_flag)
1537     {
1538     GMP_RATS_mpq_copy(&left_neighbor, &abs_norm_recip_rn);
1539     }
1540    
1541     //Take the reciprocal of both neighbors, which will put them out of order,
1542     //then swap them, which will put them back in order.
1543     GMP_RATS_mpq_swap_components(&left_neighbor);
1544     GMP_RATS_mpq_swap_components(&right_neighbor);
1545     GMP_RATS_mpq_swap(&left_neighbor, &right_neighbor);
1546     }
1547    
1548     #if 0
1549     //Print out the left neighbor and right neighbor determined, for debugging.
1550     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1551     &(left_neighbor.num),
1552     "left_neigh_num");
1553     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1554     &(left_neighbor.den),
1555     "left_neigh_den");
1556     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1557     &(right_neighbor.num),
1558     "right_neigh_num");
1559     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1560     &(right_neighbor.den),
1561     "right_neigh_den");
1562     #endif
1563    
1564    
1565     //**************************************************************************
1566     // d)While (queued_count < count)
1567     //**************************************************************************
1568     while (result->n_left_out < result->n_left_in)
1569     {
1570     //**************************************************************************
1571     // e)Queue the left neighbor, queued_count++
1572     //**************************************************************************
1573     (result->lefts + result->n_left_out)->index = -(result->n_left_out + 1);
1574     GMP_RATS_mpq_init(&((result->lefts + result->n_left_out)->neighbor));
1575     GMP_RATS_mpq_copy(&((result->lefts + result->n_left_out)->neighbor), &left_neighbor);
1576     (result->n_left_out)++;
1577    
1578     //**************************************************************************
1579     // f)If (queued_count >= count), break.
1580     //**************************************************************************
1581     //By the way, this step is to save unnecessary contortions once we've met
1582     //the quota.
1583     if (result->n_left_out >= result->n_left_in)
1584     break;
1585    
1586     //**************************************************************************
1587     // g)If (left_neighbor <= -HMAX/1), break
1588     //**************************************************************************
1589     //This breaks us when we've queued the most negative number we can--can't go
1590     //further. This only applies for cases where KMAX is also specified.
1591     if (result->hmax_supplied
1592     &&
1593     GMP_RATS_mpq_cmp(&left_neighbor, &hmax_over_one_neg, NULL) <= 0)
1594     break;
1595    
1596     //**************************************************************************
1597     // h)Advance the frame one to the left.
1598     //**************************************************************************
1599     //Advancing one frame to the left is a complicated affair, requiring several
1600     //subcases. We break it up into regions which are best visualized using
1601     //a graph of the integer lattice with dots for each rational number.
1602     if (!(result->hmax_supplied))
1603     {
1604     //This is the case where we're are looking only in the
1605     //Farey series.
1606     if (GMP_INTS_mpz_is_pos(&left_neighbor.num))
1607     {
1608     //In this case, the left neighbor and right neighbor
1609     //are both positive, and we can apply the standard
1610     //formulas.
1611     GMP_RALG_farey_predecessor(&q_temp1,
1612     &right_neighbor,
1613     &left_neighbor,
1614     &(result->kmax_in));
1615     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1616     GMP_RATS_mpq_copy(&left_neighbor, &q_temp1);
1617     }
1618     else if (GMP_INTS_mpz_is_zero(&left_neighbor.num))
1619     {
1620     //In this case, we are making the transition from positive
1621     //to negative.
1622     GMP_INTS_mpz_set_si(&(left_neighbor.num), -1);
1623     GMP_INTS_mpz_copy(&(left_neighbor.den), &(result->kmax_in));
1624     GMP_RATS_mpq_set_si(&right_neighbor, 0, 1);
1625     }
1626     else
1627     {
1628     //Here, we are purely negative and decreasing. Need to negate
1629     //the numbers, find successor, then negate.
1630     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1631     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1632     GMP_INTS_mpz_abs(&(q_temp1.num));
1633     GMP_INTS_mpz_abs(&(q_temp2.num));
1634     GMP_RALG_farey_successor(&q_temp3,
1635     &q_temp2,
1636     &q_temp1,
1637     &(result->kmax_in));
1638     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1639     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
1640     GMP_INTS_mpz_negate(&(left_neighbor.num));
1641     }
1642     }
1643     else if (GMP_RATS_mpq_cmp(&left_neighbor, &(result->corner_point), NULL) > 0)
1644     {
1645     //We are above the top corner point. In this case HMAX is the dominant
1646     //constraint, and we find our food by taking reciprocals and applying
1647     //the standard relationships in the Farey series of order HMAX.
1648     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1649     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1650     GMP_RATS_mpq_swap_components(&q_temp1);
1651     GMP_RATS_mpq_swap_components(&q_temp2);
1652     GMP_RALG_farey_successor(&q_temp3,
1653     &q_temp2,
1654     &q_temp1,
1655     &(result->hmax_in));
1656     GMP_RATS_mpq_swap_components(&q_temp3);
1657     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1658     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
1659     }
1660     else if (GMP_RATS_mpq_cmp(&left_neighbor, &(result->corner_point), NULL) == 0)
1661     {
1662     //We are precisely at the corner point. This is where we round the corner.
1663     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1664     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_minus_one));
1665     }
1666     else if (GMP_INTS_mpz_is_pos(&left_neighbor.num))
1667     {
1668     //In this region we are going straight down the Farey series.
1669     GMP_RALG_farey_predecessor(&q_temp1,
1670     &right_neighbor,
1671     &left_neighbor,
1672     &(result->kmax_in));
1673     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1674     GMP_RATS_mpq_copy(&left_neighbor, &q_temp1);
1675     }
1676     else if (GMP_INTS_mpz_is_zero(&left_neighbor.num))
1677     {
1678     //In this case, we are making the transition from positive
1679     //to negative.
1680     GMP_INTS_mpz_set_si(&(left_neighbor.num), -1);
1681     GMP_INTS_mpz_copy(&(left_neighbor.den), &(result->kmax_in));
1682     GMP_RATS_mpq_set_si(&right_neighbor, 0, 1);
1683     }
1684     else if (GMP_RATS_mpq_cmp(&left_neighbor, &corner_point_neg, NULL) > 0)
1685     {
1686     //Here, we are purely negative and decreasing. Need to negate
1687     //the numbers, find successor, then negate.
1688     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1689     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1690     GMP_INTS_mpz_abs(&(q_temp1.num));
1691     GMP_INTS_mpz_abs(&(q_temp2.num));
1692     GMP_RALG_farey_successor(&q_temp3,
1693     &q_temp2,
1694     &q_temp1,
1695     &(result->kmax_in));
1696     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1697     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
1698     GMP_INTS_mpz_negate(&(left_neighbor.num));
1699     }
1700     else if (GMP_RATS_mpq_cmp(&left_neighbor, &corner_point_neg, NULL) == 0)
1701     {
1702     //This is where we are rounding the negative corner.
1703     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1704     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point_plus_one));
1705     GMP_INTS_mpz_negate(&(left_neighbor.num));
1706     }
1707     else
1708     {
1709     //Here we're going in the negative direction along the "bottom" of the
1710     //rectangle.
1711     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1712     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1713     GMP_INTS_mpz_abs(&(q_temp1.num));
1714     GMP_INTS_mpz_abs(&(q_temp2.num));
1715     GMP_RATS_mpq_swap_components(&q_temp1);
1716     GMP_RATS_mpq_swap_components(&q_temp2);
1717     GMP_RALG_farey_predecessor(&q_temp3,
1718     &q_temp2,
1719     &q_temp1,
1720     &(result->hmax_in));
1721     GMP_RATS_mpq_swap_components(&q_temp3);
1722     GMP_INTS_mpz_negate(&(q_temp3.num));
1723     GMP_RATS_mpq_copy(&right_neighbor, &left_neighbor);
1724     GMP_RATS_mpq_copy(&left_neighbor, &q_temp3);
1725     }
1726     }
1727    
1728    
1729     done_with_left_neighbors: ;
1730    
1731     //**************************************************************************
1732     // a)If the rational number to approximate is >= HMAX/1 or there are no
1733     // right neighbors requested, terminate with no neighbors on the right.
1734     //**************************************************************************
1735     if ((result->hmax_supplied && GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->hmax_over_one), NULL) >= 0)
1736     || (n_right_in <= 0))
1737     goto done_with_right_neighbors;
1738    
1739     //**************************************************************************
1740     // b)Find the right neighbor of the rational number supplied.
1741     //**************************************************************************
1742     // c)Find the left neighbor of the rational number supplied.
1743     //**************************************************************************
1744     if (!result->hmax_supplied)
1745     {
1746     //In this case, the notion of corner point is meaningless, because
1747     //there is no constraint on H. We can just go on our merry way. Get
1748     //the two neighbors.
1749     GMP_RALG_enclosing_farey_neighbors(
1750     &(result->norm_rn),
1751     &(result->kmax_in),
1752     &(result->rn_abs_cf_rep),
1753     &equality_flag,
1754     &left_neighbor,
1755     &right_neighbor
1756     );
1757     //The enclosing Farey neighbor function is prohibited from identifying the
1758     //rational number itself as a Farey term. If the number is in the Farey
1759     //series, we must replace the left neighbor, otherwise we cannot apply
1760     //the standard recursive formulas.
1761     if (equality_flag)
1762     {
1763     GMP_RATS_mpq_copy(&left_neighbor, &(result->norm_rn));
1764     }
1765     }
1766     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) < 0)
1767     {
1768     //The rational number specified is negative and below the negative corner point.
1769     //This means that HMAX is the dominant constraint. We need to find the
1770     //neighbors in the Farey series of order HMAX, then reorder and invert, etc.
1771     GMP_RALG_enclosing_farey_neighbors(
1772     &abs_norm_recip_rn,
1773     &(result->hmax_in),
1774     &(result->rn_abs_recip_cf_rep),
1775     &equality_flag,
1776     &left_neighbor,
1777     &right_neighbor
1778     );
1779    
1780     //Again, if the number specified was already in the series of interest,
1781     //we need to swap in the left neighbor.
1782     if (equality_flag)
1783     {
1784     GMP_RATS_mpq_copy(&left_neighbor, &abs_norm_recip_rn);
1785     }
1786    
1787     //Take the reciprocal of both neighbors, which will put them out of order,
1788     //then negate them, which will put them back in order.
1789     GMP_RATS_mpq_swap_components(&left_neighbor);
1790     GMP_INTS_mpz_negate(&(left_neighbor.num));
1791     GMP_RATS_mpq_swap_components(&right_neighbor);
1792     GMP_INTS_mpz_negate(&(right_neighbor.num));
1793     }
1794     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) == 0)
1795     {
1796     //The rational number specified is the negative corner point. In this case
1797     //Because we can never return the corner point itself as a left neighbor,
1798     //we need to set the right value to be the negative of one before, and the left
1799     //to be the negative of the corner point.
1800     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point));
1801     GMP_INTS_mpz_negate(&(left_neighbor.num));
1802     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_minus_one));
1803     GMP_INTS_mpz_negate(&(right_neighbor.num));
1804     }
1805     else if ((GMP_RATS_mpq_cmp(&(result->norm_rn), &corner_point_neg, NULL) > 0)
1806     &&
1807     (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) < 0))
1808     {
1809     //The rational number specified is in the area dominated by the KMAX constraint
1810     //between -HMAX/KMAX and HMAX/KMAX. The ordinary Farey neighbor function will
1811     //handle this correctly. Again, we need to adjust the output if the number
1812     //is already formable, because the Farey neighbor function is prohibited from
1813     //returning the number itself as a neighbor.
1814     GMP_RALG_enclosing_farey_neighbors(
1815     &(result->norm_rn),
1816     &(result->kmax_in),
1817     &(result->rn_abs_cf_rep),
1818     &equality_flag,
1819     &left_neighbor,
1820     &right_neighbor
1821     );
1822     //The enclosing Farey neighbor function is prohibited from identifying the
1823     //rational number itself as a Farey term. If the number is in the Farey
1824     //series, we must replace the left neighbor, otherwise we cannot apply
1825     //the standard recursive formulas.
1826     if (equality_flag)
1827     {
1828     GMP_RATS_mpq_copy(&left_neighbor, &(result->norm_rn));
1829     }
1830     }
1831     else if (GMP_RATS_mpq_cmp(&(result->norm_rn), &(result->corner_point), NULL) == 0)
1832     {
1833     //The rational number specified is the positive corner point. In this case.
1834     //because we can never return the corner point itself as a right neighbor,
1835     //we need to set the right value to be one after, and the left
1836     //to be the corner point.
1837     GMP_RATS_mpq_copy(&left_neighbor, &(result->corner_point));
1838     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_plus_one));
1839     }
1840     else
1841     {
1842     //The only possibility left is that the number is positive and at or above the
1843     //corner point where HMAX is the dominant constraint.
1844     GMP_RALG_enclosing_farey_neighbors(
1845     &abs_norm_recip_rn,
1846     &(result->hmax_in),
1847     &(result->rn_abs_recip_cf_rep),
1848     &equality_flag,
1849     &left_neighbor,
1850     &right_neighbor
1851     );
1852    
1853     //Again, if the number specified was already in the series of interest,
1854     //we need to swap in the neighbor. This time, however, it must be the
1855     //right neighbor because taking the reciprocals will reverse the order.
1856     if (equality_flag)
1857     {
1858     GMP_RATS_mpq_copy(&right_neighbor, &abs_norm_recip_rn);
1859     }
1860    
1861     //Take the reciprocal of both neighbors, which will put them out of order,
1862     //then swap them, which will put them back in order.
1863     GMP_RATS_mpq_swap_components(&left_neighbor);
1864     GMP_RATS_mpq_swap_components(&right_neighbor);
1865     GMP_RATS_mpq_swap(&left_neighbor, &right_neighbor);
1866     }
1867    
1868     #if 0
1869     //Print out the left neighbor and right neighbor determined, for debugging.
1870     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1871     &(left_neighbor.num),
1872     "left_neigh_num");
1873     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1874     &(left_neighbor.den),
1875     "left_neigh_den");
1876     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1877     &(right_neighbor.num),
1878     "right_neigh_num");
1879     GMP_INTS_mpz_long_int_format_to_stream(stdout,
1880     &(right_neighbor.den),
1881     "right_neigh_den");
1882     #endif
1883    
1884    
1885     //**************************************************************************
1886     // d)While (queued_count < count)
1887     //**************************************************************************
1888     while (result->n_right_out < result->n_right_in)
1889     {
1890     //**************************************************************************
1891     // e)Queue the right neighbor, queued_count++
1892     //**************************************************************************
1893     (result->rights + result->n_right_out)->index = (result->n_right_out + 1);
1894     GMP_RATS_mpq_init(&((result->rights + result->n_right_out)->neighbor));
1895     GMP_RATS_mpq_copy(&((result->rights + result->n_right_out)->neighbor), &right_neighbor);
1896     (result->n_right_out)++;
1897    
1898     //**************************************************************************
1899     // f)If (queued_count >= count), break.
1900     //**************************************************************************
1901     //By the way, this step is to save unnecessary contortions once we've met
1902     //the quota.
1903     if (result->n_right_out >= result->n_right_in)
1904     break;
1905    
1906     //**************************************************************************
1907     // g)If (right_neighbor >= HMAX/1), break
1908     //**************************************************************************
1909     //This breaks us when we've queued the most positive number we can--can't go
1910     //further. This only applies for cases where KMAX is also specified.
1911     if (result->hmax_supplied
1912     &&
1913     GMP_RATS_mpq_cmp(&right_neighbor, &(result->hmax_over_one), NULL) >= 0)
1914     break;
1915    
1916     //**************************************************************************
1917     // h)Advance the frame one to the right.
1918     //**************************************************************************
1919     //Advancing one frame to the right is a complicated affair, requiring several
1920     //subcases. We break it up into regions which are best visualized using
1921     //a graph of the integer lattice with dots for each rational number.
1922     if (!(result->hmax_supplied))
1923     {
1924     //This is the case where we're are looking only in the
1925     //Farey series.
1926     if (GMP_INTS_mpz_is_neg(&right_neighbor.num))
1927     {
1928     //Neg and increasing towards zero. Can handle by symmetry.
1929     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1930     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1931     GMP_INTS_mpz_abs(&(q_temp1.num));
1932     GMP_INTS_mpz_abs(&(q_temp2.num));
1933     GMP_RALG_farey_predecessor(&q_temp3,
1934     &q_temp1,
1935     &q_temp2,
1936     &(result->kmax_in));
1937     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
1938     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
1939     GMP_INTS_mpz_negate(&(right_neighbor.num));
1940     }
1941     else if (GMP_INTS_mpz_is_zero(&right_neighbor.num))
1942     {
1943     //Right term just hit zero and are increasing.
1944     //Left will become 0/1, right becomes 1/KMAX.
1945     GMP_RATS_mpq_set_si(&left_neighbor, 0, 1);
1946     GMP_INTS_mpz_set_si(&(right_neighbor.num), 1);
1947     GMP_INTS_mpz_copy(&(right_neighbor.den), &(result->kmax_in));
1948     }
1949     else
1950     {
1951     //Are above zero and increasing. Can use standard Farey
1952     //successor formula.
1953     GMP_RALG_farey_successor(&q_temp1,
1954     &left_neighbor,
1955     &right_neighbor,
1956     &(result->kmax_in));
1957     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
1958     GMP_RATS_mpq_copy(&right_neighbor, &q_temp1);
1959     }
1960     }
1961     else if (GMP_RATS_mpq_cmp(&right_neighbor, &corner_point_neg, NULL) < 0)
1962     {
1963     //We are below the negative corner point and moving towards
1964     //zero, with HMAX the dominant constraint. We can proceed by
1965     //symmetry, producing a Farey successor and negating and inverting.
1966     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1967     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1968     GMP_INTS_mpz_abs(&(q_temp1.num));
1969     GMP_INTS_mpz_abs(&(q_temp2.num));
1970     GMP_RATS_mpq_swap_components(&q_temp1);
1971     GMP_RATS_mpq_swap_components(&q_temp2);
1972     GMP_RALG_farey_successor(&q_temp3,
1973     &q_temp1,
1974     &q_temp2,
1975     &(result->hmax_in));
1976     GMP_RATS_mpq_swap_components(&q_temp3);
1977     GMP_INTS_mpz_negate(&(q_temp3.num));
1978     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
1979     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
1980     }
1981     else if (GMP_RATS_mpq_cmp(&right_neighbor, &corner_point_neg, NULL) == 0)
1982     {
1983     //We are at the bottom negative corner point and need to "round the corner".
1984     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
1985     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_minus_one));
1986     GMP_INTS_mpz_negate(&(right_neighbor.num));
1987     }
1988     else if (GMP_INTS_mpz_is_neg(&right_neighbor.num))
1989     {
1990     //In this region we are going straight up the Farey series approaching
1991     //zero. Need to negate to use standard relationships.
1992     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
1993     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
1994     GMP_INTS_mpz_abs(&(q_temp1.num));
1995     GMP_INTS_mpz_abs(&(q_temp2.num));
1996     GMP_RALG_farey_predecessor(&q_temp3,
1997     &q_temp1,
1998     &q_temp2,
1999     &(result->kmax_in));
2000     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2001     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2002     GMP_INTS_mpz_negate(&(right_neighbor.num));
2003     }
2004     else if (GMP_INTS_mpz_is_zero(&right_neighbor.num))
2005     {
2006     //Zero crossover.
2007     GMP_RATS_mpq_set_si(&left_neighbor, 0, 1);
2008     GMP_INTS_mpz_set_si(&(right_neighbor.num), 1);
2009     GMP_INTS_mpz_copy(&(right_neighbor.den), &(result->kmax_in));
2010     }
2011     else if (GMP_RATS_mpq_cmp(&right_neighbor, &(result->corner_point), NULL) < 0)
2012     {
2013     //Below corner point. Standard relationship applies.
2014     GMP_RALG_farey_successor(&q_temp1,
2015     &left_neighbor,
2016     &right_neighbor,
2017     &(result->kmax_in));
2018     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2019     GMP_RATS_mpq_copy(&right_neighbor, &q_temp1);
2020     }
2021     else if (GMP_RATS_mpq_cmp(&right_neighbor, &(result->corner_point), NULL) == 0)
2022     {
2023     //At the positive corner point.
2024     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2025     GMP_RATS_mpq_copy(&right_neighbor, &(result->corner_point_plus_one));
2026     }
2027     else
2028     {
2029     //Above the positive corner point and heading for HMAX/1.
2030     GMP_RATS_mpq_copy(&q_temp1, &left_neighbor);
2031     GMP_RATS_mpq_copy(&q_temp2, &right_neighbor);
2032     GMP_RATS_mpq_swap_components(&q_temp1);
2033     GMP_RATS_mpq_swap_components(&q_temp2);
2034     GMP_RALG_farey_predecessor(&q_temp3,
2035     &q_temp1,
2036     &q_temp2,
2037     &(result->hmax_in));
2038     GMP_RATS_mpq_swap_components(&q_temp3);
2039     GMP_RATS_mpq_copy(&left_neighbor, &right_neighbor);
2040     GMP_RATS_mpq_copy(&right_neighbor, &q_temp3);
2041     }
2042     }
2043    
2044     done_with_right_neighbors: ;
2045    
2046     //This is a single return point so we catch all the dynamic memory
2047     //deallocation.
2048     return_point:
2049     GMP_RATS_mpq_clear(&q_temp1);
2050     GMP_RATS_mpq_clear(&q_temp2);
2051     GMP_RATS_mpq_clear(&q_temp3);
2052     GMP_RATS_mpq_clear(&q_temp4);
2053     GMP_RATS_mpq_clear(&left_neighbor);
2054     GMP_RATS_mpq_clear(&right_neighbor);
2055     GMP_RATS_mpq_clear(&left_neighbor_abs);
2056     GMP_RATS_mpq_clear(&right_neighbor_abs);
2057     GMP_RATS_mpq_clear(&hmax_over_one_neg);
2058     GMP_RATS_mpq_clear(&corner_point_neg);
2059     GMP_RATS_mpq_clear(&abs_norm_recip_rn);
2060     }
2061    
2062    
2063     //08/16/01: Visual inspection OK.
2064     void GMP_RALG_consecutive_fab_terms_result_free(
2065     GMP_RALG_fab_neighbor_collection_struct *arg
2066     )
2067     {
2068     int i;
2069    
2070     //Eyeball the input.
2071     assert(arg != NULL);
2072    
2073     //Deallocate all rational numbers and integers that MUST be allocated, i.e. they are
2074     //never conditional.
2075     GMP_RATS_mpq_clear(&(arg->rn_in));
2076     GMP_INTS_mpz_clear(&(arg->kmax_in));
2077     GMP_INTS_mpz_clear(&(arg->hmax_in));
2078     GMP_RATS_mpq_clear(&(arg->hmax_over_one));
2079     GMP_RATS_mpq_clear(&(arg->corner_point));
2080     GMP_RATS_mpq_clear(&(arg->corner_point_minus_one));
2081     GMP_RATS_mpq_clear(&(arg->corner_point_plus_one));
2082     GMP_RATS_mpq_clear(&(arg->norm_rn));
2083     GMP_RATS_mpq_clear(&(arg->abs_norm_rn));
2084    
2085     //Destroy any continued fraction representations that were
2086     //formed.
2087     if (arg->rn_abs_cf_rep_formed)
2088     {
2089     GMP_RALG_cfdecomp_destroy(&(arg->rn_abs_cf_rep));
2090     }
2091     if (arg->rn_abs_recip_cf_rep_formed)
2092     {
2093     GMP_RALG_cfdecomp_destroy(&(arg->rn_abs_recip_cf_rep));
2094     }
2095     if(arg->corner_point_cf_rep_formed)
2096     {
2097     GMP_RALG_cfdecomp_destroy(&(arg->corner_point_cf_rep));
2098     }
2099     if(arg->corner_point_recip_cf_rep_formed)
2100     {
2101     GMP_RALG_cfdecomp_destroy(&(arg->corner_point_recip_cf_rep));
2102     }
2103    
2104     //Walk through the lefts, freeing up any allocated rational numbers.
2105     for (i=0; i < arg->n_left_out; i++)
2106     {
2107     GMP_RATS_mpq_clear(&(arg->lefts[i].neighbor));
2108     }
2109    
2110     //Walk through the rights, freeing up any allocated rational numbers.
2111     for (i=0; i < arg->n_right_out; i++)
2112     {
2113     GMP_RATS_mpq_clear(&(arg->rights[i].neighbor));
2114     }
2115    
2116     //Deallocate any area assigned for lefts.
2117     if (arg->lefts)
2118     {
2119     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
2120     CCMALLOC_free(arg->lefts);
2121     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
2122     TclpFree((char *)arg->lefts);
2123     #else
2124     free(arg->lefts);
2125     #endif
2126    
2127     arg->lefts = NULL;
2128     }
2129    
2130     //Deallocate any area assigned for rights.
2131     if (arg->rights)
2132     {
2133     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
2134     CCMALLOC_free(arg->rights);
2135     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
2136     TclpFree((char *)arg->rights);
2137     #else
2138     free(arg->rights);
2139     #endif
2140    
2141     arg->rights = NULL;
2142     }
2143     }
2144    
2145    
2146     //08/16/01: Visual inspection OK.
2147     void GMP_RALG_consecutive_fab_terms_result_dump(
2148     FILE *s,
2149     GMP_RALG_fab_neighbor_collection_struct *arg
2150     )
2151     {
2152     int i;
2153     char buf[250];
2154    
2155     //Eyeball the input parameters.
2156     assert(s != NULL);
2157     assert(arg != NULL);
2158    
2159     //Announce intent.
2160     FCMIOF_stream_bannerheading(s,
2161     "Emitting Neighbor Decomp For Analysis",
2162     1);
2163    
2164     //Dump the fields, in order.
2165     GMP_INTS_mpz_long_int_format_to_stream(s,
2166     &(arg->rn_in.num),
2167     "rn_in_num");
2168     GMP_INTS_mpz_long_int_format_to_stream(s,
2169     &(arg->rn_in.den),
2170     "rn_in_den");
2171     GMP_INTS_mpz_long_int_format_to_stream(s,
2172     &(arg->kmax_in),
2173     "kmax_in");
2174     fprintf(s, " hmax_supplied: %12d\n", arg->hmax_supplied);
2175     GMP_INTS_mpz_long_int_format_to_stream(s,
2176     &(arg->hmax_in),
2177     "hmax_in");
2178     if (arg->error)
2179     {
2180     fprintf(s, " error: %s\n", arg->error);
2181     }
2182     else
2183     {
2184     fprintf(s, " error: NULL\n");
2185     }
2186    
2187     GMP_INTS_mpz_long_int_format_to_stream(s,
2188     &(arg->corner_point.num),
2189     "corner_point_num");
2190     GMP_INTS_mpz_long_int_format_to_stream(s,
2191     &(arg->corner_point.den),
2192     "corner_point_den");
2193    
2194     GMP_INTS_mpz_long_int_format_to_stream(s,
2195     &(arg->corner_point_minus_one.num),
2196     "cor_pt_minus_one_num");
2197     GMP_INTS_mpz_long_int_format_to_stream(s,
2198     &(arg->corner_point_minus_one.den),
2199     "cor_pt_minus_one_den");
2200    
2201     GMP_INTS_mpz_long_int_format_to_stream(s,
2202     &(arg->corner_point_plus_one.num),
2203     "cor_pt_plus_one_num");
2204     GMP_INTS_mpz_long_int_format_to_stream(s,
2205     &(arg->corner_point_plus_one.den),
2206     "cor_pt_plus_one_den");
2207    
2208     GMP_INTS_mpz_long_int_format_to_stream(s,
2209     &(arg->hmax_over_one.num),
2210     "hmax/1_num");
2211     GMP_INTS_mpz_long_int_format_to_stream(s,
2212     &(arg->hmax_over_one.den),
2213     "hmax/1_den");
2214    
2215     GMP_INTS_mpz_long_int_format_to_stream(s,
2216     &(arg->norm_rn.num),
2217     "norm_rn_num");
2218     GMP_INTS_mpz_long_int_format_to_stream(s,
2219     &(arg->norm_rn.den),
2220     "norm_rn_den");
2221    
2222     GMP_INTS_mpz_long_int_format_to_stream(s,
2223     &(arg->abs_norm_rn.num),
2224     "abs_norm_rn_num");
2225     GMP_INTS_mpz_long_int_format_to_stream(s,
2226     &(arg->abs_norm_rn.den),
2227     "abs_norm_rn_den");
2228    
2229     fprintf(s, " rn_is_neg: %12d\n", arg->rn_is_neg);
2230    
2231     fprintf(s, " n_left_in: %12d\n", arg->n_left_in);
2232    
2233     fprintf(s, " n_right_in: %12d\n", arg->n_right_in);
2234    
2235     fprintf(s, "rn_abs_cf_rep_formed: %12d\n", arg->rn_abs_cf_rep_formed);
2236    
2237     if (arg->rn_abs_cf_rep_formed)
2238     {
2239     GMP_RALG_cfdecomp_emit(s, "Abs Of RN In CF Decomp", &(arg->rn_abs_cf_rep), 0, 0, NULL);
2240     }
2241    
2242     fprintf(s, "rn_abs_recip_cf_rep_formed: %12d\n", arg->rn_abs_recip_cf_rep_formed);
2243    
2244     if (arg->rn_abs_recip_cf_rep_formed)
2245     {
2246     GMP_RALG_cfdecomp_emit(s, "Abs Of Recip Of RN In CF Decomp", &(arg->rn_abs_recip_cf_rep), 0, 0, NULL);
2247     }
2248    
2249     fprintf(s, "corner_point_cf_rep_formed: %12d\n", arg->corner_point_cf_rep_formed);
2250    
2251     if (arg->corner_point_cf_rep_formed)
2252     {
2253     GMP_RALG_cfdecomp_emit(s, "Corner Point CF Decomp", &(arg->corner_point_cf_rep), 0, 0, NULL);
2254     }
2255    
2256     fprintf(s, "cor_pt_recip_cf_rep_formed: %12d\n", arg->corner_point_recip_cf_rep_formed);
2257    
2258     if (arg->corner_point_recip_cf_rep_formed)
2259     {
2260     GMP_RALG_cfdecomp_emit(s, "Corner Point Recip CF Decomp", &(arg->corner_point_recip_cf_rep), 0, 0, NULL);
2261     }
2262    
2263     fprintf(s, " equality: %12d\n", arg->equality);
2264    
2265     fprintf(s, " n_left_out: %12d\n", arg->n_left_out);
2266    
2267     fprintf(s, " n_right_out: %12d\n", arg->n_right_out);
2268    
2269     for (i=0; i < arg->n_left_out; i++)
2270     {
2271     sprintf(buf, "Contents Of Left Neighbor Array [%d]", i);
2272     FCMIOF_stream_bannerheading(s,
2273     buf,
2274     0);
2275     fprintf(s, " index: %12d\n", arg->lefts[i].index);
2276     GMP_INTS_mpz_long_int_format_to_stream(s,
2277     &(arg->lefts[i].neighbor.num),
2278     "neighbor_num");
2279     GMP_INTS_mpz_long_int_format_to_stream(s,
2280     &(arg->lefts[i].neighbor.den),
2281     "neighbor_den");
2282     }
2283    
2284     for (i=0; i < arg->n_right_out; i++)
2285     {
2286     sprintf(buf, "Contents Of Right Neighbor Array [%d]", i);
2287     FCMIOF_stream_bannerheading(s,
2288     buf,
2289     0);
2290     fprintf(s, " index: %12d\n", arg->rights[i].index);
2291     GMP_INTS_mpz_long_int_format_to_stream(s,
2292     &(arg->rights[i].neighbor.num),
2293     "neighbor_num");
2294     GMP_INTS_mpz_long_int_format_to_stream(s,
2295     &(arg->rights[i].neighbor.den),
2296     "neighbor_den");
2297     }
2298    
2299     FCMIOF_stream_hline(s);
2300     }
2301    
2302    
2303     /******************************************************************/
2304     /*** VERSION CONTROL REPORTING FUNCTIONS ************************/
2305     /******************************************************************/
2306    
2307     //08/16/01: Visual inspection OK.
2308     const char *GMP_RALG_cvcinfo(void)
2309     {
2310 dashley 56 return("$Header$");
2311 dashley 25 }
2312    
2313    
2314     //08/16/01: Visual inspection OK.
2315     const char *GMP_RALG_hvcinfo(void)
2316     {
2317     return(GMP_RALG_H_VERSION);
2318     }
2319    
2320 dashley 56 //End of gmp_ralg.c.

Properties

Name Value
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25