/src/rt/bigint/bigint_ext.cpp

  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
 28#include <sys/types.h>
 29#include <signal.h>
 30#include <stdio.h>
 31#include <stdlib.h>
 32#include <unistd.h>
 33#include <time.h>
 34
 35#include "bigint.h"
 36#include "low_primes.h"
 37
 38
 39bigint bi_0, bi_1, bi_2, bi_10, bi_m1, bi_maxint, bi_minint;
 40
 41
 42/* Forwards. */
 43static void print_pos( FILE* f, bigint bi );
 44
 45
 46bigint
 47str_to_bi( char* str )
 48    {
 49    int sign;
 50    bigint biR;
 51
 52    sign = 1;
 53    if ( *str == '-' )
 54	{
 55	sign = -1;
 56	++str;
 57	}
 58    for ( biR = bi_0; *str >= '0' && *str <= '9'; ++str )
 59	biR = bi_int_add( bi_int_multiply( biR, 10 ), *str - '0' );
 60    if ( sign == -1 )
 61	biR = bi_negate( biR );
 62    return biR;
 63    }
 64
 65
 66void
 67bi_print( FILE* f, bigint bi )
 68    {
 69    if ( bi_is_negative( bi_copy( bi ) ) )
 70	{
 71	putc( '-', f );
 72	bi = bi_negate( bi );
 73	}
 74    print_pos( f, bi );
 75    }
 76
 77
 78bigint
 79bi_scan( FILE* f )
 80    {
 81    int sign;
 82    int c;
 83    bigint biR;
 84
 85    sign = 1;
 86    c = getc( f );
 87    if ( c == '-' )
 88	sign = -1;
 89    else
 90	ungetc( c, f );
 91
 92    biR = bi_0;
 93    for (;;)
 94	{
 95	c = getc( f );
 96	if ( c < '0' || c > '9' )
 97	    break;
 98	biR = bi_int_add( bi_int_multiply( biR, 10 ), c - '0' );
 99	}
100
101    if ( sign == -1 )
102	biR = bi_negate( biR );
103    return biR;
104    }
105
106
107static void
108print_pos( FILE* f, bigint bi )
109    {
110    if ( bi_compare( bi_copy( bi ), bi_10 ) >= 0 )
111	print_pos( f, bi_int_divide( bi_copy( bi ), 10 ) );
112    putc( bi_int_mod( bi, 10 ) + '0', f );
113    }
114
115
116int
117bi_int_mod( bigint bi, int m )
118    {
119    int r;
120
121    if ( m <= 0 )
122	{
123	(void) fprintf( stderr, "bi_int_mod: zero or negative modulus\n" );
124	(void) kill( getpid(), SIGFPE );
125	}
126    r = bi_int_rem( bi, m );
127    if ( r < 0 )
128	r += m;
129    return r;
130    }
131
132
133bigint
134bi_rem( bigint bia, bigint bim )
135    {
136    return bi_subtract(
137	bia, bi_multiply( bi_divide( bi_copy( bia ), bi_copy( bim ) ), bim ) );
138    }
139
140
141bigint
142bi_mod( bigint bia, bigint bim )
143    {
144    bigint biR;
145
146    if ( bi_compare( bi_copy( bim ), bi_0 ) <= 0 )
147	{
148	(void) fprintf( stderr, "bi_mod: zero or negative modulus\n" );
149	(void) kill( getpid(), SIGFPE );
150	}
151    biR = bi_rem( bia, bi_copy( bim ) );
152    if ( bi_is_negative( bi_copy( biR ) ) )
153	biR = bi_add( biR, bim );
154    else
155	bi_free( bim );
156    return biR;
157    }
158
159
160bigint
161bi_square( bigint bi )
162    {
163    bigint biR;
164
165    biR = bi_multiply( bi_copy( bi ), bi_copy( bi ) );
166    bi_free( bi );
167    return biR;
168    }
169
170
171bigint
172bi_power( bigint bi, bigint biexp )
173    {
174    bigint biR;
175
176    if ( bi_is_negative( bi_copy( biexp ) ) )
177	{
178	(void) fprintf( stderr, "bi_power: negative exponent\n" );
179	(void) kill( getpid(), SIGFPE );
180	}
181    biR = bi_1;
182    for (;;)
183	{
184	if ( bi_is_odd( bi_copy( biexp ) ) )
185	    biR = bi_multiply( biR, bi_copy( bi ) );
186	biexp = bi_half( biexp );
187	if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 )
188	    break;
189	bi = bi_multiply( bi_copy( bi ), bi );
190	}
191    bi_free( bi );
192    bi_free( biexp );
193    return biR;
194    }
195
196
197bigint
198bi_factorial( bigint bi )
199    {
200    bigint biR;
201
202    biR = bi_1;
203    while ( bi_compare( bi_copy( bi ), bi_1 ) > 0 )
204	{
205	biR = bi_multiply( biR, bi_copy( bi ) );
206	bi = bi_int_subtract( bi, 1 );
207	}
208    bi_free( bi );
209    return biR;
210    }
211
212
213int
214bi_is_even( bigint bi )
215    {
216    return ! bi_is_odd( bi );
217    }
218
219
220bigint
221bi_mod_power( bigint bi, bigint biexp, bigint bim )
222    {
223    int invert;
224    bigint biR;
225
226    invert = 0;
227    if ( bi_is_negative( bi_copy( biexp ) ) )
228	{
229	biexp = bi_negate( biexp );
230	invert = 1;
231	}
232
233    biR = bi_1;
234    for (;;)
235	{
236	if ( bi_is_odd( bi_copy( biexp ) ) )
237	    biR = bi_mod( bi_multiply( biR, bi_copy( bi ) ), bi_copy( bim ) );
238	biexp = bi_half( biexp );
239	if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 )
240	    break;
241	bi = bi_mod( bi_multiply( bi_copy( bi ), bi ), bi_copy( bim ) );
242	}
243    bi_free( bi );
244    bi_free( biexp );
245
246    if ( invert )
247	biR = bi_mod_inverse( biR, bim );
248    else
249	bi_free( bim );
250    return biR;
251    }
252
253
254bigint
255bi_mod_inverse( bigint bi, bigint bim )
256    {
257    bigint gcd, mul0, mul1;
258
259    gcd = bi_egcd( bi_copy( bim ), bi, &mul0, &mul1 );
260
261    /* Did we get gcd == 1? */
262    if ( ! bi_is_one( gcd ) )
263	{
264	(void) fprintf( stderr, "bi_mod_inverse: not relatively prime\n" );
265	(void) kill( getpid(), SIGFPE );
266	}
267
268    bi_free( mul0 );
269    return bi_mod( mul1, bim );
270    }
271
272
273/* Euclid's algorithm. */
274bigint
275bi_gcd( bigint bim, bigint bin )
276    {
277    bigint bit;
278
279    bim = bi_abs( bim );
280    bin = bi_abs( bin );
281    while ( ! bi_is_zero( bi_copy( bin ) ) )
282	{
283	bit = bi_mod( bim, bi_copy( bin ) );
284	bim = bin;
285	bin = bit;
286	}
287    bi_free( bin );
288    return bim;
289    }
290
291
292/* Extended Euclidean algorithm. */
293bigint
294bi_egcd( bigint bim, bigint bin, bigint* bim_mul, bigint* bin_mul )
295    {
296    bigint a0, b0, c0, a1, b1, c1, q, t;
297
298    if ( bi_is_negative( bi_copy( bim ) ) )
299	{
300	bigint biR;
301
302	biR = bi_egcd( bi_negate( bim ), bin, &t, bin_mul );
303	*bim_mul = bi_negate( t );
304	return biR;
305	}
306    if ( bi_is_negative( bi_copy( bin ) ) )
307	{
308	bigint biR;
309
310	biR = bi_egcd( bim, bi_negate( bin ), bim_mul, &t );
311	*bin_mul = bi_negate( t );
312	return biR;
313	}
314
315    a0 = bi_1;  b0 = bi_0;  c0 = bim;
316    a1 = bi_0;  b1 = bi_1;  c1 = bin;
317
318    while ( ! bi_is_zero( bi_copy( c1 ) ) )
319	{
320	q = bi_divide( bi_copy( c0 ), bi_copy( c1 ) );
321	t = a0;
322	a0 = bi_copy( a1 );
323	a1 = bi_subtract( t, bi_multiply( bi_copy( q ), a1 ) );
324	t = b0;
325	b0 = bi_copy( b1 );
326	b1 = bi_subtract( t, bi_multiply( bi_copy( q ), b1 ) );
327	t = c0;
328	c0 = bi_copy( c1 );
329	c1 = bi_subtract( t, bi_multiply( bi_copy( q ), c1 ) );
330	bi_free( q );
331	}
332
333    bi_free( a1 );
334    bi_free( b1 );
335    bi_free( c1 );
336    *bim_mul = a0;
337    *bin_mul = b0;
338    return c0;
339    }
340
341
342bigint
343bi_lcm( bigint bia, bigint bib )
344    {
345    bigint biR;
346
347    biR = bi_divide(
348	bi_multiply( bi_copy( bia ), bi_copy( bib ) ),
349	bi_gcd( bi_copy( bia ), bi_copy( bib ) ) );
350    bi_free( bia );
351    bi_free( bib );
352    return biR;
353    }
354
355
356/* The Jacobi symbol. */
357bigint
358bi_jacobi( bigint bia, bigint bib )
359    {
360    bigint biR;
361
362    if ( bi_is_even( bi_copy( bib ) ) )
363	{
364	(void) fprintf( stderr, "bi_jacobi: don't know how to compute Jacobi(n, even)\n" );
365	(void) kill( getpid(), SIGFPE );
366	}
367
368    if ( bi_compare( bi_copy( bia ), bi_copy( bib ) ) >= 0 )
369	return bi_jacobi( bi_mod( bia, bi_copy( bib ) ), bib );
370
371    if ( bi_is_zero( bi_copy( bia ) ) || bi_is_one( bi_copy( bia ) ) )
372	{
373	bi_free( bib );
374	return bia;
375	}
376
377    if ( bi_compare( bi_copy( bia ), bi_2 ) == 0 )
378	{
379	bi_free( bia );
380	switch ( bi_int_mod( bib, 8 ) )
381	    {
382	    case 1: case 7:
383	    return bi_1;
384	    case 3: case 5:
385	    return bi_m1;
386	    }
387	}
388
389    if ( bi_is_even( bi_copy( bia ) ) )
390	{
391	biR = bi_multiply(
392	    bi_jacobi( bi_2, bi_copy( bib ) ),
393	    bi_jacobi( bi_half( bia ), bi_copy( bib ) ) );
394	bi_free( bib );
395	return biR;
396	}
397
398    if ( bi_int_mod( bi_copy( bia ), 4 ) == 3 &&
399         bi_int_mod( bi_copy( bib ), 4 ) == 3 )
400	return bi_negate( bi_jacobi( bib, bia ) );
401    else
402	return bi_jacobi( bib, bia );
403    }
404
405
406/* Probabalistic prime checking. */
407int
408bi_is_probable_prime( bigint bi, int certainty )
409    {
410    int i, p;
411    bigint bim1;
412
413    /* First do trial division by a list of small primes.  This eliminates
414    ** many candidates.
415    */
416    for ( i = 0; i < sizeof(low_primes)/sizeof(*low_primes); ++i )
417	{
418	p = low_primes[i];
419	switch ( bi_compare( int_to_bi( p ), bi_copy( bi ) ) )
420	    {
421	    case 0:
422	    bi_free( bi );
423	    return 1;
424	    case 1:
425	    bi_free( bi );
426	    return 0;
427	    }
428	if ( bi_int_mod( bi_copy( bi ), p ) == 0 )
429	    {
430	    bi_free( bi );
431	    return 0;
432	    }
433	}
434
435    /* Now do the probabilistic tests. */
436    bim1 = bi_int_subtract( bi_copy( bi ), 1 );
437    for ( i = 0; i < certainty; ++i )
438	{
439	bigint a, j, jac;
440
441	/* Pick random test number. */
442	a = bi_random( bi_copy( bi ) );
443
444	/* Decide whether to run the Fermat test or the Solovay-Strassen
445	** test.  The Fermat test is fast but lets some composite numbers
446	** through.  Solovay-Strassen runs slower but is more certain.
447	** So the compromise here is we run the Fermat test a couple of
448	** times to quickly reject most composite numbers, and then do
449	** the rest of the iterations with Solovay-Strassen so nothing
450	** slips through.
451	*/
452	if ( i < 2 && certainty >= 5 )
453	    {
454	    /* Fermat test.  Note that this is not state of the art.  There's a
455	    ** class of numbers called Carmichael numbers which are composite
456	    ** but look prime to this test - it lets them slip through no
457	    ** matter how many reps you run.  However, it's nice and fast so
458	    ** we run it anyway to help quickly reject most of the composites.
459	    */
460	    if ( ! bi_is_one( bi_mod_power( bi_copy( a ), bi_copy( bim1 ), bi_copy( bi ) ) ) )
461		{
462		bi_free( bi );
463		bi_free( bim1 );
464		bi_free( a );
465		return 0;
466		}
467	    }
468	else
469	    {
470	    /* GCD test.  This rarely hits, but we need it for Solovay-Strassen. */
471	    if ( ! bi_is_one( bi_gcd( bi_copy( bi ), bi_copy( a ) ) ) )
472		{
473		bi_free( bi );
474		bi_free( bim1 );
475		bi_free( a );
476		return 0;
477		}
478
479	    /* Solovay-Strassen test.  First compute pseudo Jacobi. */
480	    j = bi_mod_power(
481		    bi_copy( a ), bi_half( bi_copy( bim1 ) ), bi_copy( bi ) );
482	    if ( bi_compare( bi_copy( j ), bi_copy( bim1 ) ) == 0 )
483		{
484		bi_free( j );
485		j = bi_m1;
486		}
487
488	    /* Now compute real Jacobi. */
489	    jac = bi_jacobi( bi_copy( a ), bi_copy( bi ) );
490
491	    /* If they're not equal, the number is definitely composite. */
492	    if ( bi_compare( j, jac ) != 0 )
493		{
494		bi_free( bi );
495		bi_free( bim1 );
496		bi_free( a );
497		return 0;
498		}
499	    }
500
501	bi_free( a );
502	}
503
504    bi_free( bim1 );
505
506    bi_free( bi );
507    return 1;
508    }
509
510
511bigint
512bi_generate_prime( int bits, int certainty )
513    {
514    bigint bimo2, bip;
515    int i, inc = 0;
516
517    bimo2 = bi_power( bi_2, int_to_bi( bits - 1 ) );
518    for (;;)
519	{
520	bip = bi_add( bi_random( bi_copy( bimo2 ) ), bi_copy( bimo2 ) );
521	/* By shoving the candidate numbers up to the next highest multiple
522	** of six plus or minus one, we pre-eliminate all multiples of
523	** two and/or three.
524	*/
525	switch ( bi_int_mod( bi_copy( bip ), 6 ) )
526	    {
527	    case 0: inc = 4; bip = bi_int_add( bip, 1 ); break;
528	    case 1: inc = 4;                             break;
529	    case 2: inc = 2; bip = bi_int_add( bip, 3 ); break;
530	    case 3: inc = 2; bip = bi_int_add( bip, 2 ); break;
531	    case 4: inc = 2; bip = bi_int_add( bip, 1 ); break;
532	    case 5: inc = 2;                             break;
533	    }
534	/* Starting from the generated random number, check a bunch of
535	** numbers in sequence.  This is just to avoid calls to bi_random(),
536	** which is more expensive than a simple add.
537	*/
538	for ( i = 0; i < 1000; ++i )	/* arbitrary */
539	    {
540	    if ( bi_is_probable_prime( bi_copy( bip ), certainty ) )
541		{
542		bi_free( bimo2 );
543		return bip;
544		}
545	    bip = bi_int_add( bip, inc );
546	    inc = 6 - inc;
547	    }
548	/* We ran through the whole sequence and didn't find a prime.
549	** Shrug, just try a different random starting point.
550	*/
551	bi_free( bip );
552	}
553    }