1 |
// $Header: svn://localhost/dtapublic/projs/trunk/projs/20161014_cfbrapab/cfbrapab.cpp 47 2016-10-17 00:04:56Z dashley $
|
2 |
//-------------------------------------------------------------------------------------------------
|
3 |
//This file is part of "cfbrapab", a program for determining the best rational approximation to a
|
4 |
//real number, subject to constraints on the numerator and denominator, using continued fraction
|
5 |
//algorithms.
|
6 |
//-------------------------------------------------------------------------------------------------
|
7 |
//This source code and any program in which it is compiled/used is provided under the MIT License,
|
8 |
//reproduced below.
|
9 |
//-------------------------------------------------------------------------------------------------
|
10 |
//Permission is hereby granted, free of charge, to any person obtaining a copy of
|
11 |
//this software and associated documentation files(the "Software"), to deal in the
|
12 |
//Software without restriction, including without limitation the rights to use,
|
13 |
//copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
|
14 |
//Software, and to permit persons to whom the Software is furnished to do so,
|
15 |
//subject to the following conditions :
|
16 |
//
|
17 |
//The above copyright notice and this permission notice shall be included in all
|
18 |
//copies or substantial portions of the Software.
|
19 |
//
|
20 |
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
21 |
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
22 |
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
|
23 |
//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
24 |
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
25 |
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
26 |
//SOFTWARE.
|
27 |
//-------------------------------------------------------------------------------------------------
|
28 |
#include <assert.h>
|
29 |
#include <malloc.h>
|
30 |
#include <process.h>
|
31 |
#include <stdio.h>
|
32 |
#include <string.h>
|
33 |
#include <time.h>
|
34 |
|
35 |
|
36 |
#include "bstrfunc.h"
|
37 |
#include "ccmalloc.h"
|
38 |
#include "ccmfatal.h"
|
39 |
#include "charfunc.h"
|
40 |
#include "cu_msgs.h"
|
41 |
#include "fcmiof.h"
|
42 |
#include "gmp_ints.h"
|
43 |
#include "gmp_rats.h"
|
44 |
#include "gmp_ralg.h"
|
45 |
#include "intfunc.h"
|
46 |
|
47 |
|
48 |
#define PNAME "cfbrapab"
|
49 |
#define PNAMEUC "CFBRAPAB"
|
50 |
|
51 |
|
52 |
const char *C_MAIN_cvcinfo(void)
|
53 |
{
|
54 |
return("$Header: /cvsroot/esrg/sfesrg/esrgpcpj/cfbrapab/c_main.c,v 1.6 2002/01/27 17:58:15 dtashley Exp $");
|
55 |
}
|
56 |
|
57 |
|
58 |
//This is a NULL-terminated table of pointers to functions
|
59 |
//which return version control strings for all of the files
|
60 |
//which make up the INTFAC program. This information would
|
61 |
//be helpful for debugging.
|
62 |
static const char *(*C_MAIN_vcinfoptrs[])(void) =
|
63 |
{
|
64 |
//This is the main module, should come first.
|
65 |
C_MAIN_cvcinfo,
|
66 |
|
67 |
//And now the others, in alphabetical order.
|
68 |
BSTRFUNC_hvcinfo,
|
69 |
BSTRFUNC_cvcinfo,
|
70 |
CCMALLOC_hvcinfo,
|
71 |
CCMALLOC_cvcinfo,
|
72 |
CCMFATAL_hvcinfo,
|
73 |
CCMFATAL_cvcinfo,
|
74 |
CHARFUNC_hvcinfo,
|
75 |
CHARFUNC_cvcinfo,
|
76 |
CU_MSGS_hvcinfo,
|
77 |
CU_MSGS_cvcinfo,
|
78 |
FCMIOF_hvcinfo,
|
79 |
FCMIOF_cvcinfo,
|
80 |
GMP_INTS_hvcinfo,
|
81 |
GMP_INTS_cvcinfo,
|
82 |
GMP_RALG_hvcinfo,
|
83 |
GMP_RALG_cvcinfo,
|
84 |
GMP_RATS_hvcinfo,
|
85 |
GMP_RATS_cvcinfo,
|
86 |
INTFUNC_hvcinfo,
|
87 |
INTFUNC_cvcinfo,
|
88 |
NULL
|
89 |
};
|
90 |
|
91 |
|
92 |
//This is the structure type used to hold information about all the
|
93 |
//command-line parameters.
|
94 |
//
|
95 |
struct CfbrapabCmainStruct
|
96 |
{
|
97 |
GMP_RATS_mpq_struct rn;
|
98 |
//The rational number specified on the command line.
|
99 |
//symmetry.
|
100 |
GMP_INTS_mpz_struct kmax;
|
101 |
//The value of KMAX specified on the command line. This must always
|
102 |
//be present.
|
103 |
int hmax_specified;
|
104 |
//TRUE if HMAX is specified in addition to KMAX. KMAX is mandatory
|
105 |
//in all cases.
|
106 |
GMP_INTS_mpz_struct hmax;
|
107 |
//The value of HMAX if it is specified. This is optional. This will be
|
108 |
//set to zero if it is not present on the command line.
|
109 |
int neversmaller_specified;
|
110 |
//TRUE if the -neversmaller option is specified on the command line.
|
111 |
int neverlarger_specified;
|
112 |
//TRUE if the -neverlarger option is specified on the command line.
|
113 |
int pred_specified;
|
114 |
//TRUE if the -pred option is specified on the command line.
|
115 |
int succ_specified;
|
116 |
//TRUE if the -succ option specified on the command line.
|
117 |
int n_specified;
|
118 |
//TRUE if the -n parameter is specified on the command line.
|
119 |
unsigned n;
|
120 |
//The value of n if it has been specified.
|
121 |
CU_MSGS_std_cmd_line_par_results_struct argblock;
|
122 |
//The block holding the options which are common across all
|
123 |
//of these command-line utilities.
|
124 |
};
|
125 |
|
126 |
|
127 |
//Processes the command-line parameters, and abstracts it to a
|
128 |
//the contents of a structure plus a failure flag.
|
129 |
static void process_command_line_args(struct CfbrapabCmainStruct *parblock,
|
130 |
int argc,
|
131 |
char* argv[])
|
132 |
{
|
133 |
int error_flag;
|
134 |
int first_dashed_parameter;
|
135 |
int i;
|
136 |
int recognized;
|
137 |
|
138 |
//Eyeball the input parameters.
|
139 |
assert(parblock != NULL);
|
140 |
assert(argc >= 1);
|
141 |
assert(argv != NULL);
|
142 |
|
143 |
//We have to have at least 3 total parameters. However, this is covered
|
144 |
//in main().
|
145 |
|
146 |
//Process the first parameter, which has to be the rational number we
|
147 |
//want to approximate. If there is a problem, give a helpful message
|
148 |
//and exit with an error code.
|
149 |
GMP_RATS_mpq_init(&(parblock->rn));
|
150 |
GMP_RATS_mpq_set_all_format_rat_num( argv[1],
|
151 |
&error_flag,
|
152 |
&(parblock->rn));
|
153 |
|
154 |
//If there was a parse error, announce and abort.
|
155 |
if (error_flag || GMP_RATS_mpq_is_nan(&(parblock->rn)))
|
156 |
{
|
157 |
printf("\"%s\" is not a properly formatted rational number.\n", argv[1]);
|
158 |
exit(4);
|
159 |
}
|
160 |
|
161 |
//Normalize the rational number specified as input. It is allowed to
|
162 |
//be negative.
|
163 |
GMP_RATS_mpq_normalize(&(parblock->rn));
|
164 |
|
165 |
//The next item has to be a number, it has to be
|
166 |
//an integer, it has to be positive, and it
|
167 |
//is KMAX. Parse out that. If there are any
|
168 |
//errors, abort.
|
169 |
GMP_INTS_mpz_init(&(parblock->kmax));
|
170 |
GMP_INTS_mpz_set_general_int(&(parblock->kmax), &error_flag, argv[2]);
|
171 |
if (error_flag || GMP_INTS_mpz_is_zero(&(parblock->kmax)) || GMP_INTS_mpz_is_neg(&(parblock->kmax)))
|
172 |
{
|
173 |
printf("\"%s\" is not a properly formatted positive integer.\n", argv[2]);
|
174 |
exit(4);
|
175 |
}
|
176 |
|
177 |
//Unconditionally allocate space for hmax.
|
178 |
GMP_INTS_mpz_init(&(parblock->hmax));
|
179 |
|
180 |
//If there is a third parameter, it can be two things. It can
|
181 |
//be either HMAX, or it can be the start of the parameters
|
182 |
//with dashes. First, let's decide which case applies.
|
183 |
if (argc <= 3)
|
184 |
{
|
185 |
first_dashed_parameter = 3;
|
186 |
parblock->hmax_specified = 0;
|
187 |
}
|
188 |
else
|
189 |
{
|
190 |
if (argv[3][0] == '-')
|
191 |
{
|
192 |
first_dashed_parameter = 3;
|
193 |
parblock->hmax_specified = 0;
|
194 |
}
|
195 |
else
|
196 |
{
|
197 |
first_dashed_parameter = 4;
|
198 |
parblock->hmax_specified = 1;
|
199 |
|
200 |
GMP_INTS_mpz_set_general_int(&(parblock->hmax), &error_flag, argv[3]);
|
201 |
if (error_flag || GMP_INTS_mpz_is_zero(&(parblock->hmax)) || GMP_INTS_mpz_is_neg(&(parblock->hmax)))
|
202 |
{
|
203 |
printf("\"%s\" is not a properly formatted positive integer.\n", argv[3]);
|
204 |
exit(4);
|
205 |
}
|
206 |
}
|
207 |
}
|
208 |
|
209 |
//Loop through the remaining parameters, trying to process each
|
210 |
//one either as a parameter specific to this program or else
|
211 |
//as a general parameter.
|
212 |
//
|
213 |
//Initialize the internal general parameter block.
|
214 |
CU_MSGS_cmd_line_par_results_struct_create(&(parblock->argblock));
|
215 |
parblock->neversmaller_specified = 0;
|
216 |
parblock->neverlarger_specified = 0;
|
217 |
parblock->pred_specified = 0;
|
218 |
parblock->succ_specified = 0;
|
219 |
parblock->n_specified = 0;
|
220 |
parblock->n = 0;
|
221 |
|
222 |
for (i=first_dashed_parameter; i<argc; i++)
|
223 |
{
|
224 |
if (!strcmp("-neversmaller", argv[i]))
|
225 |
{
|
226 |
parblock->neversmaller_specified = 1;
|
227 |
}
|
228 |
else if (!strcmp("-neverlarger", argv[i]))
|
229 |
{
|
230 |
parblock->neverlarger_specified = 1;
|
231 |
}
|
232 |
else if (!strcmp("-pred", argv[i]))
|
233 |
{
|
234 |
parblock->pred_specified = 1;
|
235 |
}
|
236 |
else if (!strcmp("-succ", argv[i]))
|
237 |
{
|
238 |
parblock->succ_specified = 1;
|
239 |
}
|
240 |
else if (!strcmp("-n", argv[i]))
|
241 |
{
|
242 |
parblock->n_specified = 1;
|
243 |
|
244 |
//To go along with -n, we have to have a next parameter.
|
245 |
if (i == (argc-1))
|
246 |
{
|
247 |
printf("The \"-n\" parameter must include a following count.\n");
|
248 |
exit(4);
|
249 |
}
|
250 |
|
251 |
//Bump i to index to next par.
|
252 |
i++;
|
253 |
|
254 |
//Try to parse this as a UINT24. It must be that.
|
255 |
GMP_INTS_mpz_parse_into_uint32(&(parblock->n), &error_flag, argv[i]);
|
256 |
|
257 |
//If it couldn't be parsed as an integer, flunk it.
|
258 |
if (error_flag)
|
259 |
{
|
260 |
printf("\"%s\" is not a valid unsigned integer or exceeds 24 bits.\n", argv[i]);
|
261 |
exit(4);
|
262 |
}
|
263 |
|
264 |
//If it is too large, flunk it.
|
265 |
if (parblock->n > 0x00FFFFFF)
|
266 |
{
|
267 |
printf("\"%s\" is an unsigned integer but exceeds 24 bits.\n", argv[i]);
|
268 |
exit(4);
|
269 |
}
|
270 |
|
271 |
//OK, we're cool ...
|
272 |
}
|
273 |
else
|
274 |
{
|
275 |
//Two possibilities left. Either general parameter, or else unrecognized.
|
276 |
CU_MSGS_cmd_line_par_results_struct_process_arg(&(parblock->argblock),
|
277 |
argv[i],
|
278 |
&recognized);
|
279 |
if (!recognized)
|
280 |
{
|
281 |
printf("\"%s\" is not a recognized command-line parameter.\n", argv[i]);
|
282 |
exit(4);
|
283 |
}
|
284 |
|
285 |
//Was picked up as general parameter.
|
286 |
}
|
287 |
}
|
288 |
|
289 |
//Congeal our thoughts on the "general" command-line parameters. No errors possible
|
290 |
//here.
|
291 |
CU_MSGS_cmd_line_par_results_struct_finalize(&(parblock->argblock));
|
292 |
|
293 |
//printf("Boo.\n");
|
294 |
//printf("neverlarger %d succ %d\n", parblock->neverlarger_specified, parblock->succ_specified);
|
295 |
|
296 |
//Look for mutually exclusive options among the program-specific parameters.
|
297 |
if (
|
298 |
(parblock->neversmaller_specified && (parblock->neverlarger_specified || parblock->pred_specified || parblock->succ_specified|| parblock->n_specified))
|
299 |
||
|
300 |
(parblock->neverlarger_specified && (parblock->pred_specified || parblock->succ_specified || parblock->n_specified))
|
301 |
||
|
302 |
(parblock->pred_specified && (parblock->succ_specified|| parblock->n_specified))
|
303 |
||
|
304 |
(parblock->succ_specified && parblock->n_specified)
|
305 |
)
|
306 |
{
|
307 |
printf("The \"-neversmaller\", \"-neverlarger\", \"-pred\", \"-succ\", and \"-n\" options are\nmutually exclusive.\n");
|
308 |
exit(4);
|
309 |
}
|
310 |
|
311 |
//OK, we're clean, all pars in order.
|
312 |
}
|
313 |
|
314 |
|
315 |
//Releases the dynamic memory associated with the parameter block.
|
316 |
static void release_command_line_args(struct CfbrapabCmainStruct *parblock)
|
317 |
{
|
318 |
assert(parblock != NULL);
|
319 |
|
320 |
//This function is superfluous, since in a command-line utility it doesn't really
|
321 |
//matter if everything is released. But, here goes.
|
322 |
CU_MSGS_cmd_line_par_results_struct_destroy(&(parblock->argblock));
|
323 |
GMP_RATS_mpq_clear(&(parblock->rn));
|
324 |
GMP_INTS_mpz_clear(&(parblock->kmax));
|
325 |
GMP_INTS_mpz_clear(&(parblock->hmax));
|
326 |
}
|
327 |
|
328 |
|
329 |
//Prints out a single rational number in the format endorsed
|
330 |
//by this program. This often includes DAP information
|
331 |
//and difference information. It is assumed that the
|
332 |
//previous information is terminated by a horizontal line,
|
333 |
//and this function terminates with a horizontal line.
|
334 |
static void CMAIN_print_app_in_std_form(FILE *s,
|
335 |
int index,
|
336 |
GMP_RATS_mpq_struct *rn,
|
337 |
GMP_RATS_mpq_struct *approx,
|
338 |
int nf,
|
339 |
int show_diff,
|
340 |
int show_dap,
|
341 |
GMP_INTS_mpz_struct *dap_den)
|
342 |
{
|
343 |
char sbuf[250];
|
344 |
GMP_RATS_mpq_struct diff, q_temp1;
|
345 |
GMP_INTS_mpz_struct z_temp1, quotient, remainder;
|
346 |
|
347 |
//Eyeball the input parameters.
|
348 |
assert(s != NULL);
|
349 |
assert(rn != NULL);
|
350 |
assert(approx != NULL);
|
351 |
assert(dap_den != NULL);
|
352 |
|
353 |
//Allocate.
|
354 |
GMP_RATS_mpq_init(&diff);
|
355 |
GMP_RATS_mpq_init(&q_temp1);
|
356 |
GMP_INTS_mpz_init(&z_temp1);
|
357 |
GMP_INTS_mpz_init("ient);
|
358 |
GMP_INTS_mpz_init(&remainder);
|
359 |
|
360 |
//Print out the approximation numerator.
|
361 |
sprintf(sbuf, "approx_num(%d)", index);
|
362 |
if (!nf)
|
363 |
{
|
364 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
365 |
&(approx->num),
|
366 |
sbuf);
|
367 |
}
|
368 |
else
|
369 |
{
|
370 |
int nreserved;
|
371 |
char *p;
|
372 |
|
373 |
fprintf(s, "%d\n", index);
|
374 |
|
375 |
nreserved = GMP_INTS_mpz_size_in_base_10(&(approx->num));
|
376 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
377 |
GMP_INTS_mpz_to_string(p, &(approx->num));
|
378 |
fprintf(s, "%s\n", p);
|
379 |
CCMALLOC_free(p);
|
380 |
}
|
381 |
|
382 |
if (!nf)
|
383 |
FCMIOF_hline();
|
384 |
|
385 |
//Print out the approximation denominator.
|
386 |
sprintf(sbuf, "approx_den(%d)", index);
|
387 |
if (!nf)
|
388 |
{
|
389 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
390 |
&(approx->den),
|
391 |
sbuf);
|
392 |
}
|
393 |
else
|
394 |
{
|
395 |
int nreserved;
|
396 |
char *p;
|
397 |
|
398 |
nreserved = GMP_INTS_mpz_size_in_base_10(&(approx->den));
|
399 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
400 |
GMP_INTS_mpz_to_string(p, &(approx->den));
|
401 |
fprintf(s, "%s\n", p);
|
402 |
CCMALLOC_free(p);
|
403 |
}
|
404 |
|
405 |
if (!nf)
|
406 |
FCMIOF_hline();
|
407 |
|
408 |
|
409 |
//If the "dap" flag is set, calculate and display decimal equivalent of the
|
410 |
//approximation.
|
411 |
if (show_dap)
|
412 |
{
|
413 |
//Make the calculation for decimal approximation.
|
414 |
GMP_RATS_mpq_copy(&q_temp1, approx);
|
415 |
GMP_RATS_mpq_normalize(&q_temp1);
|
416 |
GMP_INTS_mpz_mul(&z_temp1, dap_den, &q_temp1.num);
|
417 |
GMP_INTS_mpz_tdiv_qr("ient, &remainder,
|
418 |
&z_temp1, &q_temp1.den);
|
419 |
|
420 |
sprintf(sbuf, "dap_num(%d)", index);
|
421 |
if (!nf)
|
422 |
{
|
423 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
424 |
&(quotient),
|
425 |
sbuf);
|
426 |
}
|
427 |
else
|
428 |
{
|
429 |
int nreserved;
|
430 |
char *p;
|
431 |
|
432 |
nreserved = GMP_INTS_mpz_size_in_base_10(&(quotient));
|
433 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
434 |
GMP_INTS_mpz_to_string(p, &(quotient));
|
435 |
fprintf(s, "%s\n", p);
|
436 |
CCMALLOC_free(p);
|
437 |
}
|
438 |
|
439 |
if (!nf)
|
440 |
FCMIOF_hline();
|
441 |
|
442 |
//Print out the approximation denominator.
|
443 |
sprintf(sbuf, "dap_den(%d)", index);
|
444 |
if (!nf)
|
445 |
{
|
446 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
447 |
dap_den,
|
448 |
sbuf);
|
449 |
}
|
450 |
else
|
451 |
{
|
452 |
int nreserved;
|
453 |
char *p;
|
454 |
|
455 |
nreserved = GMP_INTS_mpz_size_in_base_10(dap_den);
|
456 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
457 |
GMP_INTS_mpz_to_string(p, dap_den);
|
458 |
fprintf(s, "%s\n", p);
|
459 |
CCMALLOC_free(p);
|
460 |
}
|
461 |
|
462 |
if (!nf)
|
463 |
FCMIOF_hline();
|
464 |
}
|
465 |
|
466 |
|
467 |
//If the "diff" flag is set, calculate and display the rational difference.
|
468 |
if (show_diff)
|
469 |
{
|
470 |
GMP_RATS_mpq_sub(&diff, approx, rn);
|
471 |
GMP_RATS_mpq_normalize(&diff);
|
472 |
|
473 |
sprintf(sbuf, "error_num(%d)", index);
|
474 |
if (!nf)
|
475 |
{
|
476 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
477 |
&(diff.num),
|
478 |
sbuf);
|
479 |
}
|
480 |
else
|
481 |
{
|
482 |
int nreserved;
|
483 |
char *p;
|
484 |
|
485 |
nreserved = GMP_INTS_mpz_size_in_base_10(&(diff.num));
|
486 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
487 |
GMP_INTS_mpz_to_string(p, &(diff.num));
|
488 |
fprintf(s, "%s\n", p);
|
489 |
CCMALLOC_free(p);
|
490 |
}
|
491 |
|
492 |
if (!nf)
|
493 |
FCMIOF_hline();
|
494 |
|
495 |
//Print out the approximation denominator.
|
496 |
sprintf(sbuf, "error_den(%d)", index);
|
497 |
if (!nf)
|
498 |
{
|
499 |
GMP_INTS_mpz_long_int_format_to_stream(s,
|
500 |
&(diff.den),
|
501 |
sbuf);
|
502 |
}
|
503 |
else
|
504 |
{
|
505 |
int nreserved;
|
506 |
char *p;
|
507 |
|
508 |
nreserved = GMP_INTS_mpz_size_in_base_10(&(diff.den));
|
509 |
p = CCMALLOC_malloc(sizeof(char) * nreserved);
|
510 |
GMP_INTS_mpz_to_string(p, &(diff.den));
|
511 |
fprintf(s, "%s\n", p);
|
512 |
CCMALLOC_free(p);
|
513 |
}
|
514 |
|
515 |
if (!nf)
|
516 |
FCMIOF_hline();
|
517 |
}
|
518 |
|
519 |
//Deallocate.
|
520 |
GMP_RATS_mpq_clear(&diff);
|
521 |
GMP_RATS_mpq_clear(&q_temp1);
|
522 |
GMP_INTS_mpz_clear(&z_temp1);
|
523 |
GMP_INTS_mpz_clear("ient);
|
524 |
GMP_INTS_mpz_clear(&remainder);
|
525 |
}
|
526 |
|
527 |
|
528 |
//Handles the classic case of finding the closest
|
529 |
//neighbor(s).
|
530 |
static int CMAIN_classic_closest_neighbor(struct CfbrapabCmainStruct *parblock)
|
531 |
{
|
532 |
int rv = 0;
|
533 |
GMP_RATS_mpq_struct hmax_over_one, hmax_over_kmax, rn_in_abs;
|
534 |
GMP_INTS_mpz_struct dap_denominator;
|
535 |
GMP_RALG_cf_app_struct cf_decomp;
|
536 |
GMP_RALG_fab_neighbor_collection_struct neighbor_data;
|
537 |
int error_flag;
|
538 |
|
539 |
//Allocate all dynamic memory.
|
540 |
GMP_RATS_mpq_init(&hmax_over_one);
|
541 |
GMP_RATS_mpq_init(&hmax_over_kmax);
|
542 |
GMP_INTS_mpz_init(&dap_denominator);
|
543 |
GMP_RATS_mpq_init(&rn_in_abs);
|
544 |
|
545 |
//Set the DAP denominator to 1e108.
|
546 |
GMP_INTS_mpz_set_general_int(&dap_denominator,
|
547 |
&error_flag,
|
548 |
"1e108");
|
549 |
|
550 |
//By convention, we will not mess with anything with an
|
551 |
//absolute value greater than HMAX/1. If such a condition exists, puke out.
|
552 |
//Form up the value of HMAX/1 if HMAX was specified.
|
553 |
if (parblock->hmax_specified)
|
554 |
{
|
555 |
GMP_INTS_mpz_copy(&(hmax_over_one.num), &(parblock->hmax));
|
556 |
GMP_INTS_mpz_set_ui(&(hmax_over_one.den), 1);
|
557 |
GMP_RATS_mpq_copy(&rn_in_abs, &(parblock->rn));
|
558 |
GMP_INTS_mpz_abs(&(rn_in_abs.num));
|
559 |
if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_one, NULL) > 0)
|
560 |
{
|
561 |
printf("The magnitude of the number supplied exceeds HMAX/1, and hence the number\n"
|
562 |
"has no neighbors in the series of interest. Calculation cannot continue.\n");
|
563 |
exit(4);
|
564 |
}
|
565 |
}
|
566 |
|
567 |
//If the "verbose" option is specified, we want to give the continued fraction
|
568 |
//partial quotients and convergents of either the number to approximate,
|
569 |
//its reciprocal, or none of the above, as appropriate; and give a bit more
|
570 |
//information, in addition.
|
571 |
if (parblock->argblock.verbose)
|
572 |
{
|
573 |
if (parblock->hmax_specified)
|
574 |
{
|
575 |
//Stuff HMAX/KMAX. This is necessary for comparison.
|
576 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.num), &(parblock->hmax));
|
577 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.den), &(parblock->kmax));
|
578 |
}
|
579 |
|
580 |
if (!(parblock->hmax_specified) || (GMP_RATS_mpq_cmp(&(parblock->rn), &hmax_over_kmax, NULL) < 0))
|
581 |
{
|
582 |
//Either HMAX was not specified or else we are below the corner point on the
|
583 |
//integer lattice. Get the continued fraction representation of the number
|
584 |
//rather than its reciprocal.
|
585 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.num), &(rn_in_abs.den));
|
586 |
|
587 |
//Print out the continued fraction decomposition of the rational number.
|
588 |
GMP_RALG_cfdecomp_emit(stdout,
|
589 |
"CF Rep Of Abs Value Of Number Specified",
|
590 |
&cf_decomp,
|
591 |
0,
|
592 |
1,
|
593 |
&dap_denominator);
|
594 |
|
595 |
//Destroy the decomposition--free the memory.
|
596 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
597 |
}
|
598 |
else if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_kmax, NULL) == 0)
|
599 |
{
|
600 |
//In this case, the rational number specified is exactly the same in
|
601 |
//magnitude as HMAX/KMAX. I am inclined to suppress the CF decomp.
|
602 |
printf("Rational number specified is HMAX/KMAX. CF decomp not provided.\n");
|
603 |
}
|
604 |
else
|
605 |
{
|
606 |
//The number specified is beyond the corner point. It is appropriate to
|
607 |
//provide the decomposition of the reciprocal rather than of the number
|
608 |
//itself.
|
609 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.den), &(rn_in_abs.num));
|
610 |
|
611 |
//Print out the continued fraction decomposition of the rational number.
|
612 |
GMP_RALG_cfdecomp_emit(stdout,
|
613 |
"CF Rep Of Reciprocal Of Abs Value Of Number Specified",
|
614 |
&cf_decomp,
|
615 |
0,
|
616 |
1,
|
617 |
&dap_denominator);
|
618 |
|
619 |
//Destroy the decomposition--free the memory.
|
620 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
621 |
}
|
622 |
} //End if verbose.
|
623 |
|
624 |
//Do all the work to get the neighbors of the number passed.
|
625 |
GMP_RALG_consecutive_fab_terms(
|
626 |
&(parblock->rn),
|
627 |
&(parblock->kmax),
|
628 |
(parblock->hmax_specified) ? (&(parblock->hmax)) : (NULL),
|
629 |
1,
|
630 |
1,
|
631 |
&neighbor_data);
|
632 |
|
633 |
//Print the neighbor data block for debugging.
|
634 |
#if 0
|
635 |
GMP_RALG_consecutive_fab_terms_result_dump(stdout, &neighbor_data);
|
636 |
#endif
|
637 |
|
638 |
//There are four possibilities at this point.
|
639 |
// a)Attempting to find the rational neighbors generated an error.
|
640 |
// b)The rational number specified was already in the series of interest,
|
641 |
// in which case we will use it.
|
642 |
// c)The left neighbor is closer or in a tie we want to choose it.
|
643 |
// d)The right neighbor is closer or in a tie we want to choose it.
|
644 |
if (neighbor_data.error)
|
645 |
{
|
646 |
//
|
647 |
printf("Internal error: %s\n", neighbor_data.error);
|
648 |
}
|
649 |
else if (neighbor_data.equality)
|
650 |
{
|
651 |
CMAIN_print_app_in_std_form(stdout,
|
652 |
0,
|
653 |
&(neighbor_data.rn_in),
|
654 |
&(neighbor_data.rn_in),
|
655 |
parblock->argblock.noformat,
|
656 |
0,
|
657 |
1,
|
658 |
&dap_denominator);
|
659 |
}
|
660 |
else
|
661 |
{
|
662 |
GMP_RATS_mpq_struct left_neighbor, right_neighbor,
|
663 |
left_diff, right_diff,
|
664 |
left_abs, right_abs;
|
665 |
int error_cmp;
|
666 |
int mag_cmp;
|
667 |
|
668 |
GMP_RATS_mpq_init(&left_neighbor);
|
669 |
GMP_RATS_mpq_init(&right_neighbor);
|
670 |
GMP_RATS_mpq_init(&left_diff);
|
671 |
GMP_RATS_mpq_init(&right_diff);
|
672 |
GMP_RATS_mpq_init(&left_abs);
|
673 |
GMP_RATS_mpq_init(&right_abs);
|
674 |
|
675 |
//Snatch the left neighbor.
|
676 |
if (neighbor_data.n_left_out)
|
677 |
{
|
678 |
GMP_RATS_mpq_copy(&left_neighbor, &(neighbor_data.lefts[0].neighbor));
|
679 |
}
|
680 |
|
681 |
//Snatch the right neighbor.
|
682 |
if (neighbor_data.n_right_out)
|
683 |
{
|
684 |
GMP_RATS_mpq_copy(&right_neighbor, &(neighbor_data.rights[0].neighbor));
|
685 |
}
|
686 |
|
687 |
//Calculate the differences, take their absolute
|
688 |
//values.
|
689 |
GMP_RATS_mpq_sub(&left_diff, &left_neighbor, &(neighbor_data.rn_in));
|
690 |
GMP_RATS_mpq_sub(&right_diff, &right_neighbor, &(neighbor_data.rn_in));
|
691 |
GMP_RATS_mpq_normalize(&left_diff);
|
692 |
GMP_RATS_mpq_normalize(&right_diff);
|
693 |
GMP_INTS_mpz_abs(&(left_diff.num));
|
694 |
GMP_INTS_mpz_abs(&(right_diff.num));
|
695 |
|
696 |
//Now that the differences are calculated, take the
|
697 |
//absolute values of the neighbors themselves.
|
698 |
//We will use this to break ties.
|
699 |
GMP_RATS_mpq_normalize(&left_neighbor);
|
700 |
GMP_RATS_mpq_normalize(&right_neighbor);
|
701 |
GMP_INTS_mpz_abs(&(left_neighbor.num));
|
702 |
GMP_INTS_mpz_abs(&(right_neighbor.num));
|
703 |
|
704 |
//Compare the relative differences and magnitudes.
|
705 |
error_cmp = GMP_RATS_mpq_cmp(&left_diff, &right_diff, NULL);
|
706 |
mag_cmp = GMP_RATS_mpq_cmp(&left_neighbor, &right_neighbor, NULL);
|
707 |
|
708 |
//Figure out which to present as the best approximation and
|
709 |
//do it.
|
710 |
if (!(parblock->neversmaller_specified) &&
|
711 |
((parblock->neverlarger_specified) || (error_cmp < 0) || ((error_cmp == 0) && (mag_cmp < 0))))
|
712 |
{
|
713 |
CMAIN_print_app_in_std_form(stdout,
|
714 |
-1,
|
715 |
&(neighbor_data.rn_in),
|
716 |
&(neighbor_data.lefts[0].neighbor),
|
717 |
parblock->argblock.noformat,
|
718 |
1,
|
719 |
1,
|
720 |
&dap_denominator);
|
721 |
}
|
722 |
else
|
723 |
{
|
724 |
CMAIN_print_app_in_std_form(stdout,
|
725 |
1,
|
726 |
&(neighbor_data.rn_in),
|
727 |
&(neighbor_data.rights[0].neighbor),
|
728 |
parblock->argblock.noformat,
|
729 |
1,
|
730 |
1,
|
731 |
&dap_denominator);
|
732 |
}
|
733 |
|
734 |
//Deallocate.
|
735 |
GMP_RATS_mpq_clear(&left_neighbor);
|
736 |
GMP_RATS_mpq_clear(&right_neighbor);
|
737 |
GMP_RATS_mpq_clear(&left_diff);
|
738 |
GMP_RATS_mpq_clear(&right_diff);
|
739 |
GMP_RATS_mpq_clear(&left_abs);
|
740 |
GMP_RATS_mpq_clear(&right_abs);
|
741 |
}
|
742 |
|
743 |
//Deallocate all dynamic memory.
|
744 |
GMP_RALG_consecutive_fab_terms_result_free(&neighbor_data);
|
745 |
GMP_RATS_mpq_clear(&hmax_over_one);
|
746 |
GMP_RATS_mpq_clear(&hmax_over_kmax);
|
747 |
GMP_INTS_mpz_clear(&dap_denominator);
|
748 |
GMP_RATS_mpq_clear(&rn_in_abs);
|
749 |
|
750 |
return(rv);
|
751 |
}
|
752 |
|
753 |
|
754 |
//Handles the case of finding multiple neighbors.
|
755 |
|
756 |
static int CMAIN_multiple_neighbor(struct CfbrapabCmainStruct *parblock)
|
757 |
{
|
758 |
int rv = 0;
|
759 |
GMP_RATS_mpq_struct hmax_over_one, hmax_over_kmax, rn_in_abs;
|
760 |
GMP_INTS_mpz_struct dap_denominator;
|
761 |
GMP_RALG_cf_app_struct cf_decomp;
|
762 |
GMP_RALG_fab_neighbor_collection_struct neighbor_data;
|
763 |
int error_flag;
|
764 |
int i;
|
765 |
|
766 |
//Allocate all dynamic memory.
|
767 |
GMP_RATS_mpq_init(&hmax_over_one);
|
768 |
GMP_RATS_mpq_init(&hmax_over_kmax);
|
769 |
GMP_INTS_mpz_init(&dap_denominator);
|
770 |
GMP_RATS_mpq_init(&rn_in_abs);
|
771 |
|
772 |
//Set the DAP denominator to 1e108.
|
773 |
GMP_INTS_mpz_set_general_int(&dap_denominator,
|
774 |
&error_flag,
|
775 |
"1e108");
|
776 |
|
777 |
//By convention, we will not mess with anything with an
|
778 |
//absolute value greater than HMAX/1. If such a condition exists, puke out.
|
779 |
//Form up the value of HMAX/1 if HMAX was specified.
|
780 |
if (parblock->hmax_specified)
|
781 |
{
|
782 |
GMP_INTS_mpz_copy(&(hmax_over_one.num), &(parblock->hmax));
|
783 |
GMP_INTS_mpz_set_ui(&(hmax_over_one.den), 1);
|
784 |
GMP_RATS_mpq_copy(&rn_in_abs, &(parblock->rn));
|
785 |
GMP_INTS_mpz_abs(&(rn_in_abs.num));
|
786 |
if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_one, NULL) > 0)
|
787 |
{
|
788 |
printf("The magnitude of the number supplied exceeds HMAX/1, and hence the number\n"
|
789 |
"has no neighbors in the series of interest. Calculation cannot continue.\n");
|
790 |
exit(4);
|
791 |
}
|
792 |
}
|
793 |
|
794 |
//If the "verbose" option is specified, we want to give the continued fraction
|
795 |
//partial quotients and convergents of either the number to approximate,
|
796 |
//its reciprocal, or none of the above, as appropriate; and give a bit more
|
797 |
//information, in addition.
|
798 |
if (parblock->argblock.verbose)
|
799 |
{
|
800 |
if (parblock->hmax_specified)
|
801 |
{
|
802 |
//Stuff HMAX/KMAX. This is necessary for comparison.
|
803 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.num), &(parblock->hmax));
|
804 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.den), &(parblock->kmax));
|
805 |
}
|
806 |
|
807 |
if (!(parblock->hmax_specified) || (GMP_RATS_mpq_cmp(&(parblock->rn), &hmax_over_kmax, NULL) < 0))
|
808 |
{
|
809 |
//Either HMAX was not specified or else we are below the corner point on the
|
810 |
//integer lattice. Get the continued fraction representation of the number
|
811 |
//rather than its reciprocal.
|
812 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.num), &(rn_in_abs.den));
|
813 |
|
814 |
//Print out the continued fraction decomposition of the rational number.
|
815 |
GMP_RALG_cfdecomp_emit(stdout,
|
816 |
"CF Representation Of Absolute Value Of Rational Number Specified",
|
817 |
&cf_decomp,
|
818 |
0,
|
819 |
1,
|
820 |
&dap_denominator);
|
821 |
|
822 |
//Destroy the decomposition--free the memory.
|
823 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
824 |
}
|
825 |
else if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_kmax, NULL) == 0)
|
826 |
{
|
827 |
//In this case, the rational number specified is exactly the same in
|
828 |
//magnitude as HMAX/KMAX. I am inclined to suppress the CF decomp.
|
829 |
printf("Rational number specified is HMAX/KMAX. CF decomp not provided.\n");
|
830 |
}
|
831 |
else
|
832 |
{
|
833 |
//The number specified is beyond the corner point. It is appropriate to
|
834 |
//provide the decomposition of the reciprocal rather than of the number
|
835 |
//itself.
|
836 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.den), &(rn_in_abs.num));
|
837 |
|
838 |
//Print out the continued fraction decomposition of the rational number.
|
839 |
GMP_RALG_cfdecomp_emit(stdout,
|
840 |
"CF Representation Of Reciprocal Of Absolute Value Of Rational Number Specified",
|
841 |
&cf_decomp,
|
842 |
0,
|
843 |
1,
|
844 |
&dap_denominator);
|
845 |
|
846 |
//Destroy the decomposition--free the memory.
|
847 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
848 |
}
|
849 |
} //End if verbose.
|
850 |
|
851 |
//Do all the work to get the neighbors of the number passed.
|
852 |
GMP_RALG_consecutive_fab_terms(
|
853 |
&(parblock->rn),
|
854 |
&(parblock->kmax),
|
855 |
(parblock->hmax_specified) ? (&(parblock->hmax)) : (NULL),
|
856 |
parblock->n,
|
857 |
parblock->n,
|
858 |
&neighbor_data);
|
859 |
|
860 |
//Print the neighbor data block for debugging.
|
861 |
#if 0
|
862 |
GMP_RALG_consecutive_fab_terms_result_dump(stdout, &neighbor_data);
|
863 |
#endif
|
864 |
|
865 |
//Loop through, printing out the left neighbors in order.
|
866 |
for (i = neighbor_data.n_left_out - 1; i >= 0; i--)
|
867 |
{
|
868 |
CMAIN_print_app_in_std_form(stdout,
|
869 |
-(i + 1),
|
870 |
&(neighbor_data.rn_in),
|
871 |
&(neighbor_data.lefts[i].neighbor),
|
872 |
parblock->argblock.noformat,
|
873 |
1,
|
874 |
1,
|
875 |
&dap_denominator);
|
876 |
}
|
877 |
|
878 |
//If the number itself appears in the series of interest, spit that out.
|
879 |
if (neighbor_data.equality)
|
880 |
{
|
881 |
CMAIN_print_app_in_std_form(stdout,
|
882 |
0,
|
883 |
&(neighbor_data.rn_in),
|
884 |
&(neighbor_data.norm_rn),
|
885 |
parblock->argblock.noformat,
|
886 |
1,
|
887 |
1,
|
888 |
&dap_denominator);
|
889 |
}
|
890 |
|
891 |
//Loop through, printing out the right neighbors in order.
|
892 |
for (i = 0; i < neighbor_data.n_right_out; i++)
|
893 |
{
|
894 |
CMAIN_print_app_in_std_form(stdout,
|
895 |
i+1,
|
896 |
&(neighbor_data.rn_in),
|
897 |
&(neighbor_data.rights[i].neighbor),
|
898 |
parblock->argblock.noformat,
|
899 |
1,
|
900 |
1,
|
901 |
&dap_denominator);
|
902 |
}
|
903 |
|
904 |
//Deallocate all dynamic memory.
|
905 |
GMP_RALG_consecutive_fab_terms_result_free(&neighbor_data);
|
906 |
GMP_RATS_mpq_clear(&hmax_over_one);
|
907 |
GMP_RATS_mpq_clear(&hmax_over_kmax);
|
908 |
GMP_INTS_mpz_clear(&dap_denominator);
|
909 |
GMP_RATS_mpq_clear(&rn_in_abs);
|
910 |
|
911 |
return(rv);
|
912 |
}
|
913 |
|
914 |
|
915 |
//Handles the case of finding the predecessor.
|
916 |
int CMAIN_predecessor(struct CfbrapabCmainStruct *parblock)
|
917 |
{
|
918 |
int rv = 0;
|
919 |
GMP_RATS_mpq_struct hmax_over_one, hmax_over_kmax, rn_in_abs;
|
920 |
GMP_INTS_mpz_struct dap_denominator;
|
921 |
GMP_RALG_cf_app_struct cf_decomp;
|
922 |
GMP_RALG_fab_neighbor_collection_struct neighbor_data;
|
923 |
int error_flag;
|
924 |
|
925 |
//Allocate all dynamic memory.
|
926 |
GMP_RATS_mpq_init(&hmax_over_one);
|
927 |
GMP_RATS_mpq_init(&hmax_over_kmax);
|
928 |
GMP_INTS_mpz_init(&dap_denominator);
|
929 |
GMP_RATS_mpq_init(&rn_in_abs);
|
930 |
|
931 |
//Set the DAP denominator to 1e108.
|
932 |
GMP_INTS_mpz_set_general_int(&dap_denominator,
|
933 |
&error_flag,
|
934 |
"1e108");
|
935 |
|
936 |
//By convention, we will not mess with anything with an
|
937 |
//absolute value greater than HMAX/1. If such a condition exists, puke out.
|
938 |
//Form up the value of HMAX/1 if HMAX was specified.
|
939 |
if (parblock->hmax_specified)
|
940 |
{
|
941 |
GMP_INTS_mpz_copy(&(hmax_over_one.num), &(parblock->hmax));
|
942 |
GMP_INTS_mpz_set_ui(&(hmax_over_one.den), 1);
|
943 |
GMP_RATS_mpq_copy(&rn_in_abs, &(parblock->rn));
|
944 |
GMP_INTS_mpz_abs(&(rn_in_abs.num));
|
945 |
if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_one, NULL) > 0)
|
946 |
{
|
947 |
printf("The magnitude of the number supplied exceeds HMAX/1, and hence the number\n"
|
948 |
"has no neighbors in the series of interest. Calculation cannot continue.\n");
|
949 |
exit(4);
|
950 |
}
|
951 |
}
|
952 |
|
953 |
//If the "verbose" option is specified, we want to give the continued fraction
|
954 |
//partial quotients and convergents of either the number to approximate,
|
955 |
//its reciprocal, or none of the above, as appropriate; and give a bit more
|
956 |
//information, in addition.
|
957 |
if (parblock->argblock.verbose)
|
958 |
{
|
959 |
if (parblock->hmax_specified)
|
960 |
{
|
961 |
//Stuff HMAX/KMAX. This is necessary for comparison.
|
962 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.num), &(parblock->hmax));
|
963 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.den), &(parblock->kmax));
|
964 |
}
|
965 |
|
966 |
if (!(parblock->hmax_specified) || (GMP_RATS_mpq_cmp(&(parblock->rn), &hmax_over_kmax, NULL) < 0))
|
967 |
{
|
968 |
//Either HMAX was not specified or else we are below the corner point on the
|
969 |
//integer lattice. Get the continued fraction representation of the number
|
970 |
//rather than its reciprocal.
|
971 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.num), &(rn_in_abs.den));
|
972 |
|
973 |
//Print out the continued fraction decomposition of the rational number.
|
974 |
GMP_RALG_cfdecomp_emit(stdout,
|
975 |
"CF Representation Of Absolute Value Of Rational Number Specified",
|
976 |
&cf_decomp,
|
977 |
0,
|
978 |
1,
|
979 |
&dap_denominator);
|
980 |
|
981 |
//Destroy the decomposition--free the memory.
|
982 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
983 |
}
|
984 |
else if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_kmax, NULL) == 0)
|
985 |
{
|
986 |
//In this case, the rational number specified is exactly the same in
|
987 |
//magnitude as HMAX/KMAX. I am inclined to suppress the CF decomp.
|
988 |
printf("Rational number specified is HMAX/KMAX. CF decomp not provided.\n");
|
989 |
}
|
990 |
else
|
991 |
{
|
992 |
//The number specified is beyond the corner point. It is appropriate to
|
993 |
//provide the decomposition of the reciprocal rather than of the number
|
994 |
//itself.
|
995 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.den), &(rn_in_abs.num));
|
996 |
|
997 |
//Print out the continued fraction decomposition of the rational number.
|
998 |
GMP_RALG_cfdecomp_emit(stdout,
|
999 |
"CF Representation Of Reciprocal Of Absolute Value Of Rational Number Specified",
|
1000 |
&cf_decomp,
|
1001 |
0,
|
1002 |
1,
|
1003 |
&dap_denominator);
|
1004 |
|
1005 |
//Destroy the decomposition--free the memory.
|
1006 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
1007 |
}
|
1008 |
} //End if verbose.
|
1009 |
|
1010 |
//Do all the work to get the neighbors of the number passed.
|
1011 |
GMP_RALG_consecutive_fab_terms(
|
1012 |
&(parblock->rn),
|
1013 |
&(parblock->kmax),
|
1014 |
(parblock->hmax_specified) ? (&(parblock->hmax)) : (NULL),
|
1015 |
1,
|
1016 |
0,
|
1017 |
&neighbor_data);
|
1018 |
|
1019 |
//Print the neighbor data block for debugging.
|
1020 |
#if 0
|
1021 |
GMP_RALG_consecutive_fab_terms_result_dump(stdout, &neighbor_data);
|
1022 |
#endif
|
1023 |
|
1024 |
//Print the neighbor on the left, if it exists.
|
1025 |
if (neighbor_data.n_left_out)
|
1026 |
{
|
1027 |
CMAIN_print_app_in_std_form(stdout,
|
1028 |
-1,
|
1029 |
&(neighbor_data.rn_in),
|
1030 |
&(neighbor_data.lefts[0].neighbor),
|
1031 |
parblock->argblock.noformat,
|
1032 |
0,
|
1033 |
0,
|
1034 |
&dap_denominator);
|
1035 |
}
|
1036 |
|
1037 |
//Deallocate all dynamic memory.
|
1038 |
GMP_RALG_consecutive_fab_terms_result_free(&neighbor_data);
|
1039 |
GMP_RATS_mpq_clear(&hmax_over_one);
|
1040 |
GMP_RATS_mpq_clear(&hmax_over_kmax);
|
1041 |
GMP_INTS_mpz_clear(&dap_denominator);
|
1042 |
GMP_RATS_mpq_clear(&rn_in_abs);
|
1043 |
|
1044 |
return(rv);
|
1045 |
}
|
1046 |
|
1047 |
|
1048 |
//Handles the case of finding the successor.
|
1049 |
int CMAIN_successor(struct CfbrapabCmainStruct *parblock)
|
1050 |
{
|
1051 |
int rv = 0;
|
1052 |
GMP_RATS_mpq_struct hmax_over_one, hmax_over_kmax, rn_in_abs;
|
1053 |
GMP_INTS_mpz_struct dap_denominator;
|
1054 |
GMP_RALG_cf_app_struct cf_decomp;
|
1055 |
GMP_RALG_fab_neighbor_collection_struct neighbor_data;
|
1056 |
int error_flag;
|
1057 |
|
1058 |
//Allocate all dynamic memory.
|
1059 |
GMP_RATS_mpq_init(&hmax_over_one);
|
1060 |
GMP_RATS_mpq_init(&hmax_over_kmax);
|
1061 |
GMP_INTS_mpz_init(&dap_denominator);
|
1062 |
GMP_RATS_mpq_init(&rn_in_abs);
|
1063 |
|
1064 |
//Set the DAP denominator to 1e108.
|
1065 |
GMP_INTS_mpz_set_general_int(&dap_denominator,
|
1066 |
&error_flag,
|
1067 |
"1e108");
|
1068 |
|
1069 |
//By convention, we will not mess with anything with an
|
1070 |
//absolute value greater than HMAX/1. If such a condition exists, puke out.
|
1071 |
//Form up the value of HMAX/1 if HMAX was specified.
|
1072 |
if (parblock->hmax_specified)
|
1073 |
{
|
1074 |
GMP_INTS_mpz_copy(&(hmax_over_one.num), &(parblock->hmax));
|
1075 |
GMP_INTS_mpz_set_ui(&(hmax_over_one.den), 1);
|
1076 |
GMP_RATS_mpq_copy(&rn_in_abs, &(parblock->rn));
|
1077 |
GMP_INTS_mpz_abs(&(rn_in_abs.num));
|
1078 |
if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_one, NULL) > 0)
|
1079 |
{
|
1080 |
printf("The magnitude of the number supplied exceeds HMAX/1, and hence the number\n"
|
1081 |
"has no neighbors in the series of interest. Calculation cannot continue.\n");
|
1082 |
exit(4);
|
1083 |
}
|
1084 |
}
|
1085 |
|
1086 |
//If the "verbose" option is specified, we want to give the continued fraction
|
1087 |
//partial quotients and convergents of either the number to approximate,
|
1088 |
//its reciprocal, or none of the above, as appropriate; and give a bit more
|
1089 |
//information, in addition.
|
1090 |
if (parblock->argblock.verbose)
|
1091 |
{
|
1092 |
if (parblock->hmax_specified)
|
1093 |
{
|
1094 |
//Stuff HMAX/KMAX. This is necessary for comparison.
|
1095 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.num), &(parblock->hmax));
|
1096 |
GMP_INTS_mpz_copy(&(hmax_over_kmax.den), &(parblock->kmax));
|
1097 |
}
|
1098 |
|
1099 |
if (!(parblock->hmax_specified) || (GMP_RATS_mpq_cmp(&(parblock->rn), &hmax_over_kmax, NULL) < 0))
|
1100 |
{
|
1101 |
//Either HMAX was not specified or else we are below the corner point on the
|
1102 |
//integer lattice. Get the continued fraction representation of the number
|
1103 |
//rather than its reciprocal.
|
1104 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.num), &(rn_in_abs.den));
|
1105 |
|
1106 |
//Print out the continued fraction decomposition of the rational number.
|
1107 |
GMP_RALG_cfdecomp_emit(stdout,
|
1108 |
"CF Representation Of Absolute Value Of Rational Number Specified",
|
1109 |
&cf_decomp,
|
1110 |
0,
|
1111 |
1,
|
1112 |
&dap_denominator);
|
1113 |
|
1114 |
//Destroy the decomposition--free the memory.
|
1115 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
1116 |
}
|
1117 |
else if (GMP_RATS_mpq_cmp(&(rn_in_abs), &hmax_over_kmax, NULL) == 0)
|
1118 |
{
|
1119 |
//In this case, the rational number specified is exactly the same in
|
1120 |
//magnitude as HMAX/KMAX. I am inclined to suppress the CF decomp.
|
1121 |
printf("Rational number specified is HMAX/KMAX. CF decomp not provided.\n");
|
1122 |
}
|
1123 |
else
|
1124 |
{
|
1125 |
//The number specified is beyond the corner point. It is appropriate to
|
1126 |
//provide the decomposition of the reciprocal rather than of the number
|
1127 |
//itself.
|
1128 |
GMP_RALG_cfdecomp_init(&cf_decomp, &error_flag, &(rn_in_abs.den), &(rn_in_abs.num));
|
1129 |
|
1130 |
//Print out the continued fraction decomposition of the rational number.
|
1131 |
GMP_RALG_cfdecomp_emit(stdout,
|
1132 |
"CF Representation Of Reciprocal Of Absolute Value Of Rational Number Specified",
|
1133 |
&cf_decomp,
|
1134 |
0,
|
1135 |
1,
|
1136 |
&dap_denominator);
|
1137 |
|
1138 |
//Destroy the decomposition--free the memory.
|
1139 |
GMP_RALG_cfdecomp_destroy(&cf_decomp);
|
1140 |
}
|
1141 |
} //End if verbose.
|
1142 |
|
1143 |
//Do all the work to get the neighbors of the number passed.
|
1144 |
GMP_RALG_consecutive_fab_terms(
|
1145 |
&(parblock->rn),
|
1146 |
&(parblock->kmax),
|
1147 |
(parblock->hmax_specified) ? (&(parblock->hmax)) : (NULL),
|
1148 |
0,
|
1149 |
1,
|
1150 |
&neighbor_data);
|
1151 |
|
1152 |
//Print the neighbor data block for debugging.
|
1153 |
#if 0
|
1154 |
GMP_RALG_consecutive_fab_terms_result_dump(stdout, &neighbor_data);
|
1155 |
#endif
|
1156 |
|
1157 |
//Print the neighbor on the right, if it exists.
|
1158 |
if (neighbor_data.n_right_out)
|
1159 |
{
|
1160 |
CMAIN_print_app_in_std_form(stdout,
|
1161 |
-1,
|
1162 |
&(neighbor_data.rn_in),
|
1163 |
&(neighbor_data.rights[0].neighbor),
|
1164 |
parblock->argblock.noformat,
|
1165 |
0,
|
1166 |
0,
|
1167 |
&dap_denominator);
|
1168 |
}
|
1169 |
|
1170 |
//Deallocate all dynamic memory.
|
1171 |
GMP_RALG_consecutive_fab_terms_result_free(&neighbor_data);
|
1172 |
GMP_RATS_mpq_clear(&hmax_over_one);
|
1173 |
GMP_RATS_mpq_clear(&hmax_over_kmax);
|
1174 |
GMP_INTS_mpz_clear(&dap_denominator);
|
1175 |
GMP_RATS_mpq_clear(&rn_in_abs);
|
1176 |
|
1177 |
return(rv);
|
1178 |
}
|
1179 |
|
1180 |
|
1181 |
int c_main(int argc, char* argv[])
|
1182 |
{
|
1183 |
int rv=0;
|
1184 |
struct CfbrapabCmainStruct parblock;
|
1185 |
|
1186 |
if (argc==2 && !strcmp(argv[1], "-help"))
|
1187 |
{
|
1188 |
FCMIOF_hline();
|
1189 |
printf("DESCRIPTION\n");
|
1190 |
printf(" This utility calculates best rational approximations of the form h/k\n");
|
1191 |
printf(" under the constraint k <= KMAX (i.e. in the Farey series of order KMAX),\n");
|
1192 |
printf(" or under the constraints k <= KMAX and h <= HMAX (i.e. in a rectangular\n");
|
1193 |
printf(" region of the integer lattice). This utility uses continued fraction\n");
|
1194 |
printf(" algorithms presented in the accompanying book \"A Practitioner's Guide\n");
|
1195 |
printf(" ...\", and this book (a work in progress) should be consulted both as a\n");
|
1196 |
printf(" reference to the algorithms and a reference for this utility. All\n");
|
1197 |
printf(" rational numbers calculated are in lowest terms. This utility will\n");
|
1198 |
printf(" operate on negative numbers, but all results are produced by symmetry\n");
|
1199 |
printf(" (the continued fraction representation of negative numbers is NOT\n");
|
1200 |
printf(" calculated). The default operation of this utility is to calculated the\n");
|
1201 |
printf(" closest rational number in the series of interest. If the rational\n");
|
1202 |
printf(" number supplied is equidistant between two formable rational numbers in\n");
|
1203 |
printf(" the series of interest, the neighbor smaller in magnitude is returned. If\n");
|
1204 |
printf(" the rational number supplied is already formable, it is returned in lowest\n");
|
1205 |
printf(" terms. If the rational number supplied does not have neighbors (i.e. it\n");
|
1206 |
printf(" is larger than HMAX/1), an error is generated.\n");
|
1207 |
printf("\n");
|
1208 |
printf("USAGE\n");
|
1209 |
printf(" cfbrapab srn urn_kmax [options]\n");
|
1210 |
printf(" cfbrapab srn urn_kmax urn_hmax [options]\n");
|
1211 |
printf(" cfbrapab -help\n");
|
1212 |
printf("\n");
|
1213 |
printf("OPTIONS\n");
|
1214 |
printf(" -neversmaller, -neverlarger\n");
|
1215 |
printf(" The -neversmaller option will prohibit this utility from choosing a\n");
|
1216 |
printf(" rational approximation which is smaller than the rational number\n");
|
1217 |
printf(" supplied. Thus, this option will force the utility to choose the right\n");
|
1218 |
printf(" neighbor rather than the left, regardless of relative distance. The\n");
|
1219 |
printf(" behavior if the rational number supplied is formable under the \n");
|
1220 |
printf(" constraints is unchanged. The -neverlarger option is analogous.\n");
|
1221 |
printf(" These options cannot be used with -n, -pred, or -succ.\n");
|
1222 |
printf(" -pred, -succ\n");
|
1223 |
printf(" Will cause the utility to find the predecessor or successor in the\n");
|
1224 |
printf(" series of interest to the rational number supplied (in the event the\n");
|
1225 |
printf(" number supplied is already formable under the constraints). For\n");
|
1226 |
printf(" numbers not already formable under the constraints, the left or right\n");
|
1227 |
printf(" formable neighbor will be returned. Supplying a rational number that\n");
|
1228 |
printf(" does not have a predecessor or successor (i.e. < 0/1 or > HMAX/1) will\n");
|
1229 |
printf(" generate an error. These options cannot be used with -neversmaller,\n");
|
1230 |
printf(" -neverlarger, or -n.\n");
|
1231 |
CU_MSGS_std_options(stdout, PNAME);
|
1232 |
FCMIOF_hline();
|
1233 |
CU_MSGS_toolset_info_msg(stdout, PNAME);
|
1234 |
FCMIOF_hline();
|
1235 |
}
|
1236 |
else if (argc < 3)
|
1237 |
{
|
1238 |
CU_MSGS_too_few_args_msg(stdout, PNAME);
|
1239 |
rv = 4;
|
1240 |
goto ret_pt;
|
1241 |
}
|
1242 |
else
|
1243 |
{
|
1244 |
//In this branch, we must have an invocation of the form
|
1245 |
// cfbrapab SRN KMAX <options>
|
1246 |
//or
|
1247 |
// cfbrapab SRN KMAX HMAX <options>
|
1248 |
//
|
1249 |
//Call the function to collect all the command-line parameters.
|
1250 |
//This function takes care of error processing, as well. If there
|
1251 |
//is an error of any kind, the function will simply abort and
|
1252 |
//supply the right return error code of 4.
|
1253 |
process_command_line_args(&parblock,
|
1254 |
argc,
|
1255 |
argv);
|
1256 |
|
1257 |
//If the debug option was set, emit the debugging information.
|
1258 |
if (parblock.argblock.debug)
|
1259 |
{
|
1260 |
FCMIOF_hline();
|
1261 |
CU_MSGS_emit_vcinfo_from_ptr_table(stdout,C_MAIN_vcinfoptrs,PNAMEUC);
|
1262 |
}
|
1263 |
|
1264 |
//Emit the opening horizontal line iff the -nf option isn't set.
|
1265 |
if (!(parblock.argblock.noformat))
|
1266 |
FCMIOF_hline();
|
1267 |
|
1268 |
//Print out a major mode message to indicate what we are trying to do.
|
1269 |
if (!(parblock.argblock.noformat))
|
1270 |
{
|
1271 |
if (!parblock.neversmaller_specified && !parblock.neverlarger_specified && !parblock.pred_specified && !parblock.succ_specified)
|
1272 |
{
|
1273 |
printf("MAJOR MODE: Finding closest rational number(s) under the constraints.\n");
|
1274 |
}
|
1275 |
else if (parblock.neversmaller_specified)
|
1276 |
{
|
1277 |
printf("MAJOR MODE: Finding closest rational number with magnitude not smaller under\n the constraints.\n");
|
1278 |
}
|
1279 |
else if (parblock.neverlarger_specified)
|
1280 |
{
|
1281 |
printf("MAJOR MODE: Finding closest rational number with magnitude not larger under\n the constraints.\n");
|
1282 |
}
|
1283 |
else if (parblock.pred_specified)
|
1284 |
{
|
1285 |
printf("MAJOR MODE: Finding predecessor under the constraints.\n");
|
1286 |
}
|
1287 |
else if (parblock.succ_specified)
|
1288 |
{
|
1289 |
printf("MAJOR MODE: Finding successor under the constraints.\n");
|
1290 |
}
|
1291 |
else
|
1292 |
{
|
1293 |
assert(0);
|
1294 |
}
|
1295 |
|
1296 |
FCMIOF_hline();
|
1297 |
}
|
1298 |
|
1299 |
//Echo back the command-line parameters.
|
1300 |
if (!(parblock.argblock.noformat))
|
1301 |
{
|
1302 |
GMP_INTS_mpz_long_int_format_to_stream(stdout,
|
1303 |
&(parblock.rn.num),
|
1304 |
"RI_IN Numerator");
|
1305 |
FCMIOF_hline();
|
1306 |
GMP_INTS_mpz_long_int_format_to_stream(stdout,
|
1307 |
&(parblock.rn.den),
|
1308 |
"RI_IN Denominator");
|
1309 |
FCMIOF_hline();
|
1310 |
|
1311 |
GMP_INTS_mpz_long_int_format_to_stream(stdout,
|
1312 |
&(parblock.kmax),
|
1313 |
"K_MAX");
|
1314 |
FCMIOF_hline();
|
1315 |
|
1316 |
if (parblock.hmax_specified)
|
1317 |
{
|
1318 |
GMP_INTS_mpz_long_int_format_to_stream(stdout,
|
1319 |
&(parblock.hmax),
|
1320 |
"H_MAX");
|
1321 |
FCMIOF_hline();
|
1322 |
}
|
1323 |
|
1324 |
if (parblock.n_specified)
|
1325 |
{
|
1326 |
GMP_INTS_mpz_struct temp24;
|
1327 |
|
1328 |
GMP_INTS_mpz_init(&temp24);
|
1329 |
|
1330 |
GMP_INTS_mpz_set_ui(&temp24, parblock.n);
|
1331 |
|
1332 |
GMP_INTS_mpz_long_int_format_to_stream(stdout,
|
1333 |
&temp24,
|
1334 |
"Number Of Neighbors");
|
1335 |
|
1336 |
FCMIOF_hline();
|
1337 |
|
1338 |
GMP_INTS_mpz_clear(&temp24);
|
1339 |
}
|
1340 |
}
|
1341 |
|
1342 |
//We need to split now into distinct cases
|
1343 |
//depending on the command-line parameters. We will
|
1344 |
//then hack out solutions for each case.
|
1345 |
if (!parblock.pred_specified && !parblock.succ_specified && !parblock.n_specified)
|
1346 |
{
|
1347 |
//Classic closest neighbor case.
|
1348 |
rv = CMAIN_classic_closest_neighbor(&parblock);
|
1349 |
}
|
1350 |
else if (parblock.n_specified)
|
1351 |
{
|
1352 |
//Classic multiple neighbor case.
|
1353 |
rv = CMAIN_multiple_neighbor(&parblock);
|
1354 |
}
|
1355 |
else if (parblock.pred_specified)
|
1356 |
{
|
1357 |
rv = CMAIN_predecessor(&parblock);
|
1358 |
}
|
1359 |
else if (parblock.succ_specified)
|
1360 |
{
|
1361 |
rv = CMAIN_successor(&parblock);
|
1362 |
}
|
1363 |
else
|
1364 |
{
|
1365 |
assert(0);
|
1366 |
}
|
1367 |
|
1368 |
//Emit the closing horizontal line iff the -nf option isn't set.
|
1369 |
//if (!(parblock.argblock.noformat))
|
1370 |
// FCMIOF_hline();
|
1371 |
|
1372 |
//Release all dynamic memory.
|
1373 |
release_command_line_args(&parblock);
|
1374 |
}
|
1375 |
|
1376 |
ret_pt:
|
1377 |
return(rv);
|
1378 |
}
|
1379 |
|
1380 |
// End of c_main.c.
|