1 |
//$Header$ |
2 |
//------------------------------------------------------------------------------------------------- |
3 |
//This file is part of "David T. Ashley's Shared Source Code", a set of shared components |
4 |
//integrated into many of David T. Ashley's projects. |
5 |
//------------------------------------------------------------------------------------------------- |
6 |
//This source code and any program in which it is compiled/used is provided under the MIT License, |
7 |
//reproduced below. |
8 |
//------------------------------------------------------------------------------------------------- |
9 |
//Permission is hereby granted, free of charge, to any person obtaining a copy of |
10 |
//this software and associated documentation files(the "Software"), to deal in the |
11 |
//Software without restriction, including without limitation the rights to use, |
12 |
//copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the |
13 |
//Software, and to permit persons to whom the Software is furnished to do so, |
14 |
//subject to the following conditions : |
15 |
// |
16 |
//The above copyright notice and this permission notice shall be included in all |
17 |
//copies or substantial portions of the Software. |
18 |
// |
19 |
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
20 |
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
21 |
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE |
22 |
//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
23 |
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
24 |
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE |
25 |
//SOFTWARE. |
26 |
//------------------------------------------------------------------------------------------------- |
27 |
/* GNU multiple precision rational number algorithms |
28 |
** (and some extras). |
29 |
*/ |
30 |
|
31 |
#ifndef GMP_RALG_H_INCLUDED |
32 |
#define GMP_RALG_H_INCLUDED |
33 |
|
34 |
#ifdef MODULE_GMP_RALG |
35 |
#define DECMOD_GMP_RALG |
36 |
#else |
37 |
#define DECMOD_GMP_RALG extern |
38 |
#endif |
39 |
|
40 |
typedef struct |
41 |
{ |
42 |
int n; |
43 |
//The number of continued fraction partial quotients |
44 |
//and convergents included in this structure. Indices |
45 |
//are only valid until one less than this value. |
46 |
int nallocd; |
47 |
//The number of slots allocated for partial quotients, |
48 |
//convergents, etc. This goes in blocks of the size |
49 |
//below to avoid numerous reallocation calls. |
50 |
#define GMP_RALG_CF_ALLOC_INCREMENT (10) |
51 |
GMP_INTS_mpz_struct num; |
52 |
//The numerator of the rational number whose continued |
53 |
//fraction decomposition is contained in this structure. |
54 |
//It must be non-negative. The fraction is not |
55 |
//required to be reduced. |
56 |
GMP_INTS_mpz_struct den; |
57 |
//The denominator. It must be positive. |
58 |
GMP_INTS_mpz_struct *a; |
59 |
//Pointer to array of arbitrary integers. Each |
60 |
//is one of the partial quotients. This array |
61 |
//is allocated out to nallocd items but only n |
62 |
//are used. |
63 |
GMP_INTS_mpz_struct *dividend; |
64 |
//The dividend of each division round. |
65 |
GMP_INTS_mpz_struct *divisor; |
66 |
//The divisor at each round. |
67 |
GMP_INTS_mpz_struct *remainder; |
68 |
//The remainder of each round. |
69 |
GMP_INTS_mpz_struct *p; |
70 |
GMP_INTS_mpz_struct *q; |
71 |
//The convergents. |
72 |
} GMP_RALG_cf_app_struct; |
73 |
|
74 |
|
75 |
//This structure type holds a single neighbor to a rational number we wish |
76 |
//to approximate in either the Farey series of a certain order or |
77 |
//in a rectangular region of the integer lattice. These neighbors |
78 |
//are manipulated in potentially large arrays. |
79 |
typedef struct { |
80 |
int index; |
81 |
//The relative index of the neighbor. 0 means that it is |
82 |
//the reduced and normalized form of the number to be |
83 |
//approximated (although at the present time the number to be |
84 |
//approximated is never included among the neighbors), |
85 |
//"1" means nearest neighbor to the right, |
86 |
//"2" means second nearest neighbor to the right, "-1" means |
87 |
//neighbor to the left, etc. |
88 |
GMP_RATS_mpq_struct neighbor; |
89 |
//The neighbor itself, fully reduced and normalized. |
90 |
} GMP_RALG_fab_neighbor_struct; |
91 |
|
92 |
//This structure type holds a collection of neighbors which were |
93 |
//produced on demand to a rational number, either in a Farey series |
94 |
//or in a rectangular region of the integer lattice. The collection |
95 |
//contains the neighbors plus a great deal of other information. |
96 |
typedef struct { |
97 |
GMP_RATS_mpq_struct rn_in; |
98 |
//The rational number which was supplied to have its neighbors |
99 |
//found. This number is raw as supplied by the caller and |
100 |
//not necessarily normalized. |
101 |
GMP_INTS_mpz_struct kmax_in; |
102 |
//The value of KMAX supplied. |
103 |
int hmax_supplied; |
104 |
//TRUE if the problem includes a bound on the numerator, i.e. |
105 |
//rectangular region of integer lattice rather than Farey series. |
106 |
GMP_INTS_mpz_struct hmax_in; |
107 |
//The value of HMAX if it was supplied, or zero otherwise. |
108 |
char *error; |
109 |
//NULL if no error occured in finding the neighbors. If non-NULL, |
110 |
//points to a descriptive string explaining the error. |
111 |
//If this error flag is set, contents of this structure |
112 |
//below this point should not be used--they are not |
113 |
//necessarily valid. |
114 |
GMP_RATS_mpq_struct corner_point; |
115 |
//The normalized (reduced) version of the "corner point" |
116 |
//(HMAX/KMAX). This will only be filled in if |
117 |
//HMAX is specified, or it will be 0/1 otherwise. |
118 |
GMP_RATS_mpq_struct corner_point_minus_one; |
119 |
GMP_RATS_mpq_struct corner_point_plus_one; |
120 |
//The points just before and just above the corner point |
121 |
//in the rectangular region of interest. |
122 |
GMP_RATS_mpq_struct hmax_over_one; |
123 |
//If HMAX is supplied, the rational number which is |
124 |
//HMAX/1. Used for comparisons if we are asked |
125 |
//to approximate a number out of range. |
126 |
GMP_RATS_mpq_struct norm_rn; |
127 |
//The normalized form of the rational number supplied. |
128 |
GMP_RATS_mpq_struct abs_norm_rn; |
129 |
//The absolute value of the normalized form of the rational |
130 |
//number supplied. |
131 |
int rn_is_neg; |
132 |
//Set TRUE if the rational number supplied was negative. This is |
133 |
//not a problem, because symmetry is exploited. For the purposes |
134 |
//of this structure, the Farey series extends on both sides of zero |
135 |
//forever, and a rectangular region of the integer lattice extends |
136 |
//from -HMAX/1 to HMAX/1. |
137 |
int n_left_in; |
138 |
int n_right_in; |
139 |
//The number of left neighbors to the rational number supplied |
140 |
//requested, and the number of right neighbors. Note that if |
141 |
//one of the bounds of the series of interest is encountered, |
142 |
//the actual number supplied may be less. |
143 |
GMP_RALG_cf_app_struct rn_abs_cf_rep; |
144 |
//The CF representation of the absolute value of the |
145 |
//rational number supplied. This always has to be |
146 |
//formed, because it is the only way to immediately |
147 |
//get neighbors. |
148 |
int rn_abs_cf_rep_formed; |
149 |
//Boolean flag to indicate if above has been formed. |
150 |
//This is necessary because otherwise don't have |
151 |
//any obvious clues about whether to free this. |
152 |
GMP_RALG_cf_app_struct rn_abs_recip_cf_rep; |
153 |
//The CF representation of the reciprocal of the |
154 |
//absolute value of the rational number supplied. |
155 |
//If HMAX is supplied, this may need to be formed, |
156 |
//because it will allow us to zero in on neighbors |
157 |
//where HMAX is the dominant constraint. |
158 |
int rn_abs_recip_cf_rep_formed; |
159 |
//Boolean flag to indicate if above has been formed. |
160 |
//This is necessary because otherwise don't have |
161 |
//any obvious clues about whether to free this. |
162 |
GMP_RALG_cf_app_struct corner_point_cf_rep; |
163 |
//The CF representation of the corner point. This is |
164 |
//necessary to find the predecessor to the corner point |
165 |
//in F_KMAX. This will only be calculated if HMAX |
166 |
//is specified. |
167 |
int corner_point_cf_rep_formed; |
168 |
//Boolean flag to indicate if above has been formed. |
169 |
//This is necessary because otherwise don't have |
170 |
//any obvious clues about whether to free this. |
171 |
GMP_RALG_cf_app_struct corner_point_recip_cf_rep; |
172 |
//The CF representation of the reciprocal of the |
173 |
//corner point. This is necessary to find the |
174 |
//successor to the corner point along the HMAX constaint |
175 |
//in the rectangular region of the integer lattice. |
176 |
//This will only be calculated if HMAX is specified. |
177 |
int corner_point_recip_cf_rep_formed; |
178 |
//Boolean flag to indicate if above has been formed. |
179 |
//This is necessary because otherwise don't have |
180 |
//any obvious clues about whether to free this. |
181 |
int equality; |
182 |
//TRUE if the number supplied was already in the series of interest, |
183 |
//or FALSE otherwise. |
184 |
int n_left_out; |
185 |
int n_right_out; |
186 |
//The actual number of neighbors supplied. These may be less |
187 |
//than the number requested only if a boundary of the series of interest |
188 |
//was encountered. |
189 |
GMP_RALG_fab_neighbor_struct *lefts; |
190 |
GMP_RALG_fab_neighbor_struct *rights; |
191 |
//Arrays of neighbors in the series of interest. These are |
192 |
//only stuffed as indicated by "n_left_out" and "n_right_out"--anything |
193 |
//else will hit on dynamic memory that is not initialized. The |
194 |
//left array begins -1 and proceeds leftward through increasingly |
195 |
//negative indices. The right begins with 1 and |
196 |
//proceeds through increasingly positive indices. Each neighbor is |
197 |
//irreducible and normalized. The rational number to be approximated |
198 |
//is NEVER included among the neighbors. |
199 |
} GMP_RALG_fab_neighbor_collection_struct; |
200 |
|
201 |
/******************************************************************/ |
202 |
/*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/ |
203 |
/******************************************************************/ |
204 |
//Decomposes a rational number into its partial quotients |
205 |
//and convergents. This initializes the "decomp" structure |
206 |
//by allocating space, etc. The numerator must be non-negative |
207 |
//and the denominator must be positive, or else there will be |
208 |
//a failure. To make the interface simpler, the destruction |
209 |
//call must ALWAYS be made, even in the event of a failure. |
210 |
//If a failure occurs, the rational number will be decomposed |
211 |
//as if it were 0/1. This means that if the failure flag |
212 |
//is missed, it will still seem to be a valid rational number. |
213 |
//It is not required that the rational number be reduced. |
214 |
//The two integers passed are not changed at all, but they |
215 |
//are copied into the decomp structure without modification. |
216 |
//After this operation, the destruction call must be made. |
217 |
DECMOD_GMP_RALG void GMP_RALG_cfdecomp_init( |
218 |
GMP_RALG_cf_app_struct *decomp, |
219 |
int *failure, |
220 |
GMP_INTS_mpz_struct *num, |
221 |
GMP_INTS_mpz_struct *den); |
222 |
//Destroys all the the data allocated. |
223 |
DECMOD_GMP_RALG void GMP_RALG_cfdecomp_destroy( |
224 |
GMP_RALG_cf_app_struct *decomp |
225 |
); |
226 |
|
227 |
/******************************************************************/ |
228 |
/*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/ |
229 |
/******************************************************************/ |
230 |
//Prints out a continued fraction decomposition. The banner parameter |
231 |
//may be NULL, which indicates no banner. nf flags whether to |
232 |
//print with "no format", i.e. unformatted strings, separated by |
233 |
//newlines. dap indicates whether to include decimal approximation |
234 |
//information. The dap_denominator only needs to be non-NULL if |
235 |
//"dap" is selected. |
236 |
DECMOD_GMP_RALG |
237 |
void GMP_RALG_cfdecomp_emit( |
238 |
FILE *s, |
239 |
char *banner, |
240 |
GMP_RALG_cf_app_struct *decomp, |
241 |
int nf, |
242 |
int dap, |
243 |
const GMP_INTS_mpz_struct *dap_denominator); |
244 |
|
245 |
/******************************************************************/ |
246 |
/*** FAREY SERIES PREDECESSOR AND SUCCESSOR FUNCTIONS ***********/ |
247 |
/******************************************************************/ |
248 |
//The two functions below compute Farey predecessors and |
249 |
//successors. These are straight sequences of arithmetic |
250 |
//operations--no checking for overflow is done. The behavior |
251 |
//below zero will probably not be correct, so these should |
252 |
//not be used to go below 0/1 or with operands below 0/1. |
253 |
DECMOD_GMP_RALG |
254 |
void GMP_RALG_farey_predecessor( |
255 |
GMP_RATS_mpq_struct *result, |
256 |
const GMP_RATS_mpq_struct *plus_two, |
257 |
const GMP_RATS_mpq_struct *plus_one, |
258 |
const GMP_INTS_mpz_struct *N); |
259 |
DECMOD_GMP_RALG |
260 |
void GMP_RALG_farey_successor( |
261 |
GMP_RATS_mpq_struct *result, |
262 |
const GMP_RATS_mpq_struct *minus_two, |
263 |
const GMP_RATS_mpq_struct *minus_one, |
264 |
const GMP_INTS_mpz_struct *N); |
265 |
|
266 |
//Finds the enclosing Farey neighbors to an arbitrary positive, |
267 |
//zero, or negative rational number which may or may not be in |
268 |
//the Farey series of order N. |
269 |
// |
270 |
//INPUTS |
271 |
// rn: Rational number whose neighbors to find. MUST |
272 |
// be normalized, meaning reduced and sign adjusted |
273 |
// so numerator contains the sign. |
274 |
// N: Order of series. MUST be positive. |
275 |
// cf_rep: The continued fraction representation of the |
276 |
// non-negative version of the rational number |
277 |
// whose neighbors to find. The rational number |
278 |
// need not have been reduced when the continued |
279 |
// fraction representation was formed, i.e. only |
280 |
// partial quotients and convergents are used. |
281 |
//OUTPUTS |
282 |
// equality: Set TRUE if the rational number is in the Farey |
283 |
// series, or FALSE otherwise. This affects the |
284 |
// generation of neighbors. If the number is |
285 |
// already in the Farey series its other neighbor |
286 |
// is generated, i.e. the rational number of |
287 |
// interest is never one of the neighbors returned. |
288 |
// left: The left neighbor. |
289 |
// right: The right neighbor. |
290 |
DECMOD_GMP_RALG |
291 |
void GMP_RALG_enclosing_farey_neighbors( |
292 |
const GMP_RATS_mpq_struct *rn, |
293 |
const GMP_INTS_mpz_struct *N, |
294 |
const GMP_RALG_cf_app_struct *cf_rep, |
295 |
int *equality, |
296 |
GMP_RATS_mpq_struct *left, |
297 |
GMP_RATS_mpq_struct *right); |
298 |
|
299 |
//Provides up to 2**24 - 1 (in concept--we don't recommend this!) terms in either |
300 |
//the Farey series or in a rectangular region of the integer lattice. The caller |
301 |
//must provide the space to be used, and must free it afterwards--this function |
302 |
//performs no dynamic allocation. |
303 |
DECMOD_GMP_RALG |
304 |
void GMP_RALG_consecutive_fab_terms( |
305 |
const GMP_RATS_mpq_struct *rn_in, |
306 |
//The rational number whose neighbors are to be found. |
307 |
const GMP_INTS_mpz_struct *kmax, |
308 |
//Maximum denominator. |
309 |
const GMP_INTS_mpz_struct *hmax, |
310 |
//Maximum numerator, or a NULL pointer or a non-null pointer but |
311 |
//an arbitrary integer value of 0 if only the Farey series is being |
312 |
//considered. |
313 |
int n_left_in, |
314 |
//The number of left neighbors requested. Must be non-negative. |
315 |
int n_right_in, |
316 |
//The number of right neighbors requested. Must be non-negative. |
317 |
GMP_RALG_fab_neighbor_collection_struct *result |
318 |
//The structure containing all results. Even if there is an error |
319 |
//condition, this must be freed when done by the caller. |
320 |
); |
321 |
|
322 |
//Frees the dynamic memory consumed by an earlier call to generate the |
323 |
//neighbor information. |
324 |
DECMOD_GMP_RALG |
325 |
void GMP_RALG_consecutive_fab_terms_result_free( |
326 |
GMP_RALG_fab_neighbor_collection_struct *arg |
327 |
//The structure whose members to free up. After this call, |
328 |
//the structure must be considered invalid. |
329 |
); |
330 |
|
331 |
//Debugging function. Dumps the contents of the neighbor structure to the indicated |
332 |
//stream for analysis. |
333 |
DECMOD_GMP_RALG |
334 |
void GMP_RALG_consecutive_fab_terms_result_dump( |
335 |
FILE *s, |
336 |
GMP_RALG_fab_neighbor_collection_struct *arg); |
337 |
|
338 |
/******************************************************************/ |
339 |
/*** VERSION CONTROL REPORTING FUNCTIONS ************************/ |
340 |
/******************************************************************/ |
341 |
DECMOD_GMP_RALG const char *GMP_RALG_cvcinfo(void); |
342 |
DECMOD_GMP_RALG const char *GMP_RALG_hvcinfo(void); |
343 |
#define GMP_RALG_H_VERSION ("$Header$") |
344 |
#endif |
345 |
|
346 |
//End of gmp_ralg.h. |