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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25