/src/rt/bigint/bigint_ext.cpp

http://github.com/jruderman/rust · C++ · 553 lines · 403 code · 84 blank · 66 comment · 62 complexity · 69ef421e2b423e04b66db38b55b2c8ac MD5 · raw file

  1. /* bigint_ext - external portion of large integer package
  2. **
  3. ** Copyright © 2000 by Jef Poskanzer <jef@mail.acme.com>.
  4. ** All rights reserved.
  5. **
  6. ** Redistribution and use in source and binary forms, with or without
  7. ** modification, are permitted provided that the following conditions
  8. ** are met:
  9. ** 1. Redistributions of source code must retain the above copyright
  10. ** notice, this list of conditions and the following disclaimer.
  11. ** 2. Redistributions in binary form must reproduce the above copyright
  12. ** notice, this list of conditions and the following disclaimer in the
  13. ** documentation and/or other materials provided with the distribution.
  14. **
  15. ** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  16. ** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  17. ** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  18. ** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  19. ** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  20. ** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  21. ** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  22. ** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  23. ** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  24. ** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  25. ** SUCH DAMAGE.
  26. */
  27. #include <sys/types.h>
  28. #include <signal.h>
  29. #include <stdio.h>
  30. #include <stdlib.h>
  31. #include <unistd.h>
  32. #include <time.h>
  33. #include "bigint.h"
  34. #include "low_primes.h"
  35. bigint bi_0, bi_1, bi_2, bi_10, bi_m1, bi_maxint, bi_minint;
  36. /* Forwards. */
  37. static void print_pos( FILE* f, bigint bi );
  38. bigint
  39. str_to_bi( char* str )
  40. {
  41. int sign;
  42. bigint biR;
  43. sign = 1;
  44. if ( *str == '-' )
  45. {
  46. sign = -1;
  47. ++str;
  48. }
  49. for ( biR = bi_0; *str >= '0' && *str <= '9'; ++str )
  50. biR = bi_int_add( bi_int_multiply( biR, 10 ), *str - '0' );
  51. if ( sign == -1 )
  52. biR = bi_negate( biR );
  53. return biR;
  54. }
  55. void
  56. bi_print( FILE* f, bigint bi )
  57. {
  58. if ( bi_is_negative( bi_copy( bi ) ) )
  59. {
  60. putc( '-', f );
  61. bi = bi_negate( bi );
  62. }
  63. print_pos( f, bi );
  64. }
  65. bigint
  66. bi_scan( FILE* f )
  67. {
  68. int sign;
  69. int c;
  70. bigint biR;
  71. sign = 1;
  72. c = getc( f );
  73. if ( c == '-' )
  74. sign = -1;
  75. else
  76. ungetc( c, f );
  77. biR = bi_0;
  78. for (;;)
  79. {
  80. c = getc( f );
  81. if ( c < '0' || c > '9' )
  82. break;
  83. biR = bi_int_add( bi_int_multiply( biR, 10 ), c - '0' );
  84. }
  85. if ( sign == -1 )
  86. biR = bi_negate( biR );
  87. return biR;
  88. }
  89. static void
  90. print_pos( FILE* f, bigint bi )
  91. {
  92. if ( bi_compare( bi_copy( bi ), bi_10 ) >= 0 )
  93. print_pos( f, bi_int_divide( bi_copy( bi ), 10 ) );
  94. putc( bi_int_mod( bi, 10 ) + '0', f );
  95. }
  96. int
  97. bi_int_mod( bigint bi, int m )
  98. {
  99. int r;
  100. if ( m <= 0 )
  101. {
  102. (void) fprintf( stderr, "bi_int_mod: zero or negative modulus\n" );
  103. (void) kill( getpid(), SIGFPE );
  104. }
  105. r = bi_int_rem( bi, m );
  106. if ( r < 0 )
  107. r += m;
  108. return r;
  109. }
  110. bigint
  111. bi_rem( bigint bia, bigint bim )
  112. {
  113. return bi_subtract(
  114. bia, bi_multiply( bi_divide( bi_copy( bia ), bi_copy( bim ) ), bim ) );
  115. }
  116. bigint
  117. bi_mod( bigint bia, bigint bim )
  118. {
  119. bigint biR;
  120. if ( bi_compare( bi_copy( bim ), bi_0 ) <= 0 )
  121. {
  122. (void) fprintf( stderr, "bi_mod: zero or negative modulus\n" );
  123. (void) kill( getpid(), SIGFPE );
  124. }
  125. biR = bi_rem( bia, bi_copy( bim ) );
  126. if ( bi_is_negative( bi_copy( biR ) ) )
  127. biR = bi_add( biR, bim );
  128. else
  129. bi_free( bim );
  130. return biR;
  131. }
  132. bigint
  133. bi_square( bigint bi )
  134. {
  135. bigint biR;
  136. biR = bi_multiply( bi_copy( bi ), bi_copy( bi ) );
  137. bi_free( bi );
  138. return biR;
  139. }
  140. bigint
  141. bi_power( bigint bi, bigint biexp )
  142. {
  143. bigint biR;
  144. if ( bi_is_negative( bi_copy( biexp ) ) )
  145. {
  146. (void) fprintf( stderr, "bi_power: negative exponent\n" );
  147. (void) kill( getpid(), SIGFPE );
  148. }
  149. biR = bi_1;
  150. for (;;)
  151. {
  152. if ( bi_is_odd( bi_copy( biexp ) ) )
  153. biR = bi_multiply( biR, bi_copy( bi ) );
  154. biexp = bi_half( biexp );
  155. if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 )
  156. break;
  157. bi = bi_multiply( bi_copy( bi ), bi );
  158. }
  159. bi_free( bi );
  160. bi_free( biexp );
  161. return biR;
  162. }
  163. bigint
  164. bi_factorial( bigint bi )
  165. {
  166. bigint biR;
  167. biR = bi_1;
  168. while ( bi_compare( bi_copy( bi ), bi_1 ) > 0 )
  169. {
  170. biR = bi_multiply( biR, bi_copy( bi ) );
  171. bi = bi_int_subtract( bi, 1 );
  172. }
  173. bi_free( bi );
  174. return biR;
  175. }
  176. int
  177. bi_is_even( bigint bi )
  178. {
  179. return ! bi_is_odd( bi );
  180. }
  181. bigint
  182. bi_mod_power( bigint bi, bigint biexp, bigint bim )
  183. {
  184. int invert;
  185. bigint biR;
  186. invert = 0;
  187. if ( bi_is_negative( bi_copy( biexp ) ) )
  188. {
  189. biexp = bi_negate( biexp );
  190. invert = 1;
  191. }
  192. biR = bi_1;
  193. for (;;)
  194. {
  195. if ( bi_is_odd( bi_copy( biexp ) ) )
  196. biR = bi_mod( bi_multiply( biR, bi_copy( bi ) ), bi_copy( bim ) );
  197. biexp = bi_half( biexp );
  198. if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 )
  199. break;
  200. bi = bi_mod( bi_multiply( bi_copy( bi ), bi ), bi_copy( bim ) );
  201. }
  202. bi_free( bi );
  203. bi_free( biexp );
  204. if ( invert )
  205. biR = bi_mod_inverse( biR, bim );
  206. else
  207. bi_free( bim );
  208. return biR;
  209. }
  210. bigint
  211. bi_mod_inverse( bigint bi, bigint bim )
  212. {
  213. bigint gcd, mul0, mul1;
  214. gcd = bi_egcd( bi_copy( bim ), bi, &mul0, &mul1 );
  215. /* Did we get gcd == 1? */
  216. if ( ! bi_is_one( gcd ) )
  217. {
  218. (void) fprintf( stderr, "bi_mod_inverse: not relatively prime\n" );
  219. (void) kill( getpid(), SIGFPE );
  220. }
  221. bi_free( mul0 );
  222. return bi_mod( mul1, bim );
  223. }
  224. /* Euclid's algorithm. */
  225. bigint
  226. bi_gcd( bigint bim, bigint bin )
  227. {
  228. bigint bit;
  229. bim = bi_abs( bim );
  230. bin = bi_abs( bin );
  231. while ( ! bi_is_zero( bi_copy( bin ) ) )
  232. {
  233. bit = bi_mod( bim, bi_copy( bin ) );
  234. bim = bin;
  235. bin = bit;
  236. }
  237. bi_free( bin );
  238. return bim;
  239. }
  240. /* Extended Euclidean algorithm. */
  241. bigint
  242. bi_egcd( bigint bim, bigint bin, bigint* bim_mul, bigint* bin_mul )
  243. {
  244. bigint a0, b0, c0, a1, b1, c1, q, t;
  245. if ( bi_is_negative( bi_copy( bim ) ) )
  246. {
  247. bigint biR;
  248. biR = bi_egcd( bi_negate( bim ), bin, &t, bin_mul );
  249. *bim_mul = bi_negate( t );
  250. return biR;
  251. }
  252. if ( bi_is_negative( bi_copy( bin ) ) )
  253. {
  254. bigint biR;
  255. biR = bi_egcd( bim, bi_negate( bin ), bim_mul, &t );
  256. *bin_mul = bi_negate( t );
  257. return biR;
  258. }
  259. a0 = bi_1; b0 = bi_0; c0 = bim;
  260. a1 = bi_0; b1 = bi_1; c1 = bin;
  261. while ( ! bi_is_zero( bi_copy( c1 ) ) )
  262. {
  263. q = bi_divide( bi_copy( c0 ), bi_copy( c1 ) );
  264. t = a0;
  265. a0 = bi_copy( a1 );
  266. a1 = bi_subtract( t, bi_multiply( bi_copy( q ), a1 ) );
  267. t = b0;
  268. b0 = bi_copy( b1 );
  269. b1 = bi_subtract( t, bi_multiply( bi_copy( q ), b1 ) );
  270. t = c0;
  271. c0 = bi_copy( c1 );
  272. c1 = bi_subtract( t, bi_multiply( bi_copy( q ), c1 ) );
  273. bi_free( q );
  274. }
  275. bi_free( a1 );
  276. bi_free( b1 );
  277. bi_free( c1 );
  278. *bim_mul = a0;
  279. *bin_mul = b0;
  280. return c0;
  281. }
  282. bigint
  283. bi_lcm( bigint bia, bigint bib )
  284. {
  285. bigint biR;
  286. biR = bi_divide(
  287. bi_multiply( bi_copy( bia ), bi_copy( bib ) ),
  288. bi_gcd( bi_copy( bia ), bi_copy( bib ) ) );
  289. bi_free( bia );
  290. bi_free( bib );
  291. return biR;
  292. }
  293. /* The Jacobi symbol. */
  294. bigint
  295. bi_jacobi( bigint bia, bigint bib )
  296. {
  297. bigint biR;
  298. if ( bi_is_even( bi_copy( bib ) ) )
  299. {
  300. (void) fprintf( stderr, "bi_jacobi: don't know how to compute Jacobi(n, even)\n" );
  301. (void) kill( getpid(), SIGFPE );
  302. }
  303. if ( bi_compare( bi_copy( bia ), bi_copy( bib ) ) >= 0 )
  304. return bi_jacobi( bi_mod( bia, bi_copy( bib ) ), bib );
  305. if ( bi_is_zero( bi_copy( bia ) ) || bi_is_one( bi_copy( bia ) ) )
  306. {
  307. bi_free( bib );
  308. return bia;
  309. }
  310. if ( bi_compare( bi_copy( bia ), bi_2 ) == 0 )
  311. {
  312. bi_free( bia );
  313. switch ( bi_int_mod( bib, 8 ) )
  314. {
  315. case 1: case 7:
  316. return bi_1;
  317. case 3: case 5:
  318. return bi_m1;
  319. }
  320. }
  321. if ( bi_is_even( bi_copy( bia ) ) )
  322. {
  323. biR = bi_multiply(
  324. bi_jacobi( bi_2, bi_copy( bib ) ),
  325. bi_jacobi( bi_half( bia ), bi_copy( bib ) ) );
  326. bi_free( bib );
  327. return biR;
  328. }
  329. if ( bi_int_mod( bi_copy( bia ), 4 ) == 3 &&
  330. bi_int_mod( bi_copy( bib ), 4 ) == 3 )
  331. return bi_negate( bi_jacobi( bib, bia ) );
  332. else
  333. return bi_jacobi( bib, bia );
  334. }
  335. /* Probabalistic prime checking. */
  336. int
  337. bi_is_probable_prime( bigint bi, int certainty )
  338. {
  339. int i, p;
  340. bigint bim1;
  341. /* First do trial division by a list of small primes. This eliminates
  342. ** many candidates.
  343. */
  344. for ( i = 0; i < sizeof(low_primes)/sizeof(*low_primes); ++i )
  345. {
  346. p = low_primes[i];
  347. switch ( bi_compare( int_to_bi( p ), bi_copy( bi ) ) )
  348. {
  349. case 0:
  350. bi_free( bi );
  351. return 1;
  352. case 1:
  353. bi_free( bi );
  354. return 0;
  355. }
  356. if ( bi_int_mod( bi_copy( bi ), p ) == 0 )
  357. {
  358. bi_free( bi );
  359. return 0;
  360. }
  361. }
  362. /* Now do the probabilistic tests. */
  363. bim1 = bi_int_subtract( bi_copy( bi ), 1 );
  364. for ( i = 0; i < certainty; ++i )
  365. {
  366. bigint a, j, jac;
  367. /* Pick random test number. */
  368. a = bi_random( bi_copy( bi ) );
  369. /* Decide whether to run the Fermat test or the Solovay-Strassen
  370. ** test. The Fermat test is fast but lets some composite numbers
  371. ** through. Solovay-Strassen runs slower but is more certain.
  372. ** So the compromise here is we run the Fermat test a couple of
  373. ** times to quickly reject most composite numbers, and then do
  374. ** the rest of the iterations with Solovay-Strassen so nothing
  375. ** slips through.
  376. */
  377. if ( i < 2 && certainty >= 5 )
  378. {
  379. /* Fermat test. Note that this is not state of the art. There's a
  380. ** class of numbers called Carmichael numbers which are composite
  381. ** but look prime to this test - it lets them slip through no
  382. ** matter how many reps you run. However, it's nice and fast so
  383. ** we run it anyway to help quickly reject most of the composites.
  384. */
  385. if ( ! bi_is_one( bi_mod_power( bi_copy( a ), bi_copy( bim1 ), bi_copy( bi ) ) ) )
  386. {
  387. bi_free( bi );
  388. bi_free( bim1 );
  389. bi_free( a );
  390. return 0;
  391. }
  392. }
  393. else
  394. {
  395. /* GCD test. This rarely hits, but we need it for Solovay-Strassen. */
  396. if ( ! bi_is_one( bi_gcd( bi_copy( bi ), bi_copy( a ) ) ) )
  397. {
  398. bi_free( bi );
  399. bi_free( bim1 );
  400. bi_free( a );
  401. return 0;
  402. }
  403. /* Solovay-Strassen test. First compute pseudo Jacobi. */
  404. j = bi_mod_power(
  405. bi_copy( a ), bi_half( bi_copy( bim1 ) ), bi_copy( bi ) );
  406. if ( bi_compare( bi_copy( j ), bi_copy( bim1 ) ) == 0 )
  407. {
  408. bi_free( j );
  409. j = bi_m1;
  410. }
  411. /* Now compute real Jacobi. */
  412. jac = bi_jacobi( bi_copy( a ), bi_copy( bi ) );
  413. /* If they're not equal, the number is definitely composite. */
  414. if ( bi_compare( j, jac ) != 0 )
  415. {
  416. bi_free( bi );
  417. bi_free( bim1 );
  418. bi_free( a );
  419. return 0;
  420. }
  421. }
  422. bi_free( a );
  423. }
  424. bi_free( bim1 );
  425. bi_free( bi );
  426. return 1;
  427. }
  428. bigint
  429. bi_generate_prime( int bits, int certainty )
  430. {
  431. bigint bimo2, bip;
  432. int i, inc = 0;
  433. bimo2 = bi_power( bi_2, int_to_bi( bits - 1 ) );
  434. for (;;)
  435. {
  436. bip = bi_add( bi_random( bi_copy( bimo2 ) ), bi_copy( bimo2 ) );
  437. /* By shoving the candidate numbers up to the next highest multiple
  438. ** of six plus or minus one, we pre-eliminate all multiples of
  439. ** two and/or three.
  440. */
  441. switch ( bi_int_mod( bi_copy( bip ), 6 ) )
  442. {
  443. case 0: inc = 4; bip = bi_int_add( bip, 1 ); break;
  444. case 1: inc = 4; break;
  445. case 2: inc = 2; bip = bi_int_add( bip, 3 ); break;
  446. case 3: inc = 2; bip = bi_int_add( bip, 2 ); break;
  447. case 4: inc = 2; bip = bi_int_add( bip, 1 ); break;
  448. case 5: inc = 2; break;
  449. }
  450. /* Starting from the generated random number, check a bunch of
  451. ** numbers in sequence. This is just to avoid calls to bi_random(),
  452. ** which is more expensive than a simple add.
  453. */
  454. for ( i = 0; i < 1000; ++i ) /* arbitrary */
  455. {
  456. if ( bi_is_probable_prime( bi_copy( bip ), certainty ) )
  457. {
  458. bi_free( bimo2 );
  459. return bip;
  460. }
  461. bip = bi_int_add( bip, inc );
  462. inc = 6 - inc;
  463. }
  464. /* We ran through the whole sequence and didn't find a prime.
  465. ** Shrug, just try a different random starting point.
  466. */
  467. bi_free( bip );
  468. }
  469. }