//$Header$
//-------------------------------------------------------------------------------------------------
//This file is part of "David T. Ashley's Shared Source Code", a set of shared components
//integrated into many of David T. Ashley's projects.
//-------------------------------------------------------------------------------------------------
//This source code and any program in which it is compiled/used is provided under the MIT License,
//reproduced below.
//-------------------------------------------------------------------------------------------------
//Permission is hereby granted, free of charge, to any person obtaining a copy of
//this software and associated documentation files(the "Software"), to deal in the
//Software without restriction, including without limitation the rights to use,
//copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
//Software, and to permit persons to whom the Software is furnished to do so,
//subject to the following conditions :
//
//The above copyright notice and this permission notice shall be included in all
//copies or substantial portions of the Software.
//
//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
//SOFTWARE.
//-------------------------------------------------------------------------------------------------
/* GNU multiple precision rational number algorithms
** (and some extras).
*/
#ifndef GMP_RALG_H_INCLUDED
#define GMP_RALG_H_INCLUDED
#ifdef MODULE_GMP_RALG
#define DECMOD_GMP_RALG
#else
#define DECMOD_GMP_RALG extern
#endif
typedef struct
{
int n;
//The number of continued fraction partial quotients
//and convergents included in this structure. Indices
//are only valid until one less than this value.
int nallocd;
//The number of slots allocated for partial quotients,
//convergents, etc. This goes in blocks of the size
//below to avoid numerous reallocation calls.
#define GMP_RALG_CF_ALLOC_INCREMENT (10)
GMP_INTS_mpz_struct num;
//The numerator of the rational number whose continued
//fraction decomposition is contained in this structure.
//It must be non-negative. The fraction is not
//required to be reduced.
GMP_INTS_mpz_struct den;
//The denominator. It must be positive.
GMP_INTS_mpz_struct *a;
//Pointer to array of arbitrary integers. Each
//is one of the partial quotients. This array
//is allocated out to nallocd items but only n
//are used.
GMP_INTS_mpz_struct *dividend;
//The dividend of each division round.
GMP_INTS_mpz_struct *divisor;
//The divisor at each round.
GMP_INTS_mpz_struct *remainder;
//The remainder of each round.
GMP_INTS_mpz_struct *p;
GMP_INTS_mpz_struct *q;
//The convergents.
} GMP_RALG_cf_app_struct;
//This structure type holds a single neighbor to a rational number we wish
//to approximate in either the Farey series of a certain order or
//in a rectangular region of the integer lattice. These neighbors
//are manipulated in potentially large arrays.
typedef struct {
int index;
//The relative index of the neighbor. 0 means that it is
//the reduced and normalized form of the number to be
//approximated (although at the present time the number to be
//approximated is never included among the neighbors),
//"1" means nearest neighbor to the right,
//"2" means second nearest neighbor to the right, "-1" means
//neighbor to the left, etc.
GMP_RATS_mpq_struct neighbor;
//The neighbor itself, fully reduced and normalized.
} GMP_RALG_fab_neighbor_struct;
//This structure type holds a collection of neighbors which were
//produced on demand to a rational number, either in a Farey series
//or in a rectangular region of the integer lattice. The collection
//contains the neighbors plus a great deal of other information.
typedef struct {
GMP_RATS_mpq_struct rn_in;
//The rational number which was supplied to have its neighbors
//found. This number is raw as supplied by the caller and
//not necessarily normalized.
GMP_INTS_mpz_struct kmax_in;
//The value of KMAX supplied.
int hmax_supplied;
//TRUE if the problem includes a bound on the numerator, i.e.
//rectangular region of integer lattice rather than Farey series.
GMP_INTS_mpz_struct hmax_in;
//The value of HMAX if it was supplied, or zero otherwise.
char *error;
//NULL if no error occured in finding the neighbors. If non-NULL,
//points to a descriptive string explaining the error.
//If this error flag is set, contents of this structure
//below this point should not be used--they are not
//necessarily valid.
GMP_RATS_mpq_struct corner_point;
//The normalized (reduced) version of the "corner point"
//(HMAX/KMAX). This will only be filled in if
//HMAX is specified, or it will be 0/1 otherwise.
GMP_RATS_mpq_struct corner_point_minus_one;
GMP_RATS_mpq_struct corner_point_plus_one;
//The points just before and just above the corner point
//in the rectangular region of interest.
GMP_RATS_mpq_struct hmax_over_one;
//If HMAX is supplied, the rational number which is
//HMAX/1. Used for comparisons if we are asked
//to approximate a number out of range.
GMP_RATS_mpq_struct norm_rn;
//The normalized form of the rational number supplied.
GMP_RATS_mpq_struct abs_norm_rn;
//The absolute value of the normalized form of the rational
//number supplied.
int rn_is_neg;
//Set TRUE if the rational number supplied was negative. This is
//not a problem, because symmetry is exploited. For the purposes
//of this structure, the Farey series extends on both sides of zero
//forever, and a rectangular region of the integer lattice extends
//from -HMAX/1 to HMAX/1.
int n_left_in;
int n_right_in;
//The number of left neighbors to the rational number supplied
//requested, and the number of right neighbors. Note that if
//one of the bounds of the series of interest is encountered,
//the actual number supplied may be less.
GMP_RALG_cf_app_struct rn_abs_cf_rep;
//The CF representation of the absolute value of the
//rational number supplied. This always has to be
//formed, because it is the only way to immediately
//get neighbors.
int rn_abs_cf_rep_formed;
//Boolean flag to indicate if above has been formed.
//This is necessary because otherwise don't have
//any obvious clues about whether to free this.
GMP_RALG_cf_app_struct rn_abs_recip_cf_rep;
//The CF representation of the reciprocal of the
//absolute value of the rational number supplied.
//If HMAX is supplied, this may need to be formed,
//because it will allow us to zero in on neighbors
//where HMAX is the dominant constraint.
int rn_abs_recip_cf_rep_formed;
//Boolean flag to indicate if above has been formed.
//This is necessary because otherwise don't have
//any obvious clues about whether to free this.
GMP_RALG_cf_app_struct corner_point_cf_rep;
//The CF representation of the corner point. This is
//necessary to find the predecessor to the corner point
//in F_KMAX. This will only be calculated if HMAX
//is specified.
int corner_point_cf_rep_formed;
//Boolean flag to indicate if above has been formed.
//This is necessary because otherwise don't have
//any obvious clues about whether to free this.
GMP_RALG_cf_app_struct corner_point_recip_cf_rep;
//The CF representation of the reciprocal of the
//corner point. This is necessary to find the
//successor to the corner point along the HMAX constaint
//in the rectangular region of the integer lattice.
//This will only be calculated if HMAX is specified.
int corner_point_recip_cf_rep_formed;
//Boolean flag to indicate if above has been formed.
//This is necessary because otherwise don't have
//any obvious clues about whether to free this.
int equality;
//TRUE if the number supplied was already in the series of interest,
//or FALSE otherwise.
int n_left_out;
int n_right_out;
//The actual number of neighbors supplied. These may be less
//than the number requested only if a boundary of the series of interest
//was encountered.
GMP_RALG_fab_neighbor_struct *lefts;
GMP_RALG_fab_neighbor_struct *rights;
//Arrays of neighbors in the series of interest. These are
//only stuffed as indicated by "n_left_out" and "n_right_out"--anything
//else will hit on dynamic memory that is not initialized. The
//left array begins -1 and proceeds leftward through increasingly
//negative indices. The right begins with 1 and
//proceeds through increasingly positive indices. Each neighbor is
//irreducible and normalized. The rational number to be approximated
//is NEVER included among the neighbors.
} GMP_RALG_fab_neighbor_collection_struct;
/******************************************************************/
/*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/
/******************************************************************/
//Decomposes a rational number into its partial quotients
//and convergents. This initializes the "decomp" structure
//by allocating space, etc. The numerator must be non-negative
//and the denominator must be positive, or else there will be
//a failure. To make the interface simpler, the destruction
//call must ALWAYS be made, even in the event of a failure.
//If a failure occurs, the rational number will be decomposed
//as if it were 0/1. This means that if the failure flag
//is missed, it will still seem to be a valid rational number.
//It is not required that the rational number be reduced.
//The two integers passed are not changed at all, but they
//are copied into the decomp structure without modification.
//After this operation, the destruction call must be made.
DECMOD_GMP_RALG void GMP_RALG_cfdecomp_init(
GMP_RALG_cf_app_struct *decomp,
int *failure,
GMP_INTS_mpz_struct *num,
GMP_INTS_mpz_struct *den);
//Destroys all the the data allocated.
DECMOD_GMP_RALG void GMP_RALG_cfdecomp_destroy(
GMP_RALG_cf_app_struct *decomp
);
/******************************************************************/
/*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/
/******************************************************************/
//Prints out a continued fraction decomposition. The banner parameter
//may be NULL, which indicates no banner. nf flags whether to
//print with "no format", i.e. unformatted strings, separated by
//newlines. dap indicates whether to include decimal approximation
//information. The dap_denominator only needs to be non-NULL if
//"dap" is selected.
DECMOD_GMP_RALG
void GMP_RALG_cfdecomp_emit(
FILE *s,
char *banner,
GMP_RALG_cf_app_struct *decomp,
int nf,
int dap,
const GMP_INTS_mpz_struct *dap_denominator);
/******************************************************************/
/*** FAREY SERIES PREDECESSOR AND SUCCESSOR FUNCTIONS ***********/
/******************************************************************/
//The two functions below compute Farey predecessors and
//successors. These are straight sequences of arithmetic
//operations--no checking for overflow is done. The behavior
//below zero will probably not be correct, so these should
//not be used to go below 0/1 or with operands below 0/1.
DECMOD_GMP_RALG
void GMP_RALG_farey_predecessor(
GMP_RATS_mpq_struct *result,
const GMP_RATS_mpq_struct *plus_two,
const GMP_RATS_mpq_struct *plus_one,
const GMP_INTS_mpz_struct *N);
DECMOD_GMP_RALG
void GMP_RALG_farey_successor(
GMP_RATS_mpq_struct *result,
const GMP_RATS_mpq_struct *minus_two,
const GMP_RATS_mpq_struct *minus_one,
const GMP_INTS_mpz_struct *N);
//Finds the enclosing Farey neighbors to an arbitrary positive,
//zero, or negative rational number which may or may not be in
//the Farey series of order N.
//
//INPUTS
// rn: Rational number whose neighbors to find. MUST
// be normalized, meaning reduced and sign adjusted
// so numerator contains the sign.
// N: Order of series. MUST be positive.
// cf_rep: The continued fraction representation of the
// non-negative version of the rational number
// whose neighbors to find. The rational number
// need not have been reduced when the continued
// fraction representation was formed, i.e. only
// partial quotients and convergents are used.
//OUTPUTS
// equality: Set TRUE if the rational number is in the Farey
// series, or FALSE otherwise. This affects the
// generation of neighbors. If the number is
// already in the Farey series its other neighbor
// is generated, i.e. the rational number of
// interest is never one of the neighbors returned.
// left: The left neighbor.
// right: The right neighbor.
DECMOD_GMP_RALG
void GMP_RALG_enclosing_farey_neighbors(
const GMP_RATS_mpq_struct *rn,
const GMP_INTS_mpz_struct *N,
const GMP_RALG_cf_app_struct *cf_rep,
int *equality,
GMP_RATS_mpq_struct *left,
GMP_RATS_mpq_struct *right);
//Provides up to 2**24 - 1 (in concept--we don't recommend this!) terms in either
//the Farey series or in a rectangular region of the integer lattice. The caller
//must provide the space to be used, and must free it afterwards--this function
//performs no dynamic allocation.
DECMOD_GMP_RALG
void GMP_RALG_consecutive_fab_terms(
const GMP_RATS_mpq_struct *rn_in,
//The rational number whose neighbors are to be found.
const GMP_INTS_mpz_struct *kmax,
//Maximum denominator.
const GMP_INTS_mpz_struct *hmax,
//Maximum numerator, or a NULL pointer or a non-null pointer but
//an arbitrary integer value of 0 if only the Farey series is being
//considered.
int n_left_in,
//The number of left neighbors requested. Must be non-negative.
int n_right_in,
//The number of right neighbors requested. Must be non-negative.
GMP_RALG_fab_neighbor_collection_struct *result
//The structure containing all results. Even if there is an error
//condition, this must be freed when done by the caller.
);
//Frees the dynamic memory consumed by an earlier call to generate the
//neighbor information.
DECMOD_GMP_RALG
void GMP_RALG_consecutive_fab_terms_result_free(
GMP_RALG_fab_neighbor_collection_struct *arg
//The structure whose members to free up. After this call,
//the structure must be considered invalid.
);
//Debugging function. Dumps the contents of the neighbor structure to the indicated
//stream for analysis.
DECMOD_GMP_RALG
void GMP_RALG_consecutive_fab_terms_result_dump(
FILE *s,
GMP_RALG_fab_neighbor_collection_struct *arg);
/******************************************************************/
/*** VERSION CONTROL REPORTING FUNCTIONS ************************/
/******************************************************************/
DECMOD_GMP_RALG const char *GMP_RALG_cvcinfo(void);
DECMOD_GMP_RALG const char *GMP_RALG_hvcinfo(void);
#define GMP_RALG_H_VERSION ("$Header$")
#endif
//End of gmp_ralg.h.