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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.56  
changed lines
  Added in v.71

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25