1 USING THE ARBITRARY PRECISION ROUTINES IN A C PROGRAM 2 3Part of the calc release consists of an arbitrary precision math link library. 4This link library is used by the calc program to perform its own calculations. 5If you wish, you can ignore the calc program entirely and call the arbitrary 6precision math routines from your own C programs. 7 8The link library is called libcalc.a, and provides routines to handle arbitrary 9precision arithmetic with integers, rational numbers, or complex numbers. 10There are also many numeric functions such as factorial and gcd, along 11with some transcendental functions such as sin and exp. 12 13Take a look at the sample sub-directory. It contains a few simple 14examples of how to use libcalc.a that might be helpful to look at 15after you have read this file. 16 17------------------ 18FIRST THINGS FIRST 19------------------ 20 21............................................................................... 22. . 23. You MUST call libcalc_call_me_first() prior to using libcalc lib functions! . 24. . 25............................................................................... 26 27The function libcalc_call_me_first() takes no args and returns void. You 28need call libcalc_call_me_first() only once. 29 30------------- 31INCLUDE FILES 32------------- 33 34To use any of these routines in your own programs, you need to include the 35appropriate include file. These include files are: 36 37 zmath.h (for integer arithmetic) 38 qmath.h (for rational arithmetic) 39 cmath.h (for complex number arithmetic) 40 41You never need to include more than one of the above files, even if you wish 42to use more than one type of arithmetic, since qmath.h automatically includes 43zmath.h, and cmath.h automatically includes qmath.h. 44 45The prototypes for the available routines are listed in the above include 46files. Some of these routines are meant for internal use, and so aren't 47convenient for outside use. So you should read the source for a routine 48to see if it really does what you think it does. I won't guarantee that 49obscure internal routines won't change or disappear in future releases! 50 51When calc is installed, all of libraries are installed into ${LIBDIR}. 52All of the calc header files are installed under ${INCDIRCALC}. 53 54If CALC_SRC is defined, then the calc header files will assume that 55they are in or under the current directory. However, most external 56programs most likely will not be located under calc'c source tree. 57External programs most likely want to use the installed calc header 58files under ${INCDIRCALC}. External programs most likely NOT want 59to define CALC_SRC. 60 61You need to include the following file to get the symbols and variables 62related to error handling: 63 64 lib_calc.h 65 66External programs may want to compile with: 67 68 -I${INCDIR} -L${LIBDIR} -lcalc 69 70If custom functions are also used, they may want to compile with: 71 72 -I${INCDIR} -L${LIBDIR} -lcalc -lcustcalc 73 74The CALC_SRC symbol should NOT be defined by default. However if you are 75feeling pedantic you may want to force CALC_SRC to be undefined: 76 77 -UCALC_SRC 78 79as well. 80 81------------------- 82MATH ERROR HANDLING 83------------------- 84 85The math_error() function is called by the math routines on an error 86condition, such as malloc failures, division by zero, or some form of 87an internal computation error. The routine is called in the manner of 88printf, with a format string and optional arguments: 89 90 void math_error(char *fmt, ...); 91 92Your program must handle math errors in one of three ways: 93 94 1) Print the error message and then exit 95 96 There is a math_error() function supplied with the calc library. 97 By default, this routine simply prints a message to stderr and 98 then exits. By simply linking in this link library, any calc 99 errors will result in a error message on stderr followed by 100 an exit. 101 102 2) Use setjmp and longjmp in your program 103 104 Use setjmp at some appropriate level in your program, and let 105 the longjmp in math_error() return to that level and to allow you 106 to recover from the error. This is what the calc program does. 107 108 If one sets up calc_matherr_jmpbuf, and then sets 109 calc_use_matherr_jmpbuf to non-zero then math_error() will 110 longjmp back with the return value of calc_use_matherr_jmpbuf. 111 In addition, the last calc error message will be found in 112 calc_err_msg; this error is not printed to stderr. The calc 113 error message will not have a trailing newline. 114 115 For example: 116 117 #include <setjmp.h> 118 #include "lib_calc.h" 119 120 int error; 121 122 ... 123 124 if ((error = setjmp(calc_matherr_jmpbuf)) != 0) { 125 126 /* report the error */ 127 printf("Ouch: %s\n", calc_err_msg); 128 129 /* reinitialize calc after the longjmp */ 130 reinitialize(); 131 } 132 calc_use_matherr_jmpbuf = 1; 133 134 If calc_use_matherr_jmpbuf is non-zero, then the jmp_buf value 135 calc_matherr_jmpbuf must be initialized by the setjmp() function 136 or your program will crash. 137 138 3) Supply your own math_error function: 139 140 void math_error(char *fmt, ...); 141 142 Your math_error() function may exit or transfer control to outside 143 of the calc library, but it must never return or calc will crash. 144 145External programs can obtain the appropriate calc symbols by compiling with: 146 147 -I${INCDIR} -L${LIBDIR} -lcalc 148 149------------------------- 150PARSE/SCAN ERROR HANDLING 151------------------------- 152 153The scanerror() function is called when calc encounters a parse/scan 154error. For example, scanerror() is called when calc is given code 155with a syntax error. 156 157The variable, calc_print_scanerr_msg, controls if calc prints to stderr, 158any parse/scan errors. By default, this variable it set to 1 and so 159parse/scan errors are printed to stderr. By setting this value to zero, 160parse/scan errors are not printed: 161 162 #include "lib_calc.h" 163 164 /* do not print parse/scan errors to stderr */ 165 calc_print_scanerr_msg = 0; 166 167The last calc math error or calc parse/scan error message is kept 168in the NUL terminated buffer: 169 170 char calc_err_msg[MAXERROR+1]; 171 172The value of calc_print_scanerr_msg does not change the use 173of the calc_err_msg[] buffer. Messages are stored in that 174buffer regardless of the calc_print_scanerr_msg value. 175 176The calc_print_scanerr_msg and the calc_err_msg[] buffer are declared 177lib_calc.h include file. The initialized storage for these variables 178comes from the calc library. The MAXERROR symbol is also declared in 179the lib_calc.h include file. 180 181Your program must handle parse/scan errors in one of two ways: 182 183 1) exit on error 184 185 If you do not setup the calc_scanerr_jmpbuf, then when calc 186 encounters a parse/scan error, a message will be printed to 187 stderr and calc will exit. 188 189 2) Use setjmp and longjmp in your program 190 191 Use setjmp at some appropriate level in your program, and let 192 the longjmp in scanerror() return to that level and to allow you 193 to recover from the error. This is what the calc program does. 194 195 If one sets up calc_scanerr_jmpbuf, and then sets 196 calc_use_scanerr_jmpbuf to non-zero then scanerror() will longjmp 197 back with the return with a non-zero code. In addition, the last 198 calc error message will be found in calc_err_msg[]; this error is 199 not printed to stderr. The calc error message will not have a 200 trailing newline. 201 202 For example: 203 204 #include <setjmp.h> 205 #include "lib_calc.h" 206 207 int scan_error; 208 209 ... 210 211 /* delay the printing of the parse/scan error */ 212 calc_use_scanerr_jmpbuf = 0; /* this is optional */ 213 214 if ((scan_error = setjmp(calc_scanerr_jmpbuf)) != 0) { 215 216 /* report the parse/scan */ 217 if (calc_use_scanerr_jmpbuf == 0) { 218 printf("parse error: %s\n", calc_err_msg); 219 } 220 221 /* initialize calc after the longjmp */ 222 initialize(); 223 } 224 calc_use_scanerr_jmpbuf = 1; 225 226 If calc_use_scanerr_jmpbuf is non-zero, then the jmp_buf value 227 calc_scanerr_jmpbuf must be initialized by the setjmp() function 228 or your program will crash. 229 230External programs can obtain the appropriate calc symbols by compiling with: 231 232 -I${INCDIR} -L${LIBDIR} -lcalc 233 234--------------------------- 235PARSE/SCAN WARNING HANDLING 236--------------------------- 237 238Calc parse/scan warning message are printed to stderr by the warning() 239function. The routine is called in the manner of printf, with a format 240string and optional arguments: 241 242 void warning(char *fmt, ...); 243 244The variable, calc_print_scanwarn_msg, controls if calc prints to stderr, 245any parse/scan warnings. By default, this variable it set to 1 and so 246parse/scan warnings are printed to stderr. By setting this value to zero, 247parse/scan warnings are not printed: 248 249 #include "lib_calc.h" 250 251 /* do not print parse/scan warnings to stderr */ 252 calc_print_scanwarn_msg = 0; 253 254The last calc calc parse/scan warning message is kept in the NUL 255terminated buffer: 256 257 char calc_warn_msg[MAXERROR+1]; 258 259The value of calc_print_scanwarn_msg does not change the use 260of the calc_warn_msg[] buffer. Messages are stored in that 261buffer regardless of the calc_print_scanwarn_msg value. 262 263Your program must handle parse/scan warnings in one of two ways: 264 265 1) print the warning to stderr and continue 266 267 The warning() from libcalc prints warning messages to 268 stderr and returns. The flow of execution is not changed. 269 This is what calc does by default. 270 271 2) Supply your own warning function: 272 273 void warning(char *fmt, ...); 274 275 Your warning function should simply return when it is finished. 276 277External programs can obtain the appropriate calc symbols by compiling with: 278 279 -I${INCDIR} -L${LIBDIR} -lcalc 280 281 282--------------- 283OUTPUT ROUTINES 284--------------- 285 286The output from the routines in the link library normally goes to stdout. 287You can divert that output to either another FILE handle, or else 288to a string. Read the routines in zio.c to see what is available. 289Diversions can be nested. 290 291You use math_setfp to divert output to another FILE handle. Calling 292math_setfp with stdout restores output to stdout. 293 294Use math_divertio to begin diverting output into a string. Calling 295math_getdivertedio will then return a string containing the output, and 296clears the diversion. The string is reallocated as necessary, but since 297it is in memory, there are obviously limits on the amount of data that can 298be diverted into it. The string needs freeing when you are done with it. 299 300Calling math_cleardiversions will clear all the diversions to strings, and 301is useful on an error condition to restore output to a known state. You 302should also call math_setfp on errors if you had changed that. 303 304If you wish to mix your own output with numeric output from the math routines, 305then you can call math_chr, math_str, math_fill, math_fmt, or math_flush. 306These routines output single characters, output null-terminated strings, 307output strings with space filling, output formatted strings like printf, and 308flush the output. Output from these routines is diverted as described above. 309 310You can change the default output mode by calling math_setmode, and you can 311change the default number of digits printed by calling math_setdigits. These 312routines return the previous values. The possible modes are described in 313zmath.h. 314 315-------------- 316USING INTEGERS 317-------------- 318 319The arbitrary precision integer routines define a structure called a ZVALUE. 320This is defined in zmath.h. A ZVALUE contains a pointer to an array of 321integers, the length of the array, and a sign flag. The array is allocated 322using malloc, so you need to free this array when you are done with a 323ZVALUE. To do this, you should call zfree() with the ZVALUE as an argument 324and never try to free the array yourself using free(). The reason for this 325is that sometimes the pointer points to a statically allocated arrays which 326should NOT be freed. 327 328The ZVALUE structures are passed to routines by value, and are returned 329through pointers. For example, to multiply two small integers together, 330you could do the following: 331 332 ZVALUE z1, z2, z3; 333 334 itoz(3L, &z1); 335 itoz(4L, &z2); 336 zmul(z1, z2, &z3); 337 338Use zcopy to copy one ZVALUE to another. There is no sharing of arrays 339between different ZVALUEs even if they have the same value, so you MUST 340use this routine. Simply assigning one value into another will cause 341problems when one of the copies is freed. However, the special ZVALUE 342values _zero_ and _one_ CAN be assigned to variables directly, since their 343values of 0 and 1 are so common that special checks are made for them. 344 345For initial values besides 0 or 1, you need to call itoz to convert a long 346value into a ZVALUE, as shown in the above example. Or alternatively, 347for larger numbers you can use the atoz routine to convert a string which 348represents a number into a ZVALUE. The string can be in decimal, octal, 349hex, or binary according to the leading digits. 350 351Always make sure you free a ZVALUE when you are done with it or when you 352are about to overwrite an old ZVALUE with another value by passing its 353address to a routine as a destination value, otherwise memory will be 354lost. The following shows an example of the correct way to free memory 355over a long sequence of operations. 356 357 ZVALUE z1, z2, z3; 358 359 z1 = _one_; 360 atoz("12345678987654321", &z2); 361 zadd(z1, z2, &z3); 362 zfree(z1); 363 zfree(z2); 364 zsquare(z3, &z1); 365 zfree(z3); 366 itoz(17L, &z2); 367 zsub(z1, z2, &z3); 368 zfree(z1); 369 zfree(z2); 370 zfree(z3); 371 372There are some quick checks you can make on integers. For example, whether 373or not they are zero, negative, even, and so on. These are all macros 374defined in zmath.h, and should be used instead of checking the parts of the 375ZVALUE yourself. Examples of such checks are: 376 377 ziseven(z) (number is even) 378 zisodd(z) (number is odd) 379 ziszero(z) (number is zero) 380 zisneg(z) (number is negative) 381 zispos(z) (number is positive) 382 zisunit(z) (number is 1 or -1) 383 zisone(z) (number is 1) 384 zisnegone(z) (number is -1) 385 zistwo(z) (number is 2) 386 zisabstwo(z) (number is 2 or -2) 387 zisabsleone(z) (number is -1, 0 or 1) 388 zislezero(z) (number is <= 0) 389 zisleone(z) (number is <= 1) 390 zge16b(z) (number is >= 2^16) 391 zge24b(z) (number is >= 2^24) 392 zge31b(z) (number is >= 2^31) 393 zge32b(z) (number is >= 2^32) 394 zge64b(z) (number is >= 2^64) 395 396Typically the largest unsigned long is typedefed to FULL. The following 397macros are useful in dealing with this data type: 398 399 MAXFULL (largest positive FULL value) 400 MAXUFULL (largest unsigned FULL value) 401 zgtmaxfull(z) (number is > MAXFULL) 402 zgtmaxufull(z) (number is > MAXUFULL) 403 zgtmaxlong(z) (number is > MAXLONG, largest long value) 404 zgtmaxulong(z) (number is > MAXULONG, largest unsigned long value) 405 406If zgtmaxufull(z) is false, then one may quickly convert the absolute 407value of number into a full with the macro: 408 409 ztofull(z) (convert abs(number) to FULL) 410 ztoulong(z) (convert abs(number) to an unsigned long) 411 ztolong(z) (convert abs(number) to a long) 412 413If the value is too large for ztofull(), ztoulong() or ztolong(), only 414the low order bits converted. 415 416There are two types of comparisons you can make on ZVALUEs. This is whether 417or not they are equal, or the ordering on size of the numbers. The zcmp 418function tests whether two ZVALUEs are equal, returning TRUE if they differ. 419The zrel function tests the relative sizes of two ZVALUEs, returning -1 if 420the first one is smaller, 0 if they are the same, and 1 if the first one 421is larger. 422 423--------------- 424USING FRACTIONS 425--------------- 426 427The arbitrary precision fractional routines define a structure called NUMBER. 428This is defined in qmath.h. A NUMBER contains two ZVALUEs for the numerator 429and denominator of a fraction, and a count of the number of uses there are 430for this NUMBER. The numerator and denominator are always in lowest terms, 431and the sign of the number is contained in the numerator. The denominator 432is always positive. If the NUMBER is an integer, the denominator has the 433value 1. 434 435Unlike ZVALUEs, NUMBERs are passed using pointers, and pointers to them are 436returned by functions. So the basic type for using fractions is not really 437(NUMBER), but is (NUMBER *). NUMBERs are allocated using the qalloc routine. 438This returns a pointer to a number which has the value 1. Because of the 439special property of a ZVALUE of 1, the numerator and denominator of this 440returned value can simply be overwritten with new ZVALUEs without needing 441to free them first. The following illustrates this: 442 443 NUMBER *q; 444 445 q = qalloc(); 446 itoz(55L, &q->num); 447 448A better way to create NUMBERs with particular values is to use the itoq, 449iitoq, or atoq functions. Using itoq makes a long value into a NUMBER, 450using iitoq makes a pair of longs into the numerator and denominator of a 451NUMBER (reducing them first if needed), and atoq converts a string representing 452a number into the corresponding NUMBER. The atoq function accepts input in 453integral, fractional, real, or exponential formats. Examples of allocating 454numbers are: 455 456 NUMBER *q1, *q2, *q3; 457 458 q1 = itoq(66L); 459 q2 = iitoq(2L, 3L); 460 q3 = atoq("456.78"); 461 462Also unlike ZVALUEs, NUMBERs are quickly copied. This is because they contain 463a link count, which is the number of pointers there are to the NUMBER. The 464qlink macro is used to copy a pointer to a NUMBER, and simply increments 465the link count and returns the same pointer. Since it is a macro, the 466argument should not be a function call, but a real pointer variable. The 467qcopy routine will actually make a new copy of a NUMBER, with a new link 468count of 1. This is not usually needed. 469 470NUMBERs are deleted using the qfree routine. This decrements the link count 471in the NUMBER, and if it reaches zero, then it will deallocate both of 472the ZVALUEs contained within the NUMBER, and then puts the NUMBER structure 473onto a free list for quick reuse. The following is an example of allocating 474NUMBERs, copying them, adding them, and finally deleting them again. 475 476 NUMBER *q1, *q2, *q3; 477 478 q1 = itoq(111L); 479 q2 = qlink(q1); 480 q3 = qqadd(q1, q2); 481 qfree(q1); 482 qfree(q2); 483 qfree(q3); 484 485Because of the passing of pointers and the ability to copy numbers easily, 486you might wish to use the rational number routines even for integral 487calculations. They might be slightly slower than the raw integral routines, 488but are more convenient to program with. 489 490The prototypes for the fractional routines are defined in qmath.h. 491Many of the definitions for integer functions parallel the ones defined 492in zmath.h. But there are also functions used only for fractions. 493Examples of these are qnum to return the numerator, qden to return the 494denominator, qint to return the integer part of, qfrac to return the 495fractional part of, and qinv to invert a fraction. 496 497There are some transcendental functions in the link library, such as sin 498and cos. These cannot be evaluated exactly as fractions. Therefore, 499they accept another argument which tells how accurate you want the result. 500This is an "epsilon" value, and the returned value will be within that 501quantity of the correct value. This is usually an absolute difference, 502but for some functions (such as exp), this is a relative difference. 503For example, to calculate sin(0.5) to 100 decimal places, you could do: 504 505 NUMBER *q, *ans, *epsilon; 506 507 q = atoq("0.5"); 508 epsilon = atoq("1e-100"); 509 ans = qsin(q, epsilon); 510 511There are many convenience macros similar to the ones for ZVALUEs which can 512give quick information about NUMBERs. In addition, there are some new ones 513applicable to fractions. These are all defined in qmath.h. Some of these 514macros are: 515 516 qiszero(q) (number is zero) 517 qisneg(q) (number is negative) 518 qispos(q) (number is positive) 519 qisint(q) (number is an integer) 520 qisfrac(q) (number is fractional) 521 qisunit(q) (number is 1 or -1) 522 qisone(q) (number is 1) 523 qisnegone(q) (number is -1) 524 qistwo(q) (number is 2) 525 qiseven(q) (number is an even integer) 526 qisodd(q) (number is an odd integer) 527 qistwopower(q) (number is a power of 2 >= 1) 528 529The comparisons for NUMBERs are similar to the ones for ZVALUEs. You use the 530qcmp and qrel functions. 531 532There are four predefined values for fractions. You should qlink them when 533you want to use them. These are _qzero_, _qone_, _qnegone_, and _qonehalf_. 534These have the values 0, 1, -1, and 1/2. An example of using them is: 535 536 NUMBER *q1, *q2; 537 538 q1 = qlink(&_qonehalf_); 539 q2 = qlink(&_qone_); 540 541--------------------- 542USING COMPLEX NUMBERS 543--------------------- 544 545The arbitrary precision complex arithmetic routines define a structure 546called COMPLEX. This is defined in cmath.h. This contains two NUMBERs 547for the real and imaginary parts of a complex number, and a count of the 548number of links there are to this COMPLEX number. 549 550The complex number routines work similarly to the fractional routines. 551You can allocate a COMPLEX structure using comalloc (NOT calloc!). 552You can construct a COMPLEX number with desired real and imaginary 553fractional parts using qqtoc. You can copy COMPLEX values using clink 554which increments the link count. And you free a COMPLEX value using cfree. 555The following example illustrates this: 556 557 NUMBER *q1, *q2; 558 COMPLEX *c1, *c2, *c3; 559 560 q1 = itoq(3L); 561 q2 = itoq(4L); 562 c1 = qqtoc(q1, q2); 563 qfree(q1); 564 qfree(q2); 565 c2 = clink(c1); 566 c3 = cmul(c1, c2); 567 cfree(c1); 568 cfree(c2); 569 cfree(c3); 570 571As a shortcut, when you want to manipulate a COMPLEX value by a real value, 572you can use the caddq, csubq, cmulq, and cdivq routines. These accept one 573COMPLEX value and one NUMBER value, and produce a COMPLEX value. 574 575There is no direct routine to convert a string value into a COMPLEX value. 576But you can do this yourself by converting two strings into two NUMBERS, 577and then using the qqtoc routine. 578 579COMPLEX values are always returned from these routines. To split out the 580real and imaginary parts into normal NUMBERs, you can simply qlink the 581two components, as shown in the following example: 582 583 COMPLEX *c; 584 NUMBER *rp, *ip; 585 586 c = calloc(); 587 rp = qlink(c->real); 588 ip = qlink(c->imag); 589 590There are many macros for checking quick things about complex numbers, 591similar to the ZVALUE and NUMBER macros. In addition, there are some 592only used for complex numbers. Examples of macros are: 593 594 cisreal(c) (number is real) 595 cisimag(c) (number is pure imaginary) 596 ciszero(c) (number is zero) 597 cisnegone(c) (number is -1) 598 cisone(c) (number is 1) 599 cisrunit(c) (number is 1 or -1) 600 cisiunit(c) (number is i or -i) 601 cisunit(c) (number is 1, -1, i, or -i) 602 cistwo(c) (number is 2) 603 cisint(c) (number is has integer real and imaginary parts) 604 ciseven(c) (number is has even real and imaginary parts) 605 cisodd(c) (number is has odd real and imaginary parts) 606 607There is only one comparison you can make for COMPLEX values, and that is 608for equality. The ccmp function returns TRUE if two complex numbers differ. 609 610There are three predefined values for complex numbers. You should clink 611them when you want to use them. They are _czero_, _cone_, and _conei_. 612These have the values 0, 1, and i. 613 614---------------- 615LAST THINGS LAST 616---------------- 617 618If you wish, when you are all done you can call libcalc_call_me_last() 619to free a small amount of storage associated with the libcalc_call_me_first() 620call. This is not required, but is does bring things to a closure. 621 622The function libcalc_call_me_last() takes no args and returns void. You 623need call libcalc_call_me_last() only once. 624 625## Copyright (C) 1999,2021 David I. Bell and Landon Curt Noll 626## 627## Calc is open software; you can redistribute it and/or modify it under 628## the terms of the version 2.1 of the GNU Lesser General Public License 629## as published by the Free Software Foundation. 630## 631## Calc is distributed in the hope that it will be useful, but WITHOUT 632## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 633## or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 634## Public License for more details. 635## 636## A copy of version 2.1 of the GNU Lesser General Public License is 637## distributed with calc under the filename COPYING-LGPL. You should have 638## received a copy with calc; if not, write to Free Software Foundation, Inc. 639## 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 640## 641## Under source code control: 1993/07/30 19:44:49 642## File existed as early as: 1993 643## 644## chongo <was here> /\oo/\ http://www.isthe.com/chongo/ 645## Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 646