Parent Directory | Revision Log

Revision **42** -
(**hide annotations**)
(**download**)

*Fri Oct 14 01:50:00 2016 UTC*
(7 years, 1 month ago)
by *dashley*

File MIME type: text/plain

File size: 147987 byte(s)

File MIME type: text/plain

File size: 147987 byte(s)

Move shared source code to commonize.

1 | dashley | 25 | // $Header: /cvsroot/esrg/sfesrg/esrgpcpj/shared/c_datd/gmp_ints.c,v 1.22 2002/01/27 15:18:44 dtashley Exp $ |

2 | |||

3 | //-------------------------------------------------------------------------------- | ||

4 | //Copyright 2001 David T. Ashley | ||

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

6 | //This source code and any program in which it is compiled/used is provided under the GNU GENERAL | ||

7 | //PUBLIC LICENSE, Version 3, full license text below. | ||

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

9 | // GNU GENERAL PUBLIC LICENSE | ||

10 | // Version 3, 29 June 2007 | ||

11 | // | ||

12 | // Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> | ||

13 | // Everyone is permitted to copy and distribute verbatim copies | ||

14 | // of this license document, but changing it is not allowed. | ||

15 | // | ||

16 | // Preamble | ||

17 | // | ||

18 | // The GNU General Public License is a free, copyleft license for | ||

19 | //software and other kinds of works. | ||

20 | // | ||

21 | // The licenses for most software and other practical works are designed | ||

22 | //to take away your freedom to share and change the works. By contrast, | ||

23 | //the GNU General Public License is intended to guarantee your freedom to | ||

24 | //share and change all versions of a program--to make sure it remains free | ||

25 | //software for all its users. We, the Free Software Foundation, use the | ||

26 | //GNU General Public License for most of our software; it applies also to | ||

27 | //any other work released this way by its authors. You can apply it to | ||

28 | //your programs, too. | ||

29 | // | ||

30 | // When we speak of free software, we are referring to freedom, not | ||

31 | //price. Our General Public Licenses are designed to make sure that you | ||

32 | //have the freedom to distribute copies of free software (and charge for | ||

33 | //them if you wish), that you receive source code or can get it if you | ||

34 | //want it, that you can change the software or use pieces of it in new | ||

35 | //free programs, and that you know you can do these things. | ||

36 | // | ||

37 | // To protect your rights, we need to prevent others from denying you | ||

38 | //these rights or asking you to surrender the rights. Therefore, you have | ||

39 | //certain responsibilities if you distribute copies of the software, or if | ||

40 | //you modify it: responsibilities to respect the freedom of others. | ||

41 | // | ||

42 | // For example, if you distribute copies of such a program, whether | ||

43 | //gratis or for a fee, you must pass on to the recipients the same | ||

44 | //freedoms that you received. You must make sure that they, too, receive | ||

45 | //or can get the source code. And you must show them these terms so they | ||

46 | //know their rights. | ||

47 | // | ||

48 | // Developers that use the GNU GPL protect your rights with two steps: | ||

49 | //(1) assert copyright on the software, and (2) offer you this License | ||

50 | //giving you legal permission to copy, distribute and/or modify it. | ||

51 | // | ||

52 | // For the developers' and authors' protection, the GPL clearly explains | ||

53 | //that there is no warranty for this free software. For both users' and | ||

54 | //authors' sake, the GPL requires that modified versions be marked as | ||

55 | //changed, so that their problems will not be attributed erroneously to | ||

56 | //authors of previous versions. | ||

57 | // | ||

58 | // Some devices are designed to deny users access to install or run | ||

59 | //modified versions of the software inside them, although the manufacturer | ||

60 | //can do so. This is fundamentally incompatible with the aim of | ||

61 | //protecting users' freedom to change the software. The systematic | ||

62 | //pattern of such abuse occurs in the area of products for individuals to | ||

63 | //use, which is precisely where it is most unacceptable. Therefore, we | ||

64 | //have designed this version of the GPL to prohibit the practice for those | ||

65 | //products. If such problems arise substantially in other domains, we | ||

66 | //stand ready to extend this provision to those domains in future versions | ||

67 | //of the GPL, as needed to protect the freedom of users. | ||

68 | // | ||

69 | // Finally, every program is threatened constantly by software patents. | ||

70 | //States should not allow patents to restrict development and use of | ||

71 | //software on general-purpose computers, but in those that do, we wish to | ||

72 | //avoid the special danger that patents applied to a free program could | ||

73 | //make it effectively proprietary. To prevent this, the GPL assures that | ||

74 | //patents cannot be used to render the program non-free. | ||

75 | // | ||

76 | // The precise terms and conditions for copying, distribution and | ||

77 | //modification follow. | ||

78 | // | ||

79 | // TERMS AND CONDITIONS | ||

80 | // | ||

81 | // 0. Definitions. | ||

82 | // | ||

83 | // "This License" refers to version 3 of the GNU General Public License. | ||

84 | // | ||

85 | // "Copyright" also means copyright-like laws that apply to other kinds of | ||

86 | //works, such as semiconductor masks. | ||

87 | // | ||

88 | // "The Program" refers to any copyrightable work licensed under this | ||

89 | //License. Each licensee is addressed as "you". "Licensees" and | ||

90 | //"recipients" may be individuals or organizations. | ||

91 | // | ||

92 | // To "modify" a work means to copy from or adapt all or part of the work | ||

93 | //in a fashion requiring copyright permission, other than the making of an | ||

94 | //exact copy. The resulting work is called a "modified version" of the | ||

95 | //earlier work or a work "based on" the earlier work. | ||

96 | // | ||

97 | // A "covered work" means either the unmodified Program or a work based | ||

98 | //on the Program. | ||

99 | // | ||

100 | // To "propagate" a work means to do anything with it that, without | ||

101 | //permission, would make you directly or secondarily liable for | ||

102 | //infringement under applicable copyright law, except executing it on a | ||

103 | //computer or modifying a private copy. Propagation includes copying, | ||

104 | //distribution (with or without modification), making available to the | ||

105 | //public, and in some countries other activities as well. | ||

106 | // | ||

107 | // To "convey" a work means any kind of propagation that enables other | ||

108 | //parties to make or receive copies. Mere interaction with a user through | ||

109 | //a computer network, with no transfer of a copy, is not conveying. | ||

110 | // | ||

111 | // An interactive user interface displays "Appropriate Legal Notices" | ||

112 | //to the extent that it includes a convenient and prominently visible | ||

113 | //feature that (1) displays an appropriate copyright notice, and (2) | ||

114 | //tells the user that there is no warranty for the work (except to the | ||

115 | //extent that warranties are provided), that licensees may convey the | ||

116 | //work under this License, and how to view a copy of this License. If | ||

117 | //the interface presents a list of user commands or options, such as a | ||

118 | //menu, a prominent item in the list meets this criterion. | ||

119 | // | ||

120 | // 1. Source Code. | ||

121 | // | ||

122 | // The "source code" for a work means the preferred form of the work | ||

123 | //for making modifications to it. "Object code" means any non-source | ||

124 | //form of a work. | ||

125 | // | ||

126 | // A "Standard Interface" means an interface that either is an official | ||

127 | //standard defined by a recognized standards body, or, in the case of | ||

128 | //interfaces specified for a particular programming language, one that | ||

129 | //is widely used among developers working in that language. | ||

130 | // | ||

131 | // The "System Libraries" of an executable work include anything, other | ||

132 | //than the work as a whole, that (a) is included in the normal form of | ||

133 | //packaging a Major Component, but which is not part of that Major | ||

134 | //Component, and (b) serves only to enable use of the work with that | ||

135 | //Major Component, or to implement a Standard Interface for which an | ||

136 | //implementation is available to the public in source code form. A | ||

137 | //"Major Component", in this context, means a major essential component | ||

138 | //(kernel, window system, and so on) of the specific operating system | ||

139 | //(if any) on which the executable work runs, or a compiler used to | ||

140 | //produce the work, or an object code interpreter used to run it. | ||

141 | // | ||

142 | // The "Corresponding Source" for a work in object code form means all | ||

143 | //the source code needed to generate, install, and (for an executable | ||

144 | //work) run the object code and to modify the work, including scripts to | ||

145 | //control those activities. However, it does not include the work's | ||

146 | //System Libraries, or general-purpose tools or generally available free | ||

147 | //programs which are used unmodified in performing those activities but | ||

148 | //which are not part of the work. For example, Corresponding Source | ||

149 | //includes interface definition files associated with source files for | ||

150 | //the work, and the source code for shared libraries and dynamically | ||

151 | //linked subprograms that the work is specifically designed to require, | ||

152 | //such as by intimate data communication or control flow between those | ||

153 | //subprograms and other parts of the work. | ||

154 | // | ||

155 | // The Corresponding Source need not include anything that users | ||

156 | //can regenerate automatically from other parts of the Corresponding | ||

157 | //Source. | ||

158 | // | ||

159 | // The Corresponding Source for a work in source code form is that | ||

160 | //same work. | ||

161 | // | ||

162 | // 2. Basic Permissions. | ||

163 | // | ||

164 | // All rights granted under this License are granted for the term of | ||

165 | //copyright on the Program, and are irrevocable provided the stated | ||

166 | //conditions are met. This License explicitly affirms your unlimited | ||

167 | //permission to run the unmodified Program. The output from running a | ||

168 | //covered work is covered by this License only if the output, given its | ||

169 | //content, constitutes a covered work. This License acknowledges your | ||

170 | //rights of fair use or other equivalent, as provided by copyright law. | ||

171 | // | ||

172 | // You may make, run and propagate covered works that you do not | ||

173 | //convey, without conditions so long as your license otherwise remains | ||

174 | //in force. You may convey covered works to others for the sole purpose | ||

175 | //of having them make modifications exclusively for you, or provide you | ||

176 | //with facilities for running those works, provided that you comply with | ||

177 | //the terms of this License in conveying all material for which you do | ||

178 | //not control copyright. Those thus making or running the covered works | ||

179 | //for you must do so exclusively on your behalf, under your direction | ||

180 | //and control, on terms that prohibit them from making any copies of | ||

181 | //your copyrighted material outside their relationship with you. | ||

182 | // | ||

183 | // Conveying under any other circumstances is permitted solely under | ||

184 | //the conditions stated below. Sublicensing is not allowed; section 10 | ||

185 | //makes it unnecessary. | ||

186 | // | ||

187 | // 3. Protecting Users' Legal Rights From Anti-Circumvention Law. | ||

188 | // | ||

189 | // No covered work shall be deemed part of an effective technological | ||

190 | //measure under any applicable law fulfilling obligations under article | ||

191 | //11 of the WIPO copyright treaty adopted on 20 December 1996, or | ||

192 | //similar laws prohibiting or restricting circumvention of such | ||

193 | //measures. | ||

194 | // | ||

195 | // When you convey a covered work, you waive any legal power to forbid | ||

196 | //circumvention of technological measures to the extent such circumvention | ||

197 | //is effected by exercising rights under this License with respect to | ||

198 | //the covered work, and you disclaim any intention to limit operation or | ||

199 | //modification of the work as a means of enforcing, against the work's | ||

200 | //users, your or third parties' legal rights to forbid circumvention of | ||

201 | //technological measures. | ||

202 | // | ||

203 | // 4. Conveying Verbatim Copies. | ||

204 | // | ||

205 | // You may convey verbatim copies of the Program's source code as you | ||

206 | //receive it, in any medium, provided that you conspicuously and | ||

207 | //appropriately publish on each copy an appropriate copyright notice; | ||

208 | //keep intact all notices stating that this License and any | ||

209 | //non-permissive terms added in accord with section 7 apply to the code; | ||

210 | //keep intact all notices of the absence of any warranty; and give all | ||

211 | //recipients a copy of this License along with the Program. | ||

212 | // | ||

213 | // You may charge any price or no price for each copy that you convey, | ||

214 | //and you may offer support or warranty protection for a fee. | ||

215 | // | ||

216 | // 5. Conveying Modified Source Versions. | ||

217 | // | ||

218 | // You may convey a work based on the Program, or the modifications to | ||

219 | //produce it from the Program, in the form of source code under the | ||

220 | //terms of section 4, provided that you also meet all of these conditions: | ||

221 | // | ||

222 | // a) The work must carry prominent notices stating that you modified | ||

223 | // it, and giving a relevant date. | ||

224 | // | ||

225 | // b) The work must carry prominent notices stating that it is | ||

226 | // released under this License and any conditions added under section | ||

227 | // 7. This requirement modifies the requirement in section 4 to | ||

228 | // "keep intact all notices". | ||

229 | // | ||

230 | // c) You must license the entire work, as a whole, under this | ||

231 | // License to anyone who comes into possession of a copy. This | ||

232 | // License will therefore apply, along with any applicable section 7 | ||

233 | // additional terms, to the whole of the work, and all its parts, | ||

234 | // regardless of how they are packaged. This License gives no | ||

235 | // permission to license the work in any other way, but it does not | ||

236 | // invalidate such permission if you have separately received it. | ||

237 | // | ||

238 | // d) If the work has interactive user interfaces, each must display | ||

239 | // Appropriate Legal Notices; however, if the Program has interactive | ||

240 | // interfaces that do not display Appropriate Legal Notices, your | ||

241 | // work need not make them do so. | ||

242 | // | ||

243 | // A compilation of a covered work with other separate and independent | ||

244 | //works, which are not by their nature extensions of the covered work, | ||

245 | //and which are not combined with it such as to form a larger program, | ||

246 | //in or on a volume of a storage or distribution medium, is called an | ||

247 | //"aggregate" if the compilation and its resulting copyright are not | ||

248 | //used to limit the access or legal rights of the compilation's users | ||

249 | //beyond what the individual works permit. Inclusion of a covered work | ||

250 | //in an aggregate does not cause this License to apply to the other | ||

251 | //parts of the aggregate. | ||

252 | // | ||

253 | // 6. Conveying Non-Source Forms. | ||

254 | // | ||

255 | // You may convey a covered work in object code form under the terms | ||

256 | //of sections 4 and 5, provided that you also convey the | ||

257 | //machine-readable Corresponding Source under the terms of this License, | ||

258 | //in one of these ways: | ||

259 | // | ||

260 | // a) Convey the object code in, or embodied in, a physical product | ||

261 | // (including a physical distribution medium), accompanied by the | ||

262 | // Corresponding Source fixed on a durable physical medium | ||

263 | // customarily used for software interchange. | ||

264 | // | ||

265 | // b) Convey the object code in, or embodied in, a physical product | ||

266 | // (including a physical distribution medium), accompanied by a | ||

267 | // written offer, valid for at least three years and valid for as | ||

268 | // long as you offer spare parts or customer support for that product | ||

269 | // model, to give anyone who possesses the object code either (1) a | ||

270 | // copy of the Corresponding Source for all the software in the | ||

271 | // product that is covered by this License, on a durable physical | ||

272 | // medium customarily used for software interchange, for a price no | ||

273 | // more than your reasonable cost of physically performing this | ||

274 | // conveying of source, or (2) access to copy the | ||

275 | // Corresponding Source from a network server at no charge. | ||

276 | // | ||

277 | // c) Convey individual copies of the object code with a copy of the | ||

278 | // written offer to provide the Corresponding Source. This | ||

279 | // alternative is allowed only occasionally and noncommercially, and | ||

280 | // only if you received the object code with such an offer, in accord | ||

281 | // with subsection 6b. | ||

282 | // | ||

283 | // d) Convey the object code by offering access from a designated | ||

284 | // place (gratis or for a charge), and offer equivalent access to the | ||

285 | // Corresponding Source in the same way through the same place at no | ||

286 | // further charge. You need not require recipients to copy the | ||

287 | // Corresponding Source along with the object code. If the place to | ||

288 | // copy the object code is a network server, the Corresponding Source | ||

289 | // may be on a different server (operated by you or a third party) | ||

290 | // that supports equivalent copying facilities, provided you maintain | ||

291 | // clear directions next to the object code saying where to find the | ||

292 | // Corresponding Source. Regardless of what server hosts the | ||

293 | // Corresponding Source, you remain obligated to ensure that it is | ||

294 | // available for as long as needed to satisfy these requirements. | ||

295 | // | ||

296 | // e) Convey the object code using peer-to-peer transmission, provided | ||

297 | // you inform other peers where the object code and Corresponding | ||

298 | // Source of the work are being offered to the general public at no | ||

299 | // charge under subsection 6d. | ||

300 | // | ||

301 | // A separable portion of the object code, whose source code is excluded | ||

302 | //from the Corresponding Source as a System Library, need not be | ||

303 | //included in conveying the object code work. | ||

304 | // | ||

305 | // A "User Product" is either (1) a "consumer product", which means any | ||

306 | //tangible personal property which is normally used for personal, family, | ||

307 | //or household purposes, or (2) anything designed or sold for incorporation | ||

308 | //into a dwelling. In determining whether a product is a consumer product, | ||

309 | //doubtful cases shall be resolved in favor of coverage. For a particular | ||

310 | //product received by a particular user, "normally used" refers to a | ||

311 | //typical or common use of that class of product, regardless of the status | ||

312 | //of the particular user or of the way in which the particular user | ||

313 | //actually uses, or expects or is expected to use, the product. A product | ||

314 | //is a consumer product regardless of whether the product has substantial | ||

315 | //commercial, industrial or non-consumer uses, unless such uses represent | ||

316 | //the only significant mode of use of the product. | ||

317 | // | ||

318 | // "Installation Information" for a User Product means any methods, | ||

319 | //procedures, authorization keys, or other information required to install | ||

320 | //and execute modified versions of a covered work in that User Product from | ||

321 | //a modified version of its Corresponding Source. The information must | ||

322 | //suffice to ensure that the continued functioning of the modified object | ||

323 | //code is in no case prevented or interfered with solely because | ||

324 | //modification has been made. | ||

325 | // | ||

326 | // If you convey an object code work under this section in, or with, or | ||

327 | //specifically for use in, a User Product, and the conveying occurs as | ||

328 | //part of a transaction in which the right of possession and use of the | ||

329 | //User Product is transferred to the recipient in perpetuity or for a | ||

330 | //fixed term (regardless of how the transaction is characterized), the | ||

331 | //Corresponding Source conveyed under this section must be accompanied | ||

332 | //by the Installation Information. But this requirement does not apply | ||

333 | //if neither you nor any third party retains the ability to install | ||

334 | //modified object code on the User Product (for example, the work has | ||

335 | //been installed in ROM). | ||

336 | // | ||

337 | // The requirement to provide Installation Information does not include a | ||

338 | //requirement to continue to provide support service, warranty, or updates | ||

339 | //for a work that has been modified or installed by the recipient, or for | ||

340 | //the User Product in which it has been modified or installed. Access to a | ||

341 | //network may be denied when the modification itself materially and | ||

342 | //adversely affects the operation of the network or violates the rules and | ||

343 | //protocols for communication across the network. | ||

344 | // | ||

345 | // Corresponding Source conveyed, and Installation Information provided, | ||

346 | //in accord with this section must be in a format that is publicly | ||

347 | //documented (and with an implementation available to the public in | ||

348 | //source code form), and must require no special password or key for | ||

349 | //unpacking, reading or copying. | ||

350 | // | ||

351 | // 7. Additional Terms. | ||

352 | // | ||

353 | // "Additional permissions" are terms that supplement the terms of this | ||

354 | //License by making exceptions from one or more of its conditions. | ||

355 | //Additional permissions that are applicable to the entire Program shall | ||

356 | //be treated as though they were included in this License, to the extent | ||

357 | //that they are valid under applicable law. If additional permissions | ||

358 | //apply only to part of the Program, that part may be used separately | ||

359 | //under those permissions, but the entire Program remains governed by | ||

360 | //this License without regard to the additional permissions. | ||

361 | // | ||

362 | // When you convey a copy of a covered work, you may at your option | ||

363 | //remove any additional permissions from that copy, or from any part of | ||

364 | //it. (Additional permissions may be written to require their own | ||

365 | //removal in certain cases when you modify the work.) You may place | ||

366 | //additional permissions on material, added by you to a covered work, | ||

367 | //for which you have or can give appropriate copyright permission. | ||

368 | // | ||

369 | // Notwithstanding any other provision of this License, for material you | ||

370 | //add to a covered work, you may (if authorized by the copyright holders of | ||

371 | //that material) supplement the terms of this License with terms: | ||

372 | // | ||

373 | // a) Disclaiming warranty or limiting liability differently from the | ||

374 | // terms of sections 15 and 16 of this License; or | ||

375 | // | ||

376 | // b) Requiring preservation of specified reasonable legal notices or | ||

377 | // author attributions in that material or in the Appropriate Legal | ||

378 | // Notices displayed by works containing it; or | ||

379 | // | ||

380 | // c) Prohibiting misrepresentation of the origin of that material, or | ||

381 | // requiring that modified versions of such material be marked in | ||

382 | // reasonable ways as different from the original version; or | ||

383 | // | ||

384 | // d) Limiting the use for publicity purposes of names of licensors or | ||

385 | // authors of the material; or | ||

386 | // | ||

387 | // e) Declining to grant rights under trademark law for use of some | ||

388 | // trade names, trademarks, or service marks; or | ||

389 | // | ||

390 | // f) Requiring indemnification of licensors and authors of that | ||

391 | // material by anyone who conveys the material (or modified versions of | ||

392 | // it) with contractual assumptions of liability to the recipient, for | ||

393 | // any liability that these contractual assumptions directly impose on | ||

394 | // those licensors and authors. | ||

395 | // | ||

396 | // All other non-permissive additional terms are considered "further | ||

397 | //restrictions" within the meaning of section 10. If the Program as you | ||

398 | //received it, or any part of it, contains a notice stating that it is | ||

399 | //governed by this License along with a term that is a further | ||

400 | //restriction, you may remove that term. If a license document contains | ||

401 | //a further restriction but permits relicensing or conveying under this | ||

402 | //License, you may add to a covered work material governed by the terms | ||

403 | //of that license document, provided that the further restriction does | ||

404 | //not survive such relicensing or conveying. | ||

405 | // | ||

406 | // If you add terms to a covered work in accord with this section, you | ||

407 | //must place, in the relevant source files, a statement of the | ||

408 | //additional terms that apply to those files, or a notice indicating | ||

409 | //where to find the applicable terms. | ||

410 | // | ||

411 | // Additional terms, permissive or non-permissive, may be stated in the | ||

412 | //form of a separately written license, or stated as exceptions; | ||

413 | //the above requirements apply either way. | ||

414 | // | ||

415 | // 8. Termination. | ||

416 | // | ||

417 | // You may not propagate or modify a covered work except as expressly | ||

418 | //provided under this License. Any attempt otherwise to propagate or | ||

419 | //modify it is void, and will automatically terminate your rights under | ||

420 | //this License (including any patent licenses granted under the third | ||

421 | //paragraph of section 11). | ||

422 | // | ||

423 | // However, if you cease all violation of this License, then your | ||

424 | //license from a particular copyright holder is reinstated (a) | ||

425 | //provisionally, unless and until the copyright holder explicitly and | ||

426 | //finally terminates your license, and (b) permanently, if the copyright | ||

427 | //holder fails to notify you of the violation by some reasonable means | ||

428 | //prior to 60 days after the cessation. | ||

429 | // | ||

430 | // Moreover, your license from a particular copyright holder is | ||

431 | //reinstated permanently if the copyright holder notifies you of the | ||

432 | //violation by some reasonable means, this is the first time you have | ||

433 | //received notice of violation of this License (for any work) from that | ||

434 | //copyright holder, and you cure the violation prior to 30 days after | ||

435 | //your receipt of the notice. | ||

436 | // | ||

437 | // Termination of your rights under this section does not terminate the | ||

438 | //licenses of parties who have received copies or rights from you under | ||

439 | //this License. If your rights have been terminated and not permanently | ||

440 | //reinstated, you do not qualify to receive new licenses for the same | ||

441 | //material under section 10. | ||

442 | // | ||

443 | // 9. Acceptance Not Required for Having Copies. | ||

444 | // | ||

445 | // You are not required to accept this License in order to receive or | ||

446 | //run a copy of the Program. Ancillary propagation of a covered work | ||

447 | //occurring solely as a consequence of using peer-to-peer transmission | ||

448 | //to receive a copy likewise does not require acceptance. However, | ||

449 | //nothing other than this License grants you permission to propagate or | ||

450 | //modify any covered work. These actions infringe copyright if you do | ||

451 | //not accept this License. Therefore, by modifying or propagating a | ||

452 | //covered work, you indicate your acceptance of this License to do so. | ||

453 | // | ||

454 | // 10. Automatic Licensing of Downstream Recipients. | ||

455 | // | ||

456 | // Each time you convey a covered work, the recipient automatically | ||

457 | //receives a license from the original licensors, to run, modify and | ||

458 | //propagate that work, subject to this License. You are not responsible | ||

459 | //for enforcing compliance by third parties with this License. | ||

460 | // | ||

461 | // An "entity transaction" is a transaction transferring control of an | ||

462 | //organization, or substantially all assets of one, or subdividing an | ||

463 | //organization, or merging organizations. If propagation of a covered | ||

464 | //work results from an entity transaction, each party to that | ||

465 | //transaction who receives a copy of the work also receives whatever | ||

466 | //licenses to the work the party's predecessor in interest had or could | ||

467 | //give under the previous paragraph, plus a right to possession of the | ||

468 | //Corresponding Source of the work from the predecessor in interest, if | ||

469 | //the predecessor has it or can get it with reasonable efforts. | ||

470 | // | ||

471 | // You may not impose any further restrictions on the exercise of the | ||

472 | //rights granted or affirmed under this License. For example, you may | ||

473 | //not impose a license fee, royalty, or other charge for exercise of | ||

474 | //rights granted under this License, and you may not initiate litigation | ||

475 | //(including a cross-claim or counterclaim in a lawsuit) alleging that | ||

476 | //any patent claim is infringed by making, using, selling, offering for | ||

477 | //sale, or importing the Program or any portion of it. | ||

478 | // | ||

479 | // 11. Patents. | ||

480 | // | ||

481 | // A "contributor" is a copyright holder who authorizes use under this | ||

482 | //License of the Program or a work on which the Program is based. The | ||

483 | //work thus licensed is called the contributor's "contributor version". | ||

484 | // | ||

485 | // A contributor's "essential patent claims" are all patent claims | ||

486 | //owned or controlled by the contributor, whether already acquired or | ||

487 | //hereafter acquired, that would be infringed by some manner, permitted | ||

488 | //by this License, of making, using, or selling its contributor version, | ||

489 | //but do not include claims that would be infringed only as a | ||

490 | //consequence of further modification of the contributor version. For | ||

491 | //purposes of this definition, "control" includes the right to grant | ||

492 | //patent sublicenses in a manner consistent with the requirements of | ||

493 | //this License. | ||

494 | // | ||

495 | // Each contributor grants you a non-exclusive, worldwide, royalty-free | ||

496 | //patent license under the contributor's essential patent claims, to | ||

497 | //make, use, sell, offer for sale, import and otherwise run, modify and | ||

498 | //propagate the contents of its contributor version. | ||

499 | // | ||

500 | // In the following three paragraphs, a "patent license" is any express | ||

501 | //agreement or commitment, however denominated, not to enforce a patent | ||

502 | //(such as an express permission to practice a patent or covenant not to | ||

503 | //sue for patent infringement). To "grant" such a patent license to a | ||

504 | //party means to make such an agreement or commitment not to enforce a | ||

505 | //patent against the party. | ||

506 | // | ||

507 | // If you convey a covered work, knowingly relying on a patent license, | ||

508 | //and the Corresponding Source of the work is not available for anyone | ||

509 | //to copy, free of charge and under the terms of this License, through a | ||

510 | //publicly available network server or other readily accessible means, | ||

511 | //then you must either (1) cause the Corresponding Source to be so | ||

512 | //available, or (2) arrange to deprive yourself of the benefit of the | ||

513 | //patent license for this particular work, or (3) arrange, in a manner | ||

514 | //consistent with the requirements of this License, to extend the patent | ||

515 | //license to downstream recipients. "Knowingly relying" means you have | ||

516 | //actual knowledge that, but for the patent license, your conveying the | ||

517 | //covered work in a country, or your recipient's use of the covered work | ||

518 | //in a country, would infringe one or more identifiable patents in that | ||

519 | //country that you have reason to believe are valid. | ||

520 | // | ||

521 | // If, pursuant to or in connection with a single transaction or | ||

522 | //arrangement, you convey, or propagate by procuring conveyance of, a | ||

523 | //covered work, and grant a patent license to some of the parties | ||

524 | //receiving the covered work authorizing them to use, propagate, modify | ||

525 | //or convey a specific copy of the covered work, then the patent license | ||

526 | //you grant is automatically extended to all recipients of the covered | ||

527 | //work and works based on it. | ||

528 | // | ||

529 | // A patent license is "discriminatory" if it does not include within | ||

530 | //the scope of its coverage, prohibits the exercise of, or is | ||

531 | //conditioned on the non-exercise of one or more of the rights that are | ||

532 | //specifically granted under this License. You may not convey a covered | ||

533 | //work if you are a party to an arrangement with a third party that is | ||

534 | //in the business of distributing software, under which you make payment | ||

535 | //to the third party based on the extent of your activity of conveying | ||

536 | //the work, and under which the third party grants, to any of the | ||

537 | //parties who would receive the covered work from you, a discriminatory | ||

538 | //patent license (a) in connection with copies of the covered work | ||

539 | //conveyed by you (or copies made from those copies), or (b) primarily | ||

540 | //for and in connection with specific products or compilations that | ||

541 | //contain the covered work, unless you entered into that arrangement, | ||

542 | //or that patent license was granted, prior to 28 March 2007. | ||

543 | // | ||

544 | // Nothing in this License shall be construed as excluding or limiting | ||

545 | //any implied license or other defenses to infringement that may | ||

546 | //otherwise be available to you under applicable patent law. | ||

547 | // | ||

548 | // 12. No Surrender of Others' Freedom. | ||

549 | // | ||

550 | // If conditions are imposed on you (whether by court order, agreement or | ||

551 | //otherwise) that contradict the conditions of this License, they do not | ||

552 | //excuse you from the conditions of this License. If you cannot convey a | ||

553 | //covered work so as to satisfy simultaneously your obligations under this | ||

554 | //License and any other pertinent obligations, then as a consequence you may | ||

555 | //not convey it at all. For example, if you agree to terms that obligate you | ||

556 | //to collect a royalty for further conveying from those to whom you convey | ||

557 | //the Program, the only way you could satisfy both those terms and this | ||

558 | //License would be to refrain entirely from conveying the Program. | ||

559 | // | ||

560 | // 13. Use with the GNU Affero General Public License. | ||

561 | // | ||

562 | // Notwithstanding any other provision of this License, you have | ||

563 | //permission to link or combine any covered work with a work licensed | ||

564 | //under version 3 of the GNU Affero General Public License into a single | ||

565 | //combined work, and to convey the resulting work. The terms of this | ||

566 | //License will continue to apply to the part which is the covered work, | ||

567 | //but the special requirements of the GNU Affero General Public License, | ||

568 | //section 13, concerning interaction through a network will apply to the | ||

569 | //combination as such. | ||

570 | // | ||

571 | // 14. Revised Versions of this License. | ||

572 | // | ||

573 | // The Free Software Foundation may publish revised and/or new versions of | ||

574 | //the GNU General Public License from time to time. Such new versions will | ||

575 | //be similar in spirit to the present version, but may differ in detail to | ||

576 | //address new problems or concerns. | ||

577 | // | ||

578 | // Each version is given a distinguishing version number. If the | ||

579 | //Program specifies that a certain numbered version of the GNU General | ||

580 | //Public License "or any later version" applies to it, you have the | ||

581 | //option of following the terms and conditions either of that numbered | ||

582 | //version or of any later version published by the Free Software | ||

583 | //Foundation. If the Program does not specify a version number of the | ||

584 | //GNU General Public License, you may choose any version ever published | ||

585 | //by the Free Software Foundation. | ||

586 | // | ||

587 | // If the Program specifies that a proxy can decide which future | ||

588 | //versions of the GNU General Public License can be used, that proxy's | ||

589 | //public statement of acceptance of a version permanently authorizes you | ||

590 | //to choose that version for the Program. | ||

591 | // | ||

592 | // Later license versions may give you additional or different | ||

593 | //permissions. However, no additional obligations are imposed on any | ||

594 | //author or copyright holder as a result of your choosing to follow a | ||

595 | //later version. | ||

596 | // | ||

597 | // 15. Disclaimer of Warranty. | ||

598 | // | ||

599 | // THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY | ||

600 | //APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT | ||

601 | //HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY | ||

602 | //OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, | ||

603 | //THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR | ||

604 | //PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM | ||

605 | //IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF | ||

606 | //ALL NECESSARY SERVICING, REPAIR OR CORRECTION. | ||

607 | // | ||

608 | // 16. Limitation of Liability. | ||

609 | // | ||

610 | // IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING | ||

611 | //WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS | ||

612 | //THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY | ||

613 | //GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE | ||

614 | //USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF | ||

615 | //DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD | ||

616 | //PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), | ||

617 | //EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF | ||

618 | //SUCH DAMAGES. | ||

619 | // | ||

620 | // 17. Interpretation of Sections 15 and 16. | ||

621 | // | ||

622 | // If the disclaimer of warranty and limitation of liability provided | ||

623 | //above cannot be given local legal effect according to their terms, | ||

624 | //reviewing courts shall apply local law that most closely approximates | ||

625 | //an absolute waiver of all civil liability in connection with the | ||

626 | //Program, unless a warranty or assumption of liability accompanies a | ||

627 | //copy of the Program in return for a fee. | ||

628 | // | ||

629 | // END OF TERMS AND CONDITIONS | ||

630 | // | ||

631 | // How to Apply These Terms to Your New Programs | ||

632 | // | ||

633 | // If you develop a new program, and you want it to be of the greatest | ||

634 | //possible use to the public, the best way to achieve this is to make it | ||

635 | //free software which everyone can redistribute and change under these terms. | ||

636 | // | ||

637 | // To do so, attach the following notices to the program. It is safest | ||

638 | //to attach them to the start of each source file to most effectively | ||

639 | //state the exclusion of warranty; and each file should have at least | ||

640 | //the "copyright" line and a pointer to where the full notice is found. | ||

641 | // | ||

642 | // <one line to give the program's name and a brief idea of what it does.> | ||

643 | // Copyright (C) <year> <name of author> | ||

644 | // | ||

645 | // This program is free software: you can redistribute it and/or modify | ||

646 | // it under the terms of the GNU General Public License as published by | ||

647 | // the Free Software Foundation, either version 3 of the License, or | ||

648 | // (at your option) any later version. | ||

649 | // | ||

650 | // This program is distributed in the hope that it will be useful, | ||

651 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | ||

652 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||

653 | // GNU General Public License for more details. | ||

654 | // | ||

655 | // You should have received a copy of the GNU General Public License | ||

656 | // along with this program. If not, see <http://www.gnu.org/licenses/>. | ||

657 | // | ||

658 | //Also add information on how to contact you by electronic and paper mail. | ||

659 | // | ||

660 | // If the program does terminal interaction, make it output a short | ||

661 | //notice like this when it starts in an interactive mode: | ||

662 | // | ||

663 | // <program> Copyright (C) <year> <name of author> | ||

664 | // This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. | ||

665 | // This is free software, and you are welcome to redistribute it | ||

666 | // under certain conditions; type `show c' for details. | ||

667 | // | ||

668 | //The hypothetical commands `show w' and `show c' should show the appropriate | ||

669 | //parts of the General Public License. Of course, your program's commands | ||

670 | //might be different; for a GUI interface, you would use an "about box". | ||

671 | // | ||

672 | // You should also get your employer (if you work as a programmer) or school, | ||

673 | //if any, to sign a "copyright disclaimer" for the program, if necessary. | ||

674 | //For more information on this, and how to apply and follow the GNU GPL, see | ||

675 | //<http://www.gnu.org/licenses/>. | ||

676 | // | ||

677 | // The GNU General Public License does not permit incorporating your program | ||

678 | //into proprietary programs. If your program is a subroutine library, you | ||

679 | //may consider it more useful to permit linking proprietary applications with | ||

680 | //the library. If this is what you want to do, use the GNU Lesser General | ||

681 | //Public License instead of this License. But first, please read | ||

682 | //<http://www.gnu.org/philosophy/why-not-lgpl.html>. | ||

683 | //------------------------------------------------------------------------------------------------- | ||

684 | //-------------------------------------------------------------------------------- | ||

685 | #define MODULE_GMP_INTS | ||

686 | |||

687 | #include <assert.h> | ||

688 | #include <ctype.h> | ||

689 | #include <process.h> | ||

690 | #include <stddef.h> | ||

691 | #include <stdio.h> | ||

692 | #include <string.h> | ||

693 | #include <time.h> | ||

694 | |||

695 | |||

696 | /* Only included the guarded allocation header if we are compiling | ||

697 | ** a DOS console type application. Other types of applications have | ||

698 | ** other ways of protecting for out of memory. Including the | ||

699 | ** header would do no harm in these cases, but do no good, either. | ||

700 | */ | ||

701 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

702 | #include "ccmalloc.h" | ||

703 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

704 | #include "tclalloc.h" | ||

705 | #else | ||

706 | /* Do nothing. */ | ||

707 | #endif | ||

708 | |||

709 | #include "bstrfunc.h" | ||

710 | #include "charfunc.h" | ||

711 | #include "fcmiof.h" | ||

712 | #include "gmp_ints.h" | ||

713 | #include "intfunc.h" | ||

714 | |||

715 | |||

716 | /******************************************************************/ | ||

717 | /*** CUSTOM ALLOCATION FUNCTIONS *******************************/ | ||

718 | /******************************************************************/ | ||

719 | /* We need to decide here on how memory not on the stack will be | ||

720 | ** allocated (i.e. what will become of the standard functions | ||

721 | ** like malloc, free, etc.). This is dependent on the type of | ||

722 | ** application we're making. The possible types are so far are: | ||

723 | ** APP_TYPE_SIMPLE_DOS_CONSOLE : | ||

724 | ** Simple DOS console application. | ||

725 | ** APP_TYPE_IJUSCRIPTER_IJUCONSOLE: | ||

726 | ** The Tcl tool build. | ||

727 | ** | ||

728 | ** The custom allocation functions here are a "portal" or "wrapper" | ||

729 | ** for how the integer and rational number functions should | ||

730 | ** get memory. | ||

731 | ** | ||

732 | ** The functions below are standard, except that the GNU MP team | ||

733 | ** built more generality into what allocation schemes could be | ||

734 | ** used by including size information in some calls that don't | ||

735 | ** normally get it. That is why there are some extra calls below | ||

736 | ** where the information is discarded. Other than that, these are | ||

737 | ** standard allocation calls. | ||

738 | */ | ||

739 | //07/15/01: Visual inspection only. Function deemed too | ||

740 | //simple for unit testing. | ||

741 | void *GMP_INTS_malloc( size_t size ) | ||

742 | { | ||

743 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

744 | return(CCMALLOC_malloc(size)); | ||

745 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

746 | return(TclpAlloc(size)); | ||

747 | #else | ||

748 | return(malloc(size)); | ||

749 | #endif | ||

750 | } | ||

751 | |||

752 | |||

753 | //07/15/01: Visual inspection only. Function deemed too | ||

754 | //simple for unit testing. | ||

755 | void *GMP_INTS_calloc( size_t num, size_t size ) | ||

756 | { | ||

757 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

758 | return(CCMALLOC_calloc(num, size)); | ||

759 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

760 | return(TclpCalloc(num, size)); | ||

761 | #else | ||

762 | return(calloc(num, size)); | ||

763 | #endif | ||

764 | } | ||

765 | |||

766 | |||

767 | //07/15/01: Visual inspection only. Function deemed too | ||

768 | //simple for unit testing. | ||

769 | void *GMP_INTS_realloc( void *memblock, size_t size ) | ||

770 | { | ||

771 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

772 | return(CCMALLOC_realloc(memblock, size)); | ||

773 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

774 | return(TclpRealloc(memblock, size)); | ||

775 | #else | ||

776 | return(realloc(memblock, size)); | ||

777 | #endif | ||

778 | } | ||

779 | |||

780 | |||

781 | //07/15/01: Visual inspection only. Function deemed too | ||

782 | //simple for unit testing. | ||

783 | void *GMP_INTS_realloc_w_size( void *memblock, | ||

784 | size_t old_size, | ||

785 | size_t size ) | ||

786 | { | ||

787 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

788 | return(CCMALLOC_realloc(memblock, size)); | ||

789 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

790 | return(TclpRealloc(memblock, size)); | ||

791 | #else | ||

792 | return(realloc(memblock, size)); | ||

793 | #endif | ||

794 | } | ||

795 | |||

796 | |||

797 | //07/15/01: Visual inspection only. Function deemed too | ||

798 | //simple for unit testing. | ||

799 | void GMP_INTS_free( void *memblock ) | ||

800 | { | ||

801 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

802 | CCMALLOC_free(memblock); | ||

803 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

804 | TclpFree(memblock); | ||

805 | #else | ||

806 | free(memblock); | ||

807 | #endif | ||

808 | } | ||

809 | |||

810 | |||

811 | //07/15/01: Visual inspection only. Function deemed too | ||

812 | //simple for unit testing. | ||

813 | void GMP_INTS_free_w_size( void *memblock, size_t size ) | ||

814 | { | ||

815 | #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE) | ||

816 | CCMALLOC_free(memblock); | ||

817 | #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE) | ||

818 | TclpFree(memblock); | ||

819 | #else | ||

820 | free(memblock); | ||

821 | #endif | ||

822 | } | ||

823 | |||

824 | |||

825 | /******************************************************************/ | ||

826 | /*** PORTABILITY CHECK FUNCTIONS *******************************/ | ||

827 | /******************************************************************/ | ||

828 | //Because there is the risk that Microsoft Visual C++ might | ||

829 | //change in the future, the following function can be called | ||

830 | //to see if the assumptions about data sizes are valid. This | ||

831 | //function returns TRUE if there is a problem, or FALSE | ||

832 | //otherwise. | ||

833 | |||

834 | //07/15/01: Unit testing complete. | ||

835 | int GMP_INTS_data_sizes_are_wrong(void) | ||

836 | { | ||

837 | int i; | ||

838 | GMP_INTS_limb_t tv; | ||

839 | _int64 tv64; | ||

840 | |||

841 | //Check the number of bit rolls required to get the limb | ||

842 | //to go to zero again. This had better be 32. | ||

843 | tv = 1; | ||

844 | i = 0; | ||

845 | while (tv) | ||

846 | { | ||

847 | tv <<= 1; | ||

848 | i++; | ||

849 | } | ||

850 | if (i != 32) | ||

851 | return(1); | ||

852 | |||

853 | //Check that an _int64 is really and truly 64 bits. | ||

854 | tv64 = 1; | ||

855 | i = 0; | ||

856 | while (tv64) | ||

857 | { | ||

858 | tv64 <<= 1; | ||

859 | i++; | ||

860 | } | ||

861 | if (i != 64) | ||

862 | return(1); | ||

863 | |||

864 | //Room for additional tests here if needed later. | ||

865 | |||

866 | return(0); | ||

867 | } | ||

868 | |||

869 | |||

870 | /******************************************************************/ | ||

871 | /*** ERROR STRING IDENTIFICATION AND PROCESSING FUNCTIONS *******/ | ||

872 | /******************************************************************/ | ||

873 | |||

874 | int GMP_INTS_identify_nan_string(const char *s) | ||

875 | { | ||

876 | assert(s != NULL); | ||

877 | |||

878 | if (!strcmp(s, GMP_INTS_EF_INTOVF_POS_STRING)) | ||

879 | return(0); | ||

880 | else if (!strcmp(s, GMP_INTS_EF_INTOVF_NEG_STRING)) | ||

881 | return(1); | ||

882 | else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_POS_STRING)) | ||

883 | return(2); | ||

884 | else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING)) | ||

885 | return(3); | ||

886 | else | ||

887 | return(-1); | ||

888 | } | ||

889 | |||

890 | |||

891 | const char *GMP_INTS_supply_nan_string(int idx) | ||

892 | { | ||

893 | assert((idx >= 0) && (idx <= 3)); | ||

894 | |||

895 | if (idx==0) | ||

896 | return(GMP_INTS_EF_INTOVF_POS_STRING); | ||

897 | else if (idx==1) | ||

898 | return(GMP_INTS_EF_INTOVF_NEG_STRING); | ||

899 | else if (idx==2) | ||

900 | return(GMP_INTS_EF_INTOVF_TAINT_POS_STRING); | ||

901 | else | ||

902 | return(GMP_INTS_EF_INTOVF_TAINT_NEG_STRING); | ||

903 | } | ||

904 | |||

905 | |||

906 | /******************************************************************/ | ||

907 | /*** DEBUG PRINTING FUNCTIONS **********************************/ | ||

908 | /******************************************************************/ | ||

909 | //These functions are for printing out integers and limbs | ||

910 | //and groups of limbs for unit testing and debugging. | ||

911 | |||

912 | //07/15/01: Exempt from testing because debug/development | ||

913 | //function. | ||

914 | void GMP_INTS_print_limb_group(FILE *stream, | ||

915 | GMP_INTS_limb_srcptr lg, | ||

916 | GMP_INTS_size_t n, | ||

917 | char *desc) | ||

918 | { | ||

919 | int i; | ||

920 | |||

921 | assert(stream != NULL); | ||

922 | assert(n >= 0); | ||

923 | assert(desc != NULL); | ||

924 | |||

925 | if (!lg) | ||

926 | { | ||

927 | fprintf(stream, " %s: NULL\n", desc); | ||

928 | } | ||

929 | else | ||

930 | { | ||

931 | for (i=n-1; i>=0; i--) | ||

932 | { | ||

933 | fprintf(stream, " %s[%2d]: 0x%8X\n", desc, i, lg[i]); | ||

934 | } | ||

935 | } | ||

936 | } | ||

937 | |||

938 | |||

939 | void GMP_INTS_mpz_print_int(FILE *stream, | ||

940 | const GMP_INTS_mpz_struct *arg, | ||

941 | char *desc) | ||

942 | { | ||

943 | int i; | ||

944 | |||

945 | assert(stream != NULL); | ||

946 | assert(arg != NULL); | ||

947 | assert(desc != NULL); | ||

948 | |||

949 | fprintf(stream, "Printing integer:\n %s\n", desc); | ||

950 | |||

951 | fprintf(stream, " flags: %d\n", arg->flags); | ||

952 | fprintf(stream, " ptr value to body: %p\n", arg); | ||

953 | fprintf(stream, " n_allocd: %d\n", arg->n_allocd); | ||

954 | fprintf(stream, " size: %d\n", arg->size); | ||

955 | fprintf(stream, " limbs (ptr val): %p\n", arg->limbs); | ||

956 | |||

957 | for (i=arg->n_allocd-1; i>=0; i--) | ||

958 | { | ||

959 | fprintf(stream, " limb[%3d]: %8X\n", i, arg->limbs[i]); | ||

960 | } | ||

961 | } | ||

962 | |||

963 | |||

964 | /******************************************************************/ | ||

965 | /*** LOW-LEVEL MACRO REPLACEMENTS ******************************/ | ||

966 | /******************************************************************/ | ||

967 | //The functions in this category are replacements for macros. | ||

968 | //Clarity was gained at the expense of speed. | ||

969 | |||

970 | int GMP_INTS_mpz_get_flags (const GMP_INTS_mpz_struct *arg) | ||

971 | { | ||

972 | assert(arg != NULL); | ||

973 | assert(arg->n_allocd > 0); | ||

974 | |||

975 | return(arg->flags); | ||

976 | } | ||

977 | |||

978 | |||

979 | //07/15/01: Visual inspection only. Function deemed too | ||

980 | //simple for unit testing. | ||

981 | GMP_INTS_size_t GMP_INTS_abs_of_size_t(GMP_INTS_size_t arg) | ||

982 | { | ||

983 | //Be sure that the bit pattern does not represent the maximum | ||

984 | //negative argument. Negating this would give the result of | ||

985 | //zero, which is not what is intended. | ||

986 | assert(arg != 0x80000000); | ||

987 | |||

988 | if (arg < 0) | ||

989 | return(-arg); | ||

990 | else | ||

991 | return(arg); | ||

992 | } | ||

993 | |||

994 | |||

995 | //07/15/01: Visual inspection only. Function deemed too | ||

996 | //simple for unit testing. | ||

997 | int GMP_INTS_mpz_sgn(const GMP_INTS_mpz_struct *arg) | ||

998 | { | ||

999 | assert(arg != NULL); | ||

1000 | assert(arg->n_allocd > 0); | ||

1001 | |||

1002 | if (arg->size > 0) | ||

1003 | return(1); | ||

1004 | else if (arg->size == 0) | ||

1005 | return(0); | ||

1006 | else | ||

1007 | return(-1); | ||

1008 | } | ||

1009 | |||

1010 | |||

1011 | //07/15/01: Visual inspection only. Function deemed too | ||

1012 | //simple for unit testing. | ||

1013 | int GMP_INTS_mpz_is_neg(const GMP_INTS_mpz_struct *arg) | ||

1014 | { | ||

1015 | assert(arg != NULL); | ||

1016 | assert(arg->n_allocd > 0); | ||

1017 | |||

1018 | if (GMP_INTS_mpz_sgn(arg) == -1) | ||

1019 | return(1); | ||

1020 | else | ||

1021 | return(0); | ||

1022 | } | ||

1023 | |||

1024 | |||

1025 | //07/15/01: Visual inspection only. Function deemed too | ||

1026 | //simple for unit testing. | ||

1027 | int GMP_INTS_mpz_is_zero(const GMP_INTS_mpz_struct *arg) | ||

1028 | { | ||

1029 | assert(arg != NULL); | ||

1030 | assert(arg->n_allocd > 0); | ||

1031 | |||

1032 | if (GMP_INTS_mpz_sgn(arg) == 0) | ||

1033 | return(1); | ||

1034 | else | ||

1035 | return(0); | ||

1036 | } | ||

1037 | |||

1038 | |||

1039 | //07/15/01: Visual inspection only. Function deemed too | ||

1040 | //simple for unit testing. | ||

1041 | int GMP_INTS_mpz_is_pos(const GMP_INTS_mpz_struct *arg) | ||

1042 | { | ||

1043 | assert(arg != NULL); | ||

1044 | assert(arg->n_allocd > 0); | ||

1045 | |||

1046 | if (GMP_INTS_mpz_sgn(arg) == 1) | ||

1047 | return(1); | ||

1048 | else | ||

1049 | return(0); | ||

1050 | } | ||

1051 | |||

1052 | |||

1053 | //07/15/01: Visual inspection only. Function deemed too | ||

1054 | //simple for unit testing. | ||

1055 | int GMP_INTS_mpz_is_odd(const GMP_INTS_mpz_struct *arg) | ||

1056 | { | ||

1057 | assert(arg != NULL); | ||

1058 | assert(arg->n_allocd > 0); | ||

1059 | |||

1060 | if (arg->size == 0) | ||

1061 | return 0; | ||

1062 | else if ((arg->limbs[0] & 0x1) != 0) | ||

1063 | return 1; | ||

1064 | else | ||

1065 | return 0; | ||

1066 | } | ||

1067 | |||

1068 | |||

1069 | //07/15/01: Visual inspection only. Function deemed too | ||

1070 | //simple for unit testing. | ||

1071 | int GMP_INTS_mpz_is_even(const GMP_INTS_mpz_struct *arg) | ||

1072 | { | ||

1073 | assert(arg != NULL); | ||

1074 | assert(arg->n_allocd > 0); | ||

1075 | |||

1076 | if (GMP_INTS_mpz_is_odd(arg)) | ||

1077 | return 0; | ||

1078 | else | ||

1079 | return 1; | ||

1080 | } | ||

1081 | |||

1082 | |||

1083 | void GMP_INTS_mpz_negate(GMP_INTS_mpz_struct *arg) | ||

1084 | { | ||

1085 | //Eyeball the input parameters. | ||

1086 | assert(arg != NULL); | ||

1087 | assert(arg->n_allocd > 0); | ||

1088 | assert(arg->limbs != NULL); | ||

1089 | |||

1090 | arg->size = -(arg->size); | ||

1091 | } | ||

1092 | |||

1093 | |||

1094 | //07/15/01: Visual inspection only. Function deemed too | ||

1095 | //simple for unit testing. | ||

1096 | void GMP_INTS_mpn_normalize(GMP_INTS_limb_ptr limb_array, | ||

1097 | GMP_INTS_size_t *idx) | ||

1098 | { | ||

1099 | assert(limb_array != NULL); | ||

1100 | assert(idx != NULL); | ||

1101 | assert(idx >= 0); | ||

1102 | |||

1103 | while (*idx > 0) | ||

1104 | { | ||

1105 | if (limb_array[*idx - 1] != 0) | ||

1106 | break; | ||

1107 | (*idx)--; | ||

1108 | } | ||

1109 | } | ||

1110 | |||

1111 | |||

1112 | //07/15/01: Visual inspection only. Function deemed too | ||

1113 | //simple for unit testing. | ||

1114 | void GMP_INTS_mpn_copy_limbs(GMP_INTS_limb_ptr dest, | ||

1115 | GMP_INTS_limb_srcptr src, | ||

1116 | GMP_INTS_size_t n) | ||

1117 | { | ||

1118 | GMP_INTS_size_t i; | ||

1119 | |||

1120 | assert(dest != NULL); | ||

1121 | assert(src != NULL); | ||

1122 | assert(n >= 0); | ||

1123 | |||

1124 | for (i=0; i<n; i++) | ||

1125 | dest[i] = src[i]; | ||

1126 | } | ||

1127 | |||

1128 | |||

1129 | /******************************************************************/ | ||

1130 | /*** LOW-LEVEL ARITHMETIC FUNCTIONS ****************************/ | ||

1131 | /******************************************************************/ | ||

1132 | |||

1133 | int GMP_INTS_two_op_flags_map(int flags1, int flags2) | ||

1134 | { | ||

1135 | int cf; | ||

1136 | |||

1137 | if (!flags1 && !flags2) | ||

1138 | { | ||

1139 | return(0); | ||

1140 | } | ||

1141 | else | ||

1142 | { | ||

1143 | cf = flags1 | flags2; | ||

1144 | |||

1145 | if (cf & (GMP_INTS_EF_INTOVF_POS | GMP_INTS_EF_INTOVF_TAINT_POS)) | ||

1146 | { | ||

1147 | //In either case here, the result will be tainted | ||

1148 | //positive. | ||

1149 | return(GMP_INTS_EF_INTOVF_TAINT_POS); | ||

1150 | } | ||

1151 | else if (cf & (GMP_INTS_EF_INTOVF_NEG | GMP_INTS_EF_INTOVF_TAINT_NEG)) | ||

1152 | { | ||

1153 | //In either case here, the result will be tainted | ||

1154 | //negative. | ||

1155 | return(GMP_INTS_EF_INTOVF_TAINT_NEG); | ||

1156 | } | ||

1157 | else | ||

1158 | { | ||

1159 | //This case is where the flags ints are non-zero, but | ||

1160 | //no known bits are set. This is surely some kind of | ||

1161 | //an internal software error. In debug mode, want to | ||

1162 | //alert to error. In actual operation, want to just | ||

1163 | //pretend an ordinary positive taint. | ||

1164 | assert(0); | ||

1165 | return(GMP_INTS_EF_INTOVF_TAINT_POS); | ||

1166 | } | ||

1167 | } | ||

1168 | } | ||

1169 | |||

1170 | |||

1171 | GMP_INTS_limb_t GMP_INTS_mpn_add_1 (GMP_INTS_limb_ptr res_ptr, | ||

1172 | GMP_INTS_limb_srcptr s1_ptr, | ||

1173 | GMP_INTS_size_t s1_size, | ||

1174 | GMP_INTS_limb_t s2_limb) | ||

1175 | { | ||

1176 | GMP_INTS_limb_t x; | ||

1177 | |||

1178 | assert(res_ptr != NULL); | ||

1179 | assert(s1_ptr != NULL); | ||

1180 | assert(s1_size > 0); | ||

1181 | |||

1182 | x = *s1_ptr++; | ||

1183 | s2_limb = x + s2_limb; | ||

1184 | *res_ptr++ = s2_limb; | ||

1185 | //Since limbs are unsigned, the test below tests if there | ||

1186 | //was a carry, i.e. a positive rollover. | ||

1187 | if (s2_limb < x) | ||

1188 | { | ||

1189 | while (--s1_size != 0) | ||

1190 | { | ||

1191 | x = *s1_ptr++ + 1; | ||

1192 | *res_ptr++ = x; | ||

1193 | if (x != 0) | ||

1194 | goto fin; | ||

1195 | } | ||

1196 | |||

1197 | return 1; | ||

1198 | } | ||

1199 | |||

1200 | fin: | ||

1201 | if (res_ptr != s1_ptr) | ||

1202 | { | ||

1203 | GMP_INTS_size_t i; | ||

1204 | for (i = 0; i < s1_size - 1; i++) | ||

1205 | { | ||

1206 | res_ptr[i] = s1_ptr[i]; | ||

1207 | } | ||

1208 | } | ||

1209 | |||

1210 | return 0; | ||

1211 | } | ||

1212 | |||

1213 | |||

1214 | GMP_INTS_limb_t GMP_INTS_mpn_sub_1(GMP_INTS_limb_ptr res_ptr, | ||

1215 | GMP_INTS_limb_srcptr s1_ptr, | ||

1216 | GMP_INTS_size_t s1_size, | ||

1217 | GMP_INTS_limb_t s2_limb) | ||

1218 | { | ||

1219 | GMP_INTS_limb_t x; | ||

1220 | |||

1221 | assert(res_ptr != NULL); | ||

1222 | assert(s1_ptr != NULL); | ||

1223 | assert(s1_size > 0); | ||

1224 | |||

1225 | x = *s1_ptr++; | ||

1226 | s2_limb = x - s2_limb; | ||

1227 | *res_ptr++ = s2_limb; | ||

1228 | //The test below detects a borrow. | ||

1229 | if (s2_limb > x) | ||

1230 | { | ||

1231 | while (--s1_size != 0) | ||

1232 | { | ||

1233 | x = *s1_ptr++; | ||

1234 | *res_ptr++ = x - 1; | ||

1235 | if (x != 0) | ||

1236 | goto fin; | ||

1237 | } | ||

1238 | |||

1239 | return 1; | ||

1240 | } | ||

1241 | |||

1242 | fin: | ||

1243 | if (res_ptr != s1_ptr) | ||

1244 | { | ||

1245 | GMP_INTS_size_t i; | ||

1246 | for (i = 0; i < s1_size - 1; i++) | ||

1247 | { | ||

1248 | res_ptr[i] = s1_ptr[i]; | ||

1249 | } | ||

1250 | } | ||

1251 | return 0; | ||

1252 | } | ||

1253 | |||

1254 | |||

1255 | //07/15/01: Am willing to skip unit-testing on this. | ||

1256 | //Understand the logic (i.e. passes visual inspection), | ||

1257 | //and comes from GNU-MP. Hope any defects here will be | ||

1258 | //caught in testing of GMP_INTS_mpz_mul() and other | ||

1259 | //higher-level functions. | ||

1260 | GMP_INTS_limb_t GMP_INTS_mpn_mul_1 (GMP_INTS_limb_ptr res_ptr, | ||

1261 | GMP_INTS_limb_srcptr s1_ptr, | ||

1262 | GMP_INTS_size_t s1_size, | ||

1263 | GMP_INTS_limb_t s2_limb) | ||

1264 | { | ||

1265 | GMP_INTS_limb_t cy_limb; | ||

1266 | GMP_INTS_size_t j; | ||

1267 | GMP_INTS_limb_t prod_high, prod_low; | ||

1268 | unsigned _int64 temp; | ||

1269 | |||

1270 | assert(res_ptr != NULL); | ||

1271 | assert(s1_ptr != NULL); | ||

1272 | assert(s1_size > 0); | ||

1273 | |||

1274 | /* The loop counter and index J goes from -S1_SIZE to -1. This way | ||

1275 | the loop becomes faster. */ | ||

1276 | j = -s1_size; | ||

1277 | |||

1278 | /* Offset the base pointers to compensate for the negative indices. */ | ||

1279 | s1_ptr -= j; | ||

1280 | res_ptr -= j; | ||

1281 | |||

1282 | cy_limb = 0; | ||

1283 | do | ||

1284 | { | ||

1285 | //The original code here was the following macro: | ||

1286 | //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb); | ||

1287 | //Will use the 64-bit data type of MSVC++ to achieve | ||

1288 | //the same effect. | ||

1289 | // | ||

1290 | //NOTE AS OF 07/13/01: I have looked at the assembly- | ||

1291 | //language, and the lines below are a real sore spot. | ||

1292 | //The multiply is fairly direct (although there is a | ||

1293 | //function call), but the shift does not behave as | ||

1294 | //expected--there is a function call and a loop to | ||

1295 | //go through the 32 iterations. After logical testing, | ||

1296 | //may want to clean this out--this would surely | ||

1297 | //result in a speed increase. This is a sore spot. | ||

1298 | temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb); | ||

1299 | prod_low = (GMP_INTS_limb_t)temp; | ||

1300 | prod_high = (GMP_INTS_limb_t)(temp >> 32); | ||

1301 | |||

1302 | prod_low += cy_limb; | ||

1303 | cy_limb = (prod_low < cy_limb) + prod_high; | ||

1304 | |||

1305 | res_ptr[j] = prod_low; | ||

1306 | } | ||

1307 | while (++j != 0); | ||

1308 | |||

1309 | return cy_limb; | ||

1310 | } | ||

1311 | |||

1312 | |||

1313 | //07/15/01: Am willing to skip unit-testing on this. | ||

1314 | //Understand the logic (i.e. passes visual inspection), | ||

1315 | //and comes from GNU-MP. Hope any defects here will be | ||

1316 | //caught in testing of GMP_INTS_mpz_add() and other | ||

1317 | //higher-level functions. | ||

1318 | GMP_INTS_limb_t GMP_INTS_mpn_add_n(GMP_INTS_limb_ptr res_ptr, | ||

1319 | GMP_INTS_limb_srcptr s1_ptr, | ||

1320 | GMP_INTS_limb_srcptr s2_ptr, | ||

1321 | GMP_INTS_size_t size) | ||

1322 | { | ||

1323 | GMP_INTS_limb_t x, y, cy; | ||

1324 | GMP_INTS_size_t j; | ||

1325 | |||

1326 | assert(res_ptr != NULL); | ||

1327 | assert(s1_ptr != NULL); | ||

1328 | assert(s2_ptr != NULL); | ||

1329 | |||

1330 | /* The loop counter and index J goes from -SIZE to -1. This way | ||

1331 | the loop becomes faster. */ | ||

1332 | j = -size; | ||

1333 | |||

1334 | /* Offset the base pointers to compensate for the negative indices. */ | ||

1335 | s1_ptr -= j; | ||

1336 | s2_ptr -= j; | ||

1337 | res_ptr -= j; | ||

1338 | |||

1339 | cy = 0; | ||

1340 | do | ||

1341 | { | ||

1342 | y = s2_ptr[j]; | ||

1343 | x = s1_ptr[j]; | ||

1344 | y += cy; /* add previous carry to one addend */ | ||

1345 | cy = (y < cy); /* get out carry from that addition */ | ||

1346 | y = x + y; /* add other addend */ | ||

1347 | cy = (y < x) + cy; /* get out carry from that add, combine */ | ||

1348 | res_ptr[j] = y; | ||

1349 | } | ||

1350 | while (++j != 0); | ||

1351 | |||

1352 | return cy; | ||

1353 | } | ||

1354 | |||

1355 | |||

1356 | //07/15/01: Am willing to skip unit-testing on this. | ||

1357 | //Understand the logic (i.e. passes visual inspection), | ||

1358 | //and comes from GNU-MP. Hope any defects here will be | ||

1359 | //caught in testing of GMP_INTS_mpz_mul() and other | ||

1360 | //higher-level functions. | ||

1361 | GMP_INTS_limb_t GMP_INTS_mpn_addmul_1 (GMP_INTS_limb_ptr res_ptr, | ||

1362 | GMP_INTS_limb_srcptr s1_ptr, | ||

1363 | GMP_INTS_size_t s1_size, | ||

1364 | GMP_INTS_limb_t s2_limb) | ||

1365 | { | ||

1366 | GMP_INTS_limb_t cy_limb; | ||

1367 | GMP_INTS_size_t j; | ||

1368 | GMP_INTS_limb_t prod_high, prod_low; | ||

1369 | GMP_INTS_limb_t x; | ||

1370 | unsigned _int64 temp; | ||

1371 | |||

1372 | //Eyeball the inputs carefully. | ||

1373 | assert(res_ptr != NULL); | ||

1374 | assert(s1_ptr != NULL); | ||

1375 | assert(s1_size > 0); | ||

1376 | |||

1377 | /* The loop counter and index J goes from -SIZE to -1. This way | ||

1378 | the loop becomes faster. */ | ||

1379 | j = -s1_size; | ||

1380 | |||

1381 | /* Offset the base pointers to compensate for the negative indices. */ | ||

1382 | res_ptr -= j; | ||

1383 | s1_ptr -= j; | ||

1384 | |||

1385 | cy_limb = 0; | ||

1386 | do | ||

1387 | { | ||

1388 | //The original code here was the following macro: | ||

1389 | //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb); | ||

1390 | //Will use the 64-bit data type of MSVC++ to achieve | ||

1391 | //the same effect. | ||

1392 | // | ||

1393 | //NOTE AS OF 07/14/01: I have not looked at the assembly- | ||

1394 | //language, but the assembly-language generated by what | ||

1395 | //is below is suspected to have performance problems. | ||

1396 | //May want to come back to this. | ||

1397 | temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb); | ||

1398 | prod_low = (GMP_INTS_limb_t)temp; | ||

1399 | prod_high = (GMP_INTS_limb_t)(temp >> 32); | ||

1400 | |||

1401 | prod_low += cy_limb; | ||

1402 | cy_limb = (prod_low < cy_limb) + prod_high; | ||

1403 | |||

1404 | x = res_ptr[j]; | ||

1405 | prod_low = x + prod_low; | ||

1406 | cy_limb += (prod_low < x); | ||

1407 | res_ptr[j] = prod_low; | ||

1408 | } | ||

1409 | while (++j != 0); | ||

1410 | |||

1411 | return cy_limb; | ||

1412 | } | ||

1413 | |||

1414 | |||

1415 | //07/15/01: Am willing to skip unit-testing on this. | ||

1416 | //Understand the logic (i.e. passes visual inspection), | ||

1417 | //and comes from GNU-MP. | ||

1418 | GMP_INTS_limb_t GMP_INTS_mpn_add (GMP_INTS_limb_ptr res_ptr, | ||

1419 | GMP_INTS_limb_srcptr s1_ptr, | ||

1420 | GMP_INTS_size_t s1_size, | ||

1421 | GMP_INTS_limb_srcptr s2_ptr, | ||

1422 | GMP_INTS_size_t s2_size) | ||

1423 | { | ||

1424 | GMP_INTS_limb_t cy_limb = 0; | ||

1425 | |||

1426 | assert(res_ptr != NULL); | ||

1427 | assert(s1_ptr != NULL); | ||

1428 | assert(s2_ptr != NULL); | ||

1429 | |||

1430 | //Numbers apparently must be arranged with sizes so that | ||

1431 | //LIMBS(s1) >= LIMBS(s2). | ||

1432 | //Add the parts up to the most significant limb of S2. | ||

1433 | if (s2_size != 0) | ||

1434 | cy_limb = GMP_INTS_mpn_add_n (res_ptr, | ||

1435 | s1_ptr, | ||

1436 | s2_ptr, | ||

1437 | s2_size); | ||

1438 | |||

1439 | //Process the carry result, and propagate the carries up through | ||

1440 | //the parts of S1 that don't exist in S2, i.e. propagate the | ||

1441 | //carries upward in S1. | ||

1442 | if (s1_size - s2_size != 0) | ||

1443 | cy_limb = GMP_INTS_mpn_add_1 (res_ptr + s2_size, | ||

1444 | s1_ptr + s2_size, | ||

1445 | s1_size - s2_size, | ||

1446 | cy_limb); | ||

1447 | return cy_limb; | ||

1448 | } | ||

1449 | |||

1450 | |||

1451 | //07/15/01: Am willing to skip unit-testing on this. | ||

1452 | //Understand the logic (i.e. passes visual inspection), | ||

1453 | //and comes from GNU-MP. | ||

1454 | GMP_INTS_limb_t GMP_INTS_mpn_sub_n(GMP_INTS_limb_ptr res_ptr, | ||

1455 | GMP_INTS_limb_srcptr s1_ptr, | ||

1456 | GMP_INTS_limb_srcptr s2_ptr, | ||

1457 | GMP_INTS_size_t size) | ||

1458 | { | ||

1459 | GMP_INTS_limb_t x, y, cy; | ||

1460 | GMP_INTS_size_t j; | ||

1461 | |||

1462 | assert(res_ptr != NULL); | ||

1463 | assert(s1_ptr != NULL); | ||

1464 | assert(s2_ptr != NULL); | ||

1465 | |||

1466 | /* The loop counter and index J goes from -SIZE to -1. This way | ||

1467 | the loop becomes faster. */ | ||

1468 | j = -size; | ||

1469 | |||

1470 | /* Offset the base pointers to compensate for the negative indices. */ | ||

1471 | s1_ptr -= j; | ||

1472 | s2_ptr -= j; | ||

1473 | res_ptr -= j; | ||

1474 | |||

1475 | cy = 0; | ||

1476 | do | ||

1477 | { | ||

1478 | y = s2_ptr[j]; | ||

1479 | x = s1_ptr[j]; | ||

1480 | y += cy; /* add previous carry to subtrahend */ | ||

1481 | cy = (y < cy); /* get out carry from that addition */ | ||

1482 | y = x - y; /* main subtract */ | ||

1483 | cy = (y > x) + cy; /* get out carry from the subtract, combine */ | ||

1484 | res_ptr[j] = y; | ||

1485 | } | ||

1486 | while (++j != 0); | ||

1487 | |||

1488 | return cy; | ||

1489 | } | ||

1490 | |||

1491 | |||

1492 | //07/17/01: Am willing to skip unit-testing on this. | ||

1493 | //Understand the logic (i.e. passes visual inspection), | ||

1494 | //and comes from GNU-MP. | ||

1495 | GMP_INTS_limb_t GMP_INTS_mpn_sub (GMP_INTS_limb_ptr res_ptr, | ||

1496 | GMP_INTS_limb_srcptr s1_ptr, | ||

1497 | GMP_INTS_size_t s1_size, | ||

1498 | GMP_INTS_limb_srcptr s2_ptr, | ||

1499 | GMP_INTS_size_t s2_size) | ||

1500 | { | ||

1501 | GMP_INTS_limb_t cy_limb = 0; | ||

1502 | |||

1503 | assert(res_ptr != NULL); | ||

1504 | assert(s1_ptr != NULL); | ||

1505 | assert(s2_ptr != NULL); | ||

1506 | |||

1507 | if (s2_size != 0) | ||

1508 | cy_limb = GMP_INTS_mpn_sub_n(res_ptr, | ||

1509 | s1_ptr, | ||

1510 | s2_ptr, | ||

1511 | s2_size); | ||

1512 | |||

1513 | if (s1_size - s2_size != 0) | ||

1514 | cy_limb = GMP_INTS_mpn_sub_1(res_ptr + s2_size, | ||

1515 | s1_ptr + s2_size, | ||

1516 | s1_size - s2_size, | ||

1517 | cy_limb); | ||

1518 | return cy_limb; | ||

1519 | } | ||

1520 | |||

1521 | |||

1522 | //07/17/01: Am willing to skip unit-testing on this. | ||

1523 | //Understand the logic (i.e. passes visual inspection), | ||

1524 | //and comes from GNU-MP. | ||

1525 | GMP_INTS_limb_t GMP_INTS_mpn_lshift(GMP_INTS_limb_ptr wp, | ||

1526 | GMP_INTS_limb_srcptr up, | ||

1527 | GMP_INTS_size_t usize, | ||

1528 | unsigned int cnt) | ||

1529 | { | ||

1530 | GMP_INTS_limb_t high_limb, low_limb; | ||

1531 | unsigned sh_1, sh_2; | ||

1532 | GMP_INTS_size_t i; | ||

1533 | GMP_INTS_limb_t retval; | ||

1534 | |||

1535 | assert(wp != NULL); | ||

1536 | assert(up != NULL); | ||

1537 | assert(usize > 0); | ||

1538 | assert(cnt > 0); | ||

1539 | |||

1540 | sh_1 = cnt; | ||

1541 | |||

1542 | wp += 1; | ||

1543 | sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1; | ||

1544 | //This automatically implies that can't call this function to shift more | ||

1545 | //than 31 places. | ||

1546 | i = usize - 1; | ||

1547 | low_limb = up[i]; | ||

1548 | retval = low_limb >> sh_2; //Return value is the amount shifted | ||

1549 | //off the top. | ||

1550 | high_limb = low_limb; | ||

1551 | while (--i >= 0) | ||

1552 | { | ||

1553 | low_limb = up[i]; | ||

1554 | wp[i] = (high_limb << sh_1) | (low_limb >> sh_2); | ||

1555 | high_limb = low_limb; | ||

1556 | } | ||

1557 | wp[i] = high_limb << sh_1; | ||

1558 | |||

1559 | return retval; | ||

1560 | } | ||

1561 | |||

1562 | |||

1563 | //07/17/01: Am willing to skip unit-testing on this. | ||

1564 | //Understand the logic more or less (i.e. passes visual inspection), | ||

1565 | //and comes from GNU-MP. | ||

1566 | /* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right | ||

1567 | and store the USIZE least significant limbs of the result at WP. | ||

1568 | The bits shifted out to the right are returned. | ||

1569 | |||

1570 | Argument constraints: | ||

1571 | 1. 0 < CNT < BITS_PER_MP_LIMB | ||

1572 | 2. If the result is to be written over the input, WP must be <= UP. | ||

1573 | */ | ||

1574 | GMP_INTS_limb_t GMP_INTS_mpn_rshift (GMP_INTS_limb_ptr wp, | ||

1575 | GMP_INTS_limb_srcptr up, | ||

1576 | GMP_INTS_size_t usize, | ||

1577 | unsigned int cnt) | ||

1578 | { | ||

1579 | GMP_INTS_limb_t high_limb, low_limb; | ||

1580 | unsigned sh_1, sh_2; | ||

1581 | GMP_INTS_size_t i; | ||

1582 | GMP_INTS_limb_t retval; | ||

1583 | |||

1584 | assert(wp != NULL); | ||

1585 | assert(up != NULL); | ||

1586 | assert(usize > 0); | ||

1587 | assert(cnt > 0); | ||

1588 | |||

1589 | sh_1 = cnt; | ||

1590 | |||

1591 | wp -= 1; | ||

1592 | sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1; | ||

1593 | high_limb = up[0]; | ||

1594 | retval = high_limb << sh_2; | ||

1595 | low_limb = high_limb; | ||

1596 | |||

1597 | for (i = 1; i < usize; i++) | ||

1598 | { | ||

1599 | high_limb = up[i]; | ||

1600 | wp[i] = (low_limb >> sh_1) | (high_limb << sh_2); | ||

1601 | low_limb = high_limb; | ||

1602 | } | ||

1603 | wp[i] = low_limb >> sh_1; | ||

1604 | |||

1605 | return retval; | ||

1606 | } | ||

1607 | |||

1608 | |||

1609 | //07/17/01: Am willing to skip unit-testing on this. | ||

1610 | //Understand the logic (i.e. passes visual inspection), | ||

1611 | //and comes from GNU-MP. | ||

1612 | int GMP_INTS_mpn_cmp (GMP_INTS_limb_srcptr op1_ptr, | ||

1613 | GMP_INTS_limb_srcptr op2_ptr, | ||

1614 | GMP_INTS_size_t size) | ||

1615 | { | ||

1616 | GMP_INTS_size_t i; | ||

1617 | GMP_INTS_limb_t op1_word, op2_word; | ||

1618 | |||

1619 | assert(op1_ptr != NULL); | ||

1620 | assert(op2_ptr != NULL); | ||

1621 | |||

1622 | for (i = size - 1; i >= 0; i--) | ||

1623 | { | ||

1624 | op1_word = op1_ptr[i]; | ||

1625 | op2_word = op2_ptr[i]; | ||

1626 | if (op1_word != op2_word) | ||

1627 | goto diff; | ||

1628 | } | ||

1629 | return 0; | ||

1630 | |||

1631 | diff: | ||

1632 | //This can *not* be simplified to | ||

1633 | // op2_word - op2_word | ||

1634 | //since that expression might give signed overflow. | ||

1635 | return (op1_word > op2_word) ? 1 : -1; | ||

1636 | } | ||

1637 | |||

1638 | |||

1639 | //07/15/01: Am willing to skip unit-testing on this. | ||

1640 | //Understand the logic (i.e. passes visual inspection), | ||

1641 | //and comes from GNU-MP. Hope any defects here will be | ||

1642 | //caught in testing of GMP_INTS_mpz_mul() and other | ||

1643 | //higher-level functions. | ||

1644 | void GMP_INTS_mpn_mul_basecase (GMP_INTS_limb_ptr prodp, | ||

1645 | GMP_INTS_limb_srcptr up, | ||

1646 | GMP_INTS_size_t usize, | ||

1647 | GMP_INTS_limb_srcptr vp, | ||

1648 | GMP_INTS_size_t vsize) | ||

1649 | { | ||

1650 | assert(prodp != NULL); | ||

1651 | assert(up != NULL); | ||

1652 | assert(usize > 0); | ||

1653 | assert(vp != NULL); | ||

1654 | assert(vsize > 0); | ||

1655 | |||

1656 | /* We first multiply by the low order one or two limbs, as the result can | ||

1657 | be stored, not added, to PROD. We also avoid a loop for zeroing this | ||

1658 | way. */ | ||

1659 | prodp[usize] = GMP_INTS_mpn_mul_1 (prodp, up, usize, vp[0]); | ||

1660 | prodp++; | ||

1661 | vp++; | ||

1662 | vsize--; | ||

1663 | |||

1664 | /* For each iteration in the loop, multiply U with one limb from V, and | ||

1665 | add the result to PROD. */ | ||

1666 | while (vsize != 0) | ||

1667 | { | ||

1668 | prodp[usize] = GMP_INTS_mpn_addmul_1 (prodp, up, usize, vp[0]); | ||

1669 | prodp++, | ||

1670 | vp++, | ||

1671 | vsize--; | ||

1672 | } | ||

1673 | } | ||

1674 | |||

1675 | |||

1676 | //07/15/01: No unit testing possible--this is a passthrough. | ||

1677 | //In the original GNU MP code, there were several multiplication | ||

1678 | //algorithms, and this function would select one based on the | ||

1679 | //size of the operands and other considerations. The code has been | ||

1680 | //pared so that only simple multiplication is used, which is why | ||

1681 | //this function contains only a single pass-thru function call. | ||

1682 | void GMP_INTS_mpn_mul_n (GMP_INTS_limb_ptr p, | ||

1683 | GMP_INTS_limb_srcptr a, | ||

1684 | GMP_INTS_limb_srcptr b, | ||

1685 | GMP_INTS_size_t n) | ||

1686 | { | ||

1687 | GMP_INTS_mpn_mul_basecase (p, a, n, b, n); | ||

1688 | } | ||

1689 | |||

1690 | |||

1691 | //07/16/01: Visual inspection OK. Will not perform unit testing. | ||

1692 | GMP_INTS_limb_t GMP_INTS_mpn_mul(GMP_INTS_limb_ptr prodp, | ||

1693 | GMP_INTS_limb_srcptr up, | ||

1694 | GMP_INTS_size_t un, | ||

1695 | GMP_INTS_limb_srcptr vp, | ||

1696 | GMP_INTS_size_t vn) | ||

1697 | { | ||

1698 | //This is a gutted version of the GNU MP function. The GNU | ||

1699 | //MP function considered the case of a square, and also | ||

1700 | //better algorithms that pay off with large operands. | ||

1701 | //This gutted version uses only basic multiplication | ||

1702 | //(O(N**2)). | ||

1703 | |||

1704 | //Eyeball the input parameters. | ||

1705 | assert(prodp != NULL); | ||

1706 | assert(up != NULL); | ||

1707 | assert(un >= 0); | ||

1708 | assert(vp != NULL); | ||

1709 | assert(vn >= 0); | ||

1710 | |||

1711 | /* Basic long multiplication. */ | ||

1712 | GMP_INTS_mpn_mul_basecase (prodp, up, un, vp, vn); | ||

1713 | |||

1714 | //Return the most significant limb (which might be zero). | ||

1715 | //This is different than | ||

1716 | //most other functions, which return the spillover. | ||

1717 | return prodp[un + vn - 1]; | ||

1718 | } | ||

1719 | |||

1720 | |||

1721 | /******************************************************************/ | ||

1722 | /*** LIMB SPACE REALLOCATION FUNCTIONS *************************/ | ||

1723 | /******************************************************************/ | ||

1724 | |||

1725 | void *GMP_INTS_mpz_realloc (GMP_INTS_mpz_struct *m, | ||

1726 | GMP_INTS_size_t new_size) | ||

1727 | { | ||

1728 | /* Never allocate zero space. */ | ||

1729 | if (new_size <= 0) | ||

1730 | new_size = 1; | ||

1731 | |||

1732 | m->limbs = (GMP_INTS_limb_ptr) | ||

1733 | GMP_INTS_realloc_w_size (m->limbs, | ||

1734 | m->n_allocd * sizeof(GMP_INTS_limb_t), | ||

1735 | new_size * sizeof(GMP_INTS_limb_t)); | ||

1736 | m->n_allocd = new_size; | ||

1737 | |||

1738 | return (void *) m->limbs; | ||

1739 | } | ||

1740 | |||

1741 | |||

1742 | /******************************************************************/ | ||

1743 | /*** PUBLIC INITIALIZATION AND MEMORY MANAGEMENT FUNCTIONS *****/ | ||

1744 | /******************************************************************/ | ||

1745 | |||

1746 | void GMP_INTS_mpz_init (GMP_INTS_mpz_struct *x) | ||

1747 | { | ||

1748 | assert(x != NULL); | ||

1749 | |||

1750 | //The structure (the header block) exists in the | ||

1751 | //caller's area. Most likely it is a local variable. | ||

1752 | //This is OK, because it doesn't take up much space. | ||

1753 | |||

1754 | //Start off with no errors. | ||

1755 | x->flags = 0; | ||

1756 | |||

1757 | //Allocate space for one limb, which is the most | ||

1758 | //basic amount. This will grow, almost certainly. | ||

1759 | x->limbs = GMP_INTS_malloc(sizeof(GMP_INTS_limb_t)); | ||

1760 | |||

1761 | //Indicate that one limb was allocated. | ||

1762 | x->n_allocd = 1; | ||

1763 | |||

1764 | //Set the size to 0. This signals a value of zero. | ||

1765 | x->size = 0; | ||

1766 | } | ||

1767 | |||

1768 | |||

1769 | void GMP_INTS_mpz_clear (GMP_INTS_mpz_struct *x) | ||

1770 | { | ||

1771 | //Be sure the passed pointer is not NULL. | ||

1772 | assert(x != NULL); | ||

1773 | |||

1774 | //Be sure that the amount allocated is also above zero. | ||

1775 | //Anything else represents a logical error. | ||

1776 | assert(x->n_allocd > 0); | ||

1777 | |||

1778 | //Be sure that the pointer to the allocated limbs | ||

1779 | //is not NULL. Anything else would be a logical | ||

1780 | //error. | ||

1781 | assert(x->limbs != NULL); | ||

1782 | |||

1783 | //Deallocate the space for the limbs. The pointer is | ||

1784 | //set NULL and the allocated amount set to zero | ||

1785 | // so in case clear is called again it will be | ||

1786 | //a detectable error. | ||

1787 | GMP_INTS_free_w_size(x->limbs, | ||

1788 | x->n_allocd * sizeof(GMP_INTS_limb_t)); | ||

1789 | x->limbs = NULL; | ||

1790 | x->n_allocd = 0; | ||

1791 | } | ||

1792 | |||

1793 | |||

1794 | /******************************************************************/ | ||

1795 | /*** PUBLIC ASSIGNMENT FUNCTIONS *******************************/ | ||

1796 | /******************************************************************/ | ||

1797 | |||

1798 | void GMP_INTS_mpz_copy( GMP_INTS_mpz_struct *dst, | ||

1799 | const GMP_INTS_mpz_struct *src) | ||

1800 | { | ||

1801 | GMP_INTS_size_t i, n; | ||

1802 | |||

1803 | //Eyeball the input parameters. | ||

1804 | assert(dst != NULL); | ||

1805 | assert(dst->n_allocd > 0); | ||

1806 | assert(dst->limbs != NULL); | ||

1807 | assert(src != NULL); | ||

1808 | assert(src->n_allocd > 0); | ||

1809 | assert(src->limbs != NULL); | ||

1810 | |||

1811 | //Source and destination may not be the same. | ||

1812 | assert(src != dst); | ||

1813 | |||

1814 | //Figure out the real size of the source. We need to take the absolute | ||

1815 | //value. | ||

1816 | n = GMP_INTS_abs_of_size_t(src->size); | ||

1817 | |||

1818 | //Reallocate the destination to be bigger if necessary. | ||

1819 | if (dst->n_allocd < n) | ||

1820 | { | ||

1821 | GMP_INTS_mpz_realloc (dst, n); | ||

1822 | } | ||

1823 | |||

1824 | //Copy the non-dynamic fields in the header. | ||

1825 | dst->flags = src->flags; | ||

1826 | dst->size = src->size; | ||

1827 | |||

1828 | //Copy the limbs. | ||

1829 | for (i=0; i<n; i++) | ||

1830 | dst->limbs[i] = src->limbs[i]; | ||

1831 | } | ||

1832 | |||

1833 | |||

1834 | void GMP_INTS_mpz_set_ui (GMP_INTS_mpz_struct *dest, | ||

1835 | unsigned long int val) | ||

1836 | { | ||

1837 | assert(dest != NULL); | ||

1838 | |||

1839 | /* We don't check if the allocation is enough, since the rest of the | ||

1840 | package ensures it's at least 1, which is what we need here. */ | ||

1841 | |||

1842 | dest->flags = 0; | ||

1843 | //A set operation resets any errors. | ||

1844 | |||

1845 | if (val > 0) | ||

1846 | { | ||

1847 | dest->limbs[0] = val; | ||

1848 | dest->size = 1; | ||

1849 | } | ||

1850 | else | ||

1851 | { | ||

1852 | dest->size = 0; | ||

1853 | } | ||

1854 | } | ||

1855 | |||

1856 | |||

1857 | void GMP_INTS_mpz_set_si (GMP_INTS_mpz_struct *dest, | ||

1858 | signed long int val) | ||

1859 | { | ||

1860 | assert(dest != NULL); | ||

1861 | |||

1862 | /* We don't check if the allocation is enough, since the rest of the | ||

1863 | package ensures it's at least 1, which is what we need here. */ | ||

1864 | |||

1865 | dest->flags = 0; | ||

1866 | //A set operation resets any errors. | ||

1867 | |||

1868 | if (val > 0) | ||

1869 | { | ||

1870 | dest->limbs[0] = val; | ||

1871 | dest->size = 1; | ||

1872 | } | ||

1873 | else if (val < 0) | ||

1874 | { | ||

1875 | dest->limbs[0] = (unsigned long) -val; | ||

1876 | dest->size = -1; | ||

1877 | } | ||

1878 | else | ||

1879 | { | ||

1880 | dest->size = 0; | ||

1881 | } | ||

1882 | } | ||

1883 | |||

1884 | |||

1885 | void GMP_INTS_mpz_set_simple_char_str(GMP_INTS_mpz_struct *z, | ||

1886 | const char *s) | ||

1887 | { | ||

1888 | int sign=1; | ||

1889 | int digval; | ||

1890 | GMP_INTS_mpz_struct digvalz, k10; | ||

1891 | |||

1892 | //Eyeball the arguments. | ||

1893 | assert(z != NULL); | ||

1894 | assert(z->n_allocd > 0); | ||

1895 | assert(z->limbs != NULL); | ||

1896 | assert(s != NULL); | ||

1897 | |||

1898 | //Set the arbitrary integer to zero. This will also kill | ||

1899 | //any error flags. | ||

1900 | GMP_INTS_mpz_set_ui(z, 0); | ||

1901 | |||

1902 | //Allocate an integer for our private use to hold each digit | ||

1903 | //value. | ||

1904 | GMP_INTS_mpz_init(&digvalz); | ||

1905 | |||

1906 | //Allocate the constant 10, which we will use often. | ||

1907 | GMP_INTS_mpz_init(&k10); | ||

1908 | GMP_INTS_mpz_set_ui(&k10, 10); | ||

1909 | |||

1910 | //As long as there are are digits and no flags set, keep | ||

1911 | //multiplying and adding the value of the digit. Non- | ||

1912 | //digits are simply ignored. | ||

1913 | while (!(z->flags) && (*s)) | ||

1914 | { | ||

1915 | if (*s == '-') | ||

1916 | { | ||

1917 | sign = -sign; | ||

1918 | } | ||

1919 | else | ||

1920 | { | ||

1921 | digval = CHARFUNC_digit_to_val(*s); | ||

1922 | if (digval >= 0) | ||

1923 | { | ||

1924 | GMP_INTS_mpz_set_ui(&digvalz, digval); | ||

1925 | GMP_INTS_mpz_mul(z, z, &k10); | ||

1926 | GMP_INTS_mpz_add(z, z, &digvalz); | ||

1927 | } | ||

1928 | } | ||

1929 | s++; | ||

1930 | } | ||

1931 | |||

1932 | //Adjust the final sign of the result. | ||

1933 | if (sign < 0) | ||

1934 | z->size = -(z->size); | ||

1935 | |||

1936 | //Deallocate our temporary integers. | ||

1937 | GMP_INTS_mpz_clear(&digvalz); | ||

1938 | GMP_INTS_mpz_clear(&k10); | ||

1939 | } | ||

1940 | |||

1941 | |||

1942 | void GMP_INTS_mpz_set_sci_not_num(GMP_INTS_mpz_struct *z, | ||

1943 | int *failure, | ||

1944 | const char *s) | ||

1945 | { | ||

1946 | int parse_failure; | ||

1947 | //Return code from the floating point parsing | ||

1948 | //function. | ||

1949 | char mant_sign; | ||

1950 | //Sign character, if any, from the mantissa, | ||

1951 | //or N otherwise. | ||

1952 | size_t mant_bdp; | ||

1953 | //The index to the start of the mantissa before | ||

1954 | //the decimal point. | ||

1955 | size_t mant_bdp_len; | ||

1956 | //The length of the mantissa before the decimal | ||

1957 | //point. Zero means not defined, i.e. that | ||

1958 | //no characters were parsed and interpreted as | ||

1959 | //that part of a floating point number. | ||

1960 | size_t mant_adp; | ||

1961 | size_t mant_adp_len; | ||

1962 | //Similar fields for after the decimal point. | ||

1963 | char exp_sign; | ||

1964 | //Sign of the exponent, if any, or N otherwise. | ||

1965 | size_t exp; | ||

1966 | size_t exp_len; | ||

1967 | //Similar fields as to the mantissa, but for the | ||

1968 | //exponent. | ||

1969 | size_t si; | ||

1970 | //Iteration variable. | ||

1971 | int exponent_val; | ||

1972 | //The value of the exponent. We can't accept | ||

1973 | //an exponent outside the range of a 24-bit | ||

1974 | //signed integer. The 24-bit limit is arbitrary. | ||

1975 | //For one thing, it gives room to detect overflow | ||

1976 | //as are adding and multiplying by 10. | ||

1977 | |||

1978 | //Eyeball the input parameters. | ||

1979 | assert(z != NULL); | ||

1980 | assert(z->n_allocd > 0); | ||

1981 | assert(z->limbs != NULL); | ||

1982 | assert(failure != NULL); | ||

1983 | assert(s != NULL); | ||

1984 | |||

1985 | //Start off believing no failure. | ||

1986 | *failure = 0; | ||

1987 | |||

1988 | //Set the output to zero. This is the default case for some | ||

1989 | //steps below. | ||

1990 | GMP_INTS_mpz_set_ui(z, 0); | ||

1991 | |||

1992 | //Attempt to parse the number as a general number | ||

1993 | //in scientific notation. | ||

1994 | BSTRFUNC_parse_gen_sci_not_num(s, | ||

1995 | &parse_failure, | ||

1996 | &mant_sign, | ||

1997 | &mant_bdp, | ||

1998 | &mant_bdp_len, | ||

1999 | &mant_adp, | ||

2000 | &mant_adp_len, | ||

2001 | &exp_sign, | ||

2002 | &exp, | ||

2003 | &exp_len); | ||

2004 | |||

2005 | //If it wouldn't parse as a general number, can't go further. | ||

2006 | if (parse_failure) | ||

2007 | { | ||

2008 | *failure = 1; | ||

2009 | return; | ||

2010 | } | ||

2011 | else if (!exp_len && !mant_adp_len) | ||

2012 | { | ||

2013 | //There was no exponent, and no portion after | ||

2014 | //the decimal point. Can just parse as an integer. | ||

2015 | char *temp_buf; | ||

2016 | |||

2017 | //Allocate the temporary buffer to be one character longer | ||

2018 | //than the length specified for the parsed mantissa. | ||

2019 | temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1)); | ||

2020 | |||

2021 | //Copy from the parsed area into the temporary buffer. | ||

2022 | for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++) | ||

2023 | temp_buf[si-mant_bdp] = s[si]; | ||

2024 | temp_buf[mant_bdp_len] = 0; | ||

2025 | |||

2026 | //Set the arbitrary integer to the value of the character | ||

2027 | //string. | ||

2028 | GMP_INTS_mpz_set_simple_char_str(z, temp_buf); | ||

2029 | |||

2030 | //If the number parsed as negative, invert. | ||

2031 | if (mant_sign == '-') | ||

2032 | z->size = -z->size; | ||

2033 | |||

2034 | //Deallocate the temporary buffer. | ||

2035 | GMP_INTS_free(temp_buf); | ||

2036 | } | ||

2037 | else if (!exp_len && mant_adp_len) | ||

2038 | { | ||

2039 | char *temp_buf; | ||

2040 | |||

2041 | //In this case, there are digits after the decimal point, | ||

2042 | //but no exponent specified. The only way this makes | ||

2043 | //sense is if all of the digits are zero--otherwise it | ||

2044 | //cannot be an integer. | ||

2045 | for (si=mant_adp; si<(mant_adp+mant_adp_len); si++) | ||

2046 | { | ||

2047 | if (s[si] != '0') | ||

2048 | { | ||

2049 | *failure = 1; | ||

2050 | return; | ||

2051 | } | ||

2052 | } | ||

2053 | |||

2054 | //We're clean. They are only zeros. Execute as per | ||

2055 | //integer code. | ||

2056 | |||

2057 | //Allocate the temporary buffer to be one character longer | ||

2058 | //than the length specified for the parsed mantissa. | ||

2059 | temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1)); | ||

2060 | |||

2061 | //Copy from the parsed area into the temporary buffer. | ||

2062 | for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++) | ||

2063 | temp_buf[si-mant_bdp] = s[si]; | ||

2064 | temp_buf[mant_bdp_len] = 0; | ||

2065 | |||

2066 | //Set the arbitrary integer to the value of the character | ||

2067 | //string. | ||

2068 | GMP_INTS_mpz_set_simple_char_str(z, temp_buf); | ||

2069 | |||

2070 | //If the number parsed as negative, invert. | ||

2071 | if (mant_sign == '-') | ||

2072 | z->size = -z->size; | ||

2073 | |||

2074 | //Deallocate the temporary buffer. | ||

2075 | GMP_INTS_free(temp_buf); | ||

2076 | } | ||

2077 | else if (exp_len) | ||

2078 | { | ||

2079 | //This is the most difficult case, where an exponent | ||

2080 | //is specified. There are several complex subcases, | ||

2081 | //such as: | ||

2082 | // a)If the exponent is too positive or too negative, | ||

2083 | // we can't use it. In general, we won't tackle | ||

2084 | // an exponent that won't fit in a signed 24-bit | ||

2085 | // integer. This provides a range of from | ||

2086 | // -8,388,608 to +8,388,607. This dwarfs the | ||

2087 | // 100,000 or so digit preprocessor limit, | ||

2088 | // and should be adequate for any practical | ||

2089 | // application. | ||

2090 | // b)If the exponent is zero, we ignore it. | ||

2091 | // c)If the exponent is positive, it has to | ||

2092 | // be large enough to overcome any | ||

2093 | // digits past the decimal point, otherwise | ||

2094 | // we don't end up with an integer. | ||

2095 | // d)If the exponent is negative, there have to | ||

2096 | // be enough digits so that an integer remains | ||

2097 | // after the exponent is applied. This | ||

2098 | // generally requires trailing zeros on the | ||

2099 | // part before the decimal point. | ||

2100 | |||

2101 | //First, tackle the exponent. Process the | ||

2102 | //exponent into a signed integer. We have to | ||

2103 | //balk at anything outside of 24 bits. The | ||

2104 | //procedure used automatically handles | ||

2105 | //leading zeros correctly. | ||

2106 | exponent_val = 0; | ||

2107 | for (si=exp; si<(exp+exp_len); si++) | ||

2108 | { | ||

2109 | int val; | ||

2110 | |||

2111 | val = CHARFUNC_digit_to_val(s[si]); | ||

2112 | |||

2113 | assert(val >= 0 && val <= 9); | ||

2114 | |||

2115 | exponent_val *= 10; | ||

2116 | exponent_val += val; | ||

2117 | |||

2118 | if (((exp_sign=='-') && (exponent_val>8388608)) | ||

2119 | || | ||

2120 | ((exp_sign != '-') && (exponent_val>8388607))) | ||

2121 | { | ||

2122 | *failure = 1; | ||

2123 | return; | ||

2124 | } | ||

2125 | } | ||

2126 | |||

2127 | //If we're here, the exponent has been computed and | ||

2128 | //is within 24 bits. However, we need to adjust for | ||

2129 | //the sign. | ||

2130 | if (exp_sign == '-') | ||

2131 | exponent_val = -exponent_val; | ||

2132 | |||

2133 | //We need to make accurate assertions about the | ||

2134 | //portion of the number, if any, after the decimal point. | ||

2135 | //This means that we need to effectively discard | ||

2136 | //trailing zeros. To do this, we do not need to | ||

2137 | //relocate the string, we can just back off the index | ||

2138 | //to bypass any trailing zeros. | ||

2139 | while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0')) | ||

2140 | mant_adp_len--; | ||

2141 | |||

2142 | //We also need to make accurate assertions about the | ||

2143 | //portion of the number, if any, before the decimal | ||

2144 | //point. It is known that the parsing function | ||

2145 | //isn't tolerant of multiple zeros, but zero is a | ||

2146 | //special case. Let's advance the pointer to the | ||

2147 | //part before the decimal point so that zero will | ||

2148 | //have zero length. | ||

2149 | while ((mant_bdp_len > 0) && (s[mant_bdp]=='0')) | ||

2150 | { | ||

2151 | mant_bdp++; | ||

2152 | mant_bdp_len--; | ||

2153 | } | ||

2154 | |||

2155 | //If we are dealing with zero, who cares about the | ||

2156 | //exponent? Just return the value of zero. | ||

2157 | if (!mant_bdp_len && !mant_adp_len) | ||

2158 | { | ||

2159 | *failure = 0; | ||

2160 | GMP_INTS_mpz_set_ui(z, 0); | ||

2161 | return; | ||

2162 | } | ||

2163 | |||

2164 | //Beyond this point, we have something non-zero. | ||

2165 | //If the exponent is positive, it must be at least | ||

2166 | //as large as the number of digits beyond the | ||

2167 | //decimal point in order to form an integer. If the | ||

2168 | //exponent is zero, there must be no digits after the | ||

2169 | //decimal point. If the exponent is negative, there | ||

2170 | //must be no digits after the decimal point, and the | ||

2171 | //trailing zeros on the part before the decimal point | ||

2172 | //must be adequate to handle the right decimal shift. | ||

2173 | if (exponent_val == 0) | ||

2174 | { | ||

2175 | if (mant_adp_len) | ||

2176 | { | ||

2177 | *failure = 1; | ||

2178 | return; | ||

2179 | } | ||

2180 | } | ||

2181 | else if (exponent_val > 0) | ||

2182 | { | ||

2183 | if ((int)mant_adp_len > exponent_val) | ||

2184 | { | ||

2185 | *failure = 1; | ||

2186 | return; | ||

2187 | } | ||

2188 | } | ||

2189 | else //exponent_val < 0 | ||

2190 | { | ||

2191 | if (mant_adp_len) | ||

2192 | { | ||

2193 | *failure = 1; | ||

2194 | return; | ||

2195 | } | ||

2196 | else | ||

2197 | { | ||

2198 | //Count the number of trailing zeros on the part | ||

2199 | //before the decimal point. | ||

2200 | size_t trailing_zero_count; | ||

2201 | int idx; | ||

2202 | |||

2203 | trailing_zero_count = 0; | ||

2204 | |||

2205 | for(idx = mant_bdp + mant_bdp_len - 1; | ||

2206 | (mant_bdp_len != 0) && (idx >= (int)mant_bdp); | ||

2207 | idx--) | ||

2208 | { | ||

2209 | if (s[idx] == '0') | ||

2210 | trailing_zero_count++; | ||

2211 | else | ||

2212 | break; | ||

2213 | } | ||

2214 | |||

2215 | //Check on the assertion about trailing zeros. | ||

2216 | if ((int)trailing_zero_count < -exponent_val) | ||

2217 | { | ||

2218 | *failure = 1; | ||

2219 | return; | ||

2220 | } | ||

2221 | } | ||

2222 | } | ||

2223 | |||

2224 | { | ||

2225 | //Create a string long enough to hold the digits | ||

2226 | //before the decimal point plus the ones after and | ||

2227 | //convert that to an arbitrary integer. | ||

2228 | //Form a power of 10 which is 10 exponentiated to | ||

2229 | //the absolute value of the exponent. If the | ||

2230 | //exponent was positive, multiply by it. If the | ||

2231 | //exponent was negative, divide by it. | ||

2232 | char *conv_str; | ||

2233 | size_t sidx; | ||

2234 | GMP_INTS_mpz_struct power_of_ten, k10, trash; | ||

2235 | |||

2236 | GMP_INTS_mpz_init(&power_of_ten); | ||

2237 | GMP_INTS_mpz_init(&k10); | ||

2238 | GMP_INTS_mpz_init(&trash); | ||

2239 | |||

2240 | conv_str = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + mant_adp_len + 1)); | ||

2241 | |||

2242 | sidx=0; | ||

2243 | |||

2244 | for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++) | ||

2245 | { | ||

2246 | conv_str[sidx] = s[si]; | ||

2247 | sidx++; | ||

2248 | } | ||

2249 | for (si=mant_adp; si<(mant_adp+mant_adp_len); si++) | ||

2250 | { | ||

2251 | conv_str[sidx] = s[si]; | ||

2252 | sidx++; | ||

2253 | } | ||

2254 | conv_str[sidx] = 0; | ||

2255 | |||

2256 | assert(sidx == (mant_bdp_len + mant_adp_len)); | ||

2257 | |||

2258 | GMP_INTS_mpz_set_simple_char_str(z, conv_str); | ||

2259 | |||

2260 | GMP_INTS_mpz_set_ui(&k10, 10); | ||

2261 | |||

2262 | if (exponent_val > 0) | ||

2263 | GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, exponent_val-mant_adp_len); | ||

2264 | else | ||

2265 | GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, -exponent_val); | ||

2266 | |||

2267 | if (exponent_val >= 0) | ||

2268 | { | ||

2269 | GMP_INTS_mpz_mul(z, z, &power_of_ten); | ||

2270 | } | ||

2271 | else | ||

2272 | { | ||

2273 | GMP_INTS_mpz_tdiv_qr (&k10, | ||

2274 | &trash, | ||

2275 | z, | ||

2276 | &power_of_ten); | ||

2277 | GMP_INTS_mpz_copy(z, &k10); | ||

2278 | } | ||

2279 | |||

2280 | //If the argument had a minus sign, invert. | ||

2281 | if (mant_sign == '-') | ||

2282 | z->size = -z->size; | ||

2283 | |||

2284 | GMP_INTS_free(conv_str); | ||

2285 | |||

2286 | GMP_INTS_mpz_clear(&trash); | ||

2287 | GMP_INTS_mpz_clear(&k10); | ||

2288 | GMP_INTS_mpz_clear(&power_of_ten); | ||

2289 | |||

2290 | //Finally, if the arbitrary integer has overflowed, this is | ||

2291 | //a parse failure. Must declare as such. | ||

2292 | if (z->flags) | ||

2293 | *failure = 1; | ||

2294 | } | ||

2295 | } | ||

2296 | else | ||

2297 | { | ||

2298 | *failure = 1; | ||

2299 | return; | ||

2300 | } | ||

2301 | } | ||

2302 | |||

2303 | |||

2304 | void GMP_INTS_mpz_set_general_int(GMP_INTS_mpz_struct *z, | ||

2305 | int *failure, | ||

2306 | const char *s) | ||

2307 | { | ||

2308 | //Eyeball the input parameters. | ||

2309 | assert(z != NULL); | ||

2310 | assert(z->n_allocd > 0); | ||

2311 | assert(z->limbs != NULL); | ||

2312 | assert(failure != NULL); | ||

2313 | assert(s != NULL); | ||

2314 | |||

2315 | //Try to parse it as a simple integer. | ||

2316 | if (BSTRFUNC_is_sint_wo_commas(s)) | ||

2317 | { | ||

2318 | GMP_INTS_mpz_set_simple_char_str(z, s); | ||

2319 | *failure = 0; | ||

2320 | return; | ||

2321 | } | ||

2322 | //If that didn't work, try to parse it as a simple | ||

2323 | //integer with commas. | ||

2324 | else if (BSTRFUNC_is_sint_w_commas(s)) | ||

2325 | { | ||

2326 | GMP_INTS_mpz_set_simple_char_str(z, s); | ||

2327 | *failure = 0; | ||

2328 | return; | ||

2329 | } | ||

2330 | //If neither of those worked, try to parse it as | ||

2331 | //something containing scientific notation. | ||

2332 | else | ||

2333 | { | ||

2334 | GMP_INTS_mpz_set_sci_not_num(z, failure, s); | ||

2335 | |||

2336 | if (!*failure) | ||

2337 | { | ||

2338 | //We were able to parse it that way. | ||

2339 | //Everything is set up, just return. | ||

2340 | return; | ||

2341 | } | ||

2342 | else | ||

2343 | { | ||

2344 | //We're out of options. All parsing failed. | ||

2345 | GMP_INTS_mpz_set_ui(z, 0); | ||

2346 | *failure = 1; | ||

2347 | return; | ||

2348 | } | ||

2349 | } | ||

2350 | } | ||

2351 | |||

2352 | |||

2353 | void GMP_INTS_mpz_parse_into_uint32(unsigned *result, | ||

2354 | int *failure, | ||

2355 | char *s) | ||

2356 | { | ||

2357 | GMP_INTS_mpz_struct arb_int; | ||

2358 | |||

2359 | //Eyeball the input parameters. | ||

2360 | assert(result != NULL); | ||

2361 | assert(failure != NULL); | ||

2362 | assert(s != NULL); | ||

2363 | |||

2364 | //Allocate space for the one arbitrary integer we need. | ||

2365 | GMP_INTS_mpz_init(&arb_int); | ||

2366 | |||

2367 | //Try to parse the string into an arbitrary length integer | ||

2368 | //using all methods known to man. | ||

2369 | GMP_INTS_mpz_set_general_int(&arb_int, failure, s); | ||

2370 | |||

2371 | //If the parse failed, we must declare failure and return | ||

2372 | //0. | ||

2373 | if (*failure) | ||

2374 | { | ||

2375 | *result = 0; | ||

2376 | *failure = 1; | ||

2377 | } | ||

2378 | else | ||

2379 | { | ||

2380 | //We might have success, but it might be negative or | ||

2381 | //too big. | ||

2382 | if (arb_int.size == 1) | ||

2383 | { | ||

2384 | *result = arb_int.limbs[0]; | ||

2385 | *failure = 0; | ||

2386 | } | ||

2387 | else if (arb_int.size == 0) | ||

2388 | { | ||

2389 | *result = 0; | ||

2390 | *failure = 0; | ||

2391 | } | ||

2392 | else | ||

2393 | { | ||

2394 | *result = 0; | ||

2395 | *failure = 1; | ||

2396 | } | ||

2397 | } | ||

2398 | |||

2399 | //Deallocate the arbitrary integer. | ||

2400 | GMP_INTS_mpz_clear(&arb_int); | ||

2401 | } | ||

2402 | |||

2403 | |||

2404 | void GMP_INTS_mpz_swap(GMP_INTS_mpz_struct *a, | ||

2405 | GMP_INTS_mpz_struct *b) | ||

2406 | { | ||

2407 | GMP_INTS_mpz_struct temp; | ||

2408 | |||

2409 | //Eyeball the input parameters. | ||

2410 | assert(a != NULL); | ||

2411 | assert(a->n_allocd > 0); | ||

2412 | assert(a->limbs != NULL); | ||

2413 | assert(b != NULL); | ||

2414 | assert(b->n_allocd > 0); | ||

2415 | assert(b->limbs != NULL); | ||

2416 | |||

2417 | //Make the swap via memory copy. | ||

2418 | memcpy(&temp, a, sizeof(GMP_INTS_mpz_struct)); | ||

2419 | memcpy(a, b, sizeof(GMP_INTS_mpz_struct)); | ||

2420 | memcpy(b, &temp, sizeof(GMP_INTS_mpz_struct)); | ||

2421 | } | ||

2422 | |||

2423 | |||

2424 | /******************************************************************/ | ||

2425 | /*** PUBLIC ARITHMETIC FUNCTIONS *******************************/ | ||

2426 | /******************************************************************/ | ||

2427 | |||

2428 | //07/15/01: Unit test and visual inspection passed. | ||

2429 | void GMP_INTS_mpz_add ( GMP_INTS_mpz_struct *w, | ||

2430 | const GMP_INTS_mpz_struct *u, | ||

2431 | const GMP_INTS_mpz_struct *v) | ||

2432 | { | ||

2433 | GMP_INTS_limb_srcptr up, vp; | ||

2434 | GMP_INTS_limb_ptr wp; | ||

2435 | GMP_INTS_size_t usize, vsize, wsize; | ||

2436 | GMP_INTS_size_t abs_usize; | ||

2437 | GMP_INTS_size_t abs_vsize; | ||

2438 | |||

2439 | //Look at the input parameters carefully. | ||

2440 | assert(w != NULL); | ||

2441 | assert(u != NULL); | ||

2442 | assert(v != NULL); | ||

2443 | assert(w->n_allocd > 0); | ||

2444 | assert(u->n_allocd > 0); | ||

2445 | assert(v->n_allocd > 0); | ||

2446 | assert(w->limbs != NULL); | ||

2447 | assert(u->limbs != NULL); | ||

2448 | assert(v->limbs != NULL); | ||

2449 | |||

2450 | //Handle the case of a tainted result. If either of the | ||

2451 | //two inputs are either direct overflows or tainted by | ||

2452 | //an overflow, mark the result tainted and do not perform | ||

2453 | //any arithmetic operation. | ||

2454 | { | ||

2455 | int taint; | ||

2456 | |||

2457 | taint = GMP_INTS_two_op_flags_map(u->flags, v->flags); | ||

2458 | |||

2459 | w->flags = 0; | ||

2460 | //"w" starts off with a clean slate. Must do this | ||

2461 | //after taint calculation in case locations of u or v | ||

2462 | //are the same as w. | ||

2463 | if (taint) | ||

2464 | { | ||

2465 | w->flags = taint; | ||

2466 | return; | ||

2467 | } | ||

2468 | } | ||

2469 | |||

2470 | usize = u->size; | ||

2471 | vsize = v->size; | ||

2472 | abs_usize = GMP_INTS_abs_of_size_t(usize); | ||

2473 | abs_vsize = GMP_INTS_abs_of_size_t(vsize); | ||

2474 | |||

2475 | //Arrange things so that U has at least as many | ||

2476 | //limbs as V, i.e. limbs(U) >= limbs(V); | ||

2477 | if (abs_usize < abs_vsize) | ||

2478 | { | ||

2479 | const GMP_INTS_mpz_struct *tmp_ptr; | ||

2480 | GMP_INTS_size_t tmp_size; | ||

2481 | |||

2482 | //Swap U and V. This does no harm, because we are | ||

2483 | //manipulating only local variables. This does not | ||

2484 | //affect the caller. | ||

2485 | tmp_ptr = u; | ||

2486 | u = v; | ||

2487 | v = tmp_ptr; | ||

2488 | tmp_size = usize; | ||

2489 | usize = vsize; | ||

2490 | vsize = tmp_size; | ||

2491 | tmp_size = abs_usize; | ||

2492 | abs_usize = abs_vsize; | ||

2493 | abs_vsize = tmp_size; | ||

2494 | } | ||

2495 | |||

2496 | /* True: ABS_USIZE >= ABS_VSIZE. */ | ||

2497 | |||

2498 | /* If not space for w (and possible carry), increase space. */ | ||

2499 | wsize = abs_usize + 1; | ||

2500 | if (w->n_allocd < wsize) | ||

2501 | GMP_INTS_mpz_realloc(w, wsize); | ||

2502 | |||

2503 | //These pointers must be obtained after realloc. At this point, | ||

2504 | //u or v may be the same as w. | ||

2505 | up = u->limbs; | ||

2506 | vp = v->limbs; | ||

2507 | wp = w->limbs; | ||

2508 | |||

2509 | if ((usize ^ vsize) < 0) | ||

2510 | { | ||

2511 | //U and V have different sign. Need to compare them to determine | ||

2512 | //which operand to subtract from which. | ||

2513 | |||

2514 | //This test is right since ABS_USIZE >= ABS_VSIZE. | ||

2515 | //If the equality case is ruled out, then U has more limbs | ||

2516 | //than V, which means that it is bigger in magnitude. | ||

2517 | if (abs_usize != abs_vsize) | ||

2518 | { | ||

2519 | GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize); | ||

2520 | wsize = abs_usize; | ||

2521 | |||

2522 | //Normalize the result. This was formerly a macro. | ||

2523 | //To normalize in this context means to trim the size | ||

2524 | //down to eliminate any leading zero limbs that came | ||

2525 | //about because the size of the result of an operation | ||

2526 | //was overestimated. | ||

2527 | GMP_INTS_mpn_normalize(wp, &wsize); | ||

2528 | |||

2529 | if (usize < 0) | ||

2530 | wsize = -wsize; | ||

2531 | } | ||

2532 | else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0) | ||

2533 | { | ||

2534 | GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize); | ||

2535 | wsize = abs_usize; | ||

2536 | GMP_INTS_mpn_normalize(wp, &wsize); | ||

2537 | if (usize >= 0) | ||

2538 | wsize = -wsize; | ||

2539 | } | ||

2540 | else | ||

2541 | { | ||

2542 | GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize); | ||

2543 | wsize = abs_usize; | ||

2544 | GMP_INTS_mpn_normalize(wp, &wsize); | ||

2545 | if (usize < 0) | ||

2546 | wsize = -wsize; | ||

2547 | } | ||

2548 | } | ||

2549 | else | ||

2550 | { | ||

2551 | /* U and V have same sign. Add them. */ | ||

2552 | GMP_INTS_limb_t cy_limb | ||

2553 | = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize); | ||

2554 | wp[abs_usize] = cy_limb; | ||

2555 | wsize = abs_usize + cy_limb; | ||

2556 | if (usize < 0) | ||

2557 | wsize = -wsize; | ||

2558 | } | ||

2559 | |||

2560 | w->size = wsize; | ||

2561 | |||

2562 | //Handle the case of an overflowed result. If the result | ||

2563 | //of the addition is too big or too small, mark it as | ||

2564 | //overflowed. | ||

2565 | if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2566 | { | ||

2567 | w->flags = GMP_INTS_EF_INTOVF_POS; | ||

2568 | } | ||

2569 | else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2570 | { | ||

2571 | w->flags = GMP_INTS_EF_INTOVF_NEG; | ||

2572 | } | ||

2573 | } | ||

2574 | |||

2575 | |||

2576 | //07/15/01: Unit testing skipped because of recursive | ||

2577 | //nature. Visual inspection OK. | ||

2578 | void GMP_INTS_mpz_add_ui ( GMP_INTS_mpz_struct *w, | ||

2579 | const GMP_INTS_mpz_struct *u, | ||

2580 | unsigned long int v) | ||

2581 | { | ||

2582 | //The GNU MP version of this is quite efficient, and this | ||

2583 | //makes sense since it is a common operation. However, | ||

2584 | //for simplicity just define this recursively in terms | ||

2585 | //of the ADD function. This can always be made quicker | ||

2586 | //later (by changing back to the GNU MP version). | ||

2587 | GMP_INTS_mpz_struct temp; | ||

2588 | |||

2589 | //Eyeball the inputs carefully. | ||

2590 | assert(w != NULL); | ||

2591 | assert(w->n_allocd > 0); | ||

2592 | assert(w->limbs != NULL); | ||

2593 | assert(u != NULL); | ||

2594 | assert(u->n_allocd > 0); | ||

2595 | assert(u->limbs != NULL); | ||

2596 | |||

2597 | //Create a temporary integer. | ||

2598 | GMP_INTS_mpz_init(&temp); | ||

2599 | |||

2600 | //Set the temporary integer to the value of the input | ||

2601 | //argument. | ||

2602 | GMP_INTS_mpz_set_ui(&temp, v); | ||

2603 | |||

2604 | //Do the actual addition. This recursive definition | ||

2605 | //is inherently wasteful, but I'm after clarity, not | ||

2606 | //warp speed. | ||

2607 | GMP_INTS_mpz_add(w, u, &temp); | ||

2608 | |||

2609 | //Destroy the temporary integer (this will reclaim the | ||

2610 | //memory). | ||

2611 | GMP_INTS_mpz_clear(&temp); | ||

2612 | } | ||

2613 | |||

2614 | |||

2615 | //07/15/01: Visual inspection passed. Not unit tested | ||

2616 | //because of symmetry with GMP_INTS_mpz_add(). | ||

2617 | void GMP_INTS_mpz_sub ( GMP_INTS_mpz_struct *w, | ||

2618 | const GMP_INTS_mpz_struct *u, | ||

2619 | const GMP_INTS_mpz_struct *v) | ||

2620 | { | ||

2621 | GMP_INTS_limb_srcptr up, vp; | ||

2622 | GMP_INTS_limb_ptr wp; | ||

2623 | GMP_INTS_size_t usize, vsize, wsize; | ||

2624 | GMP_INTS_size_t abs_usize; | ||

2625 | GMP_INTS_size_t abs_vsize; | ||

2626 | |||

2627 | //Look at the input parameters carefully. | ||

2628 | assert(w != NULL); | ||

2629 | assert(u != NULL); | ||

2630 | assert(v != NULL); | ||

2631 | assert(w->n_allocd > 0); | ||

2632 | assert(u->n_allocd > 0); | ||

2633 | assert(v->n_allocd > 0); | ||

2634 | assert(w->limbs != NULL); | ||

2635 | assert(u->limbs != NULL); | ||

2636 | assert(v->limbs != NULL); | ||

2637 | |||

2638 | //Handle the case of a tainted result. If either of the | ||

2639 | //two inputs are either direct overflows or tainted by | ||

2640 | //an overflow, mark the result tainted and do not perform | ||

2641 | //any arithmetic operation. | ||

2642 | { | ||

2643 | int taint; | ||

2644 | |||

2645 | taint = GMP_INTS_two_op_flags_map(u->flags, v->flags); | ||

2646 | |||

2647 | w->flags = 0; | ||

2648 | //"w" starts off with a clean slate. Must do this | ||

2649 | //after taint calculation in case locations of u or v | ||

2650 | //are the same as w. | ||

2651 | if (taint) | ||

2652 | { | ||

2653 | w->flags = taint; | ||

2654 | return; | ||

2655 | } | ||

2656 | } | ||

2657 | |||

2658 | usize = u->size; | ||

2659 | vsize = -(v->size); /* The "-" makes the difference from mpz_add */ | ||

2660 | abs_usize = GMP_INTS_abs_of_size_t(usize); | ||

2661 | abs_vsize = GMP_INTS_abs_of_size_t(vsize); | ||

2662 | |||

2663 | if (abs_usize < abs_vsize) | ||

2664 | { | ||

2665 | const GMP_INTS_mpz_struct *tmp_ptr; | ||

2666 | GMP_INTS_size_t tmp_size; | ||

2667 | |||

2668 | //Swap U and V. This does no harm, because we are | ||

2669 | //manipulating only local variables. This does not | ||

2670 | //affect the caller. | ||

2671 | tmp_ptr = u; | ||

2672 | u = v; | ||

2673 | v = tmp_ptr; | ||

2674 | tmp_size = usize; | ||

2675 | usize = vsize; | ||

2676 | vsize = tmp_size; | ||

2677 | tmp_size = abs_usize; | ||

2678 | abs_usize = abs_vsize; | ||

2679 | abs_vsize = tmp_size; | ||

2680 | } | ||

2681 | |||

2682 | /* True: ABS_USIZE >= ABS_VSIZE. */ | ||

2683 | |||

2684 | /* If not space for w (and possible carry), increase space. */ | ||

2685 | wsize = abs_usize + 1; | ||

2686 | if (w->n_allocd < wsize) | ||

2687 | GMP_INTS_mpz_realloc (w, wsize); | ||

2688 | |||

2689 | /* These must be after realloc (u or v may be the same as w). */ | ||

2690 | up = u->limbs; | ||

2691 | vp = v->limbs; | ||

2692 | wp = w->limbs; | ||

2693 | |||

2694 | if ((usize ^ vsize) < 0) | ||

2695 | { | ||

2696 | //U and V have different sign. Need to compare them to determine | ||

2697 | //which operand to subtract from which. | ||

2698 | |||

2699 | //This test is right since ABS_USIZE >= ABS_VSIZE. | ||

2700 | if (abs_usize != abs_vsize) | ||

2701 | { | ||

2702 | GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize); | ||

2703 | wsize = abs_usize; | ||

2704 | GMP_INTS_mpn_normalize(wp, &wsize); | ||

2705 | if (usize < 0) | ||

2706 | wsize = -wsize; | ||

2707 | } | ||

2708 | else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0) | ||

2709 | { | ||

2710 | GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize); | ||

2711 | wsize = abs_usize; | ||

2712 | GMP_INTS_mpn_normalize(wp, &wsize); | ||

2713 | if (usize >= 0) | ||

2714 | wsize = -wsize; | ||

2715 | } | ||

2716 | else | ||

2717 | { | ||

2718 | GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize); | ||

2719 | wsize = abs_usize; | ||

2720 | GMP_INTS_mpn_normalize (wp, &wsize); | ||

2721 | if (usize < 0) | ||

2722 | wsize = -wsize; | ||

2723 | } | ||

2724 | } | ||

2725 | else | ||

2726 | { | ||

2727 | /* U and V have same sign. Add them. */ | ||

2728 | GMP_INTS_limb_t cy_limb | ||

2729 | = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize); | ||

2730 | wp[abs_usize] = cy_limb; | ||

2731 | wsize = abs_usize + cy_limb; | ||

2732 | if (usize < 0) | ||

2733 | wsize = -wsize; | ||

2734 | } | ||

2735 | |||

2736 | w->size = wsize; | ||

2737 | |||

2738 | //Handle the case of an overflowed result. If the result | ||

2739 | //of the addition is too big or too small, mark it as | ||

2740 | //overflowed. | ||

2741 | if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2742 | { | ||

2743 | w->flags = GMP_INTS_EF_INTOVF_POS; | ||

2744 | } | ||

2745 | else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2746 | { | ||

2747 | w->flags = GMP_INTS_EF_INTOVF_NEG; | ||

2748 | } | ||

2749 | } | ||

2750 | |||

2751 | |||

2752 | //07/15/01: Unit testing skipped because of recursive | ||

2753 | //nature. Visual inspection OK. | ||

2754 | void GMP_INTS_mpz_sub_ui ( GMP_INTS_mpz_struct *w, | ||

2755 | const GMP_INTS_mpz_struct *u, | ||

2756 | unsigned long int v) | ||

2757 | { | ||

2758 | //The GNU MP version of this is quite efficient, and this | ||

2759 | //makes sense since it is a common operation. However, | ||

2760 | //for simplicity just define this recursively in terms | ||

2761 | //of the SUB function. This can always be made quicker | ||

2762 | //later (by changing back to the GNU MP version). | ||

2763 | GMP_INTS_mpz_struct temp; | ||

2764 | |||

2765 | //Eyeball the inputs carefully. | ||

2766 | assert(w != NULL); | ||

2767 | assert(w->n_allocd > 0); | ||

2768 | assert(w->limbs != NULL); | ||

2769 | assert(u != NULL); | ||

2770 | assert(u->n_allocd > 0); | ||

2771 | assert(u->limbs != NULL); | ||

2772 | |||

2773 | //Create a temporary integer. | ||

2774 | GMP_INTS_mpz_init(&temp); | ||

2775 | |||

2776 | //Set the temporary integer to the value of the input | ||

2777 | //argument. | ||

2778 | GMP_INTS_mpz_set_ui(&temp, v); | ||

2779 | |||

2780 | //Do the actual subtraction. This recursive definition | ||

2781 | //is inherently wasteful, but I'm after clarity, not | ||

2782 | //warp speed. | ||

2783 | GMP_INTS_mpz_sub(w, u, &temp); | ||

2784 | |||

2785 | //Destroy the temporary integer (this will reclaim the | ||

2786 | //memory). | ||

2787 | GMP_INTS_mpz_clear(&temp); | ||

2788 | } | ||

2789 | |||

2790 | |||

2791 | void GMP_INTS_mpz_mul ( GMP_INTS_mpz_struct *w, | ||

2792 | const GMP_INTS_mpz_struct *u, | ||

2793 | const GMP_INTS_mpz_struct *v) | ||

2794 | { | ||

2795 | GMP_INTS_size_t usize = u->size; | ||

2796 | GMP_INTS_size_t vsize = v->size; | ||

2797 | GMP_INTS_size_t wsize; | ||

2798 | GMP_INTS_size_t sign_product; | ||

2799 | GMP_INTS_limb_ptr up, vp; | ||

2800 | GMP_INTS_limb_ptr wp; | ||

2801 | GMP_INTS_limb_ptr free_me = NULL; | ||

2802 | GMP_INTS_size_t free_me_size; | ||

2803 | GMP_INTS_limb_t cy_limb; | ||

2804 | |||

2805 | //Eyeball the inputs. | ||

2806 | assert(w != NULL); | ||

2807 | assert(w->n_allocd > 0); | ||

2808 | assert(w->limbs != NULL); | ||

2809 | assert(u != NULL); | ||

2810 | assert(u->n_allocd > 0); | ||

2811 | assert(u->limbs != NULL); | ||

2812 | assert(v != NULL); | ||

2813 | assert(v->n_allocd > 0); | ||

2814 | assert(v->limbs != NULL); | ||

2815 | |||

2816 | //Handle the case of a tainted result. If either of the | ||

2817 | //two inputs are either direct overflows or tainted by | ||

2818 | //an overflow, mark the result tainted and do not perform | ||

2819 | //any arithmetic operation. | ||

2820 | { | ||

2821 | int taint; | ||

2822 | |||

2823 | taint = GMP_INTS_two_op_flags_map(u->flags, v->flags); | ||

2824 | |||

2825 | w->flags = 0; | ||

2826 | //"w" starts off with a clean slate. Must do this | ||

2827 | //after taint calculation in case locations of u or v | ||

2828 | //are the same as w. | ||

2829 | if (taint) | ||

2830 | { | ||

2831 | w->flags = taint; | ||

2832 | return; | ||

2833 | } | ||

2834 | } | ||

2835 | |||

2836 | sign_product = usize ^ vsize; | ||

2837 | usize = GMP_INTS_abs_of_size_t(usize); | ||

2838 | vsize = GMP_INTS_abs_of_size_t(vsize); | ||

2839 | |||

2840 | //Handle the case of a certain result overflow (why do the math when | ||

2841 | //the result is certain?). In general, when multiplying two inputs | ||

2842 | //whose sizes are M limbs and N limbs, the size of the result will be | ||

2843 | //either M+N or M+N-1 limbs. If M+N-1 > MAX_ALLOWED, then can declare | ||

2844 | //an early overflow. | ||

2845 | if ((usize + vsize - 1) > GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2846 | { | ||

2847 | if (sign_product < 0) | ||

2848 | w->flags = GMP_INTS_EF_INTOVF_NEG; | ||

2849 | else | ||

2850 | w->flags = GMP_INTS_EF_INTOVF_POS; | ||

2851 | |||

2852 | return; | ||

2853 | } | ||

2854 | |||

2855 | |||

2856 | if (usize < vsize) | ||

2857 | { | ||

2858 | //Temporary variables just for the swap. | ||

2859 | const GMP_INTS_mpz_struct *tmp_ptr; | ||

2860 | GMP_INTS_size_t tmp_size; | ||

2861 | |||

2862 | //Swap U and V. | ||

2863 | tmp_ptr = u; | ||

2864 | u = v; | ||

2865 | v = tmp_ptr; | ||

2866 | tmp_size = usize; | ||

2867 | usize = vsize; | ||

2868 | vsize = tmp_size; | ||

2869 | } | ||

2870 | |||

2871 | //Grab pointers to the arrays of limbs. | ||

2872 | up = u->limbs; | ||

2873 | vp = v->limbs; | ||

2874 | wp = w->limbs; | ||

2875 | |||

2876 | /* Ensure W has space enough to store the result. */ | ||

2877 | wsize = usize + vsize; | ||

2878 | if (w->n_allocd < wsize) | ||

2879 | { | ||

2880 | if (wp == up || wp == vp) | ||

2881 | { | ||

2882 | free_me = wp; | ||

2883 | free_me_size = w->n_allocd; | ||

2884 | } | ||

2885 | else | ||

2886 | { | ||

2887 | GMP_INTS_free_w_size (wp, w->n_allocd * sizeof(GMP_INTS_limb_t)); | ||

2888 | } | ||

2889 | |||

2890 | w->n_allocd = wsize; | ||

2891 | wp = (GMP_INTS_limb_ptr) | ||

2892 | GMP_INTS_malloc (wsize * sizeof(GMP_INTS_limb_t)); | ||

2893 | w->limbs = wp; | ||

2894 | } | ||

2895 | else | ||

2896 | { | ||

2897 | /* Make U and V not overlap with W. */ | ||

2898 | if (wp == up) | ||

2899 | { | ||

2900 | /* W and U are identical. Allocate temporary space for U. */ | ||

2901 | up = (GMP_INTS_limb_ptr) | ||

2902 | _alloca(usize * sizeof(GMP_INTS_limb_t)); | ||

2903 | /* Is V identical too? Keep it identical with U. */ | ||

2904 | if (wp == vp) | ||

2905 | vp = up; | ||

2906 | /* Copy to the temporary space. */ | ||

2907 | GMP_INTS_mpn_copy_limbs(up, wp, usize); | ||

2908 | } | ||

2909 | else if (wp == vp) | ||

2910 | { | ||

2911 | /* W and V are identical. Allocate temporary space for V. */ | ||

2912 | vp = (GMP_INTS_limb_ptr) | ||

2913 | _alloca(vsize * sizeof(GMP_INTS_limb_t)); | ||

2914 | /* Copy to the temporary space. */ | ||

2915 | GMP_INTS_mpn_copy_limbs(vp, wp, vsize); | ||

2916 | } | ||

2917 | } | ||

2918 | |||

2919 | if (vsize == 0) | ||

2920 | { | ||

2921 | wsize = 0; | ||

2922 | } | ||

2923 | else | ||

2924 | { | ||

2925 | cy_limb = GMP_INTS_mpn_mul (wp, up, usize, vp, vsize); | ||

2926 | wsize = usize + vsize; | ||

2927 | wsize -= cy_limb == 0; | ||

2928 | } | ||

2929 | |||

2930 | w->size = sign_product < 0 ? -wsize : wsize; | ||

2931 | |||

2932 | if (free_me != NULL) | ||

2933 | GMP_INTS_free_w_size (free_me, free_me_size * sizeof(GMP_INTS_limb_t)); | ||

2934 | |||

2935 | //Final check for overflow. | ||

2936 | if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2937 | w->flags = GMP_INTS_EF_INTOVF_POS; | ||

2938 | else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT) | ||

2939 | w->flags = GMP_INTS_EF_INTOVF_NEG; | ||

2940 | } | ||

2941 | |||

2942 | |||

2943 | //07/15/01: Unit testing skipped because of recursive | ||

2944 | //nature. Visual inspection OK. | ||

2945 | void GMP_INTS_mpz_mul_si ( GMP_INTS_mpz_struct *w, | ||

2946 | const GMP_INTS_mpz_struct *u, | ||

2947 | long int v) | ||

2948 | { | ||

2949 | GMP_INTS_mpz_struct temp; | ||

2950 | |||

2951 | //Eyeball the inputs carefully. | ||

2952 | assert(w != NULL); | ||

2953 | assert(w->n_allocd > 0); | ||

2954 | assert(w->limbs != NULL); | ||

2955 | assert(u != NULL); | ||

2956 | assert(u->n_allocd > 0); | ||

2957 | assert(u->limbs != NULL); | ||

2958 | |||

2959 | //Create a temporary integer. | ||

2960 | GMP_INTS_mpz_init(&temp); | ||

2961 | |||

2962 | //Set the temporary integer to the value of the input | ||

2963 | //argument. | ||

2964 | GMP_INTS_mpz_set_si(&temp, v); | ||

2965 | |||

2966 | //Do the actual multiplication. This recursive definition | ||

2967 | //is inherently wasteful, but I'm after clarity, not | ||

2968 | //warp speed. | ||

2969 | GMP_INTS_mpz_mul(w, u, &temp); | ||

2970 | |||

2971 | //Destroy the temporary integer (this will reclaim the | ||

2972 | //memory). | ||

2973 | GMP_INTS_mpz_clear(&temp); | ||

2974 | } | ||

2975 | |||

2976 | |||

2977 | //07/15/01: Unit testing skipped because of recursive | ||

2978 | //nature. Visual inspection OK. | ||

2979 | void GMP_INTS_mpz_mul_ui ( GMP_INTS_mpz_struct *w, | ||

2980 | const GMP_INTS_mpz_struct *u, | ||

2981 | unsigned long int v) | ||

2982 | { | ||

2983 | GMP_INTS_mpz_struct temp; | ||

2984 | |||

2985 | //Eyeball the inputs carefully. | ||

2986 | assert(w != NULL); | ||

2987 | assert(w->size >= 0); | ||

2988 | assert(w->limbs != NULL); | ||

2989 | assert(u != NULL); | ||

2990 | assert(u->size >= 0); | ||

2991 | assert(u->limbs != NULL); | ||

2992 | |||

2993 | //Create a temporary integer. | ||

2994 | GMP_INTS_mpz_init(&temp); | ||

2995 | |||

2996 | //Set the temporary integer to the value of the input | ||

2997 | //argument. | ||

2998 | GMP_INTS_mpz_set_ui(&temp, v); | ||

2999 | |||

3000 | //Do the actual multiplication. This recursive definition | ||

3001 | //is inherently wasteful, but I'm after clarity, not | ||

3002 | //warp speed. | ||

3003 | GMP_INTS_mpz_mul(w, u, &temp); | ||

3004 | |||

3005 | //Destroy the temporary integer (this will reclaim the | ||

3006 | //memory). | ||

3007 | GMP_INTS_mpz_clear(&temp); | ||

3008 | } | ||

3009 | |||

3010 | |||

3011 | void GMP_INTS_mpz_tdiv_qr ( GMP_INTS_mpz_struct *quot, | ||

3012 | GMP_INTS_mpz_struct *rem, | ||

3013 | const GMP_INTS_mpz_struct *num, | ||

3014 | const GMP_INTS_mpz_struct *den) | ||

3015 | { | ||

3016 | GMP_INTS_size_t abs_num_size, | ||

3017 | abs_den_size, | ||

3018 | quotient_sign, | ||

3019 | remainder_sign, | ||

3020 | numerator_bitsize, | ||

3021 | denominator_bitsize, | ||

3022 | division_loop_count, | ||

3023 | division_loop_count_mod_32, | ||

3024 | division_loop_count_div_32, | ||

3025 | division_counter, | ||

3026 | i; | ||

3027 | GMP_INTS_limb_t temp_limb; | ||

3028 | GMP_INTS_limb_ptr trial_divisor; | ||

3029 | |||

3030 | //Eyeball the input parameters. | ||

3031 | assert(quot != NULL); | ||

3032 | assert(quot->n_allocd > 0); | ||

3033 | assert(quot->limbs != NULL); | ||

3034 | assert(rem != NULL); | ||

3035 | assert(rem->n_allocd > 0); | ||

3036 | assert(rem->limbs != NULL); | ||

3037 | assert(num != NULL); | ||

3038 | assert(num->n_allocd > 0); | ||

3039 | assert(num->limbs != NULL); | ||

3040 | assert(den != NULL); | ||

3041 | assert(den->n_allocd > 0); | ||

3042 | assert(den->limbs != NULL); | ||

3043 | |||

3044 | //We require for this function that the numerator, denominator, quotient, and | ||

3045 | //remainder all be distinct. | ||

3046 | assert(quot != rem); | ||

3047 | assert(quot != num); | ||

3048 | assert(quot != den); | ||

3049 | assert(rem != num); | ||

3050 | assert(rem != den); | ||

3051 | assert(num != den); | ||

3052 | |||

3053 | //The GNU code was probably very efficient, but exceeded | ||

3054 | //my abilities to analyze. This is the classic | ||

3055 | //division algorithm. | ||

3056 | |||

3057 | //First, start off with the quotient and remainder having | ||

3058 | //no error flags set. These will be set if appropriate. | ||

3059 | quot->flags = 0; | ||

3060 | rem->flags = 0; | ||

3061 | |||

3062 | //First, handle tainted inputs. If the numerator or denominator | ||

3063 | //are bad or tainted, the quotient and remainder get tainted | ||

3064 | //automatically. | ||

3065 | { | ||

3066 | int taint; | ||

3067 | |||

3068 | taint = GMP_INTS_two_op_flags_map(num->flags, den->flags); | ||

3069 | |||

3070 | if (taint) | ||

3071 | { | ||

3072 | quot->flags = taint; | ||

3073 | rem->flags = taint; | ||

3074 | return; | ||

3075 | } | ||

3076 | } | ||

3077 | |||

3078 | //The second possible cause for taint is if the divisor is | ||

3079 | //zero. This will get both the value of positive overflow. | ||

3080 | if (den->size == 0) | ||

3081 | { | ||

3082 | quot->flags = GMP_INTS_EF_INTOVF_POS; | ||

3083 | rem->flags = GMP_INTS_EF_INTOVF_POS; | ||

3084 | return; | ||

3085 | } | ||

3086 | |||

3087 | //Handle the special case of a numerator of zero. If the numerator | ||

3088 | //is zero, the quotient and remainder are zero automatically. | ||

3089 | if (num->size == 0) | ||

3090 | { | ||

3091 | GMP_INTS_mpz_set_ui(quot, 0); | ||

3092 | GMP_INTS_mpz_set_ui(rem, 0); | ||

3093 | return; | ||

3094 | } | ||

3095 | |||

3096 | //Generally, nothing else can go wrong as far as taint. The | ||

3097 | //value of the quotient is confined to be no larger than the | ||

3098 | //numerator, and the value of the remainder is confined to | ||

3099 | //be no larger than denominator-1. So, generally, if the | ||

3100 | //inputs are in size bounds, the outputs will be also. | ||

3101 | |||

3102 | //Figure out how large in limbs the numerator and denominator actually | ||

3103 | //are. | ||

3104 | abs_num_size = GMP_INTS_abs_of_size_t(num->size); | ||

3105 | abs_den_size = GMP_INTS_abs_of_size_t(den->size); | ||

3106 | |||

3107 | //Figure out the sign of things. We want the following relationship | ||

3108 | //to be true: | ||

3109 | // num/den = quot + rem/den. | ||

3110 | //The way to achieve this is to assign the sign of the quotient in the traditional | ||

3111 | //way, then to assign the remainder to have the same sign as the numerator. | ||

3112 | quotient_sign = num->size ^ den->size; | ||

3113 | remainder_sign = num->size; | ||

3114 | |||

3115 | //The remainder starts off with the absolute value of the numerator, and then | ||

3116 | //we subtract from it as part of the division loop. | ||

3117 | GMP_INTS_mpz_copy(rem, num); | ||

3118 | //We know after the copy that the amount of space allocated in the remainder | ||

3119 | //MUST be at least as large as the absolute value of the numerator. So from | ||

3120 | //this point forward we use the space. | ||

3121 | assert(rem->n_allocd >= abs_num_size); | ||

3122 | |||

3123 | //Figure out the number of significant bits in the numerator and denominator. | ||

3124 | //This determines the loop count over which we do the shift division loop. | ||

3125 | numerator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_num_size; | ||

3126 | |||

3127 | i = abs_num_size - 1; | ||

3128 | |||

3129 | //We need to be extra careful here. One failure mode is that an integer | ||

3130 | //data structure is corrupted and the "size" field reflects limbs | ||

3131 | //that are zero. Need to watch that this kind of failure doesn't | ||

3132 | //cause memory access errors. | ||

3133 | assert(num->limbs[i] != 0); | ||

3134 | if (num->limbs[i] == 0) | ||

3135 | { | ||

3136 | quot->flags = GMP_INTS_EF_INTOVF_POS; | ||

3137 | rem->flags = GMP_INTS_EF_INTOVF_POS; | ||

3138 | return; | ||

3139 | } | ||

3140 | |||

3141 | temp_limb = 0x80000000; | ||

3142 | |||

3143 | while (((num->limbs[i] & temp_limb) == 0) && (temp_limb != 0)) | ||

3144 | { | ||

3145 | numerator_bitsize--; | ||

3146 | temp_limb >>= 1; | ||

3147 | } | ||

3148 | |||

3149 | denominator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_den_size; | ||

3150 | |||

3151 | i = abs_den_size - 1; | ||

3152 | |||

3153 | //We need to be extra careful here. One failure mode is that an integer | ||

3154 | //data structure is corrupted and the "size" field reflects limbs | ||

3155 | //that are zero. Need to watch that this kind of failure doesn't | ||

3156 | //cause memory access errors. | ||

3157 | assert(den->limbs[i] != 0); | ||

3158 | if (den->limbs[i] == 0) | ||

3159 | { | ||

3160 | quot->flags = GMP_INTS_EF_INTOVF_POS; | ||

3161 | rem->flags = GMP_INTS_EF_INTOVF_POS; | ||

3162 | return; | ||

3163 | } | ||

3164 | |||

3165 | temp_limb = 0x80000000; | ||

3166 | |||

3167 | while (((den->limbs[i] & temp_limb) == 0) && (temp_limb != 0)) | ||

3168 | { | ||

3169 | denominator_bitsize--; | ||

3170 | temp_limb >>= 1; | ||

3171 | } | ||

3172 | |||

3173 | //The quotient starts off with the value of zero, but we consistently may | ||

3174 | //mask 1 into it and shift left. We need to be sure that we have as much | ||

3175 | //shift space there as is in the numerator. For this purpose we need to | ||

3176 | //prepare a block of clear memory as large as the numerator's. | ||

3177 | if (quot->n_allocd < abs_num_size) | ||

3178 | { | ||

3179 | GMP_INTS_mpz_realloc(quot, abs_num_size); //Make it big enough. | ||

3180 | } | ||

3181 | //Now, zero the memory. | ||

3182 | for (i=0; i<abs_num_size; i++) | ||

3183 | quot->limbs[i] = 0; | ||

3184 | |||

3185 | //Determine the division loop count. This is the difference | ||

3186 | //in bit sizes between the numerator and denominator. It is | ||

3187 | //possible for this number to be negative, which means that the | ||

3188 | //main division loop will be executed zero times. This gives the | ||

3189 | //right results. | ||

3190 | division_loop_count = numerator_bitsize - denominator_bitsize; | ||

3191 | |||

3192 | //We need to calculate some important numbers from the division loop | ||

3193 | //count. We need to know this number MOD 32 (which tells how far to | ||

3194 | //shift the divisor bitwise to line up with the numerator), and we | ||

3195 | //also need this number DIV 32 for the limb-wise shift. | ||

3196 | division_loop_count_mod_32 = division_loop_count % 32; | ||

3197 | division_loop_count_div_32 = division_loop_count / 32; | ||

3198 | |||

3199 | //We now need a shift register in which we shift the denominator up | ||

3200 | //for repeated comparisons. We should dynamically allocate this to | ||

3201 | //be the same size as the numerator. Using _alloca() is OK, as one | ||

3202 | //of the unit tests is to be sure that _alloca() will handle integer | ||

3203 | //of the maximum allowed size. | ||

3204 | trial_divisor = _alloca(abs_num_size * sizeof(GMP_INTS_limb_t)); | ||

3205 | |||

3206 | //Our trial divisor needs to start off with the divisor shifted up | ||

3207 | //so that the most significant bit is aligned with the numerator. | ||

3208 | for (i = 0; i < abs_num_size; i++) | ||

3209 | { | ||

3210 | if ((division_loop_count < 0) || (i < division_loop_count_div_32)) | ||

3211 | { | ||

3212 | trial_divisor[i] = 0; | ||

3213 | } | ||

3214 | else | ||

3215 | { | ||

3216 | if ((i-division_loop_count_div_32) < abs_den_size) | ||

3217 | trial_divisor[i] = den->limbs[i - division_loop_count_div_32]; | ||

3218 | else | ||

3219 | trial_divisor[i] = 0; | ||

3220 | } | ||

3221 | } | ||

3222 | |||

3223 | //The code above planted the limbs in the right place. Now need to shift bits | ||

3224 | //upward by the remaining number. | ||

3225 | if ((division_loop_count > 0) && (division_loop_count_mod_32 > 0)) | ||

3226 | { | ||

3227 | //There is an existing function we can call to do the left shift. | ||

3228 | GMP_INTS_mpn_lshift(trial_divisor, | ||

3229 | trial_divisor, | ||

3230 | abs_num_size, | ||

3231 | division_loop_count_mod_32); | ||

3232 | } | ||

3233 | |||

3234 | |||

3235 | //Everything is ready to go. Now begin the division loop itself. It is possible | ||

3236 | //for the loop to execute zero times, which will happen if the denominator is longer | ||

3237 | //in bits than the numerator. In such cases, we can't execute this loop even once | ||

3238 | //because the math assumes that the numerator is at least as long as the denominator. | ||

3239 | for (division_counter = 0; division_counter < division_loop_count+1; division_counter++) | ||

3240 | { | ||

3241 | //Shift the quotient left one bit. | ||

3242 | GMP_INTS_mpn_lshift(quot->limbs, | ||

3243 | quot->limbs, | ||

3244 | abs_num_size, | ||

3245 | 1); | ||

3246 | |||

3247 | //If the remainder is at least as large as the trial divisor, subtract the trial | ||

3248 | //divisor from the remainder and mask in the quotient. | ||

3249 | if (GMP_INTS_mpn_cmp(rem->limbs, | ||

3250 | trial_divisor, | ||

3251 | abs_num_size) >= 0) | ||

3252 | { | ||

3253 | GMP_INTS_mpn_sub(rem->limbs, | ||

3254 | rem->limbs, | ||

3255 | abs_num_size, | ||

3256 | trial_divisor, | ||

3257 | abs_num_size); | ||

3258 | quot->limbs[0] |= 1; | ||

3259 | } | ||

3260 | |||

3261 | //Shift the trial divisor right one bit. | ||

3262 | GMP_INTS_mpn_rshift(trial_divisor, | ||

3263 | trial_divisor, | ||

3264 | abs_num_size, | ||

3265 | 1); | ||

3266 | } //End for each iteration of the division loop. | ||

3267 | |||

3268 | //Normalize the quotient and the remainder. The normalization | ||

3269 | //process is to bring the sizes down if we have leading | ||

3270 | //zeros. | ||

3271 | quot->size = abs_num_size; | ||

3272 | GMP_INTS_mpn_normalize(quot->limbs, &(quot->size)); | ||

3273 | rem->size = abs_num_size; | ||

3274 | GMP_INTS_mpn_normalize(rem->limbs, &(rem->size)); | ||

3275 | |||

3276 | //Adjust the signs as required. | ||

3277 | if (quotient_sign < 0) | ||

3278 | quot->size = -(quot->size); | ||

3279 | if (remainder_sign < 0) | ||

3280 | rem->size = -(rem->size); | ||

3281 | } | ||

3282 | |||

3283 | |||

3284 | void GMP_INTS_mpz_fac_ui(GMP_INTS_mpz_struct *result, | ||

3285 | unsigned long int n) | ||

3286 | { | ||

3287 | //Just multiply the numbers in ascending order. The original | ||

3288 | //GNU library contained a much more elegant algorithm, but | ||

3289 | //this is more direct. | ||

3290 | |||

3291 | unsigned long int k; | ||

3292 | |||

3293 | GMP_INTS_mpz_set_ui (result, 1L); | ||

3294 | |||

3295 | for (k = 2; (k <= n) && !(result->flags); k++) | ||

3296 | GMP_INTS_mpz_mul_ui (result, result, k); | ||

3297 | } | ||

3298 | |||

3299 | |||

3300 | /******************************************************************/ | ||

3301 | /*** PUBLIC CONVERSION AND OUTPUT FUNCTIONS ********************/ | ||

3302 | /******************************************************************/ | ||

3303 | //07/18/01: Visual inspection OK. Function returns | ||

3304 | //reasonable values even out to 100,000 digits--seems OK. | ||

3305 | int GMP_INTS_mpz_size_in_base_10(const GMP_INTS_mpz_struct *arg) | ||

3306 | { | ||

3307 | _int64 n; | ||

3308 | |||

3309 | //Eyeball the input parameter. | ||

3310 | assert(arg != NULL); | ||

3311 | assert(arg->n_allocd > 0); | ||

3312 | assert(arg->limbs != NULL); | ||

3313 | |||

3314 | //Get the number of limbs occupied by the integer. | ||

3315 | //Because even the digit zero takes some space, | ||

3316 | //don't accept zero for an answer. | ||

3317 | n = GMP_INTS_abs_of_size_t(arg->size); | ||

3318 | if (n==0) | ||

3319 | n = 1; | ||

3320 | |||

3321 | //Convert this to the number of bits. Generously | ||

3322 | //ignore any unused leading bits. | ||

3323 | n *= 32; | ||

3324 | |||

3325 | //Used a slightly high best rational approximation in F_{65535} | ||

3326 | //to go from the number of bits to the number of | ||

3327 | //digits. The division discards, so bump the result | ||

3328 | //up by 1 to compensate for possible truncation. The number | ||

3329 | //we are aproximating is ln(2)/ln(10). | ||

3330 | n *= 12655; | ||

3331 | n /= 42039; | ||

3332 | n++; | ||

3333 | |||

3334 | //Compensate for possible commas in the result. Again, | ||

3335 | //consider truncation. | ||

3336 | n *= 4; | ||

3337 | n /= 3; | ||

3338 | n++; | ||

3339 | |||

3340 | //Compensate for the minus sign, the trailing zero, | ||

3341 | //cosmic rays striking the computer from the martian | ||

3342 | //listening post camoflaged on the moon, and the | ||

3343 | //possibility that we might need to put text in the | ||

3344 | //string if any flag is set. | ||

3345 | n += 100; | ||

3346 | |||

3347 | //And that should be a good return value. | ||

3348 | return((int) n); | ||

3349 | } | ||

3350 | |||

3351 | |||

3352 | //07/19/01: Visual inspection and unit test is OK. | ||

3353 | void GMP_INTS_mpz_to_string(char *out, | ||

3354 | const GMP_INTS_mpz_struct *in) | ||

3355 | { | ||

3356 | //Eyeball the input parameters. | ||

3357 | assert(out != NULL); | ||

3358 | assert(in != NULL); | ||

3359 | assert(in->n_allocd > 0); | ||

3360 | assert(in->limbs != NULL); | ||

3361 | |||

3362 | //If any of the flags are set, stuff in the text. | ||

3363 | if (in->flags) | ||

3364 | { | ||

3365 | if (in->flags & GMP_INTS_EF_INTOVF_POS) | ||

3366 | { | ||

3367 | strcpy(out, GMP_INTS_EF_INTOVF_POS_STRING); | ||

3368 | } | ||

3369 | else if (in->flags & GMP_INTS_EF_INTOVF_NEG) | ||

3370 | { | ||

3371 | strcpy(out, GMP_INTS_EF_INTOVF_NEG_STRING); | ||

3372 | } | ||

3373 | else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_POS) | ||

3374 | { | ||

3375 | strcpy(out, GMP_INTS_EF_INTOVF_TAINT_POS_STRING); | ||

3376 | } | ||

3377 | else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_NEG) | ||

3378 | { | ||

3379 | strcpy(out, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING); | ||

3380 | } | ||

3381 | else | ||

3382 | { | ||

3383 | strcpy(out, "INTERNAL_ERROR"); | ||

3384 | } | ||

3385 | } | ||

3386 | else | ||

3387 | { | ||

3388 | //Ordinary integer conversion. | ||

3389 | GMP_INTS_mpz_struct num, den, quot, rem, k10; | ||

3390 | |||

3391 | //Allocate space for the temporary integers. | ||

3392 | GMP_INTS_mpz_init(&num); | ||

3393 | GMP_INTS_mpz_init(&den); | ||

3394 | GMP_INTS_mpz_init("); | ||

3395 | GMP_INTS_mpz_init(&rem); | ||

3396 | GMP_INTS_mpz_init(&k10); | ||

3397 | |||

3398 | //Assign the constant 10. | ||

3399 | GMP_INTS_mpz_set_ui(&k10, 10); | ||

3400 | |||

3401 | //If the integer is zero, assign that. | ||

3402 | if (in->size == 0) | ||

3403 | { | ||

3404 | strcpy(out, "0"); | ||

3405 | } | ||

3406 | else | ||

3407 | { | ||

3408 | //We have to do a full conversion. The algorithm | ||

3409 | //is division by 10, each time obtaining the least | ||

3410 | //significant digit, until finally the quotient is | ||

3411 | //zero. | ||

3412 | char *ptr; | ||

3413 | |||

3414 | ptr = out; | ||

3415 | |||

3416 | GMP_INTS_mpz_copy(&num, in); | ||

3417 | GMP_INTS_mpz_copy(&den, &k10); | ||

3418 | do | ||

3419 | { | ||

3420 | #if 0 | ||

3421 | printf("Values before division:\n"); | ||

3422 | FCMIOF_hline(); | ||

3423 | GMP_INTS_mpz_print_int(stdout, &num, "Numerator"); | ||

3424 | FCMIOF_hline(); | ||

3425 | GMP_INTS_mpz_print_int(stdout, &den, "Denominator"); | ||

3426 | FCMIOF_hline(); | ||

3427 | GMP_INTS_mpz_print_int(stdout, ", "Quotient"); | ||

3428 | FCMIOF_hline(); | ||

3429 | GMP_INTS_mpz_print_int(stdout, &rem, "Remainder"); | ||

3430 | FCMIOF_hline(); | ||

3431 | |||

3432 | if (num.size > 1) | ||

3433 | FCMIOF_hline(); | ||

3434 | #endif | ||

3435 | |||

3436 | GMP_INTS_mpz_tdiv_qr(", &rem, &num, &den); | ||

3437 | #if 0 | ||

3438 | printf("Values after division:\n"); | ||

3439 | FCMIOF_hline(); | ||

3440 | GMP_INTS_mpz_print_int(stdout, &num, "Numerator"); | ||

3441 | FCMIOF_hline(); | ||

3442 | GMP_INTS_mpz_print_int(stdout, &den, "Denominator"); | ||

3443 | FCMIOF_hline(); | ||

3444 | GMP_INTS_mpz_print_int(stdout, ", "Quotient"); | ||

3445 | FCMIOF_hline(); | ||

3446 | GMP_INTS_mpz_print_int(stdout, &rem, "Remainder"); | ||

3447 | FCMIOF_hline(); | ||

3448 | #endif | ||

3449 | |||

3450 | if (rem.size != 0) | ||

3451 | { | ||

3452 | *ptr = '0' + (char)(rem.limbs[0]); | ||

3453 | } | ||

3454 | else | ||

3455 | { | ||

3456 | *ptr = '0'; | ||

3457 | } | ||

3458 | ptr++; | ||

3459 | GMP_INTS_mpz_copy(&num, "); | ||

3460 | //printf("digit\n"); | ||

3461 | } | ||

3462 | while (!GMP_INTS_mpz_is_zero(")); | ||

3463 | |||

3464 | //Finally, if the input was negative, tack on the | ||

3465 | //minus sign. | ||

3466 | if (GMP_INTS_mpz_is_neg(in)) | ||

3467 | { | ||

3468 | *ptr = '-'; | ||

3469 | ptr++; | ||

3470 | } | ||

3471 | |||

3472 | //Finally, tack on the trailing zero terminator. | ||

3473 | *ptr = 0; | ||

3474 | ptr++; | ||

3475 | |||

3476 | //Reverse the string. | ||

3477 | BSTRFUNC_str_reverse(out); | ||

3478 | } | ||

3479 | |||

3480 | //Deallocate the integers. | ||

3481 | GMP_INTS_mpz_clear(&num); | ||

3482 | GMP_INTS_mpz_clear(&den); | ||

3483 | GMP_INTS_mpz_clear("); | ||

3484 | GMP_INTS_mpz_clear(&rem); | ||

3485 | GMP_INTS_mpz_clear(&k10); | ||

3486 | } | ||

3487 | } | ||

3488 | |||

3489 | |||

3490 | void GMP_INTS_mpz_long_int_format_to_stream(FILE *s, | ||

3491 | const GMP_INTS_mpz_struct *i, | ||

3492 | const char *desc) | ||

3493 | { | ||

3494 | int line_len; | ||

3495 | int digits_per_line; | ||

3496 | char *digits; | ||

3497 | int num_digits; | ||

3498 | int nlines; | ||

3499 | int cur_line; | ||

3500 | int number_desc_width; | ||

3501 | |||

3502 | //Eyeball the inputs, make sure the caller isn't doing | ||

3503 | //something stupid. | ||

3504 | assert(s != NULL); | ||

3505 | assert(i != NULL); | ||

3506 | assert(i->n_allocd > 0); | ||

3507 | assert(i->limbs != NULL); | ||

3508 | assert(desc != NULL); | ||

3509 | |||

3510 | //Obtain the line length assumed for formatted output. | ||

3511 | line_len = FCMIOF_get_line_len(); | ||

3512 | |||

3513 | //The description width allowed is 20. | ||

3514 | number_desc_width = 20; | ||

3515 | |||

3516 | /* The number of digits per line that we assume must be a multiple of | ||

3517 | ** three. The formula below was not examined very carefully, but it | ||

3518 | ** works fine for a line length of 78. If line length is changed, | ||

3519 | ** this formula may need to be examined very carefully and rewritten. | ||

3520 | */ | ||

3521 | digits_per_line = INTFUNC_max(3, ((((line_len-42)*3)/4)/3)*3); | ||

3522 | assert(digits_per_line >= 3); | ||

3523 | |||

3524 | /* We now need to get a digit string corresponding to this | ||

3525 | ** number. First, need to figure out how much and | ||

3526 | ** allocate the space. | ||

3527 | */ | ||

3528 | digits = GMP_INTS_malloc(GMP_INTS_mpz_size_in_base_10(i) * sizeof(char)); | ||

3529 | GMP_INTS_mpz_to_string(digits, i); | ||

3530 | |||

3531 | //If the number is negative, delete the leading minus sign. | ||

3532 | //The rest of the display algorithm needs an unsigned | ||

3533 | //series of digits. | ||

3534 | if (*digits == '-') | ||

3535 | { | ||

3536 | int i = 0; | ||

3537 | |||

3538 | do | ||

3539 | { | ||

3540 | digits[i] = digits[i+1]; | ||

3541 | i++; | ||

3542 | } | ||

3543 | while(digits[i-1]); | ||

3544 | } | ||

3545 | |||

3546 | //Figure out how many digits in the string representation. | ||

3547 | num_digits = strlen(digits); | ||

3548 | |||

3549 | /* As the first order of business, figure out how many lines the beast | ||

3550 | ** will require. | ||

3551 | */ | ||

3552 | if (i->flags) | ||

3553 | { | ||

3554 | nlines = 1; /* Only one line required for NAN verbeage. */ | ||

3555 | } | ||

3556 | else if (GMP_INTS_mpz_is_zero(i)) | ||

3557 | { | ||

3558 | nlines = 1; /* The zero value requires one line. */ | ||

3559 | } | ||

3560 | else | ||

3561 | { | ||

3562 | /* In any other case, have a formula. | ||

3563 | */ | ||

3564 | nlines = 1 + (num_digits - 1) / digits_per_line; | ||

3565 | } | ||

3566 | |||

3567 | /* Iterate through each line, spitting out whatever is appropriate. */ | ||

3568 | for (cur_line = 0; cur_line < nlines; cur_line++) | ||

3569 | { | ||

3570 | int cur_digit_on_line; | ||

3571 | |||

3572 | /* If this is the first line, spit out the description, right-aligned. | ||

3573 | ** Otherwise, spit spaces. | ||

3574 | */ | ||

3575 | if (!cur_line) | ||

3576 | { | ||

3577 | /* First line. */ | ||

3578 | int len; | ||

3579 | |||

3580 | len = strlen(desc); | ||

3581 | |||

3582 | if (len <= number_desc_width) | ||

3583 | { | ||

3584 | /* Description is shorter or equal, pad on left. */ | ||

3585 | FCMIOF_stream_repchar(s, ' ', number_desc_width - len); | ||

3586 | fprintf(s, "%s", desc); | ||

3587 | } | ||

3588 | else | ||

3589 | { | ||

3590 | /* Description is too long, truncate. */ | ||

3591 | int i; | ||

3592 | |||

3593 | for (i=0; i<number_desc_width; i++) | ||

3594 | fprintf(s, "%c", desc[i]); | ||

3595 | } | ||

3596 | |||

3597 | fprintf(s, ": "); | ||

3598 | |||

3599 | /* If the number is negative, throw in a minus sign. */ | ||

3600 | if (GMP_INTS_mpz_is_neg(i) && !(i->flags)) | ||

3601 | { | ||

3602 | fprintf(s, "- "); | ||

3603 | } | ||

3604 | else | ||

3605 | { | ||

3606 | fprintf(s, " "); | ||

3607 | } | ||

3608 | } | ||

3609 | else | ||

3610 | { | ||

3611 | /* Every line but first line. */ | ||

3612 | FCMIOF_stream_repchar(s, ' ', number_desc_width+4); | ||

3613 | } | ||

3614 | |||

3615 | for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++) | ||

3616 | { | ||

3617 | int idx_into_string; | ||

3618 | /* Index into the string which is our digit of interest. | ||

3619 | */ | ||

3620 | |||

3621 | /* Compute the index. The equation is based on the ordering | ||

3622 | ** of presentation, for example, | ||

3623 | ** | ||

3624 | ** 7 6 | ||

3625 | ** 5 4 3 | ||

3626 | ** 2 1 0. | ||

3627 | ** | ||

3628 | ** With a little thought, the equation should make sense. | ||

3629 | ** The index won't always be used to index into the string. | ||

3630 | */ | ||

3631 | idx_into_string = | ||

3632 | ((((nlines-1) - cur_line) * digits_per_line) | ||

3633 | + | ||

3634 | (digits_per_line - 1 - cur_digit_on_line)); | ||

3635 | |||

3636 | /* Print the appropriate digit or a space. The NAN case and the | ||

3637 | ** zero case need to be treated specially. | ||

3638 | */ | ||

3639 | if (i->flags) | ||

3640 | { | ||

3641 | /* Not a number. Everything is blank, except spell out | ||

3642 | ** description of condition at the end of the string of | ||

3643 | ** digits. | ||

3644 | */ | ||

3645 | int index_from_right; | ||

3646 | int virtual_index; | ||

3647 | |||

3648 | index_from_right = digits_per_line - 1 - cur_digit_on_line; | ||

3649 | //The index calculated above is calculated so that the | ||

3650 | //final position on the line has index [0]. | ||

3651 | assert(index_from_right >= 0 && index_from_right < digits_per_line); | ||

3652 | |||

3653 | //Now, calculate the "virtual index". The virtual index | ||

3654 | //is the actual number of characters from the right, taking | ||

3655 | //into account commas. | ||

3656 | virtual_index = index_from_right + index_from_right/3; | ||

3657 | |||

3658 | if (((index_from_right % 3) == 2) && cur_digit_on_line) | ||

3659 | { | ||

3660 | //We are one position past a comma. This means | ||

3661 | //that we might need a "fill" character to go | ||

3662 | //where the comma should have gone. | ||

3663 | |||

3664 | if (virtual_index + 1 < num_digits) | ||

3665 | { | ||

3666 | //The character we should print exists. | ||

3667 | fprintf(s, "%c", digits[num_digits - 2 - virtual_index]); | ||

3668 | } | ||

3669 | else | ||

3670 | { | ||

3671 | //The character doesn't exist, because the error | ||

3672 | //string is apparently too short. Must print a | ||

3673 | //space, instead. | ||

3674 | fprintf(s, " "); | ||

3675 | } | ||

3676 | } | ||

3677 | |||

3678 | //We've done the fill character, if the position we're in | ||

3679 | //is one past a comma. Now, do the ordinary character | ||

3680 | //corresponding to a digit position. | ||

3681 | if (virtual_index < num_digits) | ||

3682 | { | ||

3683 | //The character we should print exists. | ||

3684 | fprintf(s, "%c", digits[num_digits - 1 - virtual_index]); | ||

3685 | } | ||

3686 | else | ||

3687 | { | ||

3688 | //The character doesn't exist, because the error | ||

3689 | //string is apparently too short. Must print a | ||

3690 | //space, instead. | ||

3691 | fprintf(s, " "); | ||

3692 | } | ||

3693 | } | ||

3694 | else if (GMP_INTS_mpz_is_zero(i)) | ||

3695 | { | ||

3696 | /* This is the zero case. For zero, there is only one line, | ||

3697 | ** and every character except the last one is a blank. | ||

3698 | */ | ||

3699 | if (cur_digit_on_line == (digits_per_line - 1)) | ||

3700 | { | ||

3701 | fprintf(s, "0"); | ||

3702 | } | ||

3703 | else | ||

3704 | { | ||

3705 | fprintf(s, " "); | ||

3706 | } | ||

3707 | } | ||

3708 | else | ||

3709 | { | ||

3710 | /* This is a valid number which is not zero. Need to print | ||

3711 | ** the digits. | ||

3712 | */ | ||

3713 | |||

3714 | if (idx_into_string < num_digits) | ||

3715 | { | ||

3716 | int actual_index; | ||

3717 | |||

3718 | actual_index = num_digits - 1 - idx_into_string; | ||

3719 | //This is a string reversal mapping. The original | ||

3720 | //code stored strings least significant digit first, | ||

3721 | //but this code uses most significant digit first. | ||

3722 | assert((actual_index >= 0) && (actual_index < num_digits)); | ||

3723 | fprintf(s, "%c", digits[actual_index]); | ||

3724 | } | ||

3725 | else | ||

3726 | { | ||

3727 | fprintf(s, " "); | ||

3728 | } | ||

3729 | } /* End of digit case. | ||

3730 | |||

3731 | /* Now handle the commas. The rules for commas are straightforward. | ||

3732 | ** a)NAN never has a comma. | ||

3733 | ** b)Zeros never have a comma. | ||

3734 | ** c)The final line, last digit never has a comma. | ||

3735 | ** d)Everything else in multiples of three ... | ||

3736 | */ | ||

3737 | if (!(idx_into_string % 3) && (idx_into_string)) | ||

3738 | { | ||

3739 | if (i->flags) | ||

3740 | { | ||

3741 | //fprintf(s, " "); | ||

3742 | } | ||

3743 | else if (!num_digits) | ||

3744 | { | ||

3745 | fprintf(s, " "); | ||

3746 | } | ||

3747 | else | ||

3748 | { | ||

3749 | if (idx_into_string < num_digits) | ||

3750 | { | ||

3751 | fprintf(s, ","); | ||

3752 | } | ||

3753 | else | ||

3754 | { | ||

3755 | fprintf(s, " "); | ||

3756 | } | ||

3757 | } | ||

3758 | } | ||

3759 | } /* End for each digit on the current line. */ | ||

3760 | |||

3761 | /* For the first line, print out an informative message | ||

3762 | ** advising of the number of digits. For all other lines | ||

3763 | ** print nothing. | ||

3764 | */ | ||

3765 | if (!cur_line && !(i->flags)) | ||

3766 | { | ||

3767 | if (nlines == 1) | ||

3768 | fprintf(s, " "); | ||

3769 | |||

3770 | if (num_digits <= 1) | ||

3771 | { | ||

3772 | fprintf(s, " ( 1 digit )\n"); | ||

3773 | } | ||

3774 | else if (num_digits < 1000) | ||

3775 | { | ||

3776 | fprintf(s, " (%7d digits)\n", num_digits); | ||

3777 | } | ||

3778 | else | ||

3779 | { | ||

3780 | fprintf(s, " (%3d,%03d digits)\n", num_digits / 1000, num_digits % 1000); | ||

3781 | } | ||

3782 | } | ||

3783 | else | ||

3784 | { | ||

3785 | fprintf(s, "\n"); | ||

3786 | } | ||

3787 | } /* End for each line. */ | ||

3788 | |||

3789 | //Deallocate the string space. | ||

3790 | GMP_INTS_free(digits); | ||

3791 | } | ||

3792 | |||

3793 | |||

3794 | void GMP_INTS_mpz_arb_int_raw_to_stream(FILE *s, | ||

3795 | const GMP_INTS_mpz_struct *i) | ||

3796 | { | ||

3797 | int size_reqd; | ||

3798 | char *digits; | ||

3799 | |||

3800 | //Eyeball the input parameters. | ||

3801 | assert(s != NULL); | ||

3802 | assert(i != NULL); | ||

3803 | assert(i->n_allocd > 0); | ||

3804 | assert(i->limbs != NULL); | ||

3805 | |||

3806 | size_reqd = GMP_INTS_mpz_size_in_base_10(i); | ||

3807 | digits = GMP_INTS_malloc(size_reqd * sizeof(char)); | ||

3808 | GMP_INTS_mpz_to_string(digits, i); | ||

3809 | fprintf(s, "%s", digits); | ||

3810 | GMP_INTS_free(digits); | ||

3811 | } | ||

3812 | |||

3813 | |||

3814 | //07/24/01: Passed visual inspection and unit tests. | ||

3815 | void GMP_INTS_mpz_pow_ui( GMP_INTS_mpz_struct *result, | ||

3816 | const GMP_INTS_mpz_struct *base, | ||

3817 | unsigned exponent) | ||

3818 | { | ||

3819 | GMP_INTS_mpz_struct temp; | ||

3820 | //Temporary location to hold the base raised to | ||

3821 | //a binary power (repeated squaring). | ||

3822 | |||

3823 | //Eyeball the input parameters. | ||

3824 | assert(result != NULL); | ||

3825 | assert(result->n_allocd > 0); | ||

3826 | assert(result->limbs != NULL); | ||

3827 | assert(base != NULL); | ||

3828 | assert(base->n_allocd > 0); | ||

3829 | assert(base->limbs != NULL); | ||

3830 | |||

3831 | //For this function, the base and the result may not | ||

3832 | //be the same object. | ||

3833 | assert(result != base); | ||

3834 | |||

3835 | //If the base is tained, the output is tainted by association. | ||

3836 | { | ||

3837 | int taint; | ||

3838 | |||

3839 | taint = GMP_INTS_two_op_flags_map(base->flags, 0); | ||

3840 | |||

3841 | if (taint) | ||

3842 | { | ||

3843 | result->flags = taint; | ||

3844 | return; | ||

3845 | } | ||

3846 | } | ||

3847 | |||

3848 | //Allocate our temporary variable and set it to the base. | ||

3849 | GMP_INTS_mpz_init(&temp); | ||

3850 | GMP_INTS_mpz_copy(&temp, base); | ||

3851 | |||

3852 | //The result begins with the value of 1. | ||

3853 | GMP_INTS_mpz_set_ui(result, 1); | ||

3854 | |||

3855 | //Loop through, processing each bit of the exponent. This is a fairly effective | ||

3856 | //algorithm, but not the optimal one (Knuth points this out). | ||

3857 | while (exponent && !result->flags) | ||

3858 | { | ||

3859 | if (exponent & 0x1) | ||

3860 | { | ||

3861 | GMP_INTS_mpz_mul(result, result, &temp); | ||

3862 | } | ||

3863 | |||

3864 | //Square the temporary variable. Because squaring of arb integer | ||

3865 | //may be very expensive, the test against 1 (i.e. last iteration) | ||

3866 | //certainly pays for itself. | ||

3867 | if (exponent != 1) | ||

3868 | GMP_INTS_mpz_mul(&temp, &temp, &temp); | ||

3869 | |||

3870 | exponent >>= 1; | ||

3871 | } | ||

3872 | |||

3873 | //Deallocate our temporary variable. | ||

3874 | GMP_INTS_mpz_clear(&temp); | ||

3875 | } | ||

3876 | |||

3877 | |||

3878 | void GMP_INTS_mpz_abs(GMP_INTS_mpz_struct *arg) | ||

3879 | { | ||

3880 | //Eyeball the input parameter. | ||

3881 | assert(arg != NULL); | ||

3882 | assert(arg->n_allocd > 0); | ||

3883 | assert(arg->limbs != NULL); | ||

3884 | |||

3885 | //Take the absolute value. | ||

3886 | if (arg->size < 0) | ||

3887 | arg->size = -arg->size; | ||

3888 | } | ||

3889 | |||

3890 | |||

3891 | //07/29/01: Visual inspection passed. Seems to work fine--not explicitly unit-tested | ||

3892 | //directly, but was tested from Tcl. | ||

3893 | void GMP_INTS_mpz_gcd(GMP_INTS_mpz_struct *result, | ||

3894 | const GMP_INTS_mpz_struct *arg1, | ||

3895 | const GMP_INTS_mpz_struct *arg2) | ||

3896 | { | ||

3897 | GMP_INTS_mpz_struct u, v, q, r; | ||

3898 | int loop_count; | ||

3899 | |||

3900 | //Eyeball the inputs carefully. | ||

3901 | assert(result != NULL); | ||

3902 | assert(result->n_allocd > 0); | ||

3903 | assert(result->limbs != NULL); | ||

3904 | assert(arg1 != NULL); | ||

3905 | assert(arg1->n_allocd > 0); | ||

3906 | assert(arg1->limbs != NULL); | ||

3907 | assert(arg2 != NULL); | ||

3908 | assert(arg2->n_allocd > 0); | ||

3909 | assert(arg2->limbs != NULL); | ||

3910 | |||

3911 | //Args are not allowed to be same object. | ||

3912 | assert(arg1 != arg2); | ||

3913 | |||

3914 | //If either input is error or taint, taint the output. | ||

3915 | { | ||

3916 | int taint; | ||

3917 | |||

3918 | taint = GMP_INTS_two_op_flags_map(arg1->flags, arg2->flags); | ||

3919 | |||

3920 | result->flags = 0; | ||

3921 | //"result" starts off with a clean slate. Must do this | ||

3922 | //after taint calculation in case locations of arg1 or arg2 | ||

3923 | //are the same as result. | ||

3924 | if (taint) | ||

3925 | { | ||

3926 | result->flags = taint; | ||

3927 | return; | ||

3928 | } | ||

3929 | } | ||

3930 | |||

3931 | //If either input is zero, the result is 1. | ||

3932 | if (GMP_INTS_mpz_is_zero(arg1) || GMP_INTS_mpz_is_zero(arg2)) | ||

3933 | { | ||

3934 | GMP_INTS_mpz_set_ui(result, 1); | ||

3935 | return; | ||

3936 | } | ||

3937 | |||

3938 | //Allocate space for locals. | ||

3939 | GMP_INTS_mpz_init(&u); | ||

3940 | GMP_INTS_mpz_init(&v); | ||

3941 | GMP_INTS_mpz_init(&q); | ||

3942 | GMP_INTS_mpz_init(&r); | ||

3943 | |||

3944 | //We are following Knuth Vol 2, p. 337, the modern Euclidian algorithm. | ||

3945 | //Note: There are faster algorithms for GCD, but because I hacked up the | ||

3946 | //GMP multiple-precision library so badly, those aren't included. This one | ||

3947 | //is logically correct but sub-optimal. Perhaps at a later time faster | ||

3948 | //algorithms will be re-included. |