1 |
dashley |
56 |
//$Header$
|
2 |
dashley |
25 |
//-------------------------------------------------------------------------------------------------
|
3 |
dashley |
56 |
//This file is part of "David T. Ashley's Shared Source Code", a set of shared components
|
4 |
|
|
//integrated into many of David T. Ashley's projects.
|
5 |
dashley |
25 |
//-------------------------------------------------------------------------------------------------
|
6 |
dashley |
56 |
//This source code and any program in which it is compiled/used is provided under the MIT License,
|
7 |
|
|
//reproduced below.
|
8 |
|
|
//-------------------------------------------------------------------------------------------------
|
9 |
|
|
//Permission is hereby granted, free of charge, to any person obtaining a copy of
|
10 |
|
|
//this software and associated documentation files(the "Software"), to deal in the
|
11 |
|
|
//Software without restriction, including without limitation the rights to use,
|
12 |
|
|
//copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
|
13 |
|
|
//Software, and to permit persons to whom the Software is furnished to do so,
|
14 |
|
|
//subject to the following conditions :
|
15 |
dashley |
25 |
//
|
16 |
dashley |
56 |
//The above copyright notice and this permission notice shall be included in all
|
17 |
|
|
//copies or substantial portions of the Software.
|
18 |
dashley |
25 |
//
|
19 |
dashley |
56 |
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
20 |
|
|
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
21 |
|
|
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
|
22 |
|
|
//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
23 |
|
|
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
24 |
|
|
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
25 |
|
|
//SOFTWARE.
|
26 |
dashley |
25 |
//-------------------------------------------------------------------------------------------------
|
27 |
|
|
/* 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 |
dashley |
56 |
#define GMP_RALG_H_VERSION ("$Header$")
|
344 |
dashley |
25 |
#endif
|
345 |
|
|
|
346 |
dashley |
56 |
//End of gmp_ralg.h.
|