home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / bigint / bigint.doc next >
Encoding:
Text File  |  1986-04-20  |  25.0 KB  |  875 lines
  1. :gdoc
  2. :frontm
  3. :titlep
  4. :title.Big Integer Package
  5. :author.Richard G. Larson
  6. :address
  7. :aline.170 N. Scoville Ave
  8. :aline. Oak Park, IL 60302
  9. :eaddress
  10. :date.28 June 1985
  11. :copyright 1985 "Richard G. Larson".
  12. :etitlep
  13. :toc
  14. :set tocnum 2
  15. :abstract
  16. :p.Datatypes which represent arbitrary precision integers are described.
  17. A set of functions written in C
  18. to create and manipulate arbitrary precision integers is described.
  19. This package is designed to be used with the CI C86 compiler,
  20. but effort has been made to make it transportable to other C compilers.
  21. :p.Permission is hereby granted to copy this manual and the
  22. associated programs for all non commercial purposes, provided that
  23. the copyright notice and this notice are also copied.
  24. This material may not be sold,
  25. nor may it be used in any program which is sold, without written
  26. permission from the author.
  27. :body
  28. :h1.Introduction.
  29. :p
  30. An arbitrary precision integer is represented as a variable size
  31. 2's-complement number.
  32. It is stored in a structure which includes a
  33. maximum size for the structure
  34. (in digits; a digit is an unsigned short), the current size,
  35. and the digits of the current number.
  36. The only limitations on the maximum size are the fact that the
  37. number of digits must be an integer, and the fact that the C library
  38. functions may limit the size of a single memory allocation.
  39. The type BIGINT is defined to be one of these structures.
  40. The type BIGINTP is defined to be a pointer to a BIGINT.
  41. Normally, all manipulation of arbitrary precision integers
  42. is done with variables of type BIGINTP.
  43. :p
  44. Functions (some of which are implemented as macros) are provided for
  45. the allocation of storage for integers, for copying integers,
  46. for conversion to and from ASCII strings, for conversion to and
  47. from longs,
  48. and for various arithmetic operations and comparisons
  49. on integers.
  50. :p
  51. The arbitrary precision functions are provided in libraries:
  52. BIARTHSN.LIB for small model programs to be used with
  53. an 8087 numeric coprocessor,
  54. BIARTHSS.LIB for small model programs to be used without
  55. an 8087,
  56. and BIARTHBN.LIB for big model programs to be used with
  57. an 8087.
  58. The C and ASM source for the functions is provided in BIARITH.ARC.
  59. In addition, two other files are provided:
  60. BIGINT.H, which defines the arbitrary precision integer structures
  61. and gives the necessary macro definitions,
  62. and SYSSTD.H, which contains some basic constant and type definitions.
  63. These files should be #included in any program using the library.
  64. This manual is provided in both a printable form (in BIGINT.LST)
  65. and with ReadiWriter formatting commands (in BIGINT.DOC).
  66. The C source files of the sample programs from the chapter Programming
  67. Examples are provided (in FACTOR.C and PRIME.C).
  68. :h1.Library Functions.
  69. In general, the parameters of the library functions are of type BIGINTP.
  70. The comparison functions return bool (TRUE/FALSE) values.
  71. Functions which perform operations are called with a pointer
  72. to the destination operand, in which the result is to be placed,
  73. and pointers to the source operands.
  74. These functions return a pointer to the destination operand
  75. if the operation was successfully performed,
  76. or NULL if the function was unable to perform the operation
  77. (usually because the destination did not have adequate space).
  78. The functions are all designed to return NULL if any
  79. parameter is NULL:
  80. this often will allow you to safely compose functions.
  81. (See the chapter entitled Programming Examples for samples of such usage.)
  82. In some cases, if the function returns NULL, the values of the source
  83. operands may have been changed.
  84. All operands are assumed to be expressed with the minimum number of
  85. digits necessary to give a valid 2's-complement representation;
  86. the answer is always returned in this minimum form.
  87. :p
  88. All function and macro names begin with "bi".
  89. In addition, the names of some internal functions and identifiers
  90. begin with "_bi".
  91. The user should avoid using names which begin with "_bi".
  92. :h2.bialloc, Allocate Space for a Big Integer
  93. :h3.Synopsis
  94. :xmp
  95.  
  96.       #include "bigint.h"
  97.       BIGINTP bialloc(m)
  98.       int m;
  99. :exmp
  100. :h3.Function/Returns/Notes
  101. :p
  102. Allocates space for an arbitrary precision integer with maximum size m digits,
  103. and initializes the field giving the maximum size.
  104. If m < 2, allocates space for a two digit number.
  105. Does not initialize the value of the arbitrary precision integer:
  106. using the value of an uninitialized arbitrary precision
  107. integer may have disasterous effects.
  108. Returns a pointer to the allocated arbitrary precision integer, or NULL if the
  109. function was unable to allocate sufficient space.
  110. (With CI-C86, an arbitrary precision
  111. integer with m digits occupies 4+2*m bytes.)
  112. Unused arbitrary precision integers may be freed using
  113. the standard function free().
  114. :h2.bisize, Current Size of a Big Integer
  115. :h3.Synopsis
  116. :xmp
  117.  
  118.       #include "bigint.h"
  119.       int bisize(n)
  120.       BIGINTP n;
  121. :exmp
  122. :h3.Function/Returns/Notes
  123. :p
  124. Returns the number of digits currently being used in the arbitrary
  125. precision integer pointed at by n.
  126. The value returned by bisize() may be used as a parameter for bialloc()
  127. to allocate an arbitrary precision integer just large enough to hold
  128. a copy of the integer pointed at by n.
  129. :p
  130. NOTE: bisize is implemented as a macro.
  131. The function bisize must not be declared.
  132. :h2.bicopy, Copy a Big Integer
  133. :h3.Synopsis
  134. :xmp
  135.  
  136.       #include "bigint.h"
  137.       BIGINTP bicopy(p, q)
  138.       BIGINTP p, q;
  139. :exmp
  140. :h3.Function/Returns/Notes
  141. :p
  142. Copies the arbitrary precision integer pointed at by q
  143. into the arbitrary precision integer pointed at by p.
  144. Returns p, or NULL if there was not enough space
  145. to make the copy.
  146. :h2.binumstr bistrnum, Conversion To/From String
  147. :h3.Synopsis
  148. :xmp
  149.  
  150.       #include "bigint.h"
  151.       BIGINTP bistrnum(p,s)
  152.       BIGINTP p;
  153.       char *s;
  154.  
  155.       #include "bigint.h"
  156.       char *binumstr(s,p)
  157.       char *s;
  158.       BIGINTP p;
  159. :exmp
  160. :h3.Function/Returns/Notes
  161. :p
  162. bistrnum places the number represented by the ASCII string
  163. pointed at by s into the arbitrary precision integer pointed at by p.
  164. Returns the pointer p, or NULL if there was not sufficient
  165. space in the arbitrary precision integer pointed at by p.
  166. Initial non numeric characters are skipped.
  167. Conversion continues until the first non digit is found,
  168. or the end of the string is reached.
  169. If the string contains no digits, 0 is put into the the integer
  170. pointed at by p.
  171. :p
  172. binumstr puts the ASCII string representing the number
  173. contained in the arbitrary precision integer pointed at by p into the
  174. string pointed at by s.
  175. The resulting ASCII string contains only decimal digits,
  176. and if the number is negative, an initial '-'.
  177. The string contains no blanks.
  178. It is terminated with the EOS character.
  179. The function returns s, or NULL if the function was unable to
  180. allocate sufficient temporary storage for scratch work.
  181. The programmer must make certain that the string pointed
  182. at by s contains sufficient space for the resulting number:
  183. the function cannot check this.
  184. :h2.binumlng bilngnum, Convert To/From Long
  185. :h3.Synopsis
  186. :xmp
  187.  
  188.       #include "bigint.h"
  189.       long binumlng(p)
  190.       BIGINTP p;
  191.  
  192.       #include "bigint.h"
  193.       BIGINTP bilngnum(p,n)
  194.       BIGINTP p;
  195.       long n;
  196. :exmp
  197. :h3.Function/Returns/Notes
  198. :p
  199. binumlng converts the arbitrary precision integer pointed at by p to a long,
  200. and returns the long.
  201. If the number contained in the arbitrary precision integer
  202. will not fit into a long,
  203. the function returns the least significant portion of the
  204. 2's complement representation.
  205. :p
  206. bilngnum places the long n
  207. into the arbitrary precision integer pointed at by p.
  208. Returns p, or NULL if the arbitrary precision integer is not large enough.
  209. Note that the arbitrary precision integer must always be capable of holding
  210. at least two digits, even if the specific long would fit into
  211. a smaller arbitrary precision integer.
  212. Note that bialloc() always returns a pointer to a sufficiently
  213. large integer to hold a long.
  214. :h2.binumdbl, Convert To Double
  215. :h3.Synopsis
  216. :xmp
  217.  
  218.       #include "bigint.h"
  219.       double binumdbl(p)
  220.       BIGINTP p;
  221. :exmp
  222. :h3.Function/Returns/Notes
  223. :p
  224. binumdbl converts the arbitrary precision integer pointed at by p to a
  225. double,
  226. and returns the double.
  227. If the number contained in the arbitrary precision integer
  228. is to big to be represented as a double,
  229. the function returns infinity.
  230. :h2.bineg, Negate a Big Integer
  231. :h3.Synopsis
  232. :xmp
  233.  
  234.       #include "bigint.h"
  235.       BIGINTP bineg(d, s)
  236.       BIGINTP d, s;
  237. :exmp
  238. :h3.Function/Returns/Notes
  239. :p
  240. bineg puts the negative of the arbitrary precision
  241. integer pointed at by s into the arbitrary precision integer pointed at by d.
  242. Returns d, or NULL if the arbitrary precision integer
  243. pointed at by d is not large enough.
  244. :h2.biabs, Absolute Value of a Big Integer
  245. :h3.Synopsis
  246. :xmp
  247.  
  248.       #include "bigint.h"
  249.       BIGINTP biabs(d, s)
  250.       BIGINTP d, s;
  251. :exmp
  252. :h3.Function/Returns/Notes
  253. :p
  254. biabs puts the absolute value of the arbitrary precision integer
  255. pointed at by s into the arbitrary precision integer pointed at by d.
  256. Returns d, or NULL if the arbitrary precision integer pointed at by d
  257. is not large enough.
  258. :p
  259. NOTE: biabs is implemented as a macro using bineg and bicopy.
  260. The function biabs must not be declared.
  261. If it is used, the functions bineg and bicopy should be declared.
  262. :h2.biadd bisub, Add/Subtract Big Integers
  263. :h3.Synopsis
  264. :xmp
  265.  
  266.       #include "bigint.h"
  267.       BIGINTP biadd(d, s1, s2)
  268.       BIGINTP d, s1, s2;
  269.  
  270.       #include "bigint.h"
  271.       BIGINTP bisub(d, s1, s2)
  272.       BIGINTP d, s1, s2;
  273. :exmp
  274. :h3.Function/Returns/Notes
  275. :p
  276. biadd puts the sum of the arbitrary precision integers pointed at by s1 and s2
  277. into the arbitrary precision integer pointed at by d.
  278. Returns d, or NULL if there was not enough space in
  279. the arbitrary precision integer pointed at by d to do the sum.
  280. :p
  281. bisub puts the arbitrary precision integer pointed at by s2
  282. subtracted from the arbitrary precision integer
  283. pointed at by s1, into the arbitrary precision integer pointed at by d.
  284. Returns d, or NULL if there was not have enough space to do the subtraction.
  285. (NULL return can be caused by lack of space in the number
  286. pointed at by d, or by s2.)
  287. :h2.bimul, Multiply Big Integers
  288. :h3.Synopsis
  289. :xmp
  290.  
  291.       #include "bigint.h"
  292.       BIGINTP bimul(d, s1, s2)
  293.       BIGINTP d, s1, s2;
  294. :exmp
  295. :h3.Function/Returns/Notes
  296. :p
  297. bimul puts the product of the arbitrary precision integers
  298. pointed at by s1 and s2
  299. into the arbitrary precision integer pointed at by d.
  300. Returns d, or NULL if there was not enough space to do the multiplication.
  301. (NULL return can be caused by lack of space in any operand.)
  302. :h2.biquo, Quotient and Remainder of Big Integers
  303. :h3.Synopsis
  304. :xmp
  305.  
  306.       #include "bigint.h"
  307.       BIGINTP bimul(d, s1, s2)
  308.       BIGINTP d, s1, s2;
  309. :exmp
  310. :h3.Function/Returns/Notes
  311. :p
  312. biquo puts the integer quotient q of the arbitrary precision integers
  313. pointed at by s1 and s2
  314. into the arbitrary precision integer pointed at by d,
  315. and replaces the contents of the arbitrary precision integer
  316. pointed at by s1 with the remainder r.
  317. Returns d, or NULL if there was not enough space to do the division,
  318. or if division by 0 is attempted.
  319. If a is the integer pointed at by s1, and b is the integer pointed at by s2,
  320. then we have
  321. :xmp
  322.       a = q*b + r
  323. :exmp
  324. with the sign of r the same as that of a,
  325. and the absolute value of r less than the absolute value of b.
  326. The sign of q is determined by the rules of algebra.
  327. :h2.bisqrt, Square Root of Big Integer
  328. :h3.Synopsis
  329. :xmp
  330.  
  331.       #include "bigint.h"
  332.       BIGINTP bisqrt(d, s)
  333.       BIGINTP d, s;
  334. :exmp
  335. :h3.Function/Returns/Notes
  336. :p
  337. bisqrt puts the floor of the square root of the arbitrary precision
  338. integer pointed at by s
  339. into the arbitrary precision integer pointed at by d.
  340. Returns d, or NULL if there was not enough space to do the computation,
  341. or if the square root of a negative number is computed.
  342. :h2.bilt bile bigt bige bieq bine, Numeric Comparison
  343. :h3.Synopsis
  344. :xmp
  345.  
  346.       #include "sysstd.h"
  347.       #include "bigint.h"
  348.  
  349.       bool bilt(p, q)
  350.       BIGINTP p, q;
  351.  
  352.       bool bile(p, q)
  353.       BIGINTP p, q;
  354.  
  355.       bool bigt(p, q)
  356.       BIGINTP p, q;
  357.  
  358.       bool bige(p, q)
  359.       BIGINTP p, q;
  360.  
  361.       bool bieq(p, q)
  362.       BIGINTP p, q;
  363.  
  364.       bool bine(p, q)
  365.       BIGINTP p, q;
  366. :exmp
  367. :h3.Function/Returns/Notes
  368. :p
  369. bilt (bile, bigt, bige, bieq, bine) returns TRUE if
  370. the arbitrary precision integer
  371. pointed at by the first argument is less than (respectively
  372. less than or equal to,
  373. greater than, greater than or equal to, equal to, not equal to)
  374. the arbitrary precision integer pointed at by the second argument,
  375. and otherwise returns FALSE.
  376. :p
  377. NOTE: bile, bigt, bige are implemented as macros using bilt.
  378. The functions bile, bigt, bige must not be declared.
  379. If they are used, the function bilt should be declared.
  380. bine is implemented as a macro using bieq.
  381. The function bine must not be declared.
  382. If it is used, the function bieq should be declared.
  383. :h2.bilt0 bile0 bigt0 bige0 bieq0 bine0, 0-Comparison
  384. :h3.Synopsis
  385. :xmp
  386.  
  387.       #include "sysstd.h"
  388.       #include "bigint.h"
  389.  
  390.       bool bilt0(p)
  391.       BIGINTP p;
  392.  
  393.       bool bilet0(p)
  394.       BIGINTP p;
  395.  
  396.       bool bigt0(p)
  397.       BIGINTP p;
  398.  
  399.       bool bige0(p)
  400.       BIGINTP p;
  401.  
  402.       bool bieq0(p)
  403.       BIGINTP p;
  404.  
  405.       bool bine0(p)
  406.       BIGINTP p;
  407. :exmp
  408. :h3.Function/Returns/Notes
  409. :p
  410. bilt0 (bile0, bigt0, bige0, bieq0, bine0) returns TRUE if its argument
  411. is less than 0 (respectively
  412. less than or equal to 0, greater than 0, greater than or equal to 0,
  413. equal to 0, not equal to 0),
  414. and otherwise returns FALSE.
  415. :p
  416. NOTE: The functions bilt0, bile0, bigt0, bige0, bieq0, and bine0
  417. are implemented as macros, and must not be declared.
  418. :h2.bieven biodd, Parity
  419. :h3.Synopsis
  420. :xmp
  421.  
  422.       #include "sysstd.h"
  423.       #include "bigint.h"
  424.  
  425.       bool bieven(p)
  426.       BIGINTP p;
  427.  
  428.       bool biodd(p)
  429.       BIGINTP p;
  430. :exmp
  431. :h3.Function/Returns/Notes
  432. :p
  433. bieven returns TRUE if its argument is an even integer,
  434. and otherwise returns FALSE.
  435. :p
  436. biodd returns TRUE if its argument is an odd integer,
  437. and otherwise returns FALSE.
  438. :p
  439. NOTE: The functions bieven and biodd
  440. are implemented as macros, and must not be declared.
  441. :h1.Programming Examples.
  442. :p
  443. :h2.Prime Factorization
  444. :p
  445. Here is the standard small computer prime factorization program.
  446. You should note how the big integer package insulates you from
  447. specific questions of data representation.
  448. Also note the examples of function composition.
  449. :xmp
  450.  
  451. /*
  452.  * A demonstration program for the big integer package
  453.  * which reads and factors integers.
  454.  *                              R. G. Larson
  455.  *                              21 June 1984
  456.  */
  457. #include "stdio.h"
  458. #include "sysstd.h"
  459. #include "bigint.h"
  460.  
  461. #define STRSIZE 300 /* max size of ASCII string rep of number */
  462. #define NUMSIZE  65 /* max number of digits in number */
  463.  
  464. main ()
  465. {
  466.     BIGINTP bialloc(), bistrnum(), bilngnum(),
  467.             bicopy(), bineg(), bimul(), biadd();
  468.     bool    bilt(), bieq();
  469.  
  470.     BIGINTP d, n, one, two, three, five, tmp, incr[8];
  471.     char    s[STRSIZE];
  472.     int     i;
  473.  
  474.     /* misc small constants */
  475.     one   = bilngnum(bialloc(2), 1L);
  476.     two   = bilngnum(bialloc(2), 2L);
  477.     three = bilngnum(bialloc(2), 3L);
  478.     five  = bilngnum(bialloc(2), 5L);
  479.  
  480.     /* increments for the mod 30 trial divisor loop */
  481.     incr[0] = bilngnum(bialloc(2), 4L);
  482.     incr[1] = bilngnum(bialloc(2), 2L);
  483.     incr[2] = bilngnum(bialloc(2), 4L);
  484.     incr[3] = bilngnum(bialloc(2), 2L);
  485.     incr[4] = bilngnum(bialloc(2), 4L);
  486.     incr[5] = bilngnum(bialloc(2), 6L);
  487.     incr[6] = bilngnum(bialloc(2), 2L);
  488.     incr[7] = bilngnum(bialloc(2), 6L);
  489.  
  490.     tmp = bialloc(NUMSIZE);
  491.     n = bialloc(NUMSIZE);
  492.     d = bialloc(NUMSIZE);
  493.  
  494.     printf ("Enter number to factor: ");
  495.     fgets (s, STRSIZE, stdin);
  496.     /* n = abs(input) */
  497.     if (!biabs(n, bistrnum(tmp, s)))
  498.         abort("Input Error");
  499.     while (bine0(n))   {
  500.  
  501.         trydiv (n, two);
  502.         trydiv (n, three);
  503.         trydiv (n, five);
  504.         bilngnum (d, 7L);
  505.         while (bile(bimul(tmp, d, d), n))
  506.             for (i = 0; i < 8; i++)   {
  507.                 trydiv (n, d);
  508.                 bicopy(d, biadd(tmp, d, incr[i]));
  509.             }
  510.         if (bine(n, one))
  511.             putans (n, 1);
  512.         printf ("Number factored.\n\n");
  513.  
  514.         printf ("Enter number to factor: ");
  515.         fgets (s, STRSIZE, stdin);
  516.         /* n = abs(input) */
  517.         if (!biabs(n, bistrnum(tmp, s)))
  518.             abort("Input Error");
  519.     }
  520. }
  521.  
  522. void trydiv (n, d)
  523. BIGINTP n, d;
  524. {
  525.     BIGINTP bialloc(), bicopy(), biquo();
  526.     static BIGINTP x, y;
  527.     static bool first_time = TRUE;
  528.     int exp = 0;
  529.  
  530.     if (first_time)   {
  531.         x = bialloc(NUMSIZE);
  532.         y = bialloc(NUMSIZE);
  533.         first_time = FALSE;
  534.     }
  535.     bicopy(x, n);
  536.     biquo(y, x, d);
  537.     while (bieq0(x))   {
  538.         exp++;
  539.         bicopy(n, y);
  540.         bicopy(x, y);
  541.         biquo(y, x, d);
  542.     }
  543.     if (exp)
  544.         putans (d, exp);
  545. }
  546.  
  547. void putans (d, e)
  548. BIGINTP d;
  549. int e;
  550. {
  551.     char *binumstr();
  552.     char s[STRSIZE];
  553.  
  554.     binumstr(s, d);
  555.     printf ("Factor: %s ^ %d\n", s, e);
  556. }
  557. :exmp
  558. :h2.Primality Testing
  559. :p
  560. This program factors out small divisors from a number,
  561. and then applies the strong pseudoprime test to the remaining factor
  562. twenty times.
  563. If it fails the strong pseudoprime test once, the number is composite.
  564. If it passes it twenty times, it is almost certainly prime.
  565. :xmp
  566.  
  567. /*
  568.  * A program which removes small prime factors, and then
  569.  * applies the strong pseudo prime test to the remaining factor
  570.  *                              R. G. Larson
  571.  *                              23 June 1985
  572.  */
  573. #include "stdio.h"
  574. #include "sysstd.h"
  575. #include "bigint.h"
  576.  
  577. #define STRSIZE      300 /* max size of ASCII string rep of number */
  578. #define NUMSIZE       65 /* max number of digits in number */
  579. #define DIVLIM     10000 /* limit for small divisors */
  580.  
  581. extern BIGINTP bialloc(), bimul(), biquo(), bicopy(), bilngnum(),
  582.                bicopy(), bisub(), bistrnum(), bineg(), biadd();
  583. extern bool    bilt(), bieq();
  584. extern char    *binumstr();
  585.  
  586. void sqmod (l, n)
  587. BIGINTP l, n;
  588. /*
  589.  * Sets l to l**2 (mod n)
  590.  */
  591. {
  592.    BIGINTP tmp1, tmp2;
  593.    int s;
  594.  
  595.    s = bisize(n);
  596.    tmp1 = bialloc(2 * s + 2);
  597.    tmp2 = bialloc(s + 2);
  598.    bimul(tmp1, l, l);
  599. /* l = tmp1 (mod n) */
  600.    biquo(tmp2, tmp1, n);
  601.    bicopy(l, tmp1);
  602.    free(tmp1);
  603.    free(tmp2);
  604. }
  605.  
  606. void expmod (l, a, e, n)
  607. BIGINTP l, a, e, n;
  608. /*
  609.  * Sets l = a**e (mod n)
  610.  */
  611. {
  612.    static bool    first_time = TRUE;
  613.    static BIGINTP two;
  614.  
  615.    BIGINTP atmp, etmp, tmp1, tmp2;
  616.    int s;
  617.  
  618.    if (first_time)   {
  619.       first_time = FALSE;
  620.       two = bilngnum(bialloc(2), (long) 2);
  621.    }
  622.  
  623.    s = bisize(n);
  624.    atmp = bialloc(s + 2);
  625.    etmp = bialloc(s + 2);
  626.    tmp1 = bialloc(2 * s + 2);
  627.    tmp2 = bialloc(s + 2);
  628.    bicopy(atmp, a);
  629.    bicopy(etmp, e);
  630.    bilngnum(l, (long) 1);
  631.    while (bigt0(etmp))   {
  632.       if (biodd(etmp))   {
  633.          bimul(tmp1, l, atmp);
  634.          biquo(tmp2, tmp1, n);
  635.          bicopy(l, tmp1);
  636.       }
  637.       sqmod(atmp, n);
  638.       biquo(etmp, bicopy(tmp2, etmp), two);
  639.    }
  640.    free(atmp);
  641.    free(etmp);
  642.    free(tmp1);
  643.    free(tmp2);
  644. }
  645.  
  646. bool prime (n, a)
  647. BIGINTP n, a;
  648. /*
  649.  * Test if n passes strong pseudo prime test with powers of a.
  650.  */
  651. {
  652.    static bool first_time = TRUE;
  653.    static BIGINTP one, two;
  654.  
  655.    BIGINTP k, l, tmp, minus_one;
  656.    int e, s;
  657.  
  658.    if (first_time)   {
  659.       first_time = FALSE;
  660.       one       = bilngnum(bialloc(2), (long) 1);
  661.       two       = bilngnum(bialloc(2), (long) 2);
  662.    }
  663.  
  664.    s = bisize(n);
  665.    tmp       = bialloc(s + 2);
  666.    k         = bialloc(s + 2);
  667.    l         = bialloc(s + 2);
  668.    minus_one = bialloc(s + 2);
  669.    bisub(k, n, one);
  670.    bicopy(minus_one, k); /* n - 1 */
  671.    e = 0;
  672.    while (bieven(k))   {
  673.       e++;
  674.       biquo(k, bicopy(tmp, k), two);
  675.    }
  676.    expmod(l, a, k, n);
  677.    if (bieq(l, one) | bieq(l, minus_one))   {
  678.       free(tmp);
  679.       free(k);
  680.       free(l);
  681.       free(minus_one);
  682.       return (TRUE);
  683.    }
  684. /* here l = a ** (odd part of n-1) , and not +/-1 */
  685.    for (e--; e > 0; e--)   {
  686.       sqmod(l, n);
  687.       if (bieq(l, minus_one))   {
  688.          free(tmp);
  689.          free(k);
  690.          free(l);
  691.          free(minus_one);
  692.          return (TRUE);
  693.       }
  694.       if (bieq(l, one))   {
  695.          free(tmp);
  696.          free(k);
  697.          free(l);
  698.          free(minus_one);
  699.          return (FALSE);
  700.       }
  701.    }
  702.    free(tmp);
  703.    free(k);
  704.    free(l);
  705.    free(minus_one);
  706.    return (FALSE);
  707. }
  708.  
  709. bool primetest(n)
  710. BIGINTP n;
  711. /*
  712.  * Does strong pseudo prime test on n
  713.  */
  714. {
  715.    static bool    first_time = TRUE;
  716.    static BIGINTP na;
  717.    static int     a[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 27, 0};
  718.  
  719.    int i;
  720.  
  721.    if (first_time)   {
  722.       first_time = FALSE;
  723.       na = bialloc(2);
  724.    }
  725.    for (i = 0; a[i]; i++)   {
  726.       bilngnum(na, (long) a[i]);
  727.       if (!prime(n, na))
  728.          return (FALSE);
  729.    }
  730.    return (TRUE);
  731. }
  732.  
  733. void trydiv (n, d)
  734. BIGINTP n, d;
  735. {
  736.    static BIGINTP x = NULL,
  737.                   y = NULL;
  738.    static int     xymax = 0;
  739.           int     exp = 0;
  740.  
  741.    if (bisize(n) + 2 > xymax)   {
  742.       xymax = bisize(n) + 2;
  743.       if (x)
  744.          free(x);
  745.       x = bialloc(xymax);
  746.       if (y)
  747.          free(y);
  748.       y = bialloc(xymax);
  749.    }
  750.    bicopy(x, n);
  751.    biquo(y, x, d);
  752.    while (bieq0(x))   {
  753.       exp++;
  754.       bicopy(n, y);
  755.       bicopy(x, y);
  756.       biquo(y, x, d);
  757.    }
  758.    if (exp)
  759.       putans (d, exp);
  760. }
  761.  
  762. void putans (d, e)
  763. BIGINTP d;
  764. int e;
  765. {
  766.    char s[STRSIZE];
  767.  
  768.    binumstr(s, d);
  769.    printf ("Factor: %s ^ %d\n", s, e);
  770. }
  771.  
  772. main ()
  773. {
  774. /* the increments for the mod 30 trial divisors: */
  775.    static int incr[] = {4, 2, 4, 2, 4, 6, 2, 6};
  776.  
  777.    BIGINTP nd, n, one, two, three, five, tmp;
  778.    char    s[STRSIZE];
  779.    long    d;
  780.    int     i;
  781.    bool    test;
  782.  
  783.    /* misc small constants */
  784.    one   = bilngnum(bialloc(2), 1L);
  785.    two   = bilngnum(bialloc(2), 2L);
  786.    three = bilngnum(bialloc(2), 3L);
  787.    five  = bilngnum(bialloc(2), 5L);
  788.  
  789.    tmp = bialloc(NUMSIZE);
  790.    n   = bialloc(NUMSIZE);
  791.    nd  = bialloc(NUMSIZE);
  792.  
  793.    printf ("Enter number to test: ");
  794.    fgets (s, STRSIZE, stdin);
  795.    /* n = abs(input) */
  796.    if (!biabs(n, bistrnum(tmp, s)))
  797.       abort("Input Error");
  798.  
  799.    while (bine0(n))   {
  800.       trydiv (n, two);
  801.       trydiv (n, three);
  802.       trydiv (n, five);
  803.       d = 7L;
  804.       while ((test = bile(bilngnum(tmp, (long) d * d), n)) &&
  805.              (d <= DIVLIM))
  806.          for (i = 0; i < 8; i++)   {
  807.             bilngnum(nd, d);
  808.             trydiv (n, nd);
  809.             d += (long) incr[i];
  810.          }
  811.       if (bieq(n, one))
  812.          printf ("Number factored.\n\n");
  813.       else if (!test)   {
  814.          putans(n, 1);
  815.          printf ("Number factored.\n\n");
  816.       }
  817.       else   {
  818.          binumstr(s, n);
  819.          printf ("Remaining factor: %s ", s);
  820.          if (primetest(n))
  821.             printf ("is probably prime.\n\n");
  822.          else
  823.             printf ("is composite.\n\n");
  824.       }
  825.  
  826.       printf ("Enter number to test: ");
  827.       fgets (s, STRSIZE, stdin);
  828.       /* n = abs(input) */
  829.       if (!biabs(n, bistrnum(tmp, s)))
  830.          abort("Input Error");
  831.    }
  832. }
  833. :exmp
  834. :h1.Installation Notes.
  835. :p
  836. The libraries BIARTHSN.LIB, BIARTHSS.LIB, and BIARTHBN.LIB
  837. have been compiled with Computer Innovations' C86 compiler, version
  838. 2.20J.
  839. :p
  840. In addition to the code for the
  841. multi-precision integer arithmetic package, the libraries also include
  842. code which
  843. replaces two modules in the standard CI-C86 libraries.
  844. The modules ZLDIVMOD and ZLMUL will be linked into the final EXE file
  845. instead of the corresponding modules from the CI-C86 Version 2.20
  846. library.
  847. If you are using another version of the CI-C86 libraries,
  848. you should compare the source code of these modules (in BIARITH.ARC)
  849. with the source code of the corresponding CI-C86 module:
  850. these modules contain internal functions with non standard linkage
  851. conventions.
  852. These linkage conventions could differ in different versions.
  853. Note that the BIARITH version of ZLDIVMOD requires an 8087.
  854. If in doubt, you can remove ZLDIVMOD and ZLMUL from the library.
  855. However, you will pay a severe speed penalty if you use the CI-C86
  856. version of ZLDIVMOD.
  857. :p
  858. Some functions in BIARITH.ARC are provided in both a C and an ASM file.
  859. In this case, the ASM file is a faster version of the C file.
  860. You should use the ASM version, if possible.
  861. Using the two replacements for the CI-C86 library modules and the
  862. three assembly language functions provided, instead of the CI-C86
  863. modules and the C versions of the functions speeded
  864. up the program PRIME given in section 3.2 by more than a factor of 3.
  865. :p
  866. If you need another version of the library,
  867. you will need to recreate it from the source in
  868. BIARITH.ARC.
  869. All of the C files are written in "generic" C.
  870. They should be usable to create versions of the libraries for
  871. other compilers.
  872. The ASM files for the replacements for the modules
  873. in the C86 library are also provided.
  874. :egdoc
  875. i