Parent Directory | Revision Log

Revision **56** -
(**show annotations**)
(**download**)

*Sat Oct 29 01:53:01 2016 UTC*
(7 years, 1 month ago)
by *dashley*

File MIME type: text/plain

File size: 17771 byte(s)

File MIME type: text/plain

File size: 17771 byte(s)

License and property (keyword) changes.

1 | //$Header$ |

2 | //------------------------------------------------------------------------------------------------- |

3 | //This file is part of "David T. Ashley's Shared Source Code", a set of shared components |

4 | //integrated into many of David T. Ashley's projects. |

5 | //------------------------------------------------------------------------------------------------- |

6 | //This source code and any program in which it is compiled/used is provided under the MIT License, |

7 | //reproduced below. |

8 | //------------------------------------------------------------------------------------------------- |

9 | //Permission is hereby granted, free of charge, to any person obtaining a copy of |

10 | //this software and associated documentation files(the "Software"), to deal in the |

11 | //Software without restriction, including without limitation the rights to use, |

12 | //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the |

13 | //Software, and to permit persons to whom the Software is furnished to do so, |

14 | //subject to the following conditions : |

15 | // |

16 | //The above copyright notice and this permission notice shall be included in all |

17 | //copies or substantial portions of the Software. |

18 | // |

19 | //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |

20 | //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |

21 | //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE |

22 | //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |

23 | //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |

24 | //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE |

25 | //SOFTWARE. |

26 | //------------------------------------------------------------------------------------------------- |

27 | /* GNU multiple precision rational number algorithms |

28 | ** (and some extras). |

29 | */ |

30 | |

31 | #ifndef GMP_RALG_H_INCLUDED |

32 | #define GMP_RALG_H_INCLUDED |

33 | |

34 | #ifdef MODULE_GMP_RALG |

35 | #define DECMOD_GMP_RALG |

36 | #else |

37 | #define DECMOD_GMP_RALG extern |

38 | #endif |

39 | |

40 | typedef struct |

41 | { |

42 | int n; |

43 | //The number of continued fraction partial quotients |

44 | //and convergents included in this structure. Indices |

45 | //are only valid until one less than this value. |

46 | int nallocd; |

47 | //The number of slots allocated for partial quotients, |

48 | //convergents, etc. This goes in blocks of the size |

49 | //below to avoid numerous reallocation calls. |

50 | #define GMP_RALG_CF_ALLOC_INCREMENT (10) |

51 | GMP_INTS_mpz_struct num; |

52 | //The numerator of the rational number whose continued |

53 | //fraction decomposition is contained in this structure. |

54 | //It must be non-negative. The fraction is not |

55 | //required to be reduced. |

56 | GMP_INTS_mpz_struct den; |

57 | //The denominator. It must be positive. |

58 | GMP_INTS_mpz_struct *a; |

59 | //Pointer to array of arbitrary integers. Each |

60 | //is one of the partial quotients. This array |

61 | //is allocated out to nallocd items but only n |

62 | //are used. |

63 | GMP_INTS_mpz_struct *dividend; |

64 | //The dividend of each division round. |

65 | GMP_INTS_mpz_struct *divisor; |

66 | //The divisor at each round. |

67 | GMP_INTS_mpz_struct *remainder; |

68 | //The remainder of each round. |

69 | GMP_INTS_mpz_struct *p; |

70 | GMP_INTS_mpz_struct *q; |

71 | //The convergents. |

72 | } GMP_RALG_cf_app_struct; |

73 | |

74 | |

75 | //This structure type holds a single neighbor to a rational number we wish |

76 | //to approximate in either the Farey series of a certain order or |

77 | //in a rectangular region of the integer lattice. These neighbors |

78 | //are manipulated in potentially large arrays. |

79 | typedef struct { |

80 | int index; |

81 | //The relative index of the neighbor. 0 means that it is |

82 | //the reduced and normalized form of the number to be |

83 | //approximated (although at the present time the number to be |

84 | //approximated is never included among the neighbors), |

85 | //"1" means nearest neighbor to the right, |

86 | //"2" means second nearest neighbor to the right, "-1" means |

87 | //neighbor to the left, etc. |

88 | GMP_RATS_mpq_struct neighbor; |

89 | //The neighbor itself, fully reduced and normalized. |

90 | } GMP_RALG_fab_neighbor_struct; |

91 | |

92 | //This structure type holds a collection of neighbors which were |

93 | //produced on demand to a rational number, either in a Farey series |

94 | //or in a rectangular region of the integer lattice. The collection |

95 | //contains the neighbors plus a great deal of other information. |

96 | typedef struct { |

97 | GMP_RATS_mpq_struct rn_in; |

98 | //The rational number which was supplied to have its neighbors |

99 | //found. This number is raw as supplied by the caller and |

100 | //not necessarily normalized. |

101 | GMP_INTS_mpz_struct kmax_in; |

102 | //The value of KMAX supplied. |

103 | int hmax_supplied; |

104 | //TRUE if the problem includes a bound on the numerator, i.e. |

105 | //rectangular region of integer lattice rather than Farey series. |

106 | GMP_INTS_mpz_struct hmax_in; |

107 | //The value of HMAX if it was supplied, or zero otherwise. |

108 | char *error; |

109 | //NULL if no error occured in finding the neighbors. If non-NULL, |

110 | //points to a descriptive string explaining the error. |

111 | //If this error flag is set, contents of this structure |

112 | //below this point should not be used--they are not |

113 | //necessarily valid. |

114 | GMP_RATS_mpq_struct corner_point; |

115 | //The normalized (reduced) version of the "corner point" |

116 | //(HMAX/KMAX). This will only be filled in if |

117 | //HMAX is specified, or it will be 0/1 otherwise. |

118 | GMP_RATS_mpq_struct corner_point_minus_one; |

119 | GMP_RATS_mpq_struct corner_point_plus_one; |

120 | //The points just before and just above the corner point |

121 | //in the rectangular region of interest. |

122 | GMP_RATS_mpq_struct hmax_over_one; |

123 | //If HMAX is supplied, the rational number which is |

124 | //HMAX/1. Used for comparisons if we are asked |

125 | //to approximate a number out of range. |

126 | GMP_RATS_mpq_struct norm_rn; |

127 | //The normalized form of the rational number supplied. |

128 | GMP_RATS_mpq_struct abs_norm_rn; |

129 | //The absolute value of the normalized form of the rational |

130 | //number supplied. |

131 | int rn_is_neg; |

132 | //Set TRUE if the rational number supplied was negative. This is |

133 | //not a problem, because symmetry is exploited. For the purposes |

134 | //of this structure, the Farey series extends on both sides of zero |

135 | //forever, and a rectangular region of the integer lattice extends |

136 | //from -HMAX/1 to HMAX/1. |

137 | int n_left_in; |

138 | int n_right_in; |

139 | //The number of left neighbors to the rational number supplied |

140 | //requested, and the number of right neighbors. Note that if |

141 | //one of the bounds of the series of interest is encountered, |

142 | //the actual number supplied may be less. |

143 | GMP_RALG_cf_app_struct rn_abs_cf_rep; |

144 | //The CF representation of the absolute value of the |

145 | //rational number supplied. This always has to be |

146 | //formed, because it is the only way to immediately |

147 | //get neighbors. |

148 | int rn_abs_cf_rep_formed; |

149 | //Boolean flag to indicate if above has been formed. |

150 | //This is necessary because otherwise don't have |

151 | //any obvious clues about whether to free this. |

152 | GMP_RALG_cf_app_struct rn_abs_recip_cf_rep; |

153 | //The CF representation of the reciprocal of the |

154 | //absolute value of the rational number supplied. |

155 | //If HMAX is supplied, this may need to be formed, |

156 | //because it will allow us to zero in on neighbors |

157 | //where HMAX is the dominant constraint. |

158 | int rn_abs_recip_cf_rep_formed; |

159 | //Boolean flag to indicate if above has been formed. |

160 | //This is necessary because otherwise don't have |

161 | //any obvious clues about whether to free this. |

162 | GMP_RALG_cf_app_struct corner_point_cf_rep; |

163 | //The CF representation of the corner point. This is |

164 | //necessary to find the predecessor to the corner point |

165 | //in F_KMAX. This will only be calculated if HMAX |

166 | //is specified. |

167 | int corner_point_cf_rep_formed; |

168 | //Boolean flag to indicate if above has been formed. |

169 | //This is necessary because otherwise don't have |

170 | //any obvious clues about whether to free this. |

171 | GMP_RALG_cf_app_struct corner_point_recip_cf_rep; |

172 | //The CF representation of the reciprocal of the |

173 | //corner point. This is necessary to find the |

174 | //successor to the corner point along the HMAX constaint |

175 | //in the rectangular region of the integer lattice. |

176 | //This will only be calculated if HMAX is specified. |

177 | int corner_point_recip_cf_rep_formed; |

178 | //Boolean flag to indicate if above has been formed. |

179 | //This is necessary because otherwise don't have |

180 | //any obvious clues about whether to free this. |

181 | int equality; |

182 | //TRUE if the number supplied was already in the series of interest, |

183 | //or FALSE otherwise. |

184 | int n_left_out; |

185 | int n_right_out; |

186 | //The actual number of neighbors supplied. These may be less |

187 | //than the number requested only if a boundary of the series of interest |

188 | //was encountered. |

189 | GMP_RALG_fab_neighbor_struct *lefts; |

190 | GMP_RALG_fab_neighbor_struct *rights; |

191 | //Arrays of neighbors in the series of interest. These are |

192 | //only stuffed as indicated by "n_left_out" and "n_right_out"--anything |

193 | //else will hit on dynamic memory that is not initialized. The |

194 | //left array begins -1 and proceeds leftward through increasingly |

195 | //negative indices. The right begins with 1 and |

196 | //proceeds through increasingly positive indices. Each neighbor is |

197 | //irreducible and normalized. The rational number to be approximated |

198 | //is NEVER included among the neighbors. |

199 | } GMP_RALG_fab_neighbor_collection_struct; |

200 | |

201 | /******************************************************************/ |

202 | /*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/ |

203 | /******************************************************************/ |

204 | //Decomposes a rational number into its partial quotients |

205 | //and convergents. This initializes the "decomp" structure |

206 | //by allocating space, etc. The numerator must be non-negative |

207 | //and the denominator must be positive, or else there will be |

208 | //a failure. To make the interface simpler, the destruction |

209 | //call must ALWAYS be made, even in the event of a failure. |

210 | //If a failure occurs, the rational number will be decomposed |

211 | //as if it were 0/1. This means that if the failure flag |

212 | //is missed, it will still seem to be a valid rational number. |

213 | //It is not required that the rational number be reduced. |

214 | //The two integers passed are not changed at all, but they |

215 | //are copied into the decomp structure without modification. |

216 | //After this operation, the destruction call must be made. |

217 | DECMOD_GMP_RALG void GMP_RALG_cfdecomp_init( |

218 | GMP_RALG_cf_app_struct *decomp, |

219 | int *failure, |

220 | GMP_INTS_mpz_struct *num, |

221 | GMP_INTS_mpz_struct *den); |

222 | //Destroys all the the data allocated. |

223 | DECMOD_GMP_RALG void GMP_RALG_cfdecomp_destroy( |

224 | GMP_RALG_cf_app_struct *decomp |

225 | ); |

226 | |

227 | /******************************************************************/ |

228 | /*** INITIALIZATION AND DESTRUCTION FUNCTIONS *******************/ |

229 | /******************************************************************/ |

230 | //Prints out a continued fraction decomposition. The banner parameter |

231 | //may be NULL, which indicates no banner. nf flags whether to |

232 | //print with "no format", i.e. unformatted strings, separated by |

233 | //newlines. dap indicates whether to include decimal approximation |

234 | //information. The dap_denominator only needs to be non-NULL if |

235 | //"dap" is selected. |

236 | DECMOD_GMP_RALG |

237 | void GMP_RALG_cfdecomp_emit( |

238 | FILE *s, |

239 | char *banner, |

240 | GMP_RALG_cf_app_struct *decomp, |

241 | int nf, |

242 | int dap, |

243 | const GMP_INTS_mpz_struct *dap_denominator); |

244 | |

245 | /******************************************************************/ |

246 | /*** FAREY SERIES PREDECESSOR AND SUCCESSOR FUNCTIONS ***********/ |

247 | /******************************************************************/ |

248 | //The two functions below compute Farey predecessors and |

249 | //successors. These are straight sequences of arithmetic |

250 | //operations--no checking for overflow is done. The behavior |

251 | //below zero will probably not be correct, so these should |

252 | //not be used to go below 0/1 or with operands below 0/1. |

253 | DECMOD_GMP_RALG |

254 | void GMP_RALG_farey_predecessor( |

255 | GMP_RATS_mpq_struct *result, |

256 | const GMP_RATS_mpq_struct *plus_two, |

257 | const GMP_RATS_mpq_struct *plus_one, |

258 | const GMP_INTS_mpz_struct *N); |

259 | DECMOD_GMP_RALG |

260 | void GMP_RALG_farey_successor( |

261 | GMP_RATS_mpq_struct *result, |

262 | const GMP_RATS_mpq_struct *minus_two, |

263 | const GMP_RATS_mpq_struct *minus_one, |

264 | const GMP_INTS_mpz_struct *N); |

265 | |

266 | //Finds the enclosing Farey neighbors to an arbitrary positive, |

267 | //zero, or negative rational number which may or may not be in |

268 | //the Farey series of order N. |

269 | // |

270 | //INPUTS |

271 | // rn: Rational number whose neighbors to find. MUST |

272 | // be normalized, meaning reduced and sign adjusted |

273 | // so numerator contains the sign. |

274 | // N: Order of series. MUST be positive. |

275 | // cf_rep: The continued fraction representation of the |

276 | // non-negative version of the rational number |

277 | // whose neighbors to find. The rational number |

278 | // need not have been reduced when the continued |

279 | // fraction representation was formed, i.e. only |

280 | // partial quotients and convergents are used. |

281 | //OUTPUTS |

282 | // equality: Set TRUE if the rational number is in the Farey |

283 | // series, or FALSE otherwise. This affects the |

284 | // generation of neighbors. If the number is |

285 | // already in the Farey series its other neighbor |

286 | // is generated, i.e. the rational number of |

287 | // interest is never one of the neighbors returned. |

288 | // left: The left neighbor. |

289 | // right: The right neighbor. |

290 | DECMOD_GMP_RALG |

291 | void GMP_RALG_enclosing_farey_neighbors( |

292 | const GMP_RATS_mpq_struct *rn, |

293 | const GMP_INTS_mpz_struct *N, |

294 | const GMP_RALG_cf_app_struct *cf_rep, |

295 | int *equality, |

296 | GMP_RATS_mpq_struct *left, |

297 | GMP_RATS_mpq_struct *right); |

298 | |

299 | //Provides up to 2**24 - 1 (in concept--we don't recommend this!) terms in either |

300 | //the Farey series or in a rectangular region of the integer lattice. The caller |

301 | //must provide the space to be used, and must free it afterwards--this function |

302 | //performs no dynamic allocation. |

303 | DECMOD_GMP_RALG |

304 | void GMP_RALG_consecutive_fab_terms( |

305 | const GMP_RATS_mpq_struct *rn_in, |

306 | //The rational number whose neighbors are to be found. |

307 | const GMP_INTS_mpz_struct *kmax, |

308 | //Maximum denominator. |

309 | const GMP_INTS_mpz_struct *hmax, |

310 | //Maximum numerator, or a NULL pointer or a non-null pointer but |

311 | //an arbitrary integer value of 0 if only the Farey series is being |

312 | //considered. |

313 | int n_left_in, |

314 | //The number of left neighbors requested. Must be non-negative. |

315 | int n_right_in, |

316 | //The number of right neighbors requested. Must be non-negative. |

317 | GMP_RALG_fab_neighbor_collection_struct *result |

318 | //The structure containing all results. Even if there is an error |

319 | //condition, this must be freed when done by the caller. |

320 | ); |

321 | |

322 | //Frees the dynamic memory consumed by an earlier call to generate the |

323 | //neighbor information. |

324 | DECMOD_GMP_RALG |

325 | void GMP_RALG_consecutive_fab_terms_result_free( |

326 | GMP_RALG_fab_neighbor_collection_struct *arg |

327 | //The structure whose members to free up. After this call, |

328 | //the structure must be considered invalid. |

329 | ); |

330 | |

331 | //Debugging function. Dumps the contents of the neighbor structure to the indicated |

332 | //stream for analysis. |

333 | DECMOD_GMP_RALG |

334 | void GMP_RALG_consecutive_fab_terms_result_dump( |

335 | FILE *s, |

336 | GMP_RALG_fab_neighbor_collection_struct *arg); |

337 | |

338 | /******************************************************************/ |

339 | /*** VERSION CONTROL REPORTING FUNCTIONS ************************/ |

340 | /******************************************************************/ |

341 | DECMOD_GMP_RALG const char *GMP_RALG_cvcinfo(void); |

342 | DECMOD_GMP_RALG const char *GMP_RALG_hvcinfo(void); |

343 | #define GMP_RALG_H_VERSION ("$Header$") |

344 | #endif |

345 | |

346 | //End of gmp_ralg.h. |

Name | Value |
---|---|

svn:keywords |
Header |

dashley@gmail.com | ViewVC Help |

Powered by ViewVC 1.1.25 |