zprime.c

Calc Software Package for Number Calc

/*
 * zprime - rapid small prime routines
 *
 * Copyright (C) 1999  Landon Curt Noll
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.3 $
 * @(#) $Id: zprime.c,v 29.3 2006/05/20 08:43:55 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/RCS/zprime.c,v $
 *
 * Under source code control:	1994/05/29 04:34:36
 * File existed as early as:	1994
 *
 * chongo <was here> /\oo/\	http://www.isthe.com/chongo/
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */


#include "zmath.h"
#include "prime.h"
#include "jump.h"
#include "config.h"
#include "zrand.h"
#include "have_const.h"


/*
 * When performing a probabilistic primality test, check to see
 * if the number has a factor <= PTEST_PRECHECK.
 *
 * XXX - what should this value be?  Perhaps this should be a function
 *	 of the size of the text value and the number of tests?
 */
#define PTEST_PRECHECK ((FULL)101)

/*
 * product of primes that fit into a long
 */
static CONST FULL pfact_tbl[MAX_PFACT_VAL+1] = {
    1, 1, 2, 6, 6, 30, 30, 210, 210, 210, 210, 2310, 2310, 30030, 30030,
    30030, 30030, 510510, 510510, 9699690, 9699690, 9699690, 9699690,
    223092870, 223092870, 223092870, 223092870, 223092870, 223092870
#if FULL_BITS == 64
    , U(6469693230), U(6469693230), U(200560490130), U(200560490130),
    U(200560490130), U(200560490130), U(200560490130), U(200560490130),
    U(7420738134810), U(7420738134810), U(7420738134810), U(7420738134810),
    U(304250263527210), U(304250263527210), U(13082761331670030),
    U(13082761331670030), U(13082761331670030), U(13082761331670030),
    U(614889782588491410), U(614889782588491410), U(614889782588491410),
    U(614889782588491410), U(614889782588491410), U(614889782588491410)
#endif
};

/*
 * determine the top 1 bit of a 8 bit value:
 *
 *	topbit[0] == 0 by convention
 *	topbit[x] gives the highest 1 bit of x
 */
static CONST unsigned char topbit[256] = {
    0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
    5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
    6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
    7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7
};

/*
 * integer square roots of powers of 2
 *
 * isqrt_pow2[x] == (int)(sqrt(2 to the x power))	      (for 0 <= x < 64)
 *
 * We have enough table entries for a FULL that is 64 bits long.
 */
static CONST FULL isqrt_pow2[64] = {
    1, 1, 2, 2, 4, 5, 8, 11,					/*  0 ..  7 */
    16, 22, 32, 45, 64, 90, 128, 181,				/*  8 .. 15 */
    256, 362, 512, 724, 1024, 1448, 2048, 2896,			/* 16 .. 23 */
    4096, 5792, 8192, 11585, 16384, 23170, 32768, 46340,	/* 24 .. 31 */
    65536, 92681, 131072, 185363,				/* 32 .. 35 */
    262144, 370727, 524288, 741455,				/* 36 .. 39 */
    1048576, 1482910, 2097152, 2965820,				/* 40 .. 43 */
    4194304, 5931641, 8388608, 11863283,			/* 44 .. 47 */
    16777216, 23726566, 33554432, 47453132,			/* 48 .. 51 */
    67108864, 94906265, 134217728, 189812531,			/* 52 .. 55 */
    268435456, 379625062, 536870912, 759250124,			/* 56 .. 59 */
    1073741824, 1518500249, 0x80000000, 0xb504f333		/* 60 .. 63 */
};

/*
 * static functions
 */
static FULL fsqrt(FULL v);		/* quick square root of v */
static long pix(FULL x);		/* pi of x */
static FULL small_factor(ZVALUE n, FULL limit);	  /* factor or 0 */


/*
 * Determine if a value is a small (32 bit) prime
 *
 * Returns:
 *	1	z is a prime <= MAX_SM_VAL
 *	0	z is not a prime <= MAX_SM_VAL
 *	-1	z > MAX_SM_VAL
 */
FLAG
zisprime(ZVALUE z)
{
	FULL n;				/* number to test */
	FULL isqr;			/* factor limit */
	CONST unsigned short *tp;	/* pointer to a prime factor */

	z.sign = 0;
	if (zisleone(z)) {
		return 0;
	}

	/* even numbers > 2 are not prime */
	if (ziseven(z)) {
		/*
		 * "2 is the greatest odd prime because it is the least even!"
		 *					- Dr. Dan Jurca	 1978
		 */
		return zisabstwo(z);
	}

	/* ignore non-small values */
	if (zge32b(z)) {
		return -1;
	}

	/* we now know that we are dealing with a value 0 <= n < 2^32 */
	n = ztofull(z);

	/* lookup small cases in pr_map */
	if (n <= MAX_MAP_VAL) {
		return (pr_map_bit(n) ? 1 : 0);
	}

	/* ignore Saber-C warning #530 about empty for statement */
	/*	  ok to ignore in proc zisprime */
	/* a number >=2^16 and < 2^32 */
	for (isqr=fsqrt(n), tp=prime; (*tp <= isqr) && (n % *tp); ++tp) {
	}
	return ((*tp <= isqr && *tp != 1) ? 0 : 1);
}


/*
 * Determine the next small (32 bit) prime > a 32 bit value.
 *
 * given:
 *	z		search point
 *
 * Returns:
 *	0	next prime is 2^32+15
 *	1	abs(z) >= 2^32
 *	smallest prime > abs(z) otherwise
 */
FULL
znprime(ZVALUE z)
{
	FULL n;			/* search point */

	z.sign = 0;

	/* ignore large values */
	if (zge32b(z)) {
		return (FULL)1;
	}

	/* deal a search point of 0 or 1 */
	if (zisabsleone(z)) {
		return (FULL)2;
	}

	/* deal with returning a value that is beyond our reach */
	n = ztofull(z);
	if (n >= MAX_SM_PRIME) {
		return (FULL)0;
	}

	/* return the next prime */
	return next_prime(n);
}


/*
 * Compute the next prime beyond a small (32 bit) value.
 *
 * This function assumes that 2 <= n < 2^32-5.
 *
 * given:
 *	n		search point
 */
FULL
next_prime(FULL n)
{
	CONST unsigned short *tp;	/* pointer to a prime factor */
	CONST unsigned char *j;		/* current jump increment */
	int tmp;

	/* find our search point */
	n = ((n & 0x1) ? n+2 : n+1);

	/* if we can just search the bit map, then search it */
	if (n <= MAX_MAP_PRIME) {

		/* search until we find a 1 bit */
		while (pr_map_bit(n) == 0) {
			n += (FULL)2;
		}

	/* too large for our table, find the next prime the hard way */
	} else {
		FULL isqr;		/* factor limit */

		/*
		 * Our search for a prime may cause us to increment n over
		 * a perfect square, but never two perfect squares.  The largest
		 * prime gap <= 2614941711251 is 651.  Shanks conjectures that
		 * the largest gap below P is about ln(P)^2.
		 *
		 * The value fsqrt(n)^2 will always be the perfect square
		 * that is <= n.  Given the smallness of prime gaps we will
		 * deal with, we know that n could carry us across the next
		 * perfect square (fsqrt(n)+1)^2 but not the following
		 * perfect square (fsqrt(n)+2)^2.
		 *
		 * Now the factor search limit for values < (fsqrt(n)+2)^2
		 * is the same limit for (fsqrt(n)+1)^2; namely fsqrt(n)+1.
		 * Therefore setting our limit at fsqrt(n)+1 and never
		 * bothering with it after that is safe.
		 */
		isqr = fsqrt(n)+1;

		/*
		 * If our factor limit is even, then we can reduce it to
		 * the next lowest odd value.  We already tested if n
		 * was even and all of our remaining potential factors
		 * are odd.
		 */
		if ((isqr & 0x1) == 0) {
			--isqr;
		}

		/*
		 * Skip to next value not divisible by a trivial prime.
		 */
		n = firstjmp(n, tmp);
		j = jmp + jmpptr(n);

		/*
		 * Look for tiny prime factors of increasing n until we
		 * find a prime.
		 */
		do {
			/* ignore Saber-C warning #530 - empty for statement */
			/*	  ok to ignore in proc next_prime */
			/* XXX - speed up test for large n by using gcds */
			/* find a factor, or give up if not found */
			for (tp=JPRIME; (*tp <= isqr) && (n % *tp); ++tp) {
			}
		} while (*tp <= isqr && *tp != 1 && (n += nxtjmp(j)));
	}

	/* return the prime that we found */
	return n;
}


/*
 * Determine the previous small (32 bit) prime < a 32 bit value
 *
 * given:
 *	z		search point
 *
 * Returns:
 *	1	abs(z) >= 2^32
 *	0	abs(z) <= 2
 *	greatest prime < abs(z) otherwise
 */
FULL
zpprime(ZVALUE z)
{
	CONST unsigned short *tp;	/* pointer to a prime factor */
	FULL isqr;			/* isqrt(z) */
	FULL n;				/* search point */
	CONST unsigned char *j;		/* current jump increment */
	int tmp;

	z.sign	= 0;

	/* ignore large values */
	if (zge32b(z)) {
		return (FULL)1;
	}

	/* deal with special case small values */
	n = ztofull(z);
	switch ((int)n) {
	case 0:
	case 1:
	case 2:
		/* ignore values <= 2 */
		return (FULL)0;
	case 3:
		/* 3 returns the only even prime */
		return (FULL)2;
	}

	/* deal with values above the bit map */
	if (n > NXT_MAP_PRIME) {

		/* find our search point */
		n = ((n & 0x1) ? n-2 : n-1);

		/* our factor limit - see next_prime for why this works */
		isqr = fsqrt(n)+1;
		if ((isqr & 0x1) == 0) {
			--isqr;
		}

		/*
		 * Skip to previous value not divisible by a trivial prime.
		 */
		tmp = jmpindxval(n);
		if (tmp >= 0) {

			/* find next value not divisible by a trivial prime */
			n += tmp;

			/* find the previous jump index */
			j = jmp + jmpptr(n);

			/* jump back */
			n -= prevjmp(j);

		/* already not divisible by a trivial prime */
		} else {
			/* find the current jump index */
			j = jmp + jmpptr(n);
		}

		/* factor values until we find a prime */
		do {
			/* ignore Saber-C warning #530 - empty for statement */
			/*	  ok to ignore in proc zpprime */
			/* XXX - speed up test for large n by using gcds */
			/* find a factor, or give up if not found */
			for (tp=prime; (*tp <= isqr) && (n % *tp); ++tp) {
			}
		} while (*tp <= isqr && *tp != 1 && (n -= prevjmp(j)));

	/* deal with values within the bit map */
	} else if (n <= MAX_MAP_PRIME) {

		/* find our search point */
		n = ((n & 0x1) ? n-2 : n-1);

		/* search until we find a 1 bit */
		while (pr_map_bit(n) == 0) {
			n -= (FULL)2;
		}

	/* deal with values that could cross into the bit map */
	} else {
		/* MAX_MAP_PRIME < n <= NXT_MAP_PRIME returns MAX_MAP_PRIME */
		return MAX_MAP_PRIME;
	}

	/* return what we found */
	return n;
}


/*
 * Compute the number of primes <= a ZVALUE that can fit into a FULL
 *
 * given:
 *	z		compute primes <= z
 *
 * Returns:
 *	-1	error
 *	>=0	number of primes <= x
 */
long
zpix(ZVALUE z)
{
	/* pi(<0) is always 0 */
	if (zisneg(z)) {
		return (long)0;
	}

	/* firewall */
	if (zge32b(z)) {
		return (long)-1;
	}
	return pix(ztofull(z));
}


/*
 * Compute the number of primes <= a ZVALUE
 *
 * given:
 *	x		value of z
 *
 * Returns:
 *	-1	error
 *	>=0	number of primes <= x
 */
static long
pix(FULL x)
{
	long count;			/* pi(x) */
	FULL top;			/* top of the range to test */
	CONST unsigned short *tp;	/* pointer to a tiny prime */
	FULL i;

	/* compute pi(x) using the 2^8 step table */
	if (x <= MAX_PI10B) {

		/* x within the prime table, so use it */
		if (x < MAX_MAP_PRIME) {
			/* firewall - pix(x) ==0 for x < 2 */
			if (x < 2) {
				count = 0;

			} else {
				/* determine how and where we will count */
				if (x < 1024) {
					count = 1;
					tp = prime;
				} else {
					count = pi10b[x>>10];
					tp = prime+count-1;
				}
				/* count primes in the table */
				while (*tp++ <= x) {
					++count;
				}
			}

		/* x is larger than the prime table, so count the hard way */
		} else {

			/* case: count down from pi18b entry to x */
			if (x & 0x200) {
				top = (x | 0x3ff);
				count = pi10b[(top+1)>>10];
				for (i=next_prime(x); i <= top;
				     i=next_prime(i)) {
					--count;
				}

			/* case: count up from pi10b entry to x */
			} else {
				count = pi10b[x>>10];
				for (i=next_prime(x&(~0x3ff));
				     i <= x; i = next_prime(i)) {
					++count;
				}
			}
		}

	/* compute pi(x) using the 2^18 interval table */
	} else {

		/* compute sum of intervals up to our interval */
		for (count=0, i=0; i < (x>>18); ++i) {
			count += pi18b[i];
		}

		/* case: count down from pi18b entry to x */
		if (x & 0x20000) {
			top = (x | 0x3ffff);
			count += pi18b[i];
			if (top > MAX_SM_PRIME) {
				if (x < MAX_SM_PRIME) {
					for (i=next_prime(x); i < MAX_SM_PRIME;
					     i=next_prime(i)) {
						--count;
					}
					--count;
				}
			} else {
				for (i=next_prime(x); i<=top; i=next_prime(i)) {
					--count;
				}
			}

		/* case: count up from pi18b entry to x */
		} else {
			for (i=next_prime(x&(~0x3ffff));
			     i <= x; i = next_prime(i)) {
				++count;
			}
		}
	}
	return count;
}


/*
 * Compute the smallest prime factor < limit
 *
 * given:
 *	n		number to factor
 *	zlimit		ending search point
 *	res		factor, if found, or NULL
 *
 * Returns:
 *	-1	error, limit >= 2^32
 *	0	no factor found, res is not changed
 *	1	factor found, res (if non-NULL) is smallest prime factor