/[dtapublic]/sf_code/esrgpubs/acm0010/c_imp/rap/rap_c.c
ViewVC logotype

Contents of /sf_code/esrgpubs/acm0010/c_imp/rap/rap_c.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 24 - (show annotations) (download)
Sat Oct 8 06:15:00 2016 UTC (7 years, 4 months ago) by dashley
File MIME type: text/plain
File size: 583746 byte(s)
Initial commit.
1 /* $Archive:: /ACM Rational Approximation Paper And Algorithms/C-Language Implementation/rap_c.c $ */
2 /* $Revision: 1.1 $ */
3 /* $Date: 2001/09/25 21:44:55 $ */
4 /* $Author: dtashley $ */
5 /* $Workfile:: rap_c.c $ */
6 /*
7 /*
8 //--------------------------------------------------------------------------------
9 //Copyright 2008 David T. Ashley
10 //-------------------------------------------------------------------------------------------------
11 //This source code and any program in which it is compiled/used is provided under the GNU GENERAL
12 //PUBLIC LICENSE, Version 3, full license text below.
13 //-------------------------------------------------------------------------------------------------
14 // GNU GENERAL PUBLIC LICENSE
15 // Version 3, 29 June 2007
16 //
17 // Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
18 // Everyone is permitted to copy and distribute verbatim copies
19 // of this license document, but changing it is not allowed.
20 //
21 // Preamble
22 //
23 // The GNU General Public License is a free, copyleft license for
24 //software and other kinds of works.
25 //
26 // The licenses for most software and other practical works are designed
27 //to take away your freedom to share and change the works. By contrast,
28 //the GNU General Public License is intended to guarantee your freedom to
29 //share and change all versions of a program--to make sure it remains free
30 //software for all its users. We, the Free Software Foundation, use the
31 //GNU General Public License for most of our software; it applies also to
32 //any other work released this way by its authors. You can apply it to
33 //your programs, too.
34 //
35 // When we speak of free software, we are referring to freedom, not
36 //price. Our General Public Licenses are designed to make sure that you
37 //have the freedom to distribute copies of free software (and charge for
38 //them if you wish), that you receive source code or can get it if you
39 //want it, that you can change the software or use pieces of it in new
40 //free programs, and that you know you can do these things.
41 //
42 // To protect your rights, we need to prevent others from denying you
43 //these rights or asking you to surrender the rights. Therefore, you have
44 //certain responsibilities if you distribute copies of the software, or if
45 //you modify it: responsibilities to respect the freedom of others.
46 //
47 // For example, if you distribute copies of such a program, whether
48 //gratis or for a fee, you must pass on to the recipients the same
49 //freedoms that you received. You must make sure that they, too, receive
50 //or can get the source code. And you must show them these terms so they
51 //know their rights.
52 //
53 // Developers that use the GNU GPL protect your rights with two steps:
54 //(1) assert copyright on the software, and (2) offer you this License
55 //giving you legal permission to copy, distribute and/or modify it.
56 //
57 // For the developers' and authors' protection, the GPL clearly explains
58 //that there is no warranty for this free software. For both users' and
59 //authors' sake, the GPL requires that modified versions be marked as
60 //changed, so that their problems will not be attributed erroneously to
61 //authors of previous versions.
62 //
63 // Some devices are designed to deny users access to install or run
64 //modified versions of the software inside them, although the manufacturer
65 //can do so. This is fundamentally incompatible with the aim of
66 //protecting users' freedom to change the software. The systematic
67 //pattern of such abuse occurs in the area of products for individuals to
68 //use, which is precisely where it is most unacceptable. Therefore, we
69 //have designed this version of the GPL to prohibit the practice for those
70 //products. If such problems arise substantially in other domains, we
71 //stand ready to extend this provision to those domains in future versions
72 //of the GPL, as needed to protect the freedom of users.
73 //
74 // Finally, every program is threatened constantly by software patents.
75 //States should not allow patents to restrict development and use of
76 //software on general-purpose computers, but in those that do, we wish to
77 //avoid the special danger that patents applied to a free program could
78 //make it effectively proprietary. To prevent this, the GPL assures that
79 //patents cannot be used to render the program non-free.
80 //
81 // The precise terms and conditions for copying, distribution and
82 //modification follow.
83 //
84 // TERMS AND CONDITIONS
85 //
86 // 0. Definitions.
87 //
88 // "This License" refers to version 3 of the GNU General Public License.
89 //
90 // "Copyright" also means copyright-like laws that apply to other kinds of
91 //works, such as semiconductor masks.
92 //
93 // "The Program" refers to any copyrightable work licensed under this
94 //License. Each licensee is addressed as "you". "Licensees" and
95 //"recipients" may be individuals or organizations.
96 //
97 // To "modify" a work means to copy from or adapt all or part of the work
98 //in a fashion requiring copyright permission, other than the making of an
99 //exact copy. The resulting work is called a "modified version" of the
100 //earlier work or a work "based on" the earlier work.
101 //
102 // A "covered work" means either the unmodified Program or a work based
103 //on the Program.
104 //
105 // To "propagate" a work means to do anything with it that, without
106 //permission, would make you directly or secondarily liable for
107 //infringement under applicable copyright law, except executing it on a
108 //computer or modifying a private copy. Propagation includes copying,
109 //distribution (with or without modification), making available to the
110 //public, and in some countries other activities as well.
111 //
112 // To "convey" a work means any kind of propagation that enables other
113 //parties to make or receive copies. Mere interaction with a user through
114 //a computer network, with no transfer of a copy, is not conveying.
115 //
116 // An interactive user interface displays "Appropriate Legal Notices"
117 //to the extent that it includes a convenient and prominently visible
118 //feature that (1) displays an appropriate copyright notice, and (2)
119 //tells the user that there is no warranty for the work (except to the
120 //extent that warranties are provided), that licensees may convey the
121 //work under this License, and how to view a copy of this License. If
122 //the interface presents a list of user commands or options, such as a
123 //menu, a prominent item in the list meets this criterion.
124 //
125 // 1. Source Code.
126 //
127 // The "source code" for a work means the preferred form of the work
128 //for making modifications to it. "Object code" means any non-source
129 //form of a work.
130 //
131 // A "Standard Interface" means an interface that either is an official
132 //standard defined by a recognized standards body, or, in the case of
133 //interfaces specified for a particular programming language, one that
134 //is widely used among developers working in that language.
135 //
136 // The "System Libraries" of an executable work include anything, other
137 //than the work as a whole, that (a) is included in the normal form of
138 //packaging a Major Component, but which is not part of that Major
139 //Component, and (b) serves only to enable use of the work with that
140 //Major Component, or to implement a Standard Interface for which an
141 //implementation is available to the public in source code form. A
142 //"Major Component", in this context, means a major essential component
143 //(kernel, window system, and so on) of the specific operating system
144 //(if any) on which the executable work runs, or a compiler used to
145 //produce the work, or an object code interpreter used to run it.
146 //
147 // The "Corresponding Source" for a work in object code form means all
148 //the source code needed to generate, install, and (for an executable
149 //work) run the object code and to modify the work, including scripts to
150 //control those activities. However, it does not include the work's
151 //System Libraries, or general-purpose tools or generally available free
152 //programs which are used unmodified in performing those activities but
153 //which are not part of the work. For example, Corresponding Source
154 //includes interface definition files associated with source files for
155 //the work, and the source code for shared libraries and dynamically
156 //linked subprograms that the work is specifically designed to require,
157 //such as by intimate data communication or control flow between those
158 //subprograms and other parts of the work.
159 //
160 // The Corresponding Source need not include anything that users
161 //can regenerate automatically from other parts of the Corresponding
162 //Source.
163 //
164 // The Corresponding Source for a work in source code form is that
165 //same work.
166 //
167 // 2. Basic Permissions.
168 //
169 // All rights granted under this License are granted for the term of
170 //copyright on the Program, and are irrevocable provided the stated
171 //conditions are met. This License explicitly affirms your unlimited
172 //permission to run the unmodified Program. The output from running a
173 //covered work is covered by this License only if the output, given its
174 //content, constitutes a covered work. This License acknowledges your
175 //rights of fair use or other equivalent, as provided by copyright law.
176 //
177 // You may make, run and propagate covered works that you do not
178 //convey, without conditions so long as your license otherwise remains
179 //in force. You may convey covered works to others for the sole purpose
180 //of having them make modifications exclusively for you, or provide you
181 //with facilities for running those works, provided that you comply with
182 //the terms of this License in conveying all material for which you do
183 //not control copyright. Those thus making or running the covered works
184 //for you must do so exclusively on your behalf, under your direction
185 //and control, on terms that prohibit them from making any copies of
186 //your copyrighted material outside their relationship with you.
187 //
188 // Conveying under any other circumstances is permitted solely under
189 //the conditions stated below. Sublicensing is not allowed; section 10
190 //makes it unnecessary.
191 //
192 // 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
193 //
194 // No covered work shall be deemed part of an effective technological
195 //measure under any applicable law fulfilling obligations under article
196 //11 of the WIPO copyright treaty adopted on 20 December 1996, or
197 //similar laws prohibiting or restricting circumvention of such
198 //measures.
199 //
200 // When you convey a covered work, you waive any legal power to forbid
201 //circumvention of technological measures to the extent such circumvention
202 //is effected by exercising rights under this License with respect to
203 //the covered work, and you disclaim any intention to limit operation or
204 //modification of the work as a means of enforcing, against the work's
205 //users, your or third parties' legal rights to forbid circumvention of
206 //technological measures.
207 //
208 // 4. Conveying Verbatim Copies.
209 //
210 // You may convey verbatim copies of the Program's source code as you
211 //receive it, in any medium, provided that you conspicuously and
212 //appropriately publish on each copy an appropriate copyright notice;
213 //keep intact all notices stating that this License and any
214 //non-permissive terms added in accord with section 7 apply to the code;
215 //keep intact all notices of the absence of any warranty; and give all
216 //recipients a copy of this License along with the Program.
217 //
218 // You may charge any price or no price for each copy that you convey,
219 //and you may offer support or warranty protection for a fee.
220 //
221 // 5. Conveying Modified Source Versions.
222 //
223 // You may convey a work based on the Program, or the modifications to
224 //produce it from the Program, in the form of source code under the
225 //terms of section 4, provided that you also meet all of these conditions:
226 //
227 // a) The work must carry prominent notices stating that you modified
228 // it, and giving a relevant date.
229 //
230 // b) The work must carry prominent notices stating that it is
231 // released under this License and any conditions added under section
232 // 7. This requirement modifies the requirement in section 4 to
233 // "keep intact all notices".
234 //
235 // c) You must license the entire work, as a whole, under this
236 // License to anyone who comes into possession of a copy. This
237 // License will therefore apply, along with any applicable section 7
238 // additional terms, to the whole of the work, and all its parts,
239 // regardless of how they are packaged. This License gives no
240 // permission to license the work in any other way, but it does not
241 // invalidate such permission if you have separately received it.
242 //
243 // d) If the work has interactive user interfaces, each must display
244 // Appropriate Legal Notices; however, if the Program has interactive
245 // interfaces that do not display Appropriate Legal Notices, your
246 // work need not make them do so.
247 //
248 // A compilation of a covered work with other separate and independent
249 //works, which are not by their nature extensions of the covered work,
250 //and which are not combined with it such as to form a larger program,
251 //in or on a volume of a storage or distribution medium, is called an
252 //"aggregate" if the compilation and its resulting copyright are not
253 //used to limit the access or legal rights of the compilation's users
254 //beyond what the individual works permit. Inclusion of a covered work
255 //in an aggregate does not cause this License to apply to the other
256 //parts of the aggregate.
257 //
258 // 6. Conveying Non-Source Forms.
259 //
260 // You may convey a covered work in object code form under the terms
261 //of sections 4 and 5, provided that you also convey the
262 //machine-readable Corresponding Source under the terms of this License,
263 //in one of these ways:
264 //
265 // a) Convey the object code in, or embodied in, a physical product
266 // (including a physical distribution medium), accompanied by the
267 // Corresponding Source fixed on a durable physical medium
268 // customarily used for software interchange.
269 //
270 // b) Convey the object code in, or embodied in, a physical product
271 // (including a physical distribution medium), accompanied by a
272 // written offer, valid for at least three years and valid for as
273 // long as you offer spare parts or customer support for that product
274 // model, to give anyone who possesses the object code either (1) a
275 // copy of the Corresponding Source for all the software in the
276 // product that is covered by this License, on a durable physical
277 // medium customarily used for software interchange, for a price no
278 // more than your reasonable cost of physically performing this
279 // conveying of source, or (2) access to copy the
280 // Corresponding Source from a network server at no charge.
281 //
282 // c) Convey individual copies of the object code with a copy of the
283 // written offer to provide the Corresponding Source. This
284 // alternative is allowed only occasionally and noncommercially, and
285 // only if you received the object code with such an offer, in accord
286 // with subsection 6b.
287 //
288 // d) Convey the object code by offering access from a designated
289 // place (gratis or for a charge), and offer equivalent access to the
290 // Corresponding Source in the same way through the same place at no
291 // further charge. You need not require recipients to copy the
292 // Corresponding Source along with the object code. If the place to
293 // copy the object code is a network server, the Corresponding Source
294 // may be on a different server (operated by you or a third party)
295 // that supports equivalent copying facilities, provided you maintain
296 // clear directions next to the object code saying where to find the
297 // Corresponding Source. Regardless of what server hosts the
298 // Corresponding Source, you remain obligated to ensure that it is
299 // available for as long as needed to satisfy these requirements.
300 //
301 // e) Convey the object code using peer-to-peer transmission, provided
302 // you inform other peers where the object code and Corresponding
303 // Source of the work are being offered to the general public at no
304 // charge under subsection 6d.
305 //
306 // A separable portion of the object code, whose source code is excluded
307 //from the Corresponding Source as a System Library, need not be
308 //included in conveying the object code work.
309 //
310 // A "User Product" is either (1) a "consumer product", which means any
311 //tangible personal property which is normally used for personal, family,
312 //or household purposes, or (2) anything designed or sold for incorporation
313 //into a dwelling. In determining whether a product is a consumer product,
314 //doubtful cases shall be resolved in favor of coverage. For a particular
315 //product received by a particular user, "normally used" refers to a
316 //typical or common use of that class of product, regardless of the status
317 //of the particular user or of the way in which the particular user
318 //actually uses, or expects or is expected to use, the product. A product
319 //is a consumer product regardless of whether the product has substantial
320 //commercial, industrial or non-consumer uses, unless such uses represent
321 //the only significant mode of use of the product.
322 //
323 // "Installation Information" for a User Product means any methods,
324 //procedures, authorization keys, or other information required to install
325 //and execute modified versions of a covered work in that User Product from
326 //a modified version of its Corresponding Source. The information must
327 //suffice to ensure that the continued functioning of the modified object
328 //code is in no case prevented or interfered with solely because
329 //modification has been made.
330 //
331 // If you convey an object code work under this section in, or with, or
332 //specifically for use in, a User Product, and the conveying occurs as
333 //part of a transaction in which the right of possession and use of the
334 //User Product is transferred to the recipient in perpetuity or for a
335 //fixed term (regardless of how the transaction is characterized), the
336 //Corresponding Source conveyed under this section must be accompanied
337 //by the Installation Information. But this requirement does not apply
338 //if neither you nor any third party retains the ability to install
339 //modified object code on the User Product (for example, the work has
340 //been installed in ROM).
341 //
342 // The requirement to provide Installation Information does not include a
343 //requirement to continue to provide support service, warranty, or updates
344 //for a work that has been modified or installed by the recipient, or for
345 //the User Product in which it has been modified or installed. Access to a
346 //network may be denied when the modification itself materially and
347 //adversely affects the operation of the network or violates the rules and
348 //protocols for communication across the network.
349 //
350 // Corresponding Source conveyed, and Installation Information provided,
351 //in accord with this section must be in a format that is publicly
352 //documented (and with an implementation available to the public in
353 //source code form), and must require no special password or key for
354 //unpacking, reading or copying.
355 //
356 // 7. Additional Terms.
357 //
358 // "Additional permissions" are terms that supplement the terms of this
359 //License by making exceptions from one or more of its conditions.
360 //Additional permissions that are applicable to the entire Program shall
361 //be treated as though they were included in this License, to the extent
362 //that they are valid under applicable law. If additional permissions
363 //apply only to part of the Program, that part may be used separately
364 //under those permissions, but the entire Program remains governed by
365 //this License without regard to the additional permissions.
366 //
367 // When you convey a copy of a covered work, you may at your option
368 //remove any additional permissions from that copy, or from any part of
369 //it. (Additional permissions may be written to require their own
370 //removal in certain cases when you modify the work.) You may place
371 //additional permissions on material, added by you to a covered work,
372 //for which you have or can give appropriate copyright permission.
373 //
374 // Notwithstanding any other provision of this License, for material you
375 //add to a covered work, you may (if authorized by the copyright holders of
376 //that material) supplement the terms of this License with terms:
377 //
378 // a) Disclaiming warranty or limiting liability differently from the
379 // terms of sections 15 and 16 of this License; or
380 //
381 // b) Requiring preservation of specified reasonable legal notices or
382 // author attributions in that material or in the Appropriate Legal
383 // Notices displayed by works containing it; or
384 //
385 // c) Prohibiting misrepresentation of the origin of that material, or
386 // requiring that modified versions of such material be marked in
387 // reasonable ways as different from the original version; or
388 //
389 // d) Limiting the use for publicity purposes of names of licensors or
390 // authors of the material; or
391 //
392 // e) Declining to grant rights under trademark law for use of some
393 // trade names, trademarks, or service marks; or
394 //
395 // f) Requiring indemnification of licensors and authors of that
396 // material by anyone who conveys the material (or modified versions of
397 // it) with contractual assumptions of liability to the recipient, for
398 // any liability that these contractual assumptions directly impose on
399 // those licensors and authors.
400 //
401 // All other non-permissive additional terms are considered "further
402 //restrictions" within the meaning of section 10. If the Program as you
403 //received it, or any part of it, contains a notice stating that it is
404 //governed by this License along with a term that is a further
405 //restriction, you may remove that term. If a license document contains
406 //a further restriction but permits relicensing or conveying under this
407 //License, you may add to a covered work material governed by the terms
408 //of that license document, provided that the further restriction does
409 //not survive such relicensing or conveying.
410 //
411 // If you add terms to a covered work in accord with this section, you
412 //must place, in the relevant source files, a statement of the
413 //additional terms that apply to those files, or a notice indicating
414 //where to find the applicable terms.
415 //
416 // Additional terms, permissive or non-permissive, may be stated in the
417 //form of a separately written license, or stated as exceptions;
418 //the above requirements apply either way.
419 //
420 // 8. Termination.
421 //
422 // You may not propagate or modify a covered work except as expressly
423 //provided under this License. Any attempt otherwise to propagate or
424 //modify it is void, and will automatically terminate your rights under
425 //this License (including any patent licenses granted under the third
426 //paragraph of section 11).
427 //
428 // However, if you cease all violation of this License, then your
429 //license from a particular copyright holder is reinstated (a)
430 //provisionally, unless and until the copyright holder explicitly and
431 //finally terminates your license, and (b) permanently, if the copyright
432 //holder fails to notify you of the violation by some reasonable means
433 //prior to 60 days after the cessation.
434 //
435 // Moreover, your license from a particular copyright holder is
436 //reinstated permanently if the copyright holder notifies you of the
437 //violation by some reasonable means, this is the first time you have
438 //received notice of violation of this License (for any work) from that
439 //copyright holder, and you cure the violation prior to 30 days after
440 //your receipt of the notice.
441 //
442 // Termination of your rights under this section does not terminate the
443 //licenses of parties who have received copies or rights from you under
444 //this License. If your rights have been terminated and not permanently
445 //reinstated, you do not qualify to receive new licenses for the same
446 //material under section 10.
447 //
448 // 9. Acceptance Not Required for Having Copies.
449 //
450 // You are not required to accept this License in order to receive or
451 //run a copy of the Program. Ancillary propagation of a covered work
452 //occurring solely as a consequence of using peer-to-peer transmission
453 //to receive a copy likewise does not require acceptance. However,
454 //nothing other than this License grants you permission to propagate or
455 //modify any covered work. These actions infringe copyright if you do
456 //not accept this License. Therefore, by modifying or propagating a
457 //covered work, you indicate your acceptance of this License to do so.
458 //
459 // 10. Automatic Licensing of Downstream Recipients.
460 //
461 // Each time you convey a covered work, the recipient automatically
462 //receives a license from the original licensors, to run, modify and
463 //propagate that work, subject to this License. You are not responsible
464 //for enforcing compliance by third parties with this License.
465 //
466 // An "entity transaction" is a transaction transferring control of an
467 //organization, or substantially all assets of one, or subdividing an
468 //organization, or merging organizations. If propagation of a covered
469 //work results from an entity transaction, each party to that
470 //transaction who receives a copy of the work also receives whatever
471 //licenses to the work the party's predecessor in interest had or could
472 //give under the previous paragraph, plus a right to possession of the
473 //Corresponding Source of the work from the predecessor in interest, if
474 //the predecessor has it or can get it with reasonable efforts.
475 //
476 // You may not impose any further restrictions on the exercise of the
477 //rights granted or affirmed under this License. For example, you may
478 //not impose a license fee, royalty, or other charge for exercise of
479 //rights granted under this License, and you may not initiate litigation
480 //(including a cross-claim or counterclaim in a lawsuit) alleging that
481 //any patent claim is infringed by making, using, selling, offering for
482 //sale, or importing the Program or any portion of it.
483 //
484 // 11. Patents.
485 //
486 // A "contributor" is a copyright holder who authorizes use under this
487 //License of the Program or a work on which the Program is based. The
488 //work thus licensed is called the contributor's "contributor version".
489 //
490 // A contributor's "essential patent claims" are all patent claims
491 //owned or controlled by the contributor, whether already acquired or
492 //hereafter acquired, that would be infringed by some manner, permitted
493 //by this License, of making, using, or selling its contributor version,
494 //but do not include claims that would be infringed only as a
495 //consequence of further modification of the contributor version. For
496 //purposes of this definition, "control" includes the right to grant
497 //patent sublicenses in a manner consistent with the requirements of
498 //this License.
499 //
500 // Each contributor grants you a non-exclusive, worldwide, royalty-free
501 //patent license under the contributor's essential patent claims, to
502 //make, use, sell, offer for sale, import and otherwise run, modify and
503 //propagate the contents of its contributor version.
504 //
505 // In the following three paragraphs, a "patent license" is any express
506 //agreement or commitment, however denominated, not to enforce a patent
507 //(such as an express permission to practice a patent or covenant not to
508 //sue for patent infringement). To "grant" such a patent license to a
509 //party means to make such an agreement or commitment not to enforce a
510 //patent against the party.
511 //
512 // If you convey a covered work, knowingly relying on a patent license,
513 //and the Corresponding Source of the work is not available for anyone
514 //to copy, free of charge and under the terms of this License, through a
515 //publicly available network server or other readily accessible means,
516 //then you must either (1) cause the Corresponding Source to be so
517 //available, or (2) arrange to deprive yourself of the benefit of the
518 //patent license for this particular work, or (3) arrange, in a manner
519 //consistent with the requirements of this License, to extend the patent
520 //license to downstream recipients. "Knowingly relying" means you have
521 //actual knowledge that, but for the patent license, your conveying the
522 //covered work in a country, or your recipient's use of the covered work
523 //in a country, would infringe one or more identifiable patents in that
524 //country that you have reason to believe are valid.
525 //
526 // If, pursuant to or in connection with a single transaction or
527 //arrangement, you convey, or propagate by procuring conveyance of, a
528 //covered work, and grant a patent license to some of the parties
529 //receiving the covered work authorizing them to use, propagate, modify
530 //or convey a specific copy of the covered work, then the patent license
531 //you grant is automatically extended to all recipients of the covered
532 //work and works based on it.
533 //
534 // A patent license is "discriminatory" if it does not include within
535 //the scope of its coverage, prohibits the exercise of, or is
536 //conditioned on the non-exercise of one or more of the rights that are
537 //specifically granted under this License. You may not convey a covered
538 //work if you are a party to an arrangement with a third party that is
539 //in the business of distributing software, under which you make payment
540 //to the third party based on the extent of your activity of conveying
541 //the work, and under which the third party grants, to any of the
542 //parties who would receive the covered work from you, a discriminatory
543 //patent license (a) in connection with copies of the covered work
544 //conveyed by you (or copies made from those copies), or (b) primarily
545 //for and in connection with specific products or compilations that
546 //contain the covered work, unless you entered into that arrangement,
547 //or that patent license was granted, prior to 28 March 2007.
548 //
549 // Nothing in this License shall be construed as excluding or limiting
550 //any implied license or other defenses to infringement that may
551 //otherwise be available to you under applicable patent law.
552 //
553 // 12. No Surrender of Others' Freedom.
554 //
555 // If conditions are imposed on you (whether by court order, agreement or
556 //otherwise) that contradict the conditions of this License, they do not
557 //excuse you from the conditions of this License. If you cannot convey a
558 //covered work so as to satisfy simultaneously your obligations under this
559 //License and any other pertinent obligations, then as a consequence you may
560 //not convey it at all. For example, if you agree to terms that obligate you
561 //to collect a royalty for further conveying from those to whom you convey
562 //the Program, the only way you could satisfy both those terms and this
563 //License would be to refrain entirely from conveying the Program.
564 //
565 // 13. Use with the GNU Affero General Public License.
566 //
567 // Notwithstanding any other provision of this License, you have
568 //permission to link or combine any covered work with a work licensed
569 //under version 3 of the GNU Affero General Public License into a single
570 //combined work, and to convey the resulting work. The terms of this
571 //License will continue to apply to the part which is the covered work,
572 //but the special requirements of the GNU Affero General Public License,
573 //section 13, concerning interaction through a network will apply to the
574 //combination as such.
575 //
576 // 14. Revised Versions of this License.
577 //
578 // The Free Software Foundation may publish revised and/or new versions of
579 //the GNU General Public License from time to time. Such new versions will
580 //be similar in spirit to the present version, but may differ in detail to
581 //address new problems or concerns.
582 //
583 // Each version is given a distinguishing version number. If the
584 //Program specifies that a certain numbered version of the GNU General
585 //Public License "or any later version" applies to it, you have the
586 //option of following the terms and conditions either of that numbered
587 //version or of any later version published by the Free Software
588 //Foundation. If the Program does not specify a version number of the
589 //GNU General Public License, you may choose any version ever published
590 //by the Free Software Foundation.
591 //
592 // If the Program specifies that a proxy can decide which future
593 //versions of the GNU General Public License can be used, that proxy's
594 //public statement of acceptance of a version permanently authorizes you
595 //to choose that version for the Program.
596 //
597 // Later license versions may give you additional or different
598 //permissions. However, no additional obligations are imposed on any
599 //author or copyright holder as a result of your choosing to follow a
600 //later version.
601 //
602 // 15. Disclaimer of Warranty.
603 //
604 // THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
605 //APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
606 //HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
607 //OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
608 //THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
609 //PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
610 //IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
611 //ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
612 //
613 // 16. Limitation of Liability.
614 //
615 // IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
616 //WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
617 //THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
618 //GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
619 //USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
620 //DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
621 //PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
622 //EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
623 //SUCH DAMAGES.
624 //
625 // 17. Interpretation of Sections 15 and 16.
626 //
627 // If the disclaimer of warranty and limitation of liability provided
628 //above cannot be given local legal effect according to their terms,
629 //reviewing courts shall apply local law that most closely approximates
630 //an absolute waiver of all civil liability in connection with the
631 //Program, unless a warranty or assumption of liability accompanies a
632 //copy of the Program in return for a fee.
633 //
634 // END OF TERMS AND CONDITIONS
635 //
636 // How to Apply These Terms to Your New Programs
637 //
638 // If you develop a new program, and you want it to be of the greatest
639 //possible use to the public, the best way to achieve this is to make it
640 //free software which everyone can redistribute and change under these terms.
641 //
642 // To do so, attach the following notices to the program. It is safest
643 //to attach them to the start of each source file to most effectively
644 //state the exclusion of warranty; and each file should have at least
645 //the "copyright" line and a pointer to where the full notice is found.
646 //
647 // <one line to give the program's name and a brief idea of what it does.>
648 // Copyright (C) <year> <name of author>
649 //
650 // This program is free software: you can redistribute it and/or modify
651 // it under the terms of the GNU General Public License as published by
652 // the Free Software Foundation, either version 3 of the License, or
653 // (at your option) any later version.
654 //
655 // This program is distributed in the hope that it will be useful,
656 // but WITHOUT ANY WARRANTY; without even the implied warranty of
657 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
658 // GNU General Public License for more details.
659 //
660 // You should have received a copy of the GNU General Public License
661 // along with this program. If not, see <http://www.gnu.org/licenses/>.
662 //
663 //Also add information on how to contact you by electronic and paper mail.
664 //
665 // If the program does terminal interaction, make it output a short
666 //notice like this when it starts in an interactive mode:
667 //
668 // <program> Copyright (C) <year> <name of author>
669 // This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
670 // This is free software, and you are welcome to redistribute it
671 // under certain conditions; type `show c' for details.
672 //
673 //The hypothetical commands `show w' and `show c' should show the appropriate
674 //parts of the General Public License. Of course, your program's commands
675 //might be different; for a GUI interface, you would use an "about box".
676 //
677 // You should also get your employer (if you work as a programmer) or school,
678 //if any, to sign a "copyright disclaimer" for the program, if necessary.
679 //For more information on this, and how to apply and follow the GNU GPL, see
680 //<http://www.gnu.org/licenses/>.
681 //
682 // The GNU General Public License does not permit incorporating your program
683 //into proprietary programs. If your program is a subroutine library, you
684 //may consider it more useful to permit linking proprietary applications with
685 //the library. If this is what you want to do, use the GNU Lesser General
686 //Public License instead of this License. But first, please read
687 //<http://www.gnu.org/philosophy/why-not-lgpl.html>.
688 //-------------------------------------------------------------------------------------------------
689 //--------------------------------------------------------------------------------
690 ** Rational Approximation Software Submitted To CALGO To Accompany TOMS Paper
691 ** Entitled "On Best Rational Approximations Using Large Integers"
692 ** --------------------------------------------------------------------------
693 ** Please see accompanying user's guide for invocation instructions.
694 */
695
696 #include <process.h>
697 /* For the exit() function. */
698
699 #include <stdio.h>
700 /* printf() and similar functions require the stdio.h header file. */
701
702 #include <string.h>
703 /* String functions and memory manipulation functions defined there. */
704
705 #include <malloc.h>
706 /* Necessary for memory allocation functions. */
707
708
709 /* It wasn't straightforward to locate the file containing the constants TRUE
710 ** and FALSE. Must define them myself.
711 */
712 #ifndef FALSE
713 #define FALSE (0)
714 #endif
715 #ifndef TRUE
716 #define TRUE (1)
717 #endif
718
719
720 /* This source file is broken into segments, each
721 ** containing a specific family of functions
722 ** that operate on a specific data type, or at a
723 ** specific level of abstraction. In C++, these
724 ** would often tend to be made into classes (and
725 ** this could be done here), but it is a more
726 ** humane experience on anyone who would like to
727 ** recompile the program to keep everything in a
728 ** single source file.
729 **
730 ** Functions are arranged from the bottom of the calling tree upward. This
731 ** way, no function prototypes are needed to ensure the compiler will
732 ** detect definitions and invocations which are inconsistent.
733 */
734
735 /****************************************************************************/
736 /****************************************************************************/
737 /************* C O M P I L A T I O N C O N S T A N T S ***************/
738 /****************************************************************************/
739 /****************************************************************************/
740 /* Many aspects of the program's behavior can be changed by simply changing
741 ** the constants below. Many of the values chosen below are intended to
742 ** keep the software clear of platform-specific limits. Changes to the
743 ** constants may cause unforeseen difficulties.
744 */
745 #define MAX_CMDLINE_PARS (10)
746 /* The maximum number of command-line parameters that will be accepted.
747 ** This is ultimately dependent on what the subcommands take, but
748 ** 10 is a safe value, and let's me be lazy and not search to find the
749 ** command that takes the most parameters.
750 */
751
752 #define INPUT_INTEGER_MAX_DIGITS (463)
753 /* This value was chosen to allow operation in the Farey series of up to
754 ** order 2**1536.
755 */
756 #define INTERMEDIATE_CALC_MAX_DIGITS (4000)
757 /* This value wasn't chosen scientifically. It is just known that
758 ** multiplying a 463-digit integer by another will result in at most
759 ** a 926-digit result ... it would be quite difficult in the calculations
760 ** outlined in the TOMS paper to reach the 4,000 digit limit. Note that
761 ** this is also the limit for the size of the integer components of
762 ** rational results.
763 */
764 #define STDIN_MAX_CHARS (32000)
765 /* The maximum number of digits that may be read from the standard input
766 ** if batch mode is selected.
767 */
768 #define LINE_LEN (78)
769 /* The number of columns assumed to be displayable on the output.
770 ** because a standard monitor (legacy) was 80 columns wide, I've chosen
771 ** a value just a tad less.
772 */
773 #define NUMBER_DESC_WIDTH (20)
774 /* The number of characters at the left reserved for describing a particular
775 ** integer.
776 */
777 #define DIGITS_PER_LINE (27)
778 /* The number of digits to be displayed per line for very long
779 ** integers.
780 */
781 #define HORIZONTAL_BAR_SEP_CHAR ('-')
782 /* The character used to form large horizontal separators.
783 */
784
785 /****************************************************************************/
786 /****************************************************************************/
787 /********* V E R S I O N C O N T R O L F U N C T I O N S ***********/
788 /****************************************************************************/
789 /****************************************************************************/
790
791 /****************************************************************************/
792 /* vcGetVcData(): */
793 /*--------------------------------------------------------------------------*/
794 /* DESCRIPTION */
795 /* Returns the embedded version control information for this program, */
796 /* arranged as a sequence of lines to be displayed on an 80-column */
797 /* monitor. */
798 /* */
799 /* INPUTS */
800 /* n : Line number to return. The first line is number 0. */
801 /* */
802 /* OUTPUTS */
803 /* <--: Pointer to a constant character string to display on that line, */
804 /* or NULL if there are no more lines of version control */
805 /* information to display. */
806 /****************************************************************************/
807 const char *vcGetVcData(unsigned n)
808 {
809 const char *RAP_version_string_array[]
810 = {"$Revision: 1.1 $ $Date: 2001/09/25 21:44:55 $"};
811 /* Note that the info above is filled in automatically by Visual
812 ** Source Safe on every check-in--this is not manually maintained.
813 */
814
815 if (n < (sizeof(RAP_version_string_array)/sizeof(RAP_version_string_array[0])))
816 {
817 return(RAP_version_string_array[n]);
818 }
819 else
820 {
821 return(NULL);
822 }
823 }
824
825
826 /****************************************************************************/
827 /****************************************************************************/
828 /****** G E N E R A L F O R M A T T I N G F U N C T I O N S *******/
829 /****************************************************************************/
830 /****************************************************************************/
831
832 /****************************************************************************/
833 /* gfRepChar(): */
834 /*--------------------------------------------------------------------------*/
835 /* DESCRIPTION */
836 /* Repeats a character to stdout the specified number of times. */
837 /* */
838 /* INPUTS */
839 /* c : The character to repeat. */
840 /* */
841 /* n : The number of times to repeat it. */
842 /****************************************************************************/
843 void gfRepChar(char c, unsigned n)
844 {
845 while(n--)
846 printf("%c", c);
847 }
848
849
850 /****************************************************************************/
851 /* gfHline(): */
852 /*--------------------------------------------------------------------------*/
853 /* DESCRIPTION */
854 /* Dumps a horizontal line of the standard length to the standard output. */
855 /****************************************************************************/
856 void gfHline(void)
857 {
858 gfRepChar(HORIZONTAL_BAR_SEP_CHAR, LINE_LEN);
859 printf("\n");
860 }
861
862
863 /****************************************************************************/
864 /* gfBannerHeading(): */
865 /*--------------------------------------------------------------------------*/
866 /* DESCRIPTION */
867 /* Prints a banner heading bracketed by astreisks to the standard output. */
868 /* This function is useful for separating different sections of output. */
869 /****************************************************************************/
870 void gfBannerHeading(char *s, int n_extra_lines)
871 {
872 const int lr_padding = 3;
873 /* The number of spaces on each side of what is printed.
874 */
875 int i;
876 /* General iteration variable.
877 */
878 int n_asterisks;
879 int input_arg_len;
880 int n_left_spaces;
881 int n_right_spaces;
882
883 /* The order of the source deck prevents me from using the assertion
884 ** function without function prototypes, and I don't like function
885 ** prototypes, and I'm too lazy to move the functions, so I'll just
886 ** protect against naughty parameters causing bizarre behavior.
887 */
888 if (!s)
889 s = "";
890 if (n_extra_lines < 0)
891 n_extra_lines = 0;
892
893 /* Print the right number of solid lines of asterisks to the
894 ** standard output.
895 */
896 for (i=0; i<n_extra_lines; i++)
897 {
898 gfRepChar('*', LINE_LEN);
899 printf("\n");
900 }
901
902 /* Figure out how many asterisks to print on each side of the
903 ** argument, and how many spaces. We also need to figure out
904 ** how many characters of the input argument to print--if there
905 ** are too many characters, we need to truncate.
906 */
907 input_arg_len = strlen(s);
908 if(input_arg_len > (LINE_LEN - 2 * lr_padding - 2))
909 input_arg_len = LINE_LEN - 2 * lr_padding - 2;
910
911 n_asterisks = (LINE_LEN - 2*lr_padding - input_arg_len)/2;
912
913 n_left_spaces = lr_padding;
914
915 if ((LINE_LEN - 2*lr_padding - input_arg_len) % 2)
916 {
917 /* Odd, need to pad the right by one. */
918 n_right_spaces = lr_padding+1;
919 }
920 else
921 {
922 n_right_spaces = lr_padding;
923 }
924
925 /* Print the text. */
926 gfRepChar('*', n_asterisks);
927 gfRepChar(' ', n_left_spaces);
928 for (i=0; i<input_arg_len; i++)
929 printf("%c", s[i]);
930 gfRepChar(' ', n_right_spaces);
931 gfRepChar('*', n_asterisks);
932 printf("\n");
933
934 /* Print the right number of solid lines of asterisks to the
935 ** standard output.
936 */
937 for (i=0; i<n_extra_lines; i++)
938 {
939 gfRepChar('*', LINE_LEN);
940 printf("\n");
941 }
942 }
943
944
945 /****************************************************************************/
946 /****************************************************************************/
947 /******* L O G I C A L A S S E R T I O N F U N C T I O N S ********/
948 /****************************************************************************/
949 /****************************************************************************/
950
951 /****************************************************************************/
952 /* asAssert(): */
953 /*--------------------------------------------------------------------------*/
954 /* DESCRIPTION */
955 /* Terminates the program abnormally if the logical condition specified */
956 /* is false. Roughly speaking, this fills the same role as the */
957 /* C-language assert() macro. */
958 /* */
959 /* INPUTS */
960 /* c : The logical condition tested (or, more precisely, the result of */
961 /* the code inserted by the compiler to stuff this parameter). */
962 /* */
963 /* l : The line number on which the call to this function occurs. */
964 /* This is filled in by using a special preprocessor symbol. */
965 /****************************************************************************/
966 void asAssert(unsigned c, unsigned long l)
967 {
968 if (!c)
969 {
970 gfHline();
971 printf("Abnormal termination; Assertion failed, line: %lu.\n", l);
972 gfHline();
973
974 /* Now need to convince the 'C' runtime environment to kill us.
975 */
976 exit(0);
977 }
978 }
979
980
981 /****************************************************************************/
982 /* asFatal(): */
983 /*--------------------------------------------------------------------------*/
984 /* DESCRIPTION */
985 /* Terminates the program with a fatal error message. */
986 /* */
987 /* INPUTS */
988 /* *msg: The message to use. Must be fairly short. */
989 /****************************************************************************/
990 void asFatal(const char *msg)
991 {
992 asAssert(msg != NULL, __LINE__);
993
994 gfHline();
995 printf("FATAL ERROR: %s.\n", msg);
996 gfHline();
997
998 /* Now need to convince the 'C' runtime environment to kill us.
999 */
1000 exit(0);
1001 }
1002
1003
1004 /****************************************************************************/
1005 /* asInfo(): */
1006 /*--------------------------------------------------------------------------*/
1007 /* DESCRIPTION */
1008 /* Prints informational message to stdout and terminates the program. */
1009 /****************************************************************************/
1010 void asInfo(void)
1011 {
1012 static char *msglines[] =
1013 {
1014 "RAP.EXE: rational approximation calculation program.",
1015 "!",
1016 "Copyright 2001, The Association For Computing Machinery (www.acm.org),",
1017 "David T. Ashley (dtashley@aol.com), Joseph P. DeVoe (jdevoe@visteon.com),",
1018 "Cory Pratt (cory_pratt@3com.com), Karl Perttunen (kperttun@visteon.com),",
1019 "and Anatoly Zhigljavsky (zhigljavskyaa@cardiff.ac.uk).",
1020 "!",
1021 "Please submit any bug reports or inquiries to all of the e-mail addresses",
1022 "listed above.",
1023 "!",
1024 "This program accompanies a paper entitled \"On Best Rational",
1025 "Approximations Using Large Integers\" submitted to ACM TOMS in late 2000/",
1026 "early 2001. The paper fully documents the algorithms employed in this",
1027 "software. In a nutshell, this program will execute algorithms on very",
1028 "large integers (up to 463 digits) and rational numbers with very large",
1029 "integer components which are helpful in finding best rational",
1030 "approximations in Farey series of very large order and in large rectangular",
1031 "regions of the integer lattice. The source code for this program should",
1032 "be available from CALGO (collected algorithms of the ACM).",
1033 "!",
1034 "Format legend: MANDATORY_PARAMETER",
1035 " [OPTIONAL_PARAMETER]",
1036 " {a, b ...} = Exactly one must be chosen from set.",
1037 " i... = Integer",
1038 " i+... = Non-Negative Integer",
1039 " i++... = Positive Integer",
1040 " r... = Rational Number",
1041 " r+... = Non-Negative Rational Number",
1042 " r++... = Positive Rational Number",
1043 "!",
1044 "How to specify an integer:",
1045 " [-]<digits>",
1046 " <digits>.<more_digits>e[{-,+}]<exponent>",
1047 " where the exponent creates an integer. Example: 3.14e2 is an integer",
1048 " but 3.14e1 is not (the former is 314, the latter is 31.4). Note that",
1049 " an integer can be substituted for a rational number, but a rational",
1050 " number cannot always be substituted for an integer.",
1051 "!",
1052 "How to specify a rational number:",
1053 " [-]<digits>/<more_digits>",
1054 " [-]<digits>.<more_digits>[e[{-,+}]<exponent>]",
1055 "!",
1056 "RAP + i1 i2",
1057 " Adds two integers to produce an integer result. Result is i1 + i2.",
1058 "!",
1059 "RAP + r1 r2",
1060 " Adds two rational numbers to produce a rational or integral result.",
1061 " Result is r1 + r2.",
1062 "!",
1063 "RAP - i1 i2",
1064 " Subtracts two integers to produce an integer result. Result is i1 - i2.",
1065 "!",
1066 "RAP * i1 i2",
1067 " Multiplies two integers to produce an integer result. Result is i1 * i2.",
1068 "!",
1069 "RAP / i1 i2",
1070 " Divides two integers to produce an integer result. Result is i1 / i2.",
1071 "!",
1072 "RAP / r1 r2",
1073 " Divides two rational numbers to produce a rational or integral result.",
1074 " Result is r1 / r2.",
1075 "!",
1076 "RAP % i1 i2",
1077 " Divides two integers and produces the remainder from the division. The",
1078 " result is the same as i1 % i2 as defined by the 'C' programming language",
1079 " except that much longer integers are accomodated.",
1080 "!",
1081 "RAP ** i1 i2+",
1082 " Exponentiates the integer i1 to the i2+'th power.",
1083 "!",
1084 "RAP ** r1 i1+",
1085 " Exponentiates the rational number r1 to the i1+'th power.",
1086 "!",
1087 "RAP gcd i1++ i2++",
1088 " Calculates the greatest common divisor of i1 and i2 using Euclid's",
1089 " algorithm.",
1090 "!",
1091 "RAP dap r1 i1+",
1092 " Given rational number r1, forms a rational approximation with a",
1093 " different denominator i1. This functionality is very useful",
1094 " for converting a rational number to a more familiar decimal",
1095 " approximation. The accompanying user's manual gives full",
1096 " details about the formula used by the DAP command. DAP is also",
1097 " used by several other commands to provide a more familiar",
1098 " output form.",
1099 "!",
1100 "RAP cf r1+",
1101 " Forms the continued fraction partial quotients and convergents",
1102 " of a non-negative rational number r1.",
1103 "!",
1104 "RAP fn r1+ i1++",
1105 " Finds the two best neighboring rational approximations to non-negative",
1106 " rational number r1 in the Farey series of order i1.",
1107 "!",
1108 "RAP fn r1++ i1++ i2++ i3++",
1109 " Same as form above except finds i2 neighbors on either side of r1",
1110 " (rather than the default of 1), and uses i3 as the denominator for",
1111 " presenting the approximations and errors in decimal form.",
1112 "!",
1113 "RAP mind r1+ r2+",
1114 " Calculates a rational number in the interval [r1, r2] with the minimum",
1115 " denominator. If the rational number in the interval with the minimum",
1116 " denominator is not unique (i.e. there is more than one with the same ",
1117 " denominator), no assurances are made about which of the numbers with",
1118 " minimum denominator will be returned (however, it is assured that the",
1119 " result will be one of them).",
1120 "!",
1121 "RAP fab r1+ i1++ i2++",
1122 " Finds the two neighbors to r1 in a rectangular area of the integer",
1123 " lattice specified by h<=i1 and k<=i2. This is a similar concept to",
1124 " the \"fn\" command, but with numerator and denominator both constrained.",
1125 " On output, 4 lines of digits after the decimal point are assumed by",
1126 " default.",
1127 "!",
1128 "RAP fab r1+ i1++ i2++ i3++ i4++",
1129 " Same as form above, except optional parameters i3 and i4 specify the",
1130 " number of neighbors on both the left and right to generate, and the",
1131 " denominator to use for display of DAP information.",
1132 "!",
1133 "RAP fndmax r1+ r2+ i++",
1134 " Provides an upper bound on the distance between successive terms of",
1135 " the Farey series of order i in the interval [r1, r2]. It is required",
1136 " that r1 < r2 and that r1 and r2 are both in the Farey series of order",
1137 " i.",
1138 "!",
1139 "RAP fndmax r1+ r2+ i++ j++",
1140 " Same as form above except j may override the default denominator of",
1141 " 1e108 for DAP presentation.",
1142 "!",
1143 "RAP fabdmax r1+ r2+ i1+ i2+",
1144 " Provides an upper bound on the distance between successive terms of",
1145 " the \"rectangular\" Farey series (doubly-constrained) with h <= i1 and",
1146 " j <= i2 in the interval [r1, r2]. It is required that r1 < r2 and that",
1147 " both meet the constraints (i.e. are in the \"rectangular\" series).",
1148 "!",
1149 "RAP fabdmax r1+ r2+ i1+ i2+ j++",
1150 " Same as form above except j may override the default denominator of 1e108",
1151 " for DAP presentation.",
1152 "!",
1153 "If this information has scrolled off the screen and can't be read, try re-",
1154 "directing it to a file (for example, RAP >OUT.TXT) and then viewing the file",
1155 "with a text editor, or piping it through \"more\" (for example, RAP|MORE).",
1156 "!"
1157 };
1158
1159 int i;
1160
1161 /* Dump out the whole array and terminate the program. Lines beginning
1162 ** with an exclamation point should be replaced with a horizontal line.
1163 */
1164 for (i=0; i<sizeof(msglines)/sizeof(msglines[0]); i++)
1165 {
1166 if ((msglines[i])[0] == '!')
1167 {
1168 gfHline();
1169 }
1170 else
1171 {
1172 printf("%s\n", msglines[i]);
1173 }
1174 }
1175
1176 exit(0);
1177 }
1178
1179
1180
1181 /****************************************************************************/
1182 /****************************************************************************/
1183 /*************** D A T A - D R I V E N F U N C T I O N S *************/
1184 /****************************************************************************/
1185 /****************************************************************************/
1186 /* This section is reserved for data-driven functions with no I/O
1187 ** requirements and no side-effects.
1188 */
1189 /****************************************************************************/
1190 /* ddCharToLower(): */
1191 /*--------------------------------------------------------------------------*/
1192 /* DESCRIPTION */
1193 /* If a character is upper-case, translates it to lower-case. */
1194 /****************************************************************************/
1195 char ddCharToLower(char c)
1196 {
1197 switch(c)
1198 {
1199 case 'A' : return ('a');
1200 break;
1201 case 'B' : return ('b');
1202 break;
1203 case 'C' : return ('c');
1204 break;
1205 case 'D' : return ('d');
1206 break;
1207 case 'E' : return ('e');
1208 break;
1209 case 'F' : return ('f');
1210 break;
1211 case 'G' : return ('g');
1212 break;
1213 case 'H' : return ('h');
1214 break;
1215 case 'I' : return ('i');
1216 break;
1217 case 'J' : return ('j');
1218 break;
1219 case 'K' : return ('k');
1220 break;
1221 case 'L' : return ('l');
1222 break;
1223 case 'M' : return ('m');
1224 break;
1225 case 'N' : return ('n');
1226 break;
1227 case 'O' : return ('o');
1228 break;
1229 case 'P' : return ('p');
1230 break;
1231 case 'Q' : return ('q');
1232 break;
1233 case 'R' : return ('r');
1234 break;
1235 case 'S' : return ('s');
1236 break;
1237 case 'T' : return ('t');
1238 break;
1239 case 'U' : return ('u');
1240 break;
1241 case 'V' : return ('v');
1242 break;
1243 case 'W' : return ('w');
1244 break;
1245 case 'X' : return ('x');
1246 break;
1247 case 'Y' : return ('y');
1248 break;
1249 case 'Z' : return ('z');
1250 break;
1251 default: return(c);
1252 break;
1253 }
1254 }
1255
1256
1257 /****************************************************************************/
1258 /* ddIsInfoChar(): */
1259 /*--------------------------------------------------------------------------*/
1260 /* DESCRIPTION */
1261 /* Returns TRUE if the passed character is a character that should be in */
1262 /* a token, or FALSE otherwise. */
1263 /****************************************************************************/
1264 unsigned ddIsInfoChar(char c)
1265 {
1266 int i;
1267 char d;
1268 const char CONST_info_characters[] = { 'a', 'b', 'c', 'd', 'e', 'f', 'g',
1269 'h', 'i', 'j', 'k', 'l', 'm', 'n',
1270 'o', 'p', 'q', 'r', 's', 't', 'u',
1271 'v', 'w', 'x', 'y', 'z', '0', '1',
1272 '2', '3', '4', '5', '6', '7', '8',
1273 '9', '+', '-', '*', '/', '_', '.',
1274 '\\','%' };
1275
1276 d = ddCharToLower(c);
1277
1278 for (i=0; i<sizeof(CONST_info_characters)/sizeof(CONST_info_characters[0]); i++)
1279 {
1280 if (CONST_info_characters[i] == d)
1281 return(TRUE);
1282 }
1283
1284 return(FALSE);
1285 }
1286
1287
1288 /****************************************************************************/
1289 /* ddIsDiscardChar(): */
1290 /*--------------------------------------------------------------------------*/
1291 /* DESCRIPTION */
1292 /* Returns TRUE if the passed character is a character that should be */
1293 /* discarded, or FALSE otherwise. */
1294 /****************************************************************************/
1295 unsigned ddIsDiscardChar(char c)
1296 {
1297 int i;
1298 char d;
1299 const char CONST_discard_characters[] = {','};
1300
1301 d = ddCharToLower(c);
1302
1303 for (i=0; i<sizeof(CONST_discard_characters)/sizeof(CONST_discard_characters[0]); i++)
1304 {
1305 if (CONST_discard_characters[i] == d)
1306 return(TRUE);
1307 }
1308
1309 return(FALSE);
1310 }
1311
1312
1313 /****************************************************************************/
1314 /* ddIsWhitespaceChar(): */
1315 /*--------------------------------------------------------------------------*/
1316 /* DESCRIPTION */
1317 /* Returns TRUE if the passed character is a character that should be */
1318 /* treated as whitespace, or FALSE otherwise. */
1319 /****************************************************************************/
1320 unsigned ddIsWhitespaceChar(char c)
1321 {
1322 int i;
1323 char d;
1324 const char CONST_whitespace_chars[] = { ' ', '\t', '\n' };
1325
1326 d = ddCharToLower(c);
1327
1328 for (i=0; i<sizeof(CONST_whitespace_chars)/sizeof(CONST_whitespace_chars[0]); i++)
1329 {
1330 if (CONST_whitespace_chars[i] == d)
1331 return(TRUE);
1332 }
1333
1334 return(FALSE);
1335 }
1336
1337
1338 /****************************************************************************/
1339 /* ddIsDigit(): */
1340 /*--------------------------------------------------------------------------*/
1341 /* DESCRIPTION */
1342 /* Returns TRUE if the passed argument is a digit, or FALSE otherwise. */
1343 /* */
1344 /* INPUTS */
1345 /* c : Character to evaluate for digitness (or is that digiality??). */
1346 /* */
1347 /* OUTPUTS */
1348 /* <--: TRUE if a digit, FALSE otherwise. */
1349 /****************************************************************************/
1350 unsigned ddIsDigit(char c)
1351 {
1352 if ((c >= '0') && (c <= '9'))
1353 return(TRUE);
1354 else
1355 return(FALSE);
1356 }
1357
1358
1359 /****************************************************************************/
1360 /* ddDigitToValue(): */
1361 /*--------------------------------------------------------------------------*/
1362 /* DESCRIPTION */
1363 /* Converts from ASCII digit to unsigned integer value. Fatal to call */
1364 /* with non-digit. */
1365 /* */
1366 /* INPUTS */
1367 /* c : Character to convert to value. */
1368 /* */
1369 /* OUTPUTS */
1370 /* <--: Value. */
1371 /****************************************************************************/
1372 unsigned ddDigitToValue(char c)
1373 {
1374 if ((c >= '0') && (c <= '9'))
1375 {
1376 return(c - '0');
1377 }
1378 else
1379 {
1380 asAssert(0, __LINE__); /* Guaranteed fatal. */
1381 return(0); /* To satisfy compiler warning about path not
1382 ** returning a value.
1383 */
1384 }
1385 }
1386
1387
1388 /****************************************************************************/
1389 /* ddValueToDigit(): */
1390 /*--------------------------------------------------------------------------*/
1391 /* DESCRIPTION */
1392 /* Converts from unsigned value to digit '0'..'9'. Fatal to call with */
1393 /* with value out of range. */
1394 /* */
1395 /* INPUTS */
1396 /* i : Character to convert to value. */
1397 /* */
1398 /* OUTPUTS */
1399 /* <--: Value. */
1400 /****************************************************************************/
1401 char ddValueToDigit(unsigned i)
1402 {
1403 if (i < 10)
1404 {
1405 return(i + '0');
1406 }
1407 else
1408 {
1409 asAssert(0, __LINE__); /* Guaranteed fatal. */
1410 return('0'); /* To satisfy compiler warning about path not
1411 ** returning a value.
1412 */
1413 }
1414 }
1415
1416
1417 /****************************************************************************/
1418 /* ddStringContains(): */
1419 /*--------------------------------------------------------------------------*/
1420 /* DESCRIPTION */
1421 /* Determines whether STR1 contains any characters in STR2, and returns */
1422 /* TRUE if so or FALSE otherwise. */
1423 /****************************************************************************/
1424 int ddStringContains(const char *str1, const char *str2)
1425 {
1426 unsigned s1len, s2len;
1427 unsigned i, j;
1428
1429 asAssert(str1 != NULL, __LINE__);
1430 asAssert(str2 != NULL, __LINE__);
1431
1432 s1len = strlen(str1);
1433 s2len = strlen(str2);
1434
1435 for (i=0; i<s1len; i++)
1436 {
1437 for (j=0; j<s2len; j++)
1438 {
1439 if (str1[i] == str2[j])
1440 return(TRUE);
1441 }
1442 }
1443
1444 return(FALSE);
1445 }
1446
1447
1448 /****************************************************************************/
1449 /* ddStringContainsOnly(): */
1450 /*--------------------------------------------------------------------------*/
1451 /* DESCRIPTION */
1452 /* Returns TRUE if STR1 contains ONLY characters from STR2, or FALSE */
1453 /* otherwise. */
1454 /****************************************************************************/
1455 int ddStringContainsOnly(const char *str1, const char *str2)
1456 {
1457 int i, j;
1458 int l1, l2;
1459
1460 asAssert(str1 != NULL, __LINE__);
1461 asAssert(str2 != NULL, __LINE__);
1462
1463 l1 = strlen(str1);
1464 l2 = strlen(str2);
1465
1466 for (i=0; i<l1; i++)
1467 {
1468 for (j=0; j<l2; j++)
1469 {
1470 if (str1[i] == str2[j])
1471 break;
1472 }
1473
1474 if (j == l2)
1475 return(FALSE);
1476 }
1477
1478 return(TRUE);
1479 }
1480
1481
1482 /****************************************************************************/
1483 /* ddStringReverse(): */
1484 /*--------------------------------------------------------------------------*/
1485 /* DESCRIPTION */
1486 /* Reverses the order of characters in the string s. */
1487 /****************************************************************************/
1488 void ddStringReverse(char *s)
1489 {
1490 unsigned l;
1491 unsigned i;
1492 unsigned limit;
1493 char temp;
1494
1495 asAssert(s != NULL, __LINE__);
1496
1497 l = strlen(s);
1498 limit = l / 2;
1499
1500 for (i=0; i<limit; i++)
1501 {
1502 temp = s[i];
1503 s[i] = s[l-i-1];
1504 s[l-i-1] = temp;
1505 }
1506 }
1507
1508
1509 /****************************************************************************/
1510 /* ddStringDeleteLeadingChar(): */
1511 /*--------------------------------------------------------------------------*/
1512 /* DESCRIPTION */
1513 /* Deletes the leading character of a string. If there is no leading */
1514 /* character, takes no action. */
1515 /****************************************************************************/
1516 void ddStringDeleteLeadingChar(char *s)
1517 {
1518 int len;
1519 int i;
1520
1521 asAssert(s != NULL, __LINE__);
1522
1523 len = strlen(s);
1524
1525 for (i=0; i<len; i++)
1526 s[i] = s[i+1];
1527 }
1528
1529
1530 /****************************************************************************/
1531 /* ddFundamentalAdditionCell(): */
1532 /*--------------------------------------------------------------------------*/
1533 /* DESCRIPTION */
1534 /* Fundamental cell of addition, used to add one ASCII digit to another */
1535 /* ASCII digit, processing carries in and carries out. Unexpected values */
1536 /* will generate a fatal error. */
1537 /* */
1538 /* INPUTS */
1539 /* digit1: Digit 1 in. */
1540 /* digit2: Digit 2 in. */
1541 /* carry_in: Carry in. Must be '0' or '1'. */
1542 /* */
1543 /* OUTPUTS */
1544 /* *digit_out: Digit result out. */
1545 /* *carry_out: Carry out. Will be '0' or '1'. */
1546 /****************************************************************************/
1547 void ddFundamentalAdditionCell(char digit1,
1548 char digit2,
1549 char carry_in,
1550 char *digit_out,
1551 char *carry_out)
1552 {
1553 unsigned vd1;
1554 unsigned vd2;
1555 unsigned vci;
1556 unsigned total;
1557
1558 asAssert(ddIsDigit(digit1), __LINE__);
1559 asAssert(ddIsDigit(digit2), __LINE__);
1560 asAssert(ddDigitToValue(carry_in) <= 1, __LINE__);
1561 asAssert(digit_out != NULL, __LINE__);
1562 asAssert(carry_out != NULL, __LINE__);
1563
1564 vd1 = ddDigitToValue(digit1);
1565 vd2 = ddDigitToValue(digit2);
1566 vci = ddDigitToValue(carry_in);
1567
1568 total = vd1+vd2+vci;
1569
1570 *digit_out = ddValueToDigit(total % 10);
1571 *carry_out = ddValueToDigit(total / 10);
1572 }
1573
1574
1575 /****************************************************************************/
1576 /* ddFundamentalMultiplicationCell(): */
1577 /*--------------------------------------------------------------------------*/
1578 /* DESCRIPTION */
1579 /* Fundamental cell of multiplication, used to multiply one ASCII digit */
1580 /* by another, processing carries in and out. Unexpected inputs will */
1581 /* generate a fatal error. */
1582 /* */
1583 /* INPUTS */
1584 /* digit1: Digit 1 in. */
1585 /* digit2: Digit 2 in. */
1586 /* carry_in: Carry in. Because of the way multiplication works, */
1587 /* this must always be from '0' through '8'. */
1588 /* */
1589 /* OUTPUTS */
1590 /* *digit_out: Digit result out. */
1591 /* *carry_out: Carry out. Because of the way multiplication works, */
1592 /* this must always be from '0' through '8'. */
1593 /****************************************************************************/
1594 void ddFundamentalMultiplicationCell(char digit1,
1595 char digit2,
1596 char carry_in,
1597 char *digit_out,
1598 char *carry_out)
1599 {
1600 unsigned vd1;
1601 unsigned vd2;
1602 unsigned vci;
1603 unsigned total;
1604
1605 asAssert(ddIsDigit(digit1), __LINE__);
1606 asAssert(ddIsDigit(digit2), __LINE__);
1607 asAssert(ddDigitToValue(carry_in) <= 8, __LINE__);
1608 asAssert(digit_out != NULL, __LINE__);
1609 asAssert(carry_out != NULL, __LINE__);
1610
1611 vd1 = ddDigitToValue(digit1);
1612 vd2 = ddDigitToValue(digit2);
1613 vci = ddDigitToValue(carry_in);
1614
1615 total = (vd1 * vd2) + vci;
1616
1617 *digit_out = ddValueToDigit(total % 10);
1618 *carry_out = ddValueToDigit(total / 10);
1619 }
1620
1621
1622 /****************************************************************************/
1623 /* ddFundamentalSubtractionCell(): */
1624 /*--------------------------------------------------------------------------*/
1625 /* DESCRIPTION */
1626 /* Fundamental cell of subtraction, used to subtract one ASCII digit from */
1627 /* another ASCII digit, processing borrows in and borrows out. */
1628 /* Unexpected values will generate a fatal error. */
1629 /* */
1630 /* INPUTS */
1631 /* digit1: Digit 1 in. */
1632 /* digit2: Digit 2 in. */
1633 /* borrow_in: Borrow in. Must be '0' or '1'. */
1634 /* */
1635 /* OUTPUTS */
1636 /* *digit_out: Digit result out, digit1-digit2-borrow_in. */
1637 /* *borrow_out: Borrow out. Will be '0' or '1'. */
1638 /****************************************************************************/
1639 void ddFundamentalSubtractionCell(char digit1,
1640 char digit2,
1641 char borrow_in,
1642 char *digit_out,
1643 char *borrow_out)
1644 {
1645 unsigned vd1;
1646 unsigned vd2;
1647 unsigned vbi;
1648 unsigned total;
1649
1650 asAssert(ddIsDigit(digit1), __LINE__);
1651 asAssert(ddIsDigit(digit2), __LINE__);
1652 asAssert(ddDigitToValue(borrow_in) <= 1, __LINE__);
1653 asAssert(digit_out != NULL, __LINE__);
1654 asAssert(borrow_out != NULL, __LINE__);
1655
1656 vd1 = ddDigitToValue(digit1);
1657 vd2 = ddDigitToValue(digit2);
1658 vbi = ddDigitToValue(borrow_in);
1659
1660 total = 10+vd1-vd2-vbi;
1661
1662 *digit_out = ddValueToDigit(total % 10);
1663 if (total < 10)
1664 {
1665 *borrow_out = '1';
1666 /* If we went neggy neggy on the subtraction, we need a borrow.
1667 */
1668 }
1669 else
1670 {
1671 *borrow_out = '0'; /* No borrow. */
1672 }
1673 }
1674
1675
1676 /****************************************************************************/
1677 /* ddUmin(): */
1678 /*--------------------------------------------------------------------------*/
1679 /* DESCRIPTION */
1680 /* Returns the minimum of two unsigned integers. */
1681 /* */
1682 /* INPUTS */
1683 /* arg1, arg2 : Unsigned integers to compare. */
1684 /* */
1685 /* OUTPUTS */
1686 /* <-- : Samller of arg1, arg2. */
1687 /****************************************************************************/
1688 unsigned ddUmin(unsigned arg1, unsigned arg2)
1689 {
1690 if (arg1 < arg2)
1691 return(arg1);
1692 else
1693 return(arg2);
1694 }
1695
1696
1697
1698 /****************************************************************************/
1699 /* ddUmax(): */
1700 /*--------------------------------------------------------------------------*/
1701 /* DESCRIPTION */
1702 /* Returns the maximum of two unsigned integers. */
1703 /* */
1704 /* INPUTS */
1705 /* arg1, arg2 : Unsigned integers to compare. */
1706 /* */
1707 /* OUTPUTS */
1708 /* <-- : Larger of arg1, arg2. */
1709 /****************************************************************************/
1710 unsigned ddUmax(unsigned arg1, unsigned arg2)
1711 {
1712 if (arg1 > arg2)
1713 return(arg1);
1714 else
1715 return(arg2);
1716 }
1717
1718
1719 /****************************************************************************/
1720 /* ddSmin(): */
1721 /*--------------------------------------------------------------------------*/
1722 /* DESCRIPTION */
1723 /* Returns the minimum of two signed integers. */
1724 /* */
1725 /* INPUTS */
1726 /* arg1, arg2 : Signed integers to compare. */
1727 /* */
1728 /* OUTPUTS */
1729 /* <-- : Samller of arg1, arg2. */
1730 /****************************************************************************/
1731 int ddSmin(int arg1, int arg2)
1732 {
1733 if (arg1 < arg2)
1734 return(arg1);
1735 else
1736 return(arg2);
1737 }
1738
1739
1740
1741 /****************************************************************************/
1742 /* ddSmax(): */
1743 /*--------------------------------------------------------------------------*/
1744 /* DESCRIPTION */
1745 /* Returns the maximum of two signed integer. */
1746 /* */
1747 /* INPUTS */
1748 /* arg1, arg2 : Signed integers to compare. */
1749 /* */
1750 /* OUTPUTS */
1751 /* <-- : Larger of arg1, arg2. */
1752 /****************************************************************************/
1753 int ddSmax(int arg1, int arg2)
1754 {
1755 if (arg1 > arg2)
1756 return(arg1);
1757 else
1758 return(arg2);
1759 }
1760
1761
1762 /****************************************************************************/
1763 /****************************************************************************/
1764 /******* M E M O R Y A L L O C A T I O N F U N C T I O N S ********/
1765 /****************************************************************************/
1766 /****************************************************************************/
1767 /* These functions are primarily wrappers for the standard 'C' library
1768 ** functions. They will trap "out of memory" errors, and also provide a
1769 ** place to insert debugging code if it becomes necessary.
1770 */
1771 /****************************************************************************/
1772 /* maMalloc(): */
1773 /*--------------------------------------------------------------------------*/
1774 /* DESCRIPTION */
1775 /* Wrapper for malloc(). See malloc() documentation. */
1776 /****************************************************************************/
1777 void *maMalloc(size_t n)
1778 {
1779 void *rv;
1780
1781 rv = malloc(n);
1782
1783 asAssert(rv != NULL, __LINE__);
1784
1785 return(rv);
1786 }
1787
1788
1789 /****************************************************************************/
1790 /* maFree(): */
1791 /*--------------------------------------------------------------------------*/
1792 /* DESCRIPTION */
1793 /* Wrapper for free(). See free() documentation. */
1794 /****************************************************************************/
1795 void maFree(void *p)
1796 {
1797 asAssert(p != NULL, __LINE__);
1798
1799 free(p);
1800 }
1801
1802
1803 /****************************************************************************/
1804 /* maRealloc(): */
1805 /*--------------------------------------------------------------------------*/
1806 /* DESCRIPTION */
1807 /* Wrapper for realloc(). See realloc() documentation. */
1808 /****************************************************************************/
1809 void *maRealloc(void *ptr, size_t n)
1810 {
1811 void *rv;
1812
1813 rv = realloc(ptr, n);
1814
1815 asAssert(rv != NULL, __LINE__);
1816
1817 return(rv);
1818 }
1819
1820
1821 /****************************************************************************/
1822 /****************************************************************************/
1823 /******* S Y N T H E T I C I N T E G E R F U N C T I O N S ********/
1824 /****************************************************************************/
1825 /****************************************************************************/
1826 /* Definition of a synthetic integer.
1827 */
1828 struct synthetic_integer_struct
1829 {
1830 unsigned nan;
1831 /* Boolean variable used to remember errors that render the integer
1832 ** as "non a number", so that its value is invalid or unknown.
1833 ** NAN's propagate so that performing any operating involving a NAN
1834 ** also yields a NAN. A non-zero value here means that an error has
1835 ** been thrown and the synthetic integer is invalid.
1836 */
1837 unsigned neg;
1838 /* Only the absolute value of the integer is stored as a character
1839 ** string. This boolean remembers if the number is negative.
1840 ** A non-zero value here means the integer is negative. This may have
1841 ** either value when the value of zero is represented by the
1842 ** synthetic integer.
1843 */
1844 unsigned len;
1845 /* The number of characters in the string representing the integer.
1846 ** This is the same as would be returned by strlen(). This field
1847 ** is maintained to avoid the expense of counting characters before
1848 ** the terminator each time the length is needed.
1849 */
1850 unsigned char digits[INTERMEDIATE_CALC_MAX_DIGITS + 1];
1851 /* The string representing the synthetic integer. To avoid the need
1852 ** to move large blocks of characters, the integer is stored with the
1853 ** least significant digits first. The string must be \0 terminated--
1854 ** this is done primarily to make printing things for debugging
1855 ** easier, although the terminator is redundant with the "len" field
1856 ** above. The valid representation of zero is a "len" field of
1857 ** zero and a terminating \0 as element [0] here. No commas or other
1858 ** characters besides digits are allowed.
1859 */
1860 };
1861
1862
1863 typedef struct synthetic_integer_struct SYNTHETIC_INTEGER;
1864 /* Typedef for more compact reference to data type.
1865 */
1866
1867
1868 /****************************************************************************/
1869 /* siCreate(): */
1870 /*--------------------------------------------------------------------------*/
1871 /* DESCRIPTION */
1872 /* Creates a synthetic integer (allocates the memory and initializes) and */
1873 /* sets value to zero. */
1874 /* */
1875 /* INPUTS */
1876 /* *i : Pointer to pointer to synthetic integer. The caller's pointer */
1877 /* will be filled in to point to the allocated object. */
1878 /****************************************************************************/
1879 void siCreate(SYNTHETIC_INTEGER **i)
1880 {
1881 asAssert(i != NULL, __LINE__);
1882 /* Input pointer may not be NULL.
1883 */
1884
1885 *i = maMalloc(sizeof(SYNTHETIC_INTEGER));
1886 /* Allocate the memory, stuff the pointer in the caller's
1887 ** area.
1888 */
1889
1890 /* Initialize the synthetic integer to what is required for
1891 ** the number zero.
1892 */
1893 (*i)->nan = 0;
1894 (*i)->neg = 0;
1895 (*i)->len = 0;
1896 (*i)->digits[0] = 0;
1897 }
1898
1899
1900 /****************************************************************************/
1901 /* siDestroy(): */
1902 /*--------------------------------------------------------------------------*/
1903 /* DESCRIPTION */
1904 /* Destroys a synthetic integer (deallocates the memory and sets the */
1905 /* caller's pointer to NULL. */
1906 /* */
1907 /* INPUTS */
1908 /* *i : Pointer to pointer to synthetic integer. The synthetic integer */
1909 /* will be destroyed and the caller's pointer will be filled to */
1910 /* NULL. */
1911 /****************************************************************************/
1912 void siDestroy(SYNTHETIC_INTEGER **i)
1913 {
1914 asAssert(i != NULL, __LINE__);
1915 /* Input pointer to pointer may not be NULL.
1916 */
1917 asAssert(*i != NULL, __LINE__);
1918 /* Input pointer to block of memory may not be NULL.
1919 */
1920
1921 maFree(*i);
1922 /* Deallocate block of memory.
1923 */
1924
1925 *i = NULL;
1926 /* Fill in the caller's variable to NULL.
1927 */
1928 }
1929
1930
1931 /****************************************************************************/
1932 /* siCopy(): */
1933 /*--------------------------------------------------------------------------*/
1934 /* DESCRIPTION */
1935 /* Copies a synthetic integer. The integer must already be created. */
1936 /* */
1937 /* INPUTS */
1938 /* *src, *dst : Source and destination synthetic integers. */
1939 /****************************************************************************/
1940 void siCopy(SYNTHETIC_INTEGER **src, SYNTHETIC_INTEGER **dst)
1941 {
1942 asAssert(src != NULL, __LINE__);
1943 asAssert(*src != NULL, __LINE__);
1944 asAssert(dst != NULL, __LINE__);
1945 asAssert(*dst != NULL, __LINE__);
1946
1947 /* There are a lot of ways to copy, and pros and cons.
1948 ** I'll choose the one that gives the least typing
1949 ** here in the source code, but maybe I'll need to
1950 ** optimize it later.
1951 */
1952 memcpy(*dst, *src, sizeof(SYNTHETIC_INTEGER));
1953 }
1954
1955
1956 /****************************************************************************/
1957 /* siMulByTen(): */
1958 /*--------------------------------------------------------------------------*/
1959 /* DESCRIPTION */
1960 /* Multiplies a synthetic integer by 10, producing a NAN result if there */
1961 /* is an overflow in the number of digits. This is a primitive operation */
1962 /* used in multiplication of arbitrary integers. */
1963 /****************************************************************************/
1964 void siMulByTen(SYNTHETIC_INTEGER **arg)
1965 {
1966 asAssert(arg != NULL, __LINE__);
1967 asAssert(*arg != NULL, __LINE__);
1968
1969 /* If the integer is already NAN, NAN it stays.
1970 */
1971 if (!((*arg)->nan))
1972 {
1973 /* If the integer is 0, 0 it stays. */
1974 if ((*arg)->len)
1975 {
1976 /* The result will overflow iff the string is already
1977 ** at its length limit. If so, must declare a NAN. Multiplication
1978 ** by 10 is guaranteed to add exactly one digit.
1979 */
1980 asAssert((*arg)->len <= INTERMEDIATE_CALC_MAX_DIGITS, __LINE__);
1981 if ((*arg)->len == INTERMEDIATE_CALC_MAX_DIGITS)
1982 {
1983 /* Overflow. Must declare a NAN. */
1984 (*arg)->digits[0] = '\0';
1985 (*arg)->len = 0;
1986 (*arg)->nan = TRUE;
1987 (*arg)->neg = FALSE;
1988 }
1989 else
1990 {
1991 /* The digit shift will come off alright. Carry it out.
1992 */
1993 unsigned tgt_idx;
1994
1995 for (tgt_idx = (*arg)->len; tgt_idx >= 1; tgt_idx--)
1996 {
1997 (*arg)->digits[tgt_idx] = (*arg)->digits[tgt_idx-1];
1998 }
1999 (*arg)->digits[0] = '0';
2000 (*arg)->digits[(*arg)->len + 1] = '\0';
2001 ((*arg)->len)++;
2002 }
2003 }
2004 }
2005 }
2006
2007
2008 /****************************************************************************/
2009 /* siDivByTen(): */
2010 /*--------------------------------------------------------------------------*/
2011 /* DESCRIPTION */
2012 /* Divides a synthetic integer by 10 using digit-shifting (remainders */
2013 /* are discarded. This function is used in division of arbitrary */
2014 /* integers. */
2015 /****************************************************************************/
2016 void siDivByTen(SYNTHETIC_INTEGER **arg)
2017 {
2018 asAssert(arg != NULL, __LINE__);
2019 asAssert(*arg != NULL, __LINE__);
2020
2021 /* If the integer is already NAN, NAN it stays.
2022 */
2023 if (!((*arg)->nan))
2024 {
2025 /* If the integer is a single digit, it must go to zero.
2026 */
2027 if ((*arg)->len == 1)
2028 {
2029 (*arg)->len = 0;
2030 (*arg)->digits[0] = '\0';
2031 (*arg)->neg = FALSE;
2032 }
2033 else
2034 {
2035 unsigned i;
2036
2037 asAssert((*arg)->len <= INTERMEDIATE_CALC_MAX_DIGITS, __LINE__);
2038
2039 /* This is just a normal digit shift. */
2040 for (i=0; i<((*arg)->len); i++)
2041 {
2042 (*arg)->digits[i] = (*arg)->digits[i+1];
2043 }
2044
2045 ((*arg)->len)--;
2046 }
2047 }
2048 }
2049
2050
2051 /****************************************************************************/
2052 /* siMulByDigit(): */
2053 /*--------------------------------------------------------------------------*/
2054 /* DESCRIPTION */
2055 /* Multiplies a synthetic integer by an arbitrary digit, producing a */
2056 /* result that may be zero, the argument, or an integer of larger mag- */
2057 /* nitude. A NAN is declared on overflow. This is a primitive operation */
2058 /* used in multiplication of arbitrary integers. */
2059 /****************************************************************************/
2060 void siMulByDigit(SYNTHETIC_INTEGER **arg, char digit)
2061 {
2062 asAssert(arg != NULL, __LINE__);
2063 asAssert(*arg != NULL, __LINE__);
2064 asAssert(ddIsDigit(digit), __LINE__);
2065
2066 /* If the integer is already NAN, NAN it stays.
2067 */
2068 if (!((*arg)->nan))
2069 {
2070 /* If the integer is 0 or the digit is '0', force the result
2071 ** to 0.
2072 */
2073 if (!((*arg)->len) || (digit == '0'))
2074 {
2075 (*arg)->digits[0] = '\0';
2076 (*arg)->len = 0;
2077 (*arg)->nan = FALSE;
2078 (*arg)->neg = FALSE;
2079 }
2080 else
2081 {
2082 /* This is now a valid multiplication. Regrettably, we can't
2083 ** determine an overflow or lack of decisively except by carrying
2084 ** the multiplication out to its end.
2085 */
2086 char carry;
2087 char new_digit;
2088 unsigned cur_digit;
2089
2090 carry = '0';
2091
2092 asAssert((*arg)->len <= INTERMEDIATE_CALC_MAX_DIGITS, __LINE__);
2093
2094 for (cur_digit=0; cur_digit < (*arg)->len; cur_digit++)
2095 {
2096 ddFundamentalMultiplicationCell((*arg)->digits[cur_digit],
2097 digit,
2098 carry,
2099 &new_digit,
2100 &carry);
2101
2102 if (cur_digit == ((*arg)->len - 1)) /* If last digit in loop. */
2103 {
2104 if (cur_digit == (INTERMEDIATE_CALC_MAX_DIGITS-1)) /* arg was max len. */
2105 {
2106 if (carry != '0')
2107 {
2108 /* We have an SI that was at the length limit, we've processed.
2109 ** the last digit, and there is a carry out. Must declare a NAN.
2110 */
2111 (*arg)->digits[0] = '\0';
2112 (*arg)->len = 0;
2113 (*arg)->nan = TRUE;
2114 (*arg)->neg = FALSE;
2115 }
2116 else
2117 {
2118 /* We have an SI that was max length, but did not overflow.
2119 ** Close it up normally.
2120 */
2121 (*arg)->digits[cur_digit] = new_digit;
2122 /* Nothing else to do. Everything is preset. */
2123 }
2124 }
2125 else
2126 {
2127 /* We are at the last digit of multiplication, but not at the
2128 ** maximum length of a string of digits. Cleanup without fear
2129 ** of overflow.
2130 */
2131 if (carry == '0')
2132 {
2133 /* Number stayed same length. Just assign digit. Terminator
2134 ** and other info is still valid.
2135 */
2136 (*arg)->digits[cur_digit] = new_digit;
2137 }
2138 else
2139 {
2140 /* At last digit of multiplication. Number grew by one digit.
2141 ** Clean up.
2142 */
2143 (*arg)->digits[cur_digit] = new_digit;
2144 (*arg)->digits[cur_digit+1] = carry;
2145 (*arg)->digits[cur_digit+2] = '\0';
2146 ((*arg)->len)++;
2147
2148 /* Must break out of the for() loop. Otherwise, modifying
2149 ** the length above will keep us going.
2150 */
2151 break;
2152 }
2153 }
2154 }
2155 else
2156 {
2157 /* Not the last digit in the loop. Just assign digit.
2158 */
2159 (*arg)->digits[cur_digit] = new_digit;
2160 }
2161 } /* End for() each digit. */
2162 }
2163 }
2164 }
2165
2166
2167 /****************************************************************************/
2168 /* siCompareAbs(): */
2169 /*--------------------------------------------------------------------------*/
2170 /* DESCRIPTION */
2171 /* Compares the absolute value of two synthetic integers (this involves */
2172 /* comparing only the string representation without considering the */
2173 /* sign). This function is useful from within the full compare function */
2174 /* to avoid duplication of code--it isn't very useful in any other */
2175 /* context. */
2176 /* */
2177 /* INPUTS */
2178 /* *arg1, *arg2 : Two synthetic integers whose absolute values should */
2179 /* be compared. */
2180 /* */
2181 /* OUTPUT */
2182 /* <-- : (-1) : ABS(*arg1) < ABS(*arg2) */
2183 /* (0) : ABS(*arg1) = ABS(*arg2) */
2184 /* (1) : ABS(*arg1) > ABS(*arg2) */
2185 /****************************************************************************/
2186 int siCompareAbs(SYNTHETIC_INTEGER **arg1, SYNTHETIC_INTEGER **arg2)
2187 {
2188 int i;
2189
2190 asAssert(arg1 != NULL, __LINE__);
2191 asAssert(*arg1 != NULL, __LINE__);
2192 asAssert(arg2 != NULL, __LINE__);
2193 asAssert(*arg2 != NULL, __LINE__);
2194
2195 if ((*arg1)->len < (*arg2)->len)
2196 {
2197 return(-1);
2198 }
2199 else if ((*arg1)->len > (*arg2)->len)
2200 {
2201 return(1);
2202 }
2203 else if (!((*arg1)->len) && !((*arg2)->len))
2204 {
2205 /* Both are zero, return equal. */
2206 return(0);
2207 }
2208 else
2209 {
2210 /* Lengths of strings are equal--must look at the strings in
2211 ** more detail.
2212 */
2213 for (i = (*arg1)->len; i >= 0; i--)
2214 {
2215 if ((*arg1)->digits[i] < (*arg2)->digits[i])
2216 {
2217 return(-1);
2218 }
2219 else if ((*arg1)->digits[i] > (*arg2)->digits[i])
2220 {
2221 return(1);
2222 }
2223 }
2224
2225 /* Have iterated through the strings--couldn't find any differences.
2226 ** They are equal. This also covers the zero case.
2227 */
2228 return(0);
2229 }
2230 }
2231
2232
2233 /****************************************************************************/
2234 /* siCompare(): */
2235 /*--------------------------------------------------------------------------*/
2236 /* DESCRIPTION */
2237 /* Establishes the logical ordering of two synthetic integers. */
2238 /* */
2239 /* INPUTS */
2240 /* *arg1, *arg2 : Two synthetic integers to compare. */
2241 /* */
2242 /* OUTPUT */
2243 /* <-- : (-1) : *arg1 < *arg2. */
2244 /* (0) : *arg1 = *arg2 */
2245 /* (1) : *arg1 > *arg2 */
2246 /****************************************************************************/
2247 int siCompare(SYNTHETIC_INTEGER **arg1, SYNTHETIC_INTEGER **arg2)
2248 {
2249 asAssert(arg1 != NULL, __LINE__);
2250 asAssert(*arg1 != NULL, __LINE__);
2251 asAssert(arg2 != NULL, __LINE__);
2252 asAssert(*arg2 != NULL, __LINE__);
2253
2254 /* If either integer has been marked NAN, we shouldn't be
2255 ** comparing, and there is no reason to handle this case,
2256 ** as it is a minority case and harmless whatever happens.
2257 */
2258 /* Cover the two easiest cases, which are where one operand is
2259 ** positive and the other is negative.
2260 */
2261 if (((*arg1)->neg) && (!((*arg2)->neg)))
2262 {
2263 return(-1);
2264 }
2265 else if ((!((*arg1)->neg)) && ((*arg2)->neg))
2266 {
2267 return(1);
2268 }
2269 /* We know both integers are of the same sign. We can call the
2270 ** absolute value compare function to get their ordering.
2271 */
2272 else if ((*arg1)->neg)
2273 {
2274 return(siCompareAbs(arg2, arg1)); /* Reverse the order to negate the
2275 ** results provided by
2276 ** siCompareAbs().
2277 */
2278 }
2279 else if (!((*arg1)->neg))
2280 {
2281 return(siCompareAbs(arg1, arg2)); /* Normal non-neg compare.
2282 */
2283 }
2284 else
2285 {
2286 /* It should be impossible to reach this case.
2287 */
2288 asAssert(0, __LINE__); /* Fatal. */
2289 return(0); /* Suppress compiler warning. */
2290 }
2291 }
2292
2293
2294 /****************************************************************************/
2295 /* siDump(): */
2296 /*--------------------------------------------------------------------------*/
2297 /* DESCRIPTION */
2298 /* Outputs a synthetic integer to the standard output stream. The output */
2299 /* includes a title, and provisions are made for NAN, for a description, */
2300 /* and prints commas. It is the caller's responsibility to provide lines */
2301 /* before and lines after if desired. Ths function assumes that it */
2302 /* starts with the cursor in column 1. */
2303 /* */
2304 /* INPUTS */
2305 /* **si : The synthetic integer which is being output. */
2306 /* */
2307 /* *desc : The description to use for the integer. It will be */
2308 /* right-justified up against the start of the integer, */
2309 /* and this function will automatically include a */
2310 /* colon. If the description is too long, this function */
2311 /* will truncate it to fit the space available. */
2312 /****************************************************************************/
2313 void siDump(SYNTHETIC_INTEGER **si, const char *desc)
2314 {
2315 unsigned cur_line;
2316 unsigned nlines;
2317 unsigned digits_per_line;
2318
2319 /* Make sure the caller isn't doing something bad for the program's health.
2320 */
2321 asAssert(si != NULL, __LINE__);
2322 asAssert(*si != NULL, __LINE__);
2323 asAssert(desc != NULL, __LINE__);
2324
2325 /* The number of digits per line that we assume must be a multiple of
2326 ** three. Round this up in case the preprocessor constant was set
2327 ** dubiously.
2328 */
2329 digits_per_line = ddUmax(3, ((DIGITS_PER_LINE + 2) /3) * 3);
2330
2331 /* As the first order of business, figure out how many lines the beast
2332 ** will require.
2333 */
2334 if ((*si)->nan)
2335 {
2336 nlines = 1; /* Only one line required for NAN verbeage. */
2337 }
2338 else if (!((*si)->len))
2339 {
2340 nlines = 1; /* The zero value requires one line. */
2341 }
2342 else
2343 {
2344 /* In any other case, have a formula.
2345 */
2346 nlines = 1 + ((*si)->len - 1) / digits_per_line;
2347 }
2348
2349 /* Iterate through each line, spitting out whatever is appropriate. */
2350 for (cur_line = 0; cur_line < nlines; cur_line++)
2351 {
2352 unsigned cur_digit_on_line;
2353
2354 /* If this is the first line, spit out the description, right-aligned.
2355 ** Otherwise, spit spaces.
2356 */
2357 if (!cur_line)
2358 {
2359 /* First line. */
2360 unsigned len;
2361
2362 len = strlen(desc);
2363
2364 if (len <= NUMBER_DESC_WIDTH)
2365 {
2366 /* Description is shorter or equal, pad on left. */
2367 gfRepChar(' ', NUMBER_DESC_WIDTH - len);
2368 printf("%s", desc);
2369 }
2370 else
2371 {
2372 /* Description is too long, truncate. */
2373 unsigned i;
2374
2375 for (i=0; i<len; i++)
2376 printf("%c", desc[i]);
2377 }
2378
2379 printf(": ");
2380
2381 /* If the number is negative, throw in a minus sign. */
2382 if ((*si)->neg && !(*si)->nan)
2383 {
2384 printf("- ");
2385 }
2386 else
2387 {
2388 printf(" ");
2389 }
2390 }
2391 else
2392 {
2393 /* Every line but first line. */
2394 gfRepChar(' ', NUMBER_DESC_WIDTH+4);
2395 }
2396
2397 for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++)
2398 {
2399 unsigned idx_into_string;
2400 /* Index into the string which is our digit of interest.
2401 */
2402
2403 /* Compute the index. The equation is based on the ordering
2404 ** of presentation, for example,
2405 **
2406 ** 7 6
2407 ** 5 4 3
2408 ** 2 1 0.
2409 **
2410 ** With a little thought, the equation should make sense.
2411 ** The index won't always be used to index into the string.
2412 */
2413 idx_into_string =
2414 (((nlines-1) - cur_line) * digits_per_line)
2415 +
2416 (digits_per_line - 1 - cur_digit_on_line);
2417
2418 /* Print the appropriate digit or a space. The NAN case and the
2419 ** zero case need to be treated specially.
2420 */
2421 if ((*si)->nan)
2422 {
2423 /* Not a number. Everything is blank, except spell out NAN
2424 ** at the end of the string of digits.
2425 */
2426 if (cur_digit_on_line == (digits_per_line - 3))
2427 {
2428 printf("N");
2429 }
2430 else if (cur_digit_on_line == (digits_per_line - 2))
2431 {
2432 printf("A");
2433 }
2434 else if (cur_digit_on_line == (digits_per_line - 1))
2435 {
2436 printf("N");
2437 }
2438 else
2439 {
2440 printf(" ");
2441 }
2442 }
2443 else if (!(*si)->len)
2444 {
2445 /* This is the zero case. For zero, there is only one line,
2446 ** and every character except the last one is a blank.
2447 */
2448 if (cur_digit_on_line == (digits_per_line - 1))
2449 {
2450 printf("0");
2451 }
2452 else
2453 {
2454 printf(" ");
2455 }
2456 }
2457 else
2458 {
2459 /* This is a valid number which is not zero. Need to print
2460 ** the digits.
2461 */
2462
2463 if (idx_into_string < (*si)->len)
2464 {
2465 printf("%c", (*si)->digits[idx_into_string]);
2466 }
2467 else
2468 {
2469 printf(" ");
2470 }
2471 } /* End of digit case.
2472
2473 /* Now handle the commas. The rules for commas are straightforward.
2474 ** a)NAN never has a comma.
2475 ** b)Zeros never have a comma.
2476 ** c)The final line, last digit never has a comma.
2477 ** d)Everything else in multiples of three ...
2478 */
2479 if (!(idx_into_string % 3) && (idx_into_string))
2480 {
2481 if ((*si)->nan)
2482 {
2483 printf(" ");
2484 }
2485 else if (!(*si)->len)
2486 {
2487 printf(" ");
2488 }
2489 else
2490 {
2491 if (idx_into_string < (*si)->len)
2492 {
2493 printf(",");
2494 }
2495 else
2496 {
2497 printf(" ");
2498 }
2499 }
2500 }
2501 } /* End for each digit on the current line. */
2502
2503 /* For the first line, print out an informative message
2504 ** advising of the number of digits. For all other lines
2505 ** print nothing.
2506 */
2507 if (!cur_line && !((*si)->nan))
2508 {
2509 if (nlines == 1)
2510 printf(" ");
2511
2512 if ((*si)->len <= 1)
2513 {
2514 printf(" ( 1 digit)\n");
2515 }
2516 else if ((*si)->len < 1000)
2517 {
2518 printf(" (%5u digits)\n", (*si)->len);
2519 }
2520 else
2521 {
2522 printf(" (%u,%03u digits)\n", (*si)->len / 1000, (*si)->len % 1000);
2523 }
2524 }
2525 else
2526 {
2527 printf("\n");
2528 }
2529 } /* End for each line. */
2530 }
2531
2532
2533 /****************************************************************************/
2534 /* siSetToLong(): */
2535 /*--------------------------------------------------------------------------*/
2536 /* DESCRIPTION */
2537 /* Sets a synthetic integer (which must already be created) to a long */
2538 /* value. This is used primarily in preliminary unit testing. */
2539 /* */
2540 /* INPUTS */
2541 /* **si : The synthetic integer to operate on. */
2542 /* */
2543 /* val : The value to set the synthetic integer to. */
2544 /****************************************************************************/
2545 void siSetToLong(SYNTHETIC_INTEGER **si, long val)
2546 {
2547 char buf[100];
2548 char *p;
2549 unsigned len;
2550 unsigned idx;
2551
2552 asAssert(si != NULL, __LINE__);
2553 asAssert(*si != NULL, __LINE__);
2554
2555 sprintf(buf, "%ld", val);
2556
2557 if (buf[0] == '-')
2558 p = buf+1;
2559 else
2560 p = buf;
2561
2562 len = strlen(p);
2563
2564 (*si)->nan = FALSE;
2565
2566 if (val == 0)
2567 {
2568 (*si)->nan = FALSE;
2569 (*si)->neg = FALSE;
2570 (*si)->len = 0;
2571 (*si)->digits[0] = '\0';
2572 }
2573 else
2574 {
2575 if (val < 0)
2576 {
2577 (*si)->neg = TRUE;
2578 }
2579 else
2580 {
2581 (*si)->neg = FALSE;
2582 }
2583
2584 (*si)->len = len;
2585
2586 for (idx=0; idx<len; idx++)
2587 {
2588 (*si)->digits[idx] = p[len - idx - 1];
2589 }
2590
2591 (*si)->digits[len] = '\0';
2592 }
2593 }
2594
2595
2596 /****************************************************************************/
2597 /* siSetToNan(): */
2598 /*--------------------------------------------------------------------------*/
2599 /* DESCRIPTION */
2600 /* Sets a synthetic integer (which must already be created) to NAN (not */
2601 /* a number). */
2602 /* */
2603 /* INPUTS */
2604 /* **si : The synthetic integer to set to NAN. */
2605 /* */
2606 /* NOTES */
2607 /* This program "evolved", and this function evolved after many of the */
2608 /* functions that set a number to NAN. There are may places in the code */
2609 /* where an integer is set to NAN without the use of this function. */
2610 /****************************************************************************/
2611 void siSetToNan(SYNTHETIC_INTEGER **si)
2612 {
2613 asAssert(si != NULL, __LINE__);
2614 asAssert(*si != NULL, __LINE__);
2615
2616 (*si)->digits[0] = '\0';
2617 (*si)->len = 0;
2618 (*si)->nan = TRUE;
2619 (*si)->neg = FALSE;
2620 }
2621
2622
2623 /****************************************************************************/
2624 /* siSetToPowerOfTen(): */
2625 /*--------------------------------------------------------------------------*/
2626 /* DESCRIPTION */
2627 /* Sets a synthetic integer (which must already be created) to an */
2628 /* integral power of ten. This function is helpful for math in many */
2629 /* cases. */
2630 /* */
2631 /* INPUTS */
2632 /* **si : The synthetic integer to operate on. */
2633 /* */
2634 /* exp : The integral exponent value (non-negative). For */
2635 /* example, a value of 0 here will result in the */
2636 /* assignment of 1, 1 here will result in the assignment */
2637 /* of 10, 2 here will result in the assignment of 100. */
2638 /* */
2639 /* ERRORS */
2640 /* A negative exponent will result in an assertion failure. An exponent */
2641 /* too large will result in the assignment of NAN. */
2642 /* */
2643 /* NOTE */
2644 /* This program "evolved", and this function evolved after many of the */
2645 /* functions that set a number to a power of 10 "directly". There may be */
2646 /* many places in the code where an integer is set to a power of 10 */
2647 /* without the use of this function. */
2648 /****************************************************************************/
2649 void siSetToPowerOfTen(SYNTHETIC_INTEGER **si, int exp)
2650 {
2651 int i;
2652
2653 asAssert(si != NULL, __LINE__);
2654 asAssert(*si != NULL, __LINE__);
2655 asAssert(exp >= 0, __LINE__);
2656
2657 if (exp < INTERMEDIATE_CALC_MAX_DIGITS)
2658 {
2659 (*si)->len = exp+1;
2660 (*si)->nan = FALSE;
2661 (*si)->neg = FALSE;
2662 (*si)->digits[exp+1] = '\0';
2663 (*si)->digits[exp] = '1';
2664 for (i=0; i<exp; i++)
2665 (*si)->digits[i] = '0';
2666 }
2667 else
2668 {
2669 siSetToNan(si);
2670 }
2671 }
2672
2673
2674 /****************************************************************************/
2675 /* siAddTwoNonnegative(): */
2676 /*--------------------------------------------------------------------------*/
2677 /* DESCRIPTION */
2678 /* Adds two synthetic integers which are guaranteed to be non-negative. */
2679 /* The result may be set NAN if there is an overflow. This is a low- */
2680 /* level function called by other higher-level functions to do more */
2681 /* complex arithmetic. */
2682 /* */
2683 /* INPUTS */
2684 /* **arg1, **arg2 : The two operands to add. Both operands must be */
2685 /* non-negative. It is legal for both SI's to */
2686 /* use the same storage or to be the same variable. */
2687 /* */
2688 /* OUTPUTS */
2689 /* **result : The synthetic integer where the output result is */
2690 /* placed. This may be the same SI as either or both */
2691 /* of the inputs. */
2692 /****************************************************************************/
2693 void siAddTwoNonnegative(SYNTHETIC_INTEGER **arg1,
2694 SYNTHETIC_INTEGER **arg2,
2695 SYNTHETIC_INTEGER **result)
2696 {
2697 SYNTHETIC_INTEGER *callers_orig_value = NULL;
2698 unsigned i;
2699 unsigned len1, len2, lenmax;
2700 char carry;
2701 char operand1, operand2, digitresult;
2702
2703 /* Be sure the caller isn't doing anything silly. */
2704 asAssert(arg1 != NULL, __LINE__);
2705 asAssert(*arg1 != NULL, __LINE__);
2706 asAssert(!((*arg1)->neg), __LINE__);
2707 asAssert(arg2 != NULL, __LINE__);
2708 asAssert(*arg2 != NULL, __LINE__);
2709 asAssert(!((*arg2)->neg), __LINE__);
2710 asAssert(result != NULL, __LINE__);
2711 asAssert(*result != NULL, __LINE__);
2712 asAssert((arg1 != arg2) && (arg2 != result) && (result != arg1), __LINE__);
2713
2714 /* If either input is NAN, the output is NAN by definition.
2715 */
2716 if (((*arg1)->nan) || ((*arg2)->nan))
2717 {
2718 (*result)->nan = TRUE;
2719 }
2720 else
2721 {
2722 /* The two input arguments don't matter because they aren't written,
2723 ** but if the caller has allowed the result to be the same as the
2724 ** either or both of the two inputs, must use a different memory
2725 ** location. However, only do this if need to, as will slow things
2726 ** down to create and destroy.
2727 */
2728 if ((*arg1 == *result) || (*arg2 == *result))
2729 {
2730 /* Caller has allowed result to be the same as either
2731 ** or both inputs. Must replace caller's pointer temporarily.
2732 */
2733 callers_orig_value = *result;
2734 siCreate(result);
2735 }
2736
2737 /* Zero out the potential result.
2738 */
2739 (*result)->neg = FALSE;
2740 (*result)->nan = FALSE;
2741 (*result)->len = 0;
2742 (*result)->digits[0] = '\0';
2743
2744 /* Buffer important handy values.
2745 */
2746 len1 = (*arg1)->len;
2747 len2 = (*arg2)->len;
2748 lenmax = ddUmax(len1, len2);
2749
2750 /* Be absolutely sure that both lengths are permissible values. Memory
2751 ** reference errors can be hard to track down.
2752 */
2753 asAssert(lenmax <= INTERMEDIATE_CALC_MAX_DIGITS, __LINE__);
2754
2755 /* Start off with an LSD carry input of '0', as there is no carry into the
2756 ** first digit.
2757 */
2758 carry = '0';
2759
2760 /* Iterate through, performing the addition, character by character.
2761 ** characters that are unpopulated are treated as zero.
2762 */
2763 for (i=0; i<lenmax; i++)
2764 {
2765 /* Grab the right character from the arguments. */
2766 if (i < len1)
2767 operand1 = (*arg1)->digits[i];
2768 else
2769 operand1 = '0';
2770 if (i < len2)
2771 operand2 = (*arg2)->digits[i];
2772 else
2773 operand2 = '0';
2774
2775 /* Perform the actual addition. Note that the carry value is copied
2776 ** copied onto the stack for the function call, so carry can double
2777 ** as both input and output.
2778 */
2779 ddFundamentalAdditionCell(operand1, operand2, carry, &digitresult, &carry);
2780
2781 /* Stuff the result digit in the right place. */
2782 (*result)->digits[i] = digitresult;
2783 }
2784
2785 /* After we're out of the loop, there are a few possible outcomes. If the
2786 ** carry is zero, then there is no risk of any kind of overflow.
2787 */
2788 if (carry == '0')
2789 {
2790 /* No risk of overflow. */
2791 (*result)->digits[i] = '\0';
2792 (*result)->len = i;
2793 }
2794 else
2795 {
2796 /* carry == '1' */
2797 /* There are two possibilities here. Either the result has grown one
2798 ** digit more than the longest input operand, and this is legal, or the
2799 ** result has grown one digit more and there is an overflow,
2800 ** forcing declaration of a NAN.
2801 */
2802 if (i == INTERMEDIATE_CALC_MAX_DIGITS)
2803 {
2804 (*result)->nan = TRUE;
2805 (*result)->len = 0;
2806 (*result)->digits[0] = '\0';
2807 }
2808 else
2809 {
2810 (*result)->len = i+1;
2811 (*result)->digits[i] = carry;
2812 (*result)->digits[i+1] = '\0';
2813 }
2814 }
2815
2816 /* If it was necessary to evacuate the caller's result area, swap it back
2817 ** and deallocate memory.
2818 */
2819 if (callers_orig_value)
2820 {
2821 siCopy(result, &callers_orig_value);
2822 siDestroy(result);
2823 *result = callers_orig_value;
2824 }
2825 }
2826 }
2827
2828
2829 /****************************************************************************/
2830 /* siSubtractToProduceNonnegativeResult(): */
2831 /*--------------------------------------------------------------------------*/
2832 /* DESCRIPTION */
2833 /* Subtracts two synthetic integers, with the understanding that the */
2834 /* result will be non-negative. This function is used as a building */
2835 /* block for subtraction which can span zero. */
2836 /* */
2837 /* INPUTS */
2838 /* **arg1, **arg2 : The two operands to subtract. Both must be non- */
2839 /* negative. */
2840 /* */
2841 /* OUTPUTS */
2842 /* **result : **arg1 - **arg2, or a fatal error if the result */
2843 /* would be negative. */
2844 /****************************************************************************/
2845 void siSubtractToProduceNonnegativeResult(SYNTHETIC_INTEGER **arg1,
2846 SYNTHETIC_INTEGER **arg2,
2847 SYNTHETIC_INTEGER **result)
2848 {
2849 SYNTHETIC_INTEGER *callers_orig_value = NULL;
2850 unsigned i;
2851 unsigned len;
2852 char borrow;
2853 char operand1, operand2, digitresult;
2854
2855 /* Be sure the caller isn't doing anything silly. */
2856 asAssert(arg1 != NULL, __LINE__);
2857 asAssert(*arg1 != NULL, __LINE__);
2858 asAssert(!((*arg1)->neg), __LINE__);
2859 asAssert(arg2 != NULL, __LINE__);
2860 asAssert(*arg2 != NULL, __LINE__);
2861 asAssert(!((*arg2)->neg), __LINE__);
2862 asAssert(result != NULL, __LINE__);
2863 asAssert(*result != NULL, __LINE__);
2864 asAssert((arg1 != arg2) && (arg2 != result) && (result != arg1), __LINE__);
2865
2866 /* If either input is NAN, the output is NAN by definition.
2867 */
2868 if (((*arg1)->nan) || ((*arg2)->nan))
2869 {
2870 (*result)->nan = TRUE;
2871 }
2872 else
2873 {
2874 /* The two input arguments don't matter because they aren't written,
2875 ** but if the caller has allowed the result to be the same as the
2876 ** either or both of the two inputs, must use a different memory
2877 ** location. However, only do this if need to, as will slow things
2878 ** down to create and destroy.
2879 */
2880 if ((*arg1 == *result) || (*arg2 == *result))
2881 {
2882 /* Caller has allowed result to be the same as either
2883 ** or both inputs. Must replace caller's pointer temporarily.
2884 */
2885 callers_orig_value = *result;
2886 siCreate(result);
2887 }
2888
2889 /* Zero out the potential result.
2890 */
2891 (*result)->neg = FALSE;
2892 (*result)->nan = FALSE;
2893 (*result)->len = 0;
2894 (*result)->digits[0] = '\0';
2895
2896 /* Buffer important handy values.
2897 */
2898 len = (*arg1)->len; /* By definition, this one is the max. */
2899
2900 /* Be absolutely sure that the length used for iteration is a
2901 ** permissible value.
2902 */
2903 asAssert(len <= INTERMEDIATE_CALC_MAX_DIGITS, __LINE__);
2904
2905 /* Start off with an LSD borrow input of '0', as there is no borrow into the
2906 ** first digit.
2907 */
2908 borrow = '0';
2909
2910 /* Iterate through, performing the subtraction, character by character.
2911 ** Characters that are unpopulated are treated as zero.
2912 */
2913 for (i=0; i<len; i++)
2914 {
2915 /* Grab the right character from the arguments. */
2916 operand1 = (*arg1)->digits[i];
2917 if (i < (*arg2)->len)
2918 operand2 = (*arg2)->digits[i];
2919 else
2920 operand2 = '0';
2921
2922 /* Perform the actual subtraction. Note that the borrow value is
2923 ** copied onto the stack for the function call, so double can double
2924 ** as both input and output.
2925 */
2926 /* printf("Idx: %u, borrow in: %c, ", i, borrow); */
2927
2928 ddFundamentalSubtractionCell(operand1, operand2, borrow, &digitresult, &borrow);
2929
2930 /* printf("op1: %c, op2: %c, borrow out: %c, result out: %c.\n",
2931 operand1,
2932 operand2,
2933 borrow,
2934 digitresult); */
2935
2936 /* Stuff the result digit in the right place. */
2937 (*result)->digits[i] = digitresult;
2938 }
2939
2940 /* Trim the leading zeros from the subtraction result. It would be more
2941 ** elegant to just identify the ripstop during the subtraction, but
2942 ** this will work, too.
2943 */
2944 if (!i)
2945 {
2946 /* We had a 0-0 subtraction. */
2947 (*result)->digits[0] = '\0';
2948 (*result)->len = 0;
2949 }
2950 else
2951 {
2952 (*result)->digits[i] = '\0';
2953 (*result)->len = i;
2954 while ((i >= 1) && ((*result)->digits[i-1] == '0'))
2955 {
2956 (*result)->digits[i-1] = '\0';
2957 (*result)->len = i-1;
2958 i--;
2959 }
2960 }
2961
2962 /* If it was necessary to evacuate the caller's result area, swap it back
2963 ** and deallocate memory.
2964 */
2965 if (callers_orig_value)
2966 {
2967 siCopy(result, &callers_orig_value);
2968 siDestroy(result);
2969 *result = callers_orig_value;
2970 }
2971 }
2972 }
2973
2974
2975 /****************************************************************************/
2976 /* siUnrestrictedAddition(): */
2977 /*--------------------------------------------------------------------------*/
2978 /* DESCRIPTION */
2979 /* Adds two integers of any sign. NAN is produced if anything overflows. */
2980 /* */
2981 /* INPUTS */
2982 /* **arg1, **arg2 : The two synthetic integer operand to add. It is */
2983 /* permitted for both to point to the same place. */
2984 /* */
2985 /* OUTPUTS */
2986 /* **result : The result. This may point to the same place as */
2987 /* either or both of arg1, arg2. */
2988 /****************************************************************************/
2989 void siUnrestrictedAddition(SYNTHETIC_INTEGER **arg1,
2990 SYNTHETIC_INTEGER **arg2,
2991 SYNTHETIC_INTEGER **result)
2992 {
2993 SYNTHETIC_INTEGER *callers_orig_value = NULL;
2994 int arg1_rel_arg2;
2995 unsigned arg1_is_neg;
2996 unsigned arg2_is_neg;
2997 unsigned case_number;
2998
2999 /* Be sure the caller isn't doing anything silly. Because double-
3000 ** indirection is involved, another word of explanation is needed.
3001 ** The caller is passing a pointer to his pointer. This isn't necessary
3002 ** in this circumstance (we'll never need to reallocate his object, for
3003 ** example), but it is just simpler to keep everthing consistent
3004 ** throughout the program. When I say that the caller can have pointers
3005 ** the same for the args and the result, I mean that the caller can have
3006 ** separate physical pointers which point to the same thing. However, the
3007 ** caller may not have one physical pointer only and pass us three identical
3008 ** pointers to the same physical pointer. This is illegal. The caller must
3009 ** have three seperate pointers, even if they point to the same place.
3010 */
3011 asAssert(arg1 != NULL, __LINE__);
3012 asAssert(*arg1 != NULL, __LINE__);
3013 asAssert(arg2 != NULL, __LINE__);
3014 asAssert(*arg2 != NULL, __LINE__);
3015 asAssert(result != NULL, __LINE__);
3016 asAssert(*result != NULL, __LINE__);
3017 asAssert((arg1 != arg2) && (arg2 != result) && (result != arg1), __LINE__);
3018
3019 /* The two input arguments don't matter because they aren't written,
3020 ** but if the caller has allowed the result to be the same as the
3021 ** either or both of the two inputs, must use a different memory
3022 ** location. However, only do this if need to, as will slow things
3023 ** down to create and destroy.
3024 */
3025 if ((*arg1 == *result) || (*arg2 == *result))
3026 {
3027 /* Caller has allowed result to be the same as either
3028 ** or both inputs. Must replace caller's pointer temporarily.
3029 */
3030 callers_orig_value = *result;
3031 siCreate(result);
3032 }
3033
3034 /* Zero out the potential result.
3035 */
3036 (*result)->neg = FALSE;
3037 (*result)->nan = FALSE;
3038 (*result)->len = 0;
3039 (*result)->digits[0] = '\0';
3040
3041 /* If either operand is NAN, the result is NAN. That is all the
3042 ** needs to be said.
3043 */
3044 if (((*arg1)->nan) || ((*arg2)->nan))
3045 {
3046 (*result)->nan = TRUE;
3047 }
3048 else
3049 {
3050 /* Gather some vital statistics about the two numbers passed. These
3051 ** will be used to break into cases.
3052 */
3053 arg1_rel_arg2 = siCompareAbs(arg1, arg2) + 1;
3054 arg1_is_neg = (*arg1)->neg;
3055 arg2_is_neg = (*arg2)->neg;
3056
3057 /* Zero out the signs of the two numbers. They will be restored
3058 ** later.
3059 */
3060 (*arg1)->neg = FALSE;
3061 (*arg2)->neg = FALSE;
3062
3063 /* Create a case number based on the above examination. It is easier
3064 ** to do it this way so there is only one switch statement.
3065 */
3066 case_number = ((unsigned)arg1_rel_arg2 * 4)
3067 + arg1_is_neg * 2
3068 + arg2_is_neg;
3069
3070 /* printf("Case number: %u.\n", case_number); */
3071
3072 switch(case_number)
3073 {
3074 /*********************************************************************/
3075 case 0:
3076 /* abs(arg1) < abs(arg2) */
3077 /* arg1 >= 0 */
3078 /* arg2 >= 0 */
3079 /* Both arguments were non-negative. The result will be non-
3080 ** negative, perhaps NAN. Just add the two numbers, perhaps
3081 ** resulting in NAN.
3082 */
3083 siAddTwoNonnegative(arg1, arg2, result);
3084 break;
3085 /*********************************************************************/
3086 case 1:
3087 /* abs(arg1) < abs(arg2) */
3088 /* arg1 >= 0 */
3089 /* arg2 < 0 */
3090 /* arg2 is negative, and will overwhelm arg1 and produce a negative
3091 ** result. Get the positive result and flip the sign.
3092 */
3093 siSubtractToProduceNonnegativeResult(arg2, arg1, result);
3094 if (!((*result)->nan))
3095 (*result)->neg = TRUE;
3096 break;
3097 /*********************************************************************/
3098 case 2:
3099 /* abs(arg1) < abs(arg2) */
3100 /* arg1 < 0 */
3101 /* arg2 >= 0 */
3102 /* Result will be positive, dominated by arg2. */
3103 siSubtractToProduceNonnegativeResult(arg2, arg1, result);
3104 break;
3105 /*********************************************************************/
3106 case 3:
3107 /* abs(arg1) < abs(arg2) */
3108 /* arg1 < 0 */
3109 /* arg2 < 0 */
3110 /* arg1 and arg2 both negative. Add them in a positive sense and
3111 ** negate the result.
3112 */
3113 siAddTwoNonnegative(arg1, arg2, result);
3114 if (!((*result)->nan))
3115 (*result)->neg = TRUE;
3116 break;
3117 /*********************************************************************/
3118 case 4:
3119 /* abs(arg1) == abs(arg2) */
3120 /* arg1 >= 0 */
3121 /* arg2 >= 0 */
3122 /* Two are equal and non-negative. Just add.
3123 */
3124 siAddTwoNonnegative(arg1, arg2, result);
3125 break;
3126 /*********************************************************************/
3127 case 5:
3128 /* abs(arg1) == abs(arg2) */
3129 /* arg1 >= 0 */
3130 /* arg2 < 0 */
3131 /* Result is zero. Just leave the result the way it is, as we
3132 ** set the result to zero higher up in the code.
3133 */
3134 break;
3135 /*********************************************************************/
3136 case 6:
3137 /* abs(arg1) == abs(arg2) */
3138 /* arg1 < 0 */
3139 /* arg2 >= 0 */
3140 /* Result is zero. Just leave the result the way it is, as we
3141 ** set the result to zero higher up in the code.
3142 */
3143 break;
3144 /*********************************************************************/
3145 case 7:
3146 /* abs(arg1) == abs(arg2) */
3147 /* arg1 < 0 */
3148 /* arg2 < 0 */
3149 /* Two are equal and negative. Add and negate.
3150 */
3151 siAddTwoNonnegative(arg1, arg2, result);
3152 if (!((*result)->nan))
3153 (*result)->neg = TRUE;
3154 break;
3155 /*********************************************************************/
3156 case 8:
3157 /* abs(arg1) > abs(arg2) */
3158 /* arg1 >= 0 */
3159 /* arg2 >= 0 */
3160 /* Both non-negative. Just add.
3161 */
3162 siAddTwoNonnegative(arg1, arg2, result);
3163 break;
3164 /*********************************************************************/
3165 case 9:
3166 /* abs(arg1) > abs(arg2) */
3167 /* arg1 >= 0 */
3168 /* arg2 < 0 */
3169 /* Result is positive and driven by arg1. Subtract. */
3170 siSubtractToProduceNonnegativeResult(arg1, arg2, result);
3171 break;
3172 /*********************************************************************/
3173 case 10:
3174 /* abs(arg1) > abs(arg2) */
3175 /* arg1 < 0 */
3176 /* arg2 >= 0 */
3177 /* Result is negative and driven by arg1. Subtract and negate.
3178 */
3179 siSubtractToProduceNonnegativeResult(arg1, arg2, result);
3180 if (!((*result)->nan))
3181 (*result)->neg = TRUE;
3182 break;
3183 /*********************************************************************/
3184 case 11:
3185 /* abs(arg1) > abs(arg2) */
3186 /* arg1 < 0 */
3187 /* arg2 < 0 */
3188 /* Result is negative. Add and negate.
3189 */
3190 siAddTwoNonnegative(arg1, arg2, result);
3191 if (!((*result)->nan))
3192 (*result)->neg = TRUE;
3193 break;
3194 /*********************************************************************/
3195 default:
3196 /* Should never be here. Fatal out. */
3197 asAssert(0, __LINE__);
3198 break;
3199 /*********************************************************************/
3200 }
3201
3202 /* Restore the results of the input arguments. We will really mess up
3203 ** the caller if we don't.
3204 */
3205 (*arg1)->neg = arg1_is_neg;
3206 (*arg2)->neg = arg2_is_neg;
3207 }
3208
3209
3210 /* If it was necessary to evacuate the caller's result area, swap it back
3211 ** and deallocate memory.
3212 */
3213 if (callers_orig_value)
3214 {
3215 siCopy(result, &callers_orig_value);
3216 siDestroy(result);
3217 *result = callers_orig_value;
3218 }
3219 }
3220
3221 /****************************************************************************/
3222 /* siUnrestrictedSubtraction(): */
3223 /*--------------------------------------------------------------------------*/
3224 /* DESCRIPTION */
3225 /* Subtracts two integers of any sign. NAN is produced if anything */
3226 /* overflows. */
3227 /* */
3228 /* INPUTS */
3229 /* **arg1, **arg2 : The two synthetic integer operand to subtract. */
3230 /* Both may point to the same place. */
3231 /* */
3232 /* OUTPUTS */
3233 /* **result : The result, arg1-arg2. May point to the same */
3234 /* place as either or both of arg1, arg2. */
3235 /****************************************************************************/
3236 void siUnrestrictedSubtraction(SYNTHETIC_INTEGER **arg1,
3237 SYNTHETIC_INTEGER **arg2,
3238 SYNTHETIC_INTEGER **result)
3239 {
3240 unsigned neg;
3241
3242 /* Be sure the caller is doing nothing silly. */
3243 asAssert(arg1 != NULL, __LINE__);
3244 asAssert(*arg1 != NULL, __LINE__);
3245 asAssert(arg2 != NULL, __LINE__);
3246 asAssert(*arg2 != NULL, __LINE__);
3247 asAssert(result != NULL, __LINE__);
3248 asAssert(*result != NULL, __LINE__);
3249 asAssert((arg1 != arg2) && (arg2 != result) && (result != arg1), __LINE__);
3250
3251 /* If either operand is NAN, the result is NAN. */
3252 if ((*arg1)->nan || (*arg2)->nan)
3253 {
3254 (*result)->digits[0] = '\0';
3255 (*result)->len = 0;
3256 (*result)->nan = TRUE;
3257 (*result)->neg = FALSE;
3258 }
3259 else
3260 {
3261 /* Subtraction is easy to express in terms of addition, since the
3262 ** addition function already handles both operands of both signs.
3263 */
3264 if ((*arg2)->len) /* Can't negate zero in this way. */
3265 {
3266 neg = (*arg2)->neg;
3267 (*arg2)->neg = !((*arg2)->neg);
3268 }
3269
3270 siUnrestrictedAddition(arg1, arg2, result);
3271
3272 if ((*arg2)->len)
3273 (*arg2)->neg = neg;
3274 }
3275 }
3276
3277
3278 /****************************************************************************/
3279 /* siUnrestrictedMultiplication(): */
3280 /*--------------------------------------------------------------------------*/
3281 /* DESCRIPTION */
3282 /* Muiltiplies two integers of any sign. NAN is produced if anything */
3283 /* overflows. */
3284 /* */
3285 /* INPUTS */
3286 /* **arg1, **arg2 : The two synthetic integer operands to multiply. */
3287 /* Both may point to the same place. */
3288 /* */
3289 /* OUTPUTS */
3290 /* **result : The result, arg1*arg2. May point to the same */
3291 /* place as either or both of arg1, arg2. */
3292 /****************************************************************************/
3293 void siUnrestrictedMultiplication(SYNTHETIC_INTEGER **arg1,
3294 SYNTHETIC_INTEGER **arg2,
3295 SYNTHETIC_INTEGER **result)
3296 {
3297 /* Be sure the caller is doing nothing silly. */
3298 asAssert(arg1 != NULL, __LINE__);
3299 asAssert(*arg1 != NULL, __LINE__);
3300 asAssert(arg2 != NULL, __LINE__);
3301 asAssert(*arg2 != NULL, __LINE__);
3302 asAssert(result != NULL, __LINE__);
3303 asAssert(*result != NULL, __LINE__);
3304 asAssert((arg1 != arg2) && (arg2 != result) && (result != arg1), __LINE__);
3305
3306 /* If either argument is NAN, the result is NAN.
3307 */
3308 if ((*arg1)->nan || (*arg2)->nan)
3309 {
3310 (*result)->digits[0] = '\0';
3311 (*result)->len = 0;
3312 (*result)->nan = TRUE;
3313 (*result)->neg = FALSE;
3314 }
3315 /* If either argument is zero, result is zero. */
3316 else if (!((*arg1)->len) || !((*arg2)->len))
3317 {
3318 (*result)->digits[0] = '\0';
3319 (*result)->len = 0;
3320 (*result)->nan = FALSE;
3321 (*result)->neg = FALSE;
3322 }
3323 else
3324 {
3325 /* We have two non-zero operands, fit to be multiplied.
3326 */
3327 SYNTHETIC_INTEGER *rt, *temp, *rt2;
3328 int cidx, len;
3329 char cur_digit;
3330 unsigned negarg1, negarg2;
3331
3332 /* Must save the signs of the arguments and restore. */
3333 negarg1 = (*arg1)->neg;
3334 negarg2 = (*arg2)->neg;
3335 (*arg1)->neg = FALSE;
3336 (*arg2)->neg = FALSE;
3337
3338 /* Create (allocate) a synthetic integer to contain the result
3339 ** that will be returned to the caller.
3340 */
3341 siCreate(&rt);
3342
3343 /* Must create an alias. The reasons for this are complex. I coded
3344 ** myself into a corner and have to live with it.
3345 */
3346 rt2 = rt;
3347
3348 /* Create (allocate) a synthetic integer that will be used to hold
3349 ** the result of multiplying arg2 by a digit.
3350 */
3351 siCreate(&temp);
3352
3353 /* Go through, multiplying in much the same way as one would by hand.
3354 */
3355 len = (*arg1)->len;
3356
3357 for (cidx = len-1; cidx >= 0; cidx--)
3358 {
3359 #if 0
3360 printf("cidx: %d\n", cidx);
3361 printf("rt before * 10 multiply\n");
3362 siDump(&rt, "rt before");
3363 /* Multiply the pending result by 10 */
3364 #endif
3365 siMulByTen(&rt);
3366 /* siDump(&rt, "rt after"); */
3367
3368 /* Grab the current digit of interest and multiply arg2 by it. */
3369 cur_digit = (*arg1)->digits[cidx];
3370 /* printf("cur digit: %c\n", cur_digit); */
3371
3372 siCopy(arg2, &temp);
3373 /* siDump(&temp, "copy of arg2"); */
3374 siMulByDigit(&temp, cur_digit);
3375 /* siDump(&temp, "arg2 * digit"); */
3376
3377 /* Add the result of the multiplication to the pending
3378 ** result.
3379 */
3380 siUnrestrictedAddition(&rt2, &temp, &rt);
3381
3382 #if 0
3383 siDump(&rt2, "rt2");
3384 siDump(&rt, "rt");
3385 siDump(&temp, "temp");
3386 #endif
3387 }
3388
3389 /* We're done with the mathematical operation. Signs should be a self-
3390 ** healing matter.
3391 */
3392 /* Copy the result back to caller's area. */
3393 siCopy(&rt, result);
3394
3395 /* Destroy temporary areas. */
3396 siDestroy(&rt);
3397 siDestroy(&temp);
3398
3399 /* Restore the signs of the integers that were arguments, then the
3400 ** sign of the result.
3401 */
3402 (*arg1)->neg = negarg1;
3403 (*arg2)->neg = negarg2;
3404
3405 if ((negarg1 && negarg2) || (!negarg1 && !negarg2))
3406 (*result)->neg = FALSE;
3407 else
3408 (*result)->neg = TRUE;
3409
3410 /* siDump(result, "result"); */
3411 }
3412 }
3413
3414
3415 /****************************************************************************/
3416 /* siUnrestrictedDivision(): */
3417 /*--------------------------------------------------------------------------*/
3418 /* DESCRIPTION */
3419 /* Divides two synthetic integers of any sign, producing an integer */
3420 /* quotient and integer remainder. NAN is assigned to both the quotient */
3421 /* and remainder on invalid inputs. NAN due to overflow is not possible */
3422 /* in this function, as any integer divided by 1 or -1 yields an integer */
3423 /* of the same magnitude. */
3424 /* */
3425 /* INPUTS */
3426 /* **dividend : The dividend of the division. May be any value. */
3427 /* */
3428 /* **divisor : The divisor of the division. May be any value */
3429 /* except zero. */
3430 /* */
3431 /* OUTPUTS */
3432 /* **quotient : The integer quotient of the division. This should */
3433 /* behave identically to the C "/" operator when */
3434 /* applied to integers. */
3435 /* */
3436 /* **remainder : The remainder of the division. This should behave */
3437 /* identically to the C "%" operator. There was some */
3438 /* question in my mind about the best way to define a */
3439 /* remainder with a negative dividend or divisor, but */
3440 /* I simply used the behavior of the "%" operator in */
3441 /* C as a guide. */
3442 /* */
3443 /* USAGE NOTES */
3444 /* The functions in this program are inconsistent about whether they */
3445 /* allow the inputs and the outputs to be pointers to the same memory */
3446 /* locations, etc. This program was constructed impromptu, and I made */
3447 /* up the rules as I went along. I've now decided that it isn't that */
3448 /* useful to allow this, so for this function I'll require uniqueness */
3449 /* between the four parameters. This doesn't affect the correctness */
3450 /* of the program--just its internal elegance. These inconsistent */
3451 /* decisions are not visible by using the program. */
3452 /****************************************************************************/
3453 void siUnrestrictedDivision(SYNTHETIC_INTEGER **dividend,
3454 SYNTHETIC_INTEGER **divisor,
3455 SYNTHETIC_INTEGER **quotient,
3456 SYNTHETIC_INTEGER **remainder)
3457 {
3458 /* Be sure the caller is doing nothing silly. */
3459 /* No first- or second-level pointers may be NULL. */
3460 asAssert(dividend != NULL, __LINE__);
3461 asAssert(*dividend != NULL, __LINE__);
3462 asAssert(divisor != NULL, __LINE__);
3463 asAssert(*divisor != NULL, __LINE__);
3464 asAssert(quotient != NULL, __LINE__);
3465 asAssert(*quotient != NULL, __LINE__);
3466 asAssert(remainder != NULL, __LINE__);
3467 asAssert(*remainder != NULL, __LINE__);
3468 /* All pointers must be mutually exclusive. */
3469 asAssert((dividend != divisor) && (dividend != quotient) && (dividend != remainder), __LINE__);
3470 asAssert((*dividend != *divisor) && (*dividend != *quotient) && (*dividend != *remainder), __LINE__);
3471 asAssert((divisor != quotient) && (divisor != remainder), __LINE__);
3472 asAssert((*divisor != *quotient) && (*divisor != *remainder), __LINE__);
3473 asAssert(quotient != remainder, __LINE__);
3474 asAssert(*quotient != *remainder, __LINE__);
3475
3476 /* If either argument is NAN, the results are NAN. Also, if the divisor is zero,
3477 ** that causes two NAN's, as well.
3478 */
3479 if ((*dividend)->nan || (*divisor)->nan || !((*divisor)->len))
3480 {
3481 (*quotient)->digits[0] = '\0';
3482 (*quotient)->len = 0;
3483 (*quotient)->nan = TRUE;
3484 (*quotient)->neg = FALSE;
3485 (*remainder)->digits[0] = '\0';
3486 (*remainder)->len = 0;
3487 (*remainder)->nan = TRUE;
3488 (*remainder)->neg = FALSE;
3489 }
3490 /* Zero in the numerator always gets us zero on both outputs.
3491 */
3492 else if (!((*dividend)->len))
3493 {
3494 siSetToLong(quotient, 0);
3495 siSetToLong(remainder, 0);
3496 }
3497 /* Must do a legitimate long division. We know for sure that
3498 ** both dividend and divisor are non-zero.
3499 */
3500 else
3501 {
3502 unsigned dividend_is_neg;
3503 unsigned divisor_is_neg;
3504 /* Signs of dividend and divisor, saved so we only
3505 ** operate on positive numbers.
3506 */
3507 int quotient_ndigits;
3508 /* The number of digits in the quotient.
3509 */
3510 int i;
3511 /* General-purpose iteration variable.
3512 */
3513 char j;
3514 /* Trial digit for trial multiplication.
3515 */
3516 SYNTHETIC_INTEGER *si_temp1 = NULL,
3517 *si_temp2 = NULL,
3518 *si_temp3 = NULL;
3519 /* Temporary large integers. These need to be set
3520 ** to NULL here, because in the event of an
3521 ** error or unexpected condition, need to know
3522 ** which ones to release.
3523 */
3524
3525 /* Remember and reset the signs of the passed dividend and
3526 ** divisor. It makes it simpler just to consider positive
3527 ** numbers.
3528 */
3529 dividend_is_neg = (*dividend)->neg;
3530 divisor_is_neg = (*divisor)->neg;
3531 (*dividend)->neg = 0;
3532 (*divisor)->neg = 0;
3533
3534 /* If the dividend is less (in magnitude) than the divisor,
3535 ** the quotient is zero and the remainder is the
3536 ** dividend.
3537 */
3538 if (siCompare(dividend, divisor) == -1)
3539 {
3540 siSetToLong(quotient, 0);
3541 siCopy(dividend, remainder);
3542 (*remainder)->neg = dividend_is_neg;
3543 }
3544 else
3545 {
3546 /* Allocate all of the temporary integers.
3547 ** Done all in one place to eliminate paths
3548 ** through the code where may allocate something
3549 ** already allocated and cause a memory leak.
3550 */
3551 siCreate(&si_temp1);
3552 siCreate(&si_temp2);
3553 siCreate(&si_temp3);
3554
3555 /* Establish the number of digits in the
3556 ** quotient. There is a window near the maximum
3557 ** number of digits where we won't be able to
3558 ** successfully start the division because we
3559 ** can't form an upper bound which is a power
3560 ** of 10 times the divisor. In this case,
3561 ** just return NAN. This window won't be a problem
3562 ** in practice. It could be worked around, but
3563 ** no need to close it.
3564 */
3565 siCopy(divisor, &si_temp1);
3566 siMulByTen(&si_temp1);
3567 quotient_ndigits = 1;
3568 while(siCompare(&si_temp1, dividend) <= 0)
3569 {
3570 siMulByTen(&si_temp1);
3571 if (si_temp1->nan)
3572 goto nan_return;
3573 quotient_ndigits++;
3574 }
3575
3576 siDivByTen(&si_temp1);
3577
3578 /* printf("quotient number of digits: %d\n",
3579 quotient_ndigits);
3580 siDump(&si_temp1, "trial divisor"); */
3581
3582 /* We now know the number of digits that the
3583 ** quotient will contain. We also know the
3584 ** starting point to use for trial subtraction.
3585 ** Count down from the highest order digit to the
3586 ** lowest-order digit, subtracting the maximum
3587 ** amount each time.
3588 */
3589 /* si_temp1 is the trial divisor reference
3590 ** copy.
3591 */
3592 /* si_temp2 is a copy of the trial divisor,
3593 ** used to see what happens as we try multiplying
3594 ** by various trial digits.
3595 */
3596 /* remainder (in the caller's area) is used for
3597 ** staging the dividend after subtractions made.
3598 ** This will end up with the remainder in the
3599 ** end.
3600 */
3601 /* quotient (in the caller's area) contains the
3602 ** quotient. We'll stage it directly there.
3603 */
3604 (*quotient)->len = (unsigned)quotient_ndigits;
3605 (*quotient)->digits[(*quotient)->len] = '\0';
3606 siCopy(dividend, remainder);
3607
3608 /* For each digit in the result, going from most
3609 ** significant to least significant.
3610 */
3611 for (i=quotient_ndigits-1; i>=0; i--)
3612 {
3613 /* Try to find the largest digit for the spot
3614 ** that will not overflow the remainder.
3615 ** This is the right digit.
3616 */
3617 /* For each digit we might try.
3618 */
3619 for (j='0'; j<='9'; j++)
3620 {
3621 siCopy(&si_temp1, &si_temp2);
3622 /* Buffer the value so we can do a trial
3623 ** multiplication.
3624 */
3625 siMulByDigit(&si_temp2, j);
3626 /* Perform the trial multiplication
3627 ** by a digit.
3628 */
3629 siUnrestrictedSubtraction(remainder,
3630 &si_temp2,
3631 &si_temp3);
3632 /* Perform a trial subtraction involving the trial
3633 ** digit.
3634 */
3635
3636 /* We have success if what is left over after
3637 ** the trial subtraction is less than the
3638 ** trial divisor.
3639 */
3640 if (siCompare(&si_temp3, &si_temp1) < 0)
3641 break;
3642
3643 /* It is a fatal internal software error if
3644 ** we are at the digit '9' and the constraint
3645 ** still is not met.
3646 */
3647 asAssert(j != '9', __LINE__);
3648 } /* End for each possible digit. */
3649
3650 /* OK, we have the digit, and it is in the
3651 ** variable 'j'. We also have the value
3652 ** we should subtract and the result.
3653 ** Stuff the digit, and the result.
3654 */
3655 (*quotient)->digits[i] = j;
3656 siCopy(&si_temp3, remainder);
3657
3658 /* Divide the trial divisor by 10 for the
3659 ** next iteration.
3660 */
3661 siDivByTen(&si_temp1);
3662 } /* End for each digit position. */
3663
3664 /* We now need to figure out what the sign
3665 ** of the quotient and the remainder are.
3666 */
3667 if ((dividend_is_neg && !divisor_is_neg) ||
3668 (!dividend_is_neg && divisor_is_neg))
3669 (*quotient)->neg = TRUE;
3670 else
3671 (*quotient)->neg = FALSE;
3672
3673 if (dividend_is_neg && (*remainder)->len)
3674 (*remainder)->neg = TRUE;
3675 else
3676 (*remainder)->neg = FALSE;
3677 }
3678
3679 goto normal_return;
3680 nan_return: ;
3681 /* This is the return point where the quotient
3682 ** and remainder are marked NAN. Thank goodness
3683 ** they didn't forget about the "goto" in 'C'.
3684 */
3685
3686 normal_return:
3687 /* Restore the signs of the input arguments.
3688 */
3689 (*dividend)->neg = dividend_is_neg;
3690 (*divisor)->neg = divisor_is_neg;
3691
3692 /* Deallocate all allocated memory.
3693 */
3694 if (si_temp1)
3695 siDestroy(&si_temp1);
3696 if (si_temp2)
3697 siDestroy(&si_temp2);
3698 if (si_temp3)
3699 siDestroy(&si_temp3);
3700 }
3701 }
3702
3703
3704 /****************************************************************************/
3705 /* siGcd(): */
3706 /*--------------------------------------------------------------------------*/
3707 /* DESCRIPTION */
3708 /* Applies Euclid's algorithm to obtain the GCD of two positive integers. */
3709 /* Any unexpected input operands will generate a NAN for the GCD. */
3710 /* */
3711 /* INPUTS */
3712 /* **arg1, arg2 : The two positive integers whose GCD is to be found. */
3713 /* */
3714 /* OUTPUTS */
3715 /* **result : gcd(**arg1, **arg2). */
3716 /****************************************************************************/
3717 void siGcd(SYNTHETIC_INTEGER **arg1,
3718 SYNTHETIC_INTEGER **arg2,
3719 SYNTHETIC_INTEGER **result)
3720 {
3721
3722 /* Be sure the caller isn't doing anything silly. */
3723 asAssert(arg1 != NULL, __LINE__);
3724 asAssert(*arg1 != NULL, __LINE__);
3725 asAssert(arg2 != NULL, __LINE__);
3726 asAssert(*arg2 != NULL, __LINE__);
3727 asAssert(result != NULL, __LINE__);
3728 asAssert(*result != NULL, __LINE__);
3729
3730 /* A NAN input or zero or negative integers gets a NAN result.
3731 */
3732 if ((*arg1)->nan || (*arg2)->nan || !(*arg1)->len || !(*arg2)->len
3733 || (*arg1)->neg || (*arg2)->neg)
3734 {
3735 (*result)->digits[0] = '\0';
3736 (*result)->len = 0;
3737 (*result)->nan = TRUE;
3738 (*result)->neg = FALSE;
3739 }
3740 else
3741 {
3742 /* This algorithm is so well known and documented that I won't insult
3743 ** the reader with comments.
3744 */
3745 SYNTHETIC_INTEGER *h, *k, *q, *r;
3746
3747 siCreate(&h);
3748 siCreate(&k);
3749 siCreate(&q);
3750 siCreate(&r);
3751
3752 if (siCompare(arg1, arg2) == -1)
3753 {
3754 siCopy(arg1, &h);
3755 siCopy(arg2, &k);
3756 }
3757 else
3758 {
3759 siCopy(arg1, &k);
3760 siCopy(arg2, &h);
3761 }
3762
3763 siCopy(&h, &r);
3764
3765 do
3766 {
3767 siCopy(&r, result);
3768 siUnrestrictedDivision(&h,&k,&q,&r);
3769 siCopy(&k, &h);
3770 siCopy(&r, &k);
3771 } while (r->len);
3772
3773 siDestroy(&h);
3774 siDestroy(&k);
3775 siDestroy(&q);
3776 siDestroy(&r);
3777 }
3778 }
3779
3780
3781 /****************************************************************************/
3782 /* siIntegerExponentiation(): */
3783 /*--------------------------------------------------------------------------*/
3784 /* DESCRIPTION */
3785 /* Exponentiates one integer to the integral power of another. The */
3786 /* integer to be exponentiated may be of any sign, but the exponent must */
3787 /* be non-negative. Invalid or overflow conditions result in the */
3788 /* assignment of NAN. */
3789 /* */
3790 /* INPUTS */
3791 /* **arg : Integer to be exponentiated. */
3792 /* */
3793 /* **exponent : Exponent to use. Must be non-negative. */
3794 /* */
3795 /* OUTPUTS */
3796 /* **result : The result, arg1**arg2. */
3797 /* */
3798 /* NOTE */
3799 /* In other functions (written chronologically earlier), it was allowed */
3800 /* for the arguments and the result to be the same physical object. This */
3801 /* didn't prove to be a useful feature, so it is disallowed in this */
3802 /* function. */
3803 /****************************************************************************/
3804 void siIntegerExponentiation(SYNTHETIC_INTEGER **arg,
3805 SYNTHETIC_INTEGER **exponent,
3806 SYNTHETIC_INTEGER **result)
3807 {
3808 /* Be sure the caller is doing nothing silly. */
3809 asAssert(arg != NULL, __LINE__);
3810 asAssert(*arg != NULL, __LINE__);
3811 asAssert(exponent != NULL, __LINE__);
3812 asAssert(*exponent != NULL, __LINE__);
3813 asAssert(result != NULL, __LINE__);
3814 asAssert(*result != NULL, __LINE__);
3815 asAssert((arg != exponent) && (exponent != result) && (result != arg), __LINE__);
3816 asAssert((*arg != *exponent) && (*exponent != *result) && (*result != *arg), __LINE__);
3817
3818 /* If either argument is NAN, the result is NAN. */
3819 if (((*arg)->nan) || ((*exponent)->nan))
3820 {
3821 (*result)->digits[0] = '\0';
3822 (*result)->len = 0;
3823 (*result)->nan = TRUE;
3824 (*result)->neg = FALSE;
3825 }
3826 /* If we're trying to raise 0 to the 0th power, this is a no-no,
3827 ** as well.
3828 */
3829 else if (!((*arg)->len) && !((*exponent)->len))
3830 {
3831 (*result)->digits[0] = '\0';
3832 (*result)->len = 0;
3833 (*result)->nan = TRUE;
3834 (*result)->neg = FALSE;
3835 }
3836 /* A negative exponent is a no-no, as well.
3837 */
3838 else if ((*exponent)->neg)
3839 {
3840 (*result)->digits[0] = '\0';
3841 (*result)->len = 0;
3842 (*result)->nan = TRUE;
3843 (*result)->neg = FALSE;
3844 }
3845 /* Zero to any power but zero is zero.
3846 */
3847 else if (!((*arg)->len))
3848 {
3849 siSetToLong(result, 0);
3850 }
3851 /* Any number but zero to the zeroth power is one. */
3852 else if (!((*exponent)->len))
3853 {
3854 siSetToLong(result, 1);
3855 }
3856 else
3857 {
3858 unsigned binary_exponent;
3859 SYNTHETIC_INTEGER *temp_si1, *result_copy1, *temp_si2, *temp_si3;
3860
3861 /* Here we have a non-zero integer raised to a non-zero power.
3862 ** The best algorithm to do this is to examine the bit pattern
3863 ** corresponding to the exponent--we can generate arg**2,
3864 ** arg**4, arg**8, etc., just by repeatedly squaring it. NAN's
3865 ** are picked up automatically.
3866 */
3867 /* First, it is convenient convert the exponent to an
3868 ** unsigned integer, so we can examine the
3869 ** bit pattern directly. Will agree that anything over
3870 ** 32000 is treated as 32000, to avoid platform issues.
3871 */
3872 siCreate(&temp_si1);
3873 siCreate(&temp_si2);
3874 siCreate(&temp_si3);
3875 siCreate(&result_copy1);
3876
3877 siSetToLong(&temp_si1, 32000);
3878
3879 if (siCompare(exponent, &temp_si1) > 0)
3880 {
3881 asFatal("exponent too large (max is 32,000)");
3882 }
3883 else
3884 {
3885 /* Must convert from the string representation
3886 ** buried in the synthetic integer to a binary
3887 ** unsigned integer. This is a bit tricky because the
3888 ** string is in reverse order from the traditional.
3889 */
3890 {
3891 int sidx;
3892
3893 sidx = (*exponent)->len - 1;
3894 binary_exponent = 0;
3895
3896 while (sidx >= 0)
3897 {
3898 binary_exponent *= 10;
3899 binary_exponent += ddDigitToValue((*exponent)->digits[sidx]);
3900 sidx--;
3901 }
3902 }
3903 }
3904
3905 /* Start off with a result of 1 */
3906 siSetToLong(result, 1);
3907 siSetToLong(&result_copy1, 1);
3908
3909 /* Start off with arg**1 */
3910 siCopy(arg, &temp_si1);
3911 siCopy(arg, &temp_si2);
3912 siCopy(arg, &temp_si3);
3913
3914 while (binary_exponent)
3915 {
3916 /* Multiply result by the working register iff we have a 1 as the
3917 ** LSB of the exponent.
3918 */
3919 if (binary_exponent & 1)
3920 {
3921 siUnrestrictedMultiplication(
3922 &result_copy1,
3923 &temp_si1,
3924 result);
3925
3926 /* Clone out the result, because we're keeping an additional copy. */
3927 siCopy(result, &result_copy1);
3928 }
3929
3930
3931 /* Square the quantity cloned from the argument, and clone the copies
3932 ** back.
3933 */
3934 siUnrestrictedMultiplication(
3935 &temp_si1,
3936 &temp_si2,
3937 &temp_si3);
3938 siCopy(&temp_si3, &temp_si1);
3939 siCopy(&temp_si3, &temp_si2);
3940
3941 /* Shift the exponent down one bit. */
3942 binary_exponent >>= 1;
3943 }
3944
3945 /* Destroy all of the data structures created on the fly. */
3946 siDestroy(&temp_si1);
3947 siDestroy(&temp_si2);
3948 siDestroy(&temp_si3);
3949 siDestroy(&result_copy1);
3950
3951 /* At this point the result should contain the right value and
3952 ** we should be ready to return.
3953 */
3954 }
3955 }
3956
3957
3958 /****************************************************************************/
3959 /****************************************************************************/
3960 /********** R A T I O N A L N U M B E R F U N C T I O N S *********/
3961 /****************************************************************************/
3962 /****************************************************************************/
3963 /* This section is reserved for functions which operate on rational numbers.
3964 ** For directness, rational numbers are represented as integer pairs (no
3965 ** sense adding more data structures to this program.
3966 */
3967
3968 /****************************************************************************/
3969 /* rnSum(): */
3970 /*--------------------------------------------------------------------------*/
3971 /* DESCRIPTION */
3972 /* Calculates the sum of two rational numbers, supplying the results in */
3973 /* lowest terms. */
3974 /* */
3975 /* INPUTS */
3976 /* **h1_formal_par_in, */
3977 /* **k1_formal_par_in, */
3978 /* **h2_formal_par_in, */
3979 /* **k2_formal_par_in : The numerators and denominators of the two */
3980 /* rational numbers to add. Zero denominators */
3981 /* are not allowed, but duplicate pointers of all */
3982 /* types are. */
3983 /* **h_result, */
3984 /* **k_result : The result, in lowest terms. Zero is repre- */
3985 /* sented canonically as 0/1. The memory must be */
3986 /* allocated in the caller's area, and the two */
3987 /* data items must be distinct. */
3988 /****************************************************************************/
3989 void rnSum(SYNTHETIC_INTEGER **h1_formal_par_in,
3990 SYNTHETIC_INTEGER **k1_formal_par_in,
3991 SYNTHETIC_INTEGER **h2_formal_par_in,
3992 SYNTHETIC_INTEGER **k2_formal_par_in,
3993 SYNTHETIC_INTEGER **h_result,
3994 SYNTHETIC_INTEGER **k_result)
3995 {
3996 SYNTHETIC_INTEGER *arg1_h,
3997 *arg1_k,
3998 *arg2_h,
3999 *arg2_k,
4000 *common_denominator,
4001 *arg1_numerator,
4002 *arg2_numerator,
4003 *cumulative_numerator,
4004 *result_gcd,
4005 *trash_remainder;
4006
4007 unsigned is_neg = FALSE;
4008
4009 /* Be sure the caller isn't doing anything silly with pointers.
4010 */
4011 asAssert(h1_formal_par_in != NULL, __LINE__);
4012 asAssert(*h1_formal_par_in != NULL, __LINE__);
4013 asAssert(k1_formal_par_in != NULL, __LINE__);
4014 asAssert(*k1_formal_par_in != NULL, __LINE__);
4015 asAssert(h2_formal_par_in != NULL, __LINE__);
4016 asAssert(*h2_formal_par_in != NULL, __LINE__);
4017 asAssert(k2_formal_par_in != NULL, __LINE__);
4018 asAssert(*k2_formal_par_in != NULL, __LINE__);
4019 asAssert(h_result != NULL, __LINE__);
4020 asAssert(*h_result != NULL, __LINE__);
4021 asAssert(k_result != NULL, __LINE__);
4022 asAssert(*k_result != NULL, __LINE__);
4023 asAssert(h_result != k_result, __LINE__);
4024 asAssert(*h_result != *k_result, __LINE__);
4025
4026 /* Allocate space for temporary variables.
4027 */
4028 siCreate(&arg1_h);
4029 siCreate(&arg1_k);
4030 siCreate(&arg2_h);
4031 siCreate(&arg2_k);
4032 siCreate(&common_denominator);
4033 siCreate(&arg1_numerator);
4034 siCreate(&arg2_numerator);
4035 siCreate(&cumulative_numerator);
4036 siCreate(&result_gcd);
4037 siCreate(&trash_remainder);
4038
4039 /* Copy over the formal parameters to local variables. This
4040 ** means we can modify the inputs without worry.
4041 */
4042 siCopy(h1_formal_par_in, &arg1_h);
4043 siCopy(k1_formal_par_in, &arg1_k);
4044 siCopy(h2_formal_par_in, &arg2_h);
4045 siCopy(k2_formal_par_in, &arg2_k);
4046
4047 /* If any of the four inputs are NAN, the results must be
4048 ** NAN.
4049 */
4050 if (arg1_h->nan || arg1_k->nan || arg2_h->nan || arg2_k->nan)
4051 {
4052 goto nan_return;
4053 }
4054
4055 /* If either denominator is zero, the result is NAN.
4056 */
4057 if (!(arg1_k->len) || !(arg2_k->len))
4058 {
4059 goto nan_return;
4060 }
4061
4062 /* Adjust the signs of the the arguments so that the denominators are
4063 ** always positive.
4064 */
4065 if (arg1_h->neg && arg1_k->neg)
4066 {
4067 arg1_h->neg = FALSE;
4068 arg1_k->neg = FALSE;
4069 }
4070 else if (!(arg1_h->neg) && arg1_k->neg)
4071 {
4072 arg1_h->neg = TRUE;
4073 arg1_k->neg = FALSE;
4074 }
4075
4076 if (arg2_h->neg && arg2_k->neg)
4077 {
4078 arg2_h->neg = FALSE;
4079 arg2_k->neg = FALSE;
4080 }
4081 else if (!(arg2_h->neg) && arg2_k->neg)
4082 {
4083 arg2_h->neg = TRUE;
4084 arg2_k->neg = FALSE;
4085 }
4086
4087 /* Calculate the common denominator.
4088 */
4089 siUnrestrictedMultiplication(&arg1_k, &arg2_k, &common_denominator);
4090
4091 /* Calculate the arg1 numerator.
4092 */
4093 siUnrestrictedMultiplication(&arg1_h, &arg2_k, &arg1_numerator);
4094
4095 /* Calculate the arg2 numerator.
4096 */
4097 siUnrestrictedMultiplication(&arg2_h, &arg1_k, &arg2_numerator);
4098
4099 /* Add the numerators to get the additive result. */
4100 siUnrestrictedAddition(&arg1_numerator,
4101 &arg2_numerator,
4102 &cumulative_numerator);
4103
4104 /* If we NAN'd out ...
4105 */
4106 if (cumulative_numerator->nan)
4107 goto nan_return;
4108
4109 /* If we've reached a result of zero, it won't be legal to form
4110 ** a gcd(). Assign the canonical form of zero.
4111 */
4112 if (!(cumulative_numerator->len))
4113 {
4114 siSetToLong(h_result, 0);
4115 siSetToLong(k_result, 1);
4116 goto normal_return;
4117 }
4118
4119 /* Snatch the sign temporarily.
4120 */
4121 if (cumulative_numerator->neg)
4122 {
4123 is_neg = TRUE;
4124 cumulative_numerator->neg = FALSE;
4125 }
4126
4127 /* Form the gcd. */
4128 siGcd(&cumulative_numerator,
4129 &common_denominator,
4130 &result_gcd);
4131
4132 /* Set the sign back. */
4133 cumulative_numerator->neg = is_neg;
4134
4135 /* Divide both numerator and denominator by
4136 ** the gcd().
4137 */
4138 siUnrestrictedDivision(&cumulative_numerator,
4139 &result_gcd,
4140 h_result,
4141 &trash_remainder);
4142 siUnrestrictedDivision(&common_denominator,
4143 &result_gcd,
4144 k_result,
4145 &trash_remainder);
4146
4147 goto normal_return;
4148 nan_return: ;
4149 siSetToNan(h_result);
4150 siSetToNan(k_result);
4151 normal_return: ;
4152
4153 /* Destroy temporary variables.
4154 */
4155 siDestroy(&arg1_h);
4156 siDestroy(&arg1_k);
4157 siDestroy(&arg2_h);
4158 siDestroy(&arg2_k);
4159 siDestroy(&common_denominator);
4160 siDestroy(&arg1_numerator);
4161 siDestroy(&arg2_numerator);
4162 siDestroy(&cumulative_numerator);
4163 siDestroy(&result_gcd);
4164 siDestroy(&trash_remainder);
4165 }
4166
4167
4168 /****************************************************************************/
4169 /* rnDifference(): */
4170 /*--------------------------------------------------------------------------*/
4171 /* DESCRIPTION */
4172 /* Calculates the difference of two rational numbers, returning the */
4173 /* results in lowest terms. */
4174 /* */
4175 /* INPUTS */
4176 /* **h1_formal_par_in, */
4177 /* **k1_formal_par_in, */
4178 /* **h2_formal_par_in, */
4179 /* **k2_formal_par_in : The numerators and denominators of the two */
4180 /* rational numbers to subtract. Zero denominators */
4181 /* are not allowed, but duplicate pointers of all */
4182 /* types are. */
4183 /* **h_result, */
4184 /* **k_result : The result, in lowest terms. Zero is repre- */
4185 /* sented canonically as 0/1. The memory must be */
4186 /* allocated in the caller's area, and the two */
4187 /* data items must be distinct. */
4188 /****************************************************************************/
4189 void rnDifference(SYNTHETIC_INTEGER **h1_formal_par_in,
4190 SYNTHETIC_INTEGER **k1_formal_par_in,
4191 SYNTHETIC_INTEGER **h2_formal_par_in,
4192 SYNTHETIC_INTEGER **k2_formal_par_in,
4193 SYNTHETIC_INTEGER **h_result,
4194 SYNTHETIC_INTEGER **k_result)
4195 {
4196 SYNTHETIC_INTEGER *arg1_h,
4197 *arg1_k,
4198 *arg2_h,
4199 *arg2_k,
4200 *common_denominator,
4201 *arg1_numerator,
4202 *arg2_numerator,
4203 *cumulative_numerator,
4204 *result_gcd,
4205 *trash_remainder;
4206
4207 unsigned is_neg = FALSE;
4208
4209 /* Be sure the caller isn't doing anything silly with pointers.
4210 */
4211 asAssert(h1_formal_par_in != NULL, __LINE__);
4212 asAssert(*h1_formal_par_in != NULL, __LINE__);
4213 asAssert(k1_formal_par_in != NULL, __LINE__);
4214 asAssert(*k1_formal_par_in != NULL, __LINE__);
4215 asAssert(h2_formal_par_in != NULL, __LINE__);
4216 asAssert(*h2_formal_par_in != NULL, __LINE__);
4217 asAssert(k2_formal_par_in != NULL, __LINE__);
4218 asAssert(*k2_formal_par_in != NULL, __LINE__);
4219 asAssert(h_result != NULL, __LINE__);
4220 asAssert(*h_result != NULL, __LINE__);
4221 asAssert(k_result != NULL, __LINE__);
4222 asAssert(*k_result != NULL, __LINE__);
4223 asAssert(h_result != k_result, __LINE__);
4224 asAssert(*h_result != *k_result, __LINE__);
4225
4226 /* Allocate space for temporary variables.
4227 */
4228 siCreate(&arg1_h);
4229 siCreate(&arg1_k);
4230 siCreate(&arg2_h);
4231 siCreate(&arg2_k);
4232 siCreate(&common_denominator);
4233 siCreate(&arg1_numerator);
4234 siCreate(&arg2_numerator);
4235 siCreate(&cumulative_numerator);
4236 siCreate(&result_gcd);
4237 siCreate(&trash_remainder);
4238
4239 /* Copy over the formal parameters to local variables. This
4240 ** means we can modify the inputs without worry.
4241 */
4242 siCopy(h1_formal_par_in, &arg1_h);
4243 siCopy(k1_formal_par_in, &arg1_k);
4244 siCopy(h2_formal_par_in, &arg2_h);
4245 siCopy(k2_formal_par_in, &arg2_k);
4246
4247 /* If any of the four inputs are NAN, the results must be
4248 ** NAN.
4249 */
4250 if (arg1_h->nan || arg1_k->nan || arg2_h->nan || arg2_k->nan)
4251 {
4252 goto nan_return;
4253 }
4254
4255 /* If either denominator is zero, the result is NAN.
4256 */
4257 if (!(arg1_k->len) || !(arg2_k->len))
4258 {
4259 goto nan_return;
4260 }
4261
4262 /* Adjust the signs of the the arguments so that the denominators are
4263 ** always positive.
4264 */
4265 if (arg1_h->neg && arg1_k->neg)
4266 {
4267 arg1_h->neg = FALSE;
4268 arg1_k->neg = FALSE;
4269 }
4270 else if (!(arg1_h->neg) && arg1_k->neg)
4271 {
4272 arg1_h->neg = TRUE;
4273 arg1_k->neg = FALSE;
4274 }
4275
4276 if (arg2_h->neg && arg2_k->neg)
4277 {
4278 arg2_h->neg = FALSE;
4279 arg2_k->neg = FALSE;
4280 }
4281 else if (!(arg2_h->neg) && arg2_k->neg)
4282 {
4283 arg2_h->neg = TRUE;
4284 arg2_k->neg = FALSE;
4285 }
4286
4287 /* Calculate the common denominator.
4288 */
4289 siUnrestrictedMultiplication(&arg1_k, &arg2_k, &common_denominator);
4290
4291 /* Calculate the arg1 numerator.
4292 */
4293 siUnrestrictedMultiplication(&arg1_h, &arg2_k, &arg1_numerator);
4294
4295 /* Calculate the arg2 numerator.
4296 */
4297 siUnrestrictedMultiplication(&arg2_h, &arg1_k, &arg2_numerator);
4298
4299 /* Subtract the numerators to get the result. */
4300 siUnrestrictedSubtraction(&arg1_numerator,
4301 &arg2_numerator,
4302 &cumulative_numerator);
4303
4304 /* If we NAN'd out ...
4305 */
4306 if (cumulative_numerator->nan)
4307 goto nan_return;
4308
4309 /* If we've reached a result of zero, it won't be legal to form
4310 ** a gcd(). Assign the canonical form of zero.
4311 */
4312 if (!(cumulative_numerator->len))
4313 {
4314 siSetToLong(h_result, 0);
4315 siSetToLong(k_result, 1);
4316 goto normal_return;
4317 }
4318
4319 /* Snatch the sign temporarily.
4320 */
4321 if (cumulative_numerator->neg)
4322 {
4323 is_neg = TRUE;
4324 cumulative_numerator->neg = FALSE;
4325 }
4326
4327 /* Form the gcd. */
4328 siGcd(&cumulative_numerator,
4329 &common_denominator,
4330 &result_gcd);
4331
4332 /* Set the sign back. */
4333 cumulative_numerator->neg = is_neg;
4334
4335 /* Divide both numerator and denominator by
4336 ** the gcd().
4337 */
4338 siUnrestrictedDivision(&cumulative_numerator,
4339 &result_gcd,
4340 h_result,
4341 &trash_remainder);
4342 siUnrestrictedDivision(&common_denominator,
4343 &result_gcd,
4344 k_result,
4345 &trash_remainder);
4346
4347 goto normal_return;
4348 nan_return: ;
4349 siSetToNan(h_result);
4350 siSetToNan(k_result);
4351 normal_return: ;
4352
4353 /* Destroy temporary variables.
4354 */
4355 siDestroy(&arg1_h);
4356 siDestroy(&arg1_k);
4357 siDestroy(&arg2_h);
4358 siDestroy(&arg2_k);
4359 siDestroy(&common_denominator);
4360 siDestroy(&arg1_numerator);
4361 siDestroy(&arg2_numerator);
4362 siDestroy(&cumulative_numerator);
4363 siDestroy(&result_gcd);
4364 siDestroy(&trash_remainder);
4365 }
4366
4367
4368 /****************************************************************************/
4369 /* rnProduct(): */
4370 /*--------------------------------------------------------------------------*/
4371 /* DESCRIPTION */
4372 /* Calculates the product of two rational numbers, returning the */
4373 /* results in lowest terms. */
4374 /* */
4375 /* INPUTS */
4376 /* **h1_formal_par_in, */
4377 /* **k1_formal_par_in, */
4378 /* **h2_formal_par_in, */
4379 /* **k2_formal_par_in : The numerators and denominators of the two */
4380 /* rational numbers to multiply. Zero denominators */
4381 /* are not allowed, but duplicate pointers of all */
4382 /* types are. */
4383 /* **h_result, */
4384 /* **k_result : The result, in lowest terms. Zero is repre- */
4385 /* sented canonically as 0/1. The memory must be */
4386 /* allocated in the caller's area, and the two */
4387 /* data items must be distinct. */
4388 /****************************************************************************/
4389 void rnProduct(SYNTHETIC_INTEGER **h1_formal_par_in,
4390 SYNTHETIC_INTEGER **k1_formal_par_in,
4391 SYNTHETIC_INTEGER **h2_formal_par_in,
4392 SYNTHETIC_INTEGER **k2_formal_par_in,
4393 SYNTHETIC_INTEGER **h_result,
4394 SYNTHETIC_INTEGER **k_result)
4395 {
4396 SYNTHETIC_INTEGER *arg1_h,
4397 *arg1_k,
4398 *arg2_h,
4399 *arg2_k,
4400 *result_numerator,
4401 *result_denominator,
4402 *gcd,
4403 *trash_remainder;
4404
4405 unsigned is_neg = FALSE;
4406
4407 /* Be sure the caller isn't doing anything silly with pointers.
4408 */
4409 asAssert(h1_formal_par_in != NULL, __LINE__);
4410 asAssert(*h1_formal_par_in != NULL, __LINE__);
4411 asAssert(k1_formal_par_in != NULL, __LINE__);
4412 asAssert(*k1_formal_par_in != NULL, __LINE__);
4413 asAssert(h2_formal_par_in != NULL, __LINE__);
4414 asAssert(*h2_formal_par_in != NULL, __LINE__);
4415 asAssert(k2_formal_par_in != NULL, __LINE__);
4416 asAssert(*k2_formal_par_in != NULL, __LINE__);
4417 asAssert(h_result != NULL, __LINE__);
4418 asAssert(*h_result != NULL, __LINE__);
4419 asAssert(k_result != NULL, __LINE__);
4420 asAssert(*k_result != NULL, __LINE__);
4421 asAssert(h_result != k_result, __LINE__);
4422 asAssert(*h_result != *k_result, __LINE__);
4423
4424 /* Allocate space for temporary variables.
4425 */
4426 siCreate(&arg1_h);
4427 siCreate(&arg1_k);
4428 siCreate(&arg2_h);
4429 siCreate(&arg2_k);
4430 siCreate(&result_numerator);
4431 siCreate(&result_denominator);
4432 siCreate(&gcd);
4433 siCreate(&trash_remainder);
4434
4435 /* Copy over the formal parameters to local variables. This
4436 ** means we can modify the inputs without worry.
4437 */
4438 siCopy(h1_formal_par_in, &arg1_h);
4439 siCopy(k1_formal_par_in, &arg1_k);
4440 siCopy(h2_formal_par_in, &arg2_h);
4441 siCopy(k2_formal_par_in, &arg2_k);
4442
4443 /* If any of the four inputs are NAN, the results must be
4444 ** NAN.
4445 */
4446 if (arg1_h->nan || arg1_k->nan || arg2_h->nan || arg2_k->nan)
4447 {
4448 goto nan_return;
4449 }
4450
4451 /* If either denominator is zero, the result is NAN.
4452 */
4453 if (!(arg1_k->len) || !(arg2_k->len))
4454 {
4455 goto nan_return;
4456 }
4457
4458 /* Carry out the multiplication of numerators and denominators to get
4459 ** the two products.
4460 */
4461 siUnrestrictedMultiplication(&arg1_h, &arg2_h, &result_numerator);
4462 siUnrestrictedMultiplication(&arg1_k, &arg2_k, &result_denominator);
4463
4464 /* If either result is NAN, the results are NAN.
4465 */
4466 if (result_numerator->nan || result_denominator->nan)
4467 {
4468 goto nan_return;
4469 }
4470
4471 /* If the numerator is zero, we must return canonical zero.
4472 */
4473 if (!(result_numerator->len))
4474 {
4475 siSetToLong(h_result, 0);
4476 siSetToLong(k_result, 1);
4477 goto normal_return;
4478 }
4479
4480 /* We need to acquire the sign and make both results positive otherwise
4481 ** can't obtain the GCD.
4482 */
4483 if (result_numerator->neg)
4484 {
4485 is_neg = TRUE;
4486 result_numerator->neg = FALSE;
4487 }
4488 if (result_denominator->neg)
4489 {
4490 is_neg = !is_neg;
4491 result_denominator->neg = FALSE;
4492 }
4493
4494 /* Form the gcd. */
4495 siGcd(&result_numerator,
4496 &result_denominator,
4497 &gcd);
4498
4499 /* Set the sign back. */
4500 result_numerator->neg = is_neg;
4501
4502 /* Divide both numerator and denominator by
4503 ** the gcd().
4504 */
4505 siUnrestrictedDivision(&result_numerator,
4506 &gcd,
4507 h_result,
4508 &trash_remainder);
4509 siUnrestrictedDivision(&result_denominator,
4510 &gcd,
4511 k_result,
4512 &trash_remainder);
4513
4514 goto normal_return;
4515 nan_return: ;
4516 siSetToNan(h_result);
4517 siSetToNan(k_result);
4518 normal_return: ;
4519
4520 /* Destroy temporary variables.
4521 */
4522 siDestroy(&arg1_h);
4523 siDestroy(&arg1_k);
4524 siDestroy(&arg2_h);
4525 siDestroy(&arg2_k);
4526 siDestroy(&result_numerator);
4527 siDestroy(&result_denominator);
4528 siDestroy(&gcd);
4529 siDestroy(&trash_remainder);
4530 }
4531
4532
4533 /****************************************************************************/
4534 /* rnQuotient(): */
4535 /*--------------------------------------------------------------------------*/
4536 /* DESCRIPTION */
4537 /* Calculates the quotient of two rational numbers, returning the */
4538 /* results in lowest terms. */
4539 /* */
4540 /* INPUTS */
4541 /* **h1_formal_par_in, */
4542 /* **k1_formal_par_in, */
4543 /* **h2_formal_par_in, */
4544 /* **k2_formal_par_in : The numerators and denominators of the two */
4545 /* rational numbers to divide. Zero denominators */
4546 /* and zero divisors are not allowed, but duplicate */
4547 /* pointers of all types are. */
4548 /* **h_result, */
4549 /* **k_result : The result, in lowest terms. Zero is repre- */
4550 /* sented canonically as 0/1. The memory must be */
4551 /* allocated in the caller's area, and the two */
4552 /* data items must be distinct. */
4553 /****************************************************************************/
4554 void rnQuotient(SYNTHETIC_INTEGER **h1_formal_par_in,
4555 SYNTHETIC_INTEGER **k1_formal_par_in,
4556 SYNTHETIC_INTEGER **h2_formal_par_in,
4557 SYNTHETIC_INTEGER **k2_formal_par_in,
4558 SYNTHETIC_INTEGER **h_result,
4559 SYNTHETIC_INTEGER **k_result)
4560 {
4561 SYNTHETIC_INTEGER *arg1_h,
4562 *arg1_k,
4563 *arg2_h,
4564 *arg2_k,
4565 *result_numerator,
4566 *result_denominator,
4567 *gcd,
4568 *trash_remainder;
4569
4570 unsigned is_neg = FALSE;
4571
4572 /* Be sure the caller isn't doing anything silly with pointers.
4573 */
4574 asAssert(h1_formal_par_in != NULL, __LINE__);
4575 asAssert(*h1_formal_par_in != NULL, __LINE__);
4576 asAssert(k1_formal_par_in != NULL, __LINE__);
4577 asAssert(*k1_formal_par_in != NULL, __LINE__);
4578 asAssert(h2_formal_par_in != NULL, __LINE__);
4579 asAssert(*h2_formal_par_in != NULL, __LINE__);
4580 asAssert(k2_formal_par_in != NULL, __LINE__);
4581 asAssert(*k2_formal_par_in != NULL, __LINE__);
4582 asAssert(h_result != NULL, __LINE__);
4583 asAssert(*h_result != NULL, __LINE__);
4584 asAssert(k_result != NULL, __LINE__);
4585 asAssert(*k_result != NULL, __LINE__);
4586 asAssert(h_result != k_result, __LINE__);
4587 asAssert(*h_result != *k_result, __LINE__);
4588
4589 /* Allocate space for temporary variables.
4590 */
4591 siCreate(&arg1_h);
4592 siCreate(&arg1_k);
4593 siCreate(&arg2_h);
4594 siCreate(&arg2_k);
4595 siCreate(&result_numerator);
4596 siCreate(&result_denominator);
4597 siCreate(&gcd);
4598 siCreate(&trash_remainder);
4599
4600 /* Copy over the formal parameters to local variables. This
4601 ** means we can modify the inputs without worry.
4602 */
4603 siCopy(h1_formal_par_in, &arg1_h);
4604 siCopy(k1_formal_par_in, &arg1_k);
4605 siCopy(h2_formal_par_in, &arg2_h);
4606 siCopy(k2_formal_par_in, &arg2_k);
4607
4608 /* If any of the four inputs are NAN, the results must be
4609 ** NAN.
4610 */
4611 if (arg1_h->nan || arg1_k->nan || arg2_h->nan || arg2_k->nan)
4612 {
4613 goto nan_return;
4614 }
4615
4616 /* If either denominator is zero, the result is NAN.
4617 */
4618 if (!(arg1_k->len) || !(arg2_k->len))
4619 {
4620 goto nan_return;
4621 }
4622
4623 /* Carry out the multiplication of numerators and denominators to get
4624 ** the two products.
4625 */
4626 siUnrestrictedMultiplication(&arg1_h, &arg2_k, &result_numerator);
4627 siUnrestrictedMultiplication(&arg1_k, &arg2_h, &result_denominator);
4628
4629 /* If either result is NAN, the results are NAN.
4630 */
4631 if (result_numerator->nan || result_denominator->nan)
4632 {
4633 goto nan_return;
4634 }
4635
4636 /* If the denominator of the result is zero, this is a no-no-NAN, too.
4637 */
4638 if (!(result_denominator->len))
4639 {
4640 goto nan_return;
4641 }
4642
4643 /* If the numerator is zero, we must return canonical zero.
4644 */
4645 if (!(result_numerator->len))
4646 {
4647 siSetToLong(h_result, 0);
4648 siSetToLong(k_result, 1);
4649 goto normal_return;
4650 }
4651
4652 /* We need to acquire the sign and make both results positive otherwise
4653 ** can't obtain the GCD.
4654 */
4655 if (result_numerator->neg)
4656 {
4657 is_neg = TRUE;
4658 result_numerator->neg = FALSE;
4659 }
4660 if (result_denominator->neg)
4661 {
4662 is_neg = !is_neg;
4663 result_denominator->neg = FALSE;
4664 }
4665
4666 /* Form the gcd. */
4667 siGcd(&result_numerator,
4668 &result_denominator,
4669 &gcd);
4670
4671 /* Set the sign back. */
4672 result_numerator->neg = is_neg;
4673
4674 /* Divide both numerator and denominator by
4675 ** the gcd().
4676 */
4677 siUnrestrictedDivision(&result_numerator,
4678 &gcd,
4679 h_result,
4680 &trash_remainder);
4681 siUnrestrictedDivision(&result_denominator,
4682 &gcd,
4683 k_result,
4684 &trash_remainder);
4685
4686 goto normal_return;
4687 nan_return: ;
4688 siSetToNan(h_result);
4689 siSetToNan(k_result);
4690 normal_return: ;
4691
4692 /* Destroy temporary variables.
4693 */
4694 siDestroy(&arg1_h);
4695 siDestroy(&arg1_k);
4696 siDestroy(&arg2_h);
4697 siDestroy(&arg2_k);
4698 siDestroy(&result_numerator);
4699 siDestroy(&result_denominator);
4700 siDestroy(&gcd);
4701 siDestroy(&trash_remainder);
4702 }
4703
4704
4705 /****************************************************************************/
4706 /* rnCanonize(): */
4707 /*--------------------------------------------------------------------------*/
4708 /* DESCRIPTION */
4709 /* Puts the rational number passed into a canonical form, which means: */
4710 /* a)If either synthetic integer is NAN, both are marked NAN. */
4711 /* b)If the denominator is zero, both synthetic integers are marked */
4712 /* NAN. */
4713 /* c)If the number has a value of zero, it is set to the canonical */
4714 /* form of 0/1. */
4715 /* d)If the number represents a positive value, the signs of both */
4716 /* synthetic integers are set +. If the number represents a */
4717 /* negative value, the numerator will be set - and the denominator */
4718 /* -. */
4719 /* e)Any g.c.d. will be removed so that the number is in lowest */
4720 /* terms. */
4721 /* */
4722 /* INPUTS (AND OUTPUTS) */
4723 /* **h_formal_par_in, */
4724 /* **k_formal_par_in : The numerator and denominator of the number to */
4725 /* operate on. */
4726 /****************************************************************************/
4727 void rnCanonize(SYNTHETIC_INTEGER **h_formal_par_in,
4728 SYNTHETIC_INTEGER **k_formal_par_in)
4729 {
4730 int is_neg = FALSE;
4731 SYNTHETIC_INTEGER *gcd, *trash_remainder, *si_temp;
4732
4733 /* Be sure the caller isn't doing anything silly with pointers.
4734 */
4735 asAssert(h_formal_par_in != NULL, __LINE__);
4736 asAssert(*h_formal_par_in != NULL, __LINE__);
4737 asAssert(k_formal_par_in != NULL, __LINE__);
4738 asAssert(*k_formal_par_in != NULL, __LINE__);
4739 asAssert(h_formal_par_in != k_formal_par_in, __LINE__);
4740 asAssert(*h_formal_par_in != *k_formal_par_in, __LINE__);
4741
4742 /* Create the local integers.
4743 */
4744 siCreate(&gcd);
4745 siCreate(&trash_remainder);
4746 siCreate(&si_temp);
4747
4748 /* If either input is NAN, must mark both NAN.
4749 */
4750 if ((*h_formal_par_in)->nan || (*k_formal_par_in)->nan)
4751 goto nan_return;
4752
4753 /* If the denominator is zero, must also mark both NAN.
4754 */
4755 if (!((*k_formal_par_in)->len))
4756 goto nan_return;
4757
4758 /* If the numerator is zero, the right answer is canonical zero, and
4759 ** must return canonical 0/1.
4760 */
4761 if (!((*h_formal_par_in)->len))
4762 goto zero_return;
4763
4764 /* If we made it to this point, we have two non-zero integers, and the
4765 ** result will be non-zero. Must canonize the sign and factor out
4766 ** any g.c.d.
4767 */
4768 if ((((*h_formal_par_in)->neg) && !((*k_formal_par_in)->neg))
4769 ||
4770 (!((*h_formal_par_in)->neg) && ((*k_formal_par_in)->neg)))
4771 {
4772 is_neg = TRUE;
4773 }
4774
4775 (*h_formal_par_in)->neg = FALSE;
4776 (*k_formal_par_in)->neg = FALSE;
4777
4778 siGcd(h_formal_par_in, k_formal_par_in, &gcd);
4779
4780 siUnrestrictedDivision(h_formal_par_in, &gcd,
4781 &si_temp, &trash_remainder);
4782 siCopy(&si_temp, h_formal_par_in);
4783 siUnrestrictedDivision(k_formal_par_in, &gcd,
4784 &si_temp, &trash_remainder);
4785 siCopy(&si_temp, k_formal_par_in);
4786
4787 if (is_neg)
4788 (*h_formal_par_in)->neg = TRUE;
4789
4790 goto normal_return;
4791 nan_return: ;
4792 siSetToNan(h_formal_par_in);
4793 siSetToNan(k_formal_par_in);
4794 goto normal_return;
4795 zero_return:
4796 siSetToLong(h_formal_par_in, 0);
4797 siSetToLong(k_formal_par_in, 1);
4798 goto normal_return;
4799 normal_return: ;
4800
4801 /* Destroy the local integers.
4802 */
4803 siDestroy(&gcd);
4804 siDestroy(&trash_remainder);
4805 siDestroy(&si_temp);
4806 }
4807
4808
4809 /****************************************************************************/
4810 /* rnCompare(): */
4811 /*--------------------------------------------------------------------------*/
4812 /* DESCRIPTION */
4813 /* Compares two rational numbers to determine thier ranking on the number */
4814 /* line. */
4815 /* */
4816 /* INPUTS . */
4817 /* **h_1, **k_1, */
4818 /* **h_2, **k_2 : The two rational numbers to compare. They are not */
4819 /* required to be in canonical form, but the */
4820 /* denominators may not be zero. */
4821 /* */
4822 /* OUTPUT */
4823 /* <-- : -1 if h1/k1 < h2/k2 */
4824 /* 0 if h1/k1 = h2/k2 */
4825 /* 1 if h1/k1 > h2/k2 */
4826 /****************************************************************************/
4827 int rnCompare(SYNTHETIC_INTEGER **h_1,
4828 SYNTHETIC_INTEGER **k_1,
4829 SYNTHETIC_INTEGER **h_2,
4830 SYNTHETIC_INTEGER **k_2)
4831 {
4832 int rv=0;
4833
4834 SYNTHETIC_INTEGER *h1,
4835 *k1,
4836 *h2,
4837 *k2,
4838 *left_cross_product,
4839 *right_cross_product;
4840
4841 /* Be sure the caller isn't doing anything silly with pointers.
4842 */
4843 asAssert(h_1 != NULL, __LINE__);
4844 asAssert(*h_1 != NULL, __LINE__);
4845 asAssert(k_1 != NULL, __LINE__);
4846 asAssert(*k_1 != NULL, __LINE__);
4847 asAssert(h_2 != NULL, __LINE__);
4848 asAssert(*h_2 != NULL, __LINE__);
4849 asAssert(k_2 != NULL, __LINE__);
4850 asAssert(*k_2 != NULL, __LINE__);
4851
4852 /* Allocate all of our temps.
4853 */
4854 siCreate(&h1);
4855 siCreate(&k1);
4856 siCreate(&h2);
4857 siCreate(&k2);
4858 siCreate(&left_cross_product);
4859 siCreate(&right_cross_product);
4860
4861 /* Copy the formal parameters to the temps so we can mess
4862 ** around with them.
4863 */
4864 siCopy(h_1, &h1);
4865 siCopy(k_1, &k1);
4866 siCopy(h_2, &h2);
4867 siCopy(k_2, &k2);
4868
4869 /* The denominators may not be zero.
4870 */
4871 asAssert((k1->len != 0) && (k2->len != 0) ,__LINE__);
4872
4873 /* We need to normalize the signs on the two fractions.
4874 ** The cross-product inequality utilized breaks down
4875 ** if either denominator is negative.
4876 */
4877 if ((!(h1->neg) && (k1->neg)) || ((h1->neg) && !(k1->neg)))
4878 {
4879 h1->neg = TRUE;
4880 k1->neg = FALSE;
4881 }
4882 else
4883 {
4884 h1->neg = FALSE;
4885 k1->neg = FALSE;
4886 }
4887
4888 if ((!(h2->neg) && (k2->neg)) || ((h2->neg) && !(k2->neg)))
4889 {
4890 h2->neg = TRUE;
4891 k2->neg = FALSE;
4892 }
4893 else
4894 {
4895 h2->neg = FALSE;
4896 k2->neg = FALSE;
4897 }
4898
4899 /* Form the left cross-product and right cross-product. If we get
4900 ** a NAN, we have to kill this function, as it would compromise
4901 ** the integrity of the program (there is no notion of NAN for a
4902 ** comparison result).
4903 */
4904 siUnrestrictedMultiplication(&h1, &k2, &left_cross_product);
4905 siUnrestrictedMultiplication(&h2, &k1, &right_cross_product);
4906 asAssert(!(left_cross_product->nan) && !(right_cross_product->nan), __LINE__);
4907
4908 /* The comparison result is the same as the relationship between the
4909 ** two cross-products. A few words should be said about this.
4910 **
4911 ** Let a/b and c/d be the two rational numbers.
4912 ** Assume that b>0 and d>0.
4913 ** ad < bc
4914 ** -> ad/b < c
4915 ** ->a/b < c/d.
4916 ** and the same logic can be applied to the other relational
4917 ** operators.
4918 */
4919 rv = siCompare(&left_cross_product, &right_cross_product);
4920
4921 #if 0
4922 siDump(&left_cross_product, "left");
4923 siDump(&right_cross_product, "right");
4924 #endif
4925
4926 /* Destroy all of our temps.
4927 */
4928 siDestroy(&h1);
4929 siDestroy(&k1);
4930 siDestroy(&h2);
4931 siDestroy(&k2);
4932 siDestroy(&left_cross_product);
4933 siDestroy(&right_cross_product);
4934
4935 /* Return the return value.
4936 */
4937 return(rv);
4938 }
4939
4940
4941 /****************************************************************************/
4942 /* rnDap(): */
4943 /*--------------------------------------------------------------------------*/
4944 /* DESCRIPTION */
4945 /* Converts a rational number to an approximately equivalent value */
4946 /* with a different denominator. This functionality is the same as */
4947 /* implemented by the DAP command. */
4948 /* */
4949 /* INPUTS . */
4950 /* **h, **k : The rational number to represent differently. */
4951 /* The denominator may not be zero. */
4952 /* */
4953 /* **D : The new denominator. May not be zero. */
4954 /* */
4955 /* OUTPUTS */
4956 /* **N : The new numerator. */
4957 /****************************************************************************/
4958 void rnDap(SYNTHETIC_INTEGER **h,
4959 SYNTHETIC_INTEGER **k,
4960 SYNTHETIC_INTEGER **N,
4961 SYNTHETIC_INTEGER **D)
4962 {
4963 SYNTHETIC_INTEGER *numerator_product,
4964 *quotient,
4965 *trash_remainder;
4966
4967 /* Be sure the caller is doing nothing silly with pointers.
4968 */
4969 asAssert(h != NULL, __LINE__);
4970 asAssert(*h != NULL, __LINE__);
4971 asAssert(k != NULL, __LINE__);
4972 asAssert(*k != NULL, __LINE__);
4973 asAssert(N != NULL, __LINE__);
4974 asAssert(*N != NULL, __LINE__);
4975 asAssert(D != NULL, __LINE__);
4976 asAssert(*D != NULL, __LINE__);
4977
4978 /* Allocate our locals.
4979 */
4980 siCreate(&numerator_product);
4981 siCreate(&quotient);
4982 siCreate(&trash_remainder);
4983
4984 /* If anything is NAN, the result is necessarily NAN.
4985 */
4986 if ((*h)->nan || (*k)->nan || (*D)->nan)
4987 goto nan_return;
4988
4989 /* If the old denominator is zero, the result is necessarily NAN.
4990 */
4991 if (!((*k)->len))
4992 goto nan_return;
4993
4994 /* If the new denominator requested is zero, the result is necessarily NAN.
4995 */
4996 if (!((*D)->len))
4997 goto nan_return;
4998
4999 /* Calculate the numerator as specified in the manual under the DAP
5000 ** command.
5001 */
5002 siUnrestrictedMultiplication(h, D, &numerator_product);
5003
5004 /* Divide to get the value.
5005 */
5006 siUnrestrictedDivision(&numerator_product, k, &quotient, &trash_remainder);
5007
5008 /* Copy to caller's area.
5009 */
5010 siCopy(&quotient, N);
5011
5012 goto normal_return;
5013 nan_return:
5014 siSetToNan(N);
5015 normal_return: ;
5016
5017 /* Destroy our locals.
5018 */
5019 siDestroy(&numerator_product);
5020 siDestroy(&quotient);
5021 siDestroy(&trash_remainder);
5022 }
5023
5024
5025 /****************************************************************************/
5026 /* rnFareyTraverse(): */
5027 /*--------------------------------------------------------------------------*/
5028 /* DESCRIPTION */
5029 /* Uses the standard recursive formulas to traverse the Farey series */
5030 /* in either direction. Two successive terms, in lowest terms, are */
5031 /* needed. */
5032 /* */
5033 /* INPUTS . */
5034 /* **h1in, **k1in: The first rational number in the series of interest, */
5035 /* which must have a positive denominator and a non- */
5036 /* negative numerator. */
5037 /* */
5038 /* **h2in, **k2in: The second rational number in the series of interest, */
5039 /* same rules as above and must be greater than first */
5040 /* term h1in/k1in. */
5041 /* */
5042 /* **N : The order of the series to form. */
5043 /* */
5044 /* dir : (-1): The series is advanced backwards. r(2)=r(1), */
5045 /* and r(1) = new term < r(2). */
5046 /* (+1): The series is advanced forwards. r(1)=r(2), */
5047 /* and r(2) = new term > r(1). */
5048 /****************************************************************************/
5049 void rnFareyTraverse(SYNTHETIC_INTEGER **h1in,
5050 SYNTHETIC_INTEGER **k1in,
5051 SYNTHETIC_INTEGER **h2in,
5052 SYNTHETIC_INTEGER **k2in,
5053 SYNTHETIC_INTEGER **N,
5054 int dir)
5055 {
5056 SYNTHETIC_INTEGER *t1, *t2, *t3, *new_h, *new_k;
5057
5058 /* Be sure the caller isn't doing anything silly with pointers.
5059 */
5060 asAssert(h1in != NULL, __LINE__);
5061 asAssert(*h1in != NULL, __LINE__);
5062 asAssert(k1in != NULL, __LINE__);
5063 asAssert(*k1in != NULL, __LINE__);
5064 asAssert(h2in != NULL, __LINE__);
5065 asAssert(*h2in != NULL, __LINE__);
5066 asAssert(k2in != NULL, __LINE__);
5067 asAssert(*k2in != NULL, __LINE__);
5068 asAssert(N != NULL, __LINE__);
5069 asAssert(*N != NULL, __LINE__);
5070
5071 /* It is also critical that none of the pointers be duplicates,
5072 ** but this isn't checked.
5073 */
5074
5075 /* Allocate local variables.
5076 */
5077 siCreate(&t1);
5078 siCreate(&t2);
5079 siCreate(&t3);
5080 siCreate(&new_h);
5081 siCreate(&new_k);
5082
5083 /* Split into cases based on the direction.
5084 */
5085 if (dir > 0)
5086 {
5087 /* Forward
5088 */
5089 /* Calculate the term floor((k[j-2] + N)/k[j-1]), which we'll need
5090 ** twice. The result is left in "t2".
5091 */
5092 siUnrestrictedAddition(k1in, N, &t1);
5093 siUnrestrictedDivision(&t1, k2in, &t2, &t3);
5094
5095 /* Multiply by h[j-1] and place in "t1".
5096 */
5097 siUnrestrictedMultiplication(&t2, h2in, &t1);
5098
5099 /* Subtract off h[j-2] and this is our new h.
5100 */
5101 siUnrestrictedSubtraction(&t1, h1in, &new_h);
5102
5103 /* Multiply intermediate term by k[j-1] and place in "t1".
5104 */
5105 siUnrestrictedMultiplication(&t2, k2in, &t1);
5106
5107 /* Subtract off k[j-2] and this is our new k.
5108 */
5109 siUnrestrictedSubtraction(&t1, k1in, &new_k);
5110
5111 /* Place the new results.
5112 */
5113 siCopy(h2in, h1in);
5114 siCopy(k2in, k1in);
5115 siCopy(&new_h, h2in);
5116 siCopy(&new_k, k2in);
5117 }
5118 else
5119 {
5120 /* Reverse: steps symmetrical with forward.
5121 */
5122 siUnrestrictedAddition(k2in, N, &t1);
5123 siUnrestrictedDivision(&t1, k1in, &t2, &t3);
5124 siUnrestrictedMultiplication(&t2, h1in, &t1);
5125 siUnrestrictedSubtraction(&t1, h2in, &new_h);
5126 siUnrestrictedMultiplication(&t2, k1in, &t1);
5127 siUnrestrictedSubtraction(&t1, k2in, &new_k);
5128
5129 /* Place the new results.
5130 */
5131 siCopy(h1in, h2in);
5132 siCopy(k1in, k2in);
5133 siCopy(&new_h, h1in);
5134 siCopy(&new_k, k1in);
5135 }
5136
5137
5138 /* Destroy local variables.
5139 */
5140 siDestroy(&t1);
5141 siDestroy(&t2);
5142 siDestroy(&t3);
5143 siDestroy(&new_h);
5144 siDestroy(&new_k);
5145 }
5146
5147
5148 /****************************************************************************/
5149 /****************************************************************************/
5150 /******* C O N T I N U E D F R A C T I O N F U N C T I O N S ******/
5151 /****************************************************************************/
5152 /****************************************************************************/
5153 /* This section is reserved for functions which form and manipulate
5154 ** continued fractions and convergents. Note that every number has at
5155 ** least a 0th order CF partial quotient and at least a 0th order
5156 ** convergent (the integer part).
5157 */
5158 typedef struct
5159 {
5160 SYNTHETIC_INTEGER *raw_numerator;
5161 SYNTHETIC_INTEGER *raw_denominator;
5162 /* The raw numerator and denominator passed to the function which
5163 ** calculates partial quotients and convergents. This might
5164 ** not be in lowest terms. No negative integers are allowed,
5165 ** and zero must be represented canonically as 0/(D>0).
5166 */
5167 SYNTHETIC_INTEGER *numerator;
5168 SYNTHETIC_INTEGER *denominator;
5169 /* The numerator and denominator, in lowest terms, of the
5170 ** rational number whose continued fraction expansion we
5171 ** are forming. No negative integers are allowed, and
5172 ** zero must be represented canonically as 0/1.
5173 */
5174 int n;
5175 /* The number of elements in the parallel arrays of partial
5176 ** quotients, etc., which are maintained. This must be
5177 ** at least 1, which would mean that only the 0th-order
5178 ** items are filled in.
5179 */
5180 SYNTHETIC_INTEGER **a;
5181 /* Pointer to a dynamically allocated array of n pointers,
5182 ** each of which points to a synthetic integer. This
5183 ** pointer may be shifted on any operation on this data
5184 ** structure (due to the behavior of "realloc()"), so
5185 ** a caller should never retain internal pointers when
5186 ** making function calls which might resize anything
5187 ** internally in this data structure. The same sizing arguments
5188 ** here apply to all the other parallel elements, so it won't
5189 ** be described again.
5190 */
5191 SYNTHETIC_INTEGER **p;
5192 SYNTHETIC_INTEGER **q;
5193 /* Convergents of continued fraction expansion of the
5194 ** rational number.
5195 */
5196 } CF_EXPANSION;
5197
5198
5199 /****************************************************************************/
5200 /* pqCreate(): */
5201 /*--------------------------------------------------------------------------*/
5202 /* DESCRIPTION */
5203 /* Creates the continued fraction partial quotient expansion of a non- */
5204 /* negative rational number, and also creates the convergents. */
5205 /* */
5206 /* INPUTS */
5207 /* **h_in */
5208 /* **k_in : The numerator and denominator of the rational */
5209 /* number whose continued fraction expansion and */
5210 /* convergents to form. Both must be positive. */
5211 /* */
5212 /* **expansion : The continued fraction expansion and the */
5213 /* convergents. */
5214 /****************************************************************************/
5215 void pqCreate(SYNTHETIC_INTEGER **h_in,
5216 SYNTHETIC_INTEGER **k_in,
5217 CF_EXPANSION **expansion)
5218 {
5219 SYNTHETIC_INTEGER *dividend,
5220 *divisor,
5221 *quotient,
5222 *remainder,
5223 *p_k_minus_1,
5224 *p_k_minus_2,
5225 *q_k_minus_1,
5226 *q_k_minus_2,
5227 *si_temp1,
5228 *si_temp2,
5229 *si_temp3,
5230 *si_temp4;
5231
5232 /* Be sure that the caller isn't doing anything silly.
5233 */
5234 asAssert(h_in != NULL, __LINE__);
5235 asAssert(*h_in != NULL, __LINE__);
5236 asAssert(k_in != NULL, __LINE__);
5237 asAssert(*k_in != NULL, __LINE__);
5238 asAssert(expansion != NULL, __LINE__);
5239
5240 /* The numerator and denominator in cannot be negative,
5241 ** and the denominator cannot be zero.
5242 */
5243 asAssert(!((*h_in)->neg), __LINE__);
5244 asAssert(!((*k_in)->neg), __LINE__);
5245 asAssert((*k_in)->len != 0, __LINE__);
5246
5247 /* Allocate space for all of the temporary integers
5248 ** we use during the process of forming partial
5249 ** quotients and convergents.
5250 */
5251 siCreate(&dividend);
5252 siCreate(&divisor);
5253 siCreate(&quotient);
5254 siCreate(&remainder);
5255 siCreate(&p_k_minus_1);
5256 siCreate(&p_k_minus_2);
5257 siCreate(&q_k_minus_1);
5258 siCreate(&q_k_minus_2);
5259 siCreate(&si_temp1);
5260 siCreate(&si_temp2);
5261 siCreate(&si_temp3);
5262 siCreate(&si_temp4);
5263
5264 /* printf("Entering function.\n"); */
5265
5266 /* Allocate the memory for the head data block.
5267 */
5268 *expansion = maMalloc(sizeof(CF_EXPANSION));
5269
5270 /* Set all of the data elements of the head data block
5271 ** to known default values.
5272 */
5273 (*expansion)->raw_numerator = NULL;
5274 (*expansion)->raw_denominator = NULL;
5275 (*expansion)->numerator = NULL;
5276 (*expansion)->denominator = NULL;
5277 (*expansion)->n = 0;
5278 (*expansion)->a = NULL;
5279 (*expansion)->p = NULL;
5280 (*expansion)->q = NULL;
5281
5282 /* Assign in the original raw numerator and
5283 ** denominator.
5284 */
5285 siCreate(&((*expansion)->raw_numerator));
5286 siCreate(&((*expansion)->raw_denominator));
5287 siCopy(h_in, &((*expansion)->raw_numerator));
5288 siCopy(k_in, &((*expansion)->raw_denominator));
5289
5290 /* Allocate the space for the lowest-terms numerator
5291 ** and denominator, but let them be zero for now.
5292 */
5293 siCreate(&((*expansion)->numerator));
5294 siCreate(&((*expansion)->denominator));
5295
5296 /* Begin with the dividend and divisor as the
5297 ** numerator and denominator.
5298 */
5299 siCopy(&((*expansion)->raw_numerator), &dividend);
5300 siCopy(&((*expansion)->raw_denominator), &divisor);
5301
5302 /* Enter a do ... while() loop to compute the continued
5303 ** fraction partial quotients and convergents. Because
5304 ** the convergents don't require any "look-ahead" into
5305 ** the partial quotients, they can be computed
5306 ** in parallel
5307 */
5308
5309 do
5310 {
5311 int curidx;
5312 /* Current index.
5313 */
5314
5315 /* Buffer the current index value to the local variable (less
5316 ** typing). The current index is what we're currently
5317 ** operating on.
5318 */
5319 curidx = (*expansion)->n;
5320
5321 /* Increment the number of elements in the
5322 ** three parallel arrays.
5323 */
5324 ((*expansion)->n)++;
5325
5326 /* Allocate the memory for the memory block in the
5327 ** four parallel arrays. There are two cases to
5328 ** consider, either this is the first element or
5329 ** else not the first.
5330 */
5331 if (!((*expansion)->a))
5332 {
5333 /* First time, first element.
5334 */
5335 (*expansion)->a =
5336 maMalloc(sizeof(SYNTHETIC_INTEGER *));
5337 (*expansion)->p =
5338 maMalloc(sizeof(SYNTHETIC_INTEGER *));
5339 (*expansion)->q =
5340 maMalloc(sizeof(SYNTHETIC_INTEGER *));
5341 }
5342 else
5343 {
5344 /* Not the first time. Blow up the arrays of
5345 ** pointers to a larger size to accomodate one more
5346 ** pointer.
5347 */
5348 (*expansion)->a =
5349 maRealloc((*expansion)->a,
5350 (sizeof(SYNTHETIC_INTEGER *)) * (curidx+1));
5351 (*expansion)->p =
5352 maRealloc((*expansion)->p,
5353 (sizeof(SYNTHETIC_INTEGER *)) * (curidx+1));
5354 (*expansion)->q =
5355 maRealloc((*expansion)->q,
5356 (sizeof(SYNTHETIC_INTEGER *)) * (curidx+1));
5357 }
5358
5359 /* Allocate the synthetic integers to go along with
5360 ** blown-up arrays.
5361 */
5362 siCreate((*expansion)->a + curidx);
5363 siCreate((*expansion)->p + curidx);
5364 siCreate((*expansion)->q + curidx);
5365
5366 /* Calculate the current continued fraction partial quotient,
5367 ** and bump the dividend and divisor for the next round.
5368 */
5369 siUnrestrictedDivision(&dividend,
5370 &divisor,
5371 (*expansion)->a + curidx,
5372 &remainder);
5373 /* siDump((*expansion)->a + curidx, "a_k"); */
5374 siCopy(&divisor, &dividend);
5375 siCopy(&remainder, &divisor);
5376
5377 /* Calculate the convergents using the standard
5378 ** recursive formulas.
5379 */
5380 if (curidx == 0)
5381 {
5382 /* p(0) = a(0) */
5383 siCopy((*expansion)->a + 0, (*expansion)->p + 0);
5384 /* q(0) = 1 */
5385 siSetToLong((*expansion)->q + 0, 1);
5386 }
5387 else if (curidx == 1)
5388 {
5389 /* p(1) = a(1)p(0) + 1 */
5390 siSetToLong(&si_temp1, 1);
5391 siUnrestrictedMultiplication((*expansion)->a + 1,
5392 (*expansion)->p + 0,
5393 &si_temp2);
5394 siUnrestrictedAddition(&si_temp2,
5395 &si_temp1,
5396 (*expansion)->p + 1);
5397 /* q(1) = a(1) */
5398 siCopy((*expansion)->a + 1, (*expansion)->q + 1);
5399 }
5400 else /* curidx >= 2 */
5401 {
5402 /* In this case, apply the full recursive formulas. */
5403 /* p(k) = a(k)p(k-1) + p(k-2) */
5404 siUnrestrictedMultiplication((*expansion)->a + curidx,
5405 (*expansion)->p + (curidx-1),
5406 &si_temp2);
5407 siUnrestrictedAddition(&si_temp2,
5408 (*expansion)->p + (curidx-2),
5409 (*expansion)->p + curidx);
5410
5411 /* q(k) = a(k)q(k-1) + q(k-2) */
5412 siUnrestrictedMultiplication((*expansion)->a + curidx,
5413 (*expansion)->q + (curidx-1),
5414 &si_temp2);
5415 siUnrestrictedAddition(&si_temp2,
5416 (*expansion)->q + (curidx-2),
5417 (*expansion)->q + curidx);
5418 }
5419
5420 /* siDump((*expansion)->p + curidx, "p_k"); */
5421 /* siDump((*expansion)->q + curidx, "q_k"); */
5422
5423 /* Put in a safety test for the loop. I've never had this loop fail to
5424 ** terminate. Mathematically, it *can't* fail to terminate, but if there
5425 ** were a bug somewhere in the large integer math or a NAN condition,
5426 ** I wouldn't want to have a machine lockup without a diagnostic message.
5427 ** The value of 32000 is used because that approaches the limits of the
5428 ** minimum that an ANSI 'C' int is required to hold. If we've gone
5429 ** around this loop 32000 times, something is very wrong.
5430 */
5431 asAssert(curidx < 32000, __LINE__);
5432 }
5433 while (remainder->len);
5434
5435 /* The lowest terms representation of the rational number supplied will
5436 ** be the final convergent. This doesn't apply very much to this program,
5437 ** because the parsing functions automatically strip out the gcd() before
5438 ** passing data on, so the "raw" and "final" will be the same. Will assign
5439 ** it anyway.
5440 */
5441 siCopy((*expansion)->p + ((*expansion)->n - 1), &((*expansion)->numerator));
5442 siCopy((*expansion)->q + ((*expansion)->n - 1), &((*expansion)->denominator));
5443
5444 /* Destroy the temporary integers.
5445 */
5446 siDestroy(&dividend);
5447 siDestroy(&divisor);
5448 siDestroy(&quotient);
5449 siDestroy(&remainder);
5450 siDestroy(&p_k_minus_1);
5451 siDestroy(&p_k_minus_2);
5452 siDestroy(&q_k_minus_1);
5453 siDestroy(&q_k_minus_2);
5454 siDestroy(&si_temp1);
5455 siDestroy(&si_temp2);
5456 siDestroy(&si_temp3);
5457 siDestroy(&si_temp4);
5458 }
5459
5460
5461 /****************************************************************************/
5462 /* pqDestroy(): */
5463 /*--------------------------------------------------------------------------*/
5464 /* DESCRIPTION */
5465 /* Deallocates the dynamic data structures associated with the CF */
5466 /* expansion and convergents, and sets the caller's pointer to NULL. */
5467 /* */
5468 /* **expansion : The continued fraction expansion and the */
5469 /* convergents. */
5470 /****************************************************************************/
5471 void pqDestroy(CF_EXPANSION **expansion)
5472 {
5473 int idx;
5474
5475 /* Be sure the caller isn't doing anything silly. with pointers.
5476 */
5477 asAssert(expansion != NULL, __LINE__);
5478 asAssert(*expansion != NULL, __LINE__);
5479
5480 /* The general strategy at this point is to deallocate the data structure
5481 ** from the bottom up. It must be done this way because once a higher
5482 ** data structure is deallocated, the memory ain't yours, so any access
5483 ** even to deallocate "lower" pointers is a violation. Checks all along
5484 ** the way are performed to be sure that nothing looks suspicious.
5485 */
5486
5487 asAssert((*expansion)->n > 0, __LINE__); /* Even a 0th order expansion (an integer)
5488 ** has the "n" at 1.
5489 */
5490 /* Check for suspicious conditions in the base data
5491 ** structure.
5492 */
5493 asAssert((*expansion)->raw_numerator != NULL, __LINE__);
5494 asAssert((*expansion)->raw_denominator != NULL, __LINE__);
5495 asAssert((*expansion)->numerator != NULL, __LINE__);
5496 asAssert((*expansion)->denominator != NULL, __LINE__);
5497 asAssert((*expansion)->a != NULL, __LINE__);
5498 asAssert((*expansion)->p != NULL, __LINE__);
5499 asAssert((*expansion)->q != NULL, __LINE__);
5500
5501
5502 /* Deallocate the directly linked synthetic integers.
5503 */
5504 siDestroy(&((*expansion)->raw_numerator));
5505 siDestroy(&((*expansion)->raw_denominator));
5506 siDestroy(&((*expansion)->numerator));
5507 siDestroy(&((*expansion)->denominator));
5508
5509 /* Deallocate each of the synthetic integers that are the partial
5510 ** quotients and convergents.
5511 */
5512 for (idx = 0; idx < ((*expansion)->n); idx++)
5513 {
5514 asAssert((((*expansion)->a)[idx]) != NULL, __LINE__);
5515 asAssert((((*expansion)->p)[idx]) != NULL, __LINE__);
5516 asAssert((((*expansion)->q)[idx]) != NULL, __LINE__);
5517
5518 siDestroy((*expansion)->a + idx);
5519 siDestroy((*expansion)->p + idx);
5520 siDestroy((*expansion)->q + idx);
5521 }
5522
5523 /* Deallocate the arrays of pointers.
5524 */
5525 maFree((*expansion)->a);
5526 maFree((*expansion)->p);
5527 maFree((*expansion)->q);
5528
5529 /* Deallocate the base data structure and set the caller's pointer
5530 ** to NULL.
5531 */
5532 maFree(*expansion);
5533 *expansion = NULL;
5534 }
5535
5536
5537 /****************************************************************************/
5538 /* pqDump(): */
5539 /*--------------------------------------------------------------------------*/
5540 /* DESCRIPTION */
5541 /* Prints the entire CF expansion and list of convergents to the standard */
5542 /* output stream, with an optional description. */
5543 /* */
5544 /* **expansion : The continued fraction expansion and the */
5545 /* convergents. */
5546 /* *desc : Description to use. If this is a zero-length */
5547 /* string, omits the description. */
5548 /* p_raw : Boolean flag to indicate if the raw numerator and */
5549 /* denominator should be printed distinctly from the */
5550 /* final CF convergent. For most of the application */
5551 /* in this program, the answer is no, because the */
5552 /* rational number is already reduced at the time it */
5553 /* is CF'd, and the raw rational number will be */
5554 /* identical to the final convergent. To print it */
5555 /* would just waste space. */
5556 /****************************************************************************/
5557 void pqDump(CF_EXPANSION **expansion, char *desc, int p_raw)
5558 {
5559 int i;
5560 char buf[100];
5561
5562 /* Be sure the caller isn't doing any pointer suicide.
5563 */
5564 asAssert(expansion != NULL, __LINE__);
5565 asAssert(*expansion != NULL, __LINE__);
5566 asAssert(desc != NULL, __LINE__);
5567
5568 /* If a description was passed print it out as a banner headline.
5569 */
5570 if (strlen(desc))
5571 {
5572 gfBannerHeading(desc, 0);
5573 }
5574
5575 /* Banner announcing inputs. */
5576 gfBannerHeading("Inputs To CF Calculation", 0);
5577 gfHline();
5578
5579 /* Print out the raw numerator and denominator. */
5580 siDump(&((*expansion)->raw_numerator), "h_in");
5581 gfHline();
5582 siDump(&((*expansion)->raw_denominator), "k_in");
5583 gfHline();
5584
5585 /* Print out each of the continued-fraction partial
5586 ** quotients.
5587 */
5588 /* Banner announcing partial quotients. */
5589 gfBannerHeading("CF Partial Quotients", 0);
5590 gfHline();
5591
5592 for (i=0; i<((*expansion)->n); i++)
5593 {
5594 sprintf(buf, "a(%d)", i);
5595 siDump((*expansion)->a + i, buf);
5596 gfHline();
5597 }
5598
5599 /* Priint out the convergents. */
5600 /* Banner announcing partial convergents. */
5601 gfBannerHeading("CF Convergents", 0);
5602 gfHline();
5603
5604 for (i=0; i<((*expansion)->n); i++)
5605 {
5606 sprintf(buf, "p(%d)", i);
5607 siDump((*expansion)->p + i, buf);
5608 sprintf(buf, "q(%d)", i);
5609 siDump((*expansion)->q + i, buf);
5610 gfHline();
5611 }
5612 }
5613
5614
5615 /****************************************************************************/
5616