/* @(#) $Revision: 1.11 $ */
/* @(#) RCS control in //prime.corp/usr/local/src/cmd/prime/src/map1.c */
/*
 * map1 - sieve numbers of the form k*2^n+/-1 for factors 2^31 < q < 2^39
 *
 * usage:
 *	map1 setname
 *
 *	setname		base state filename
 *
 * The state of the sieve comes in several files:
 *
 *	setname.parm	version, low_k, k_len, n global parameters
 *	setname.neg	array of bits bit x represents (low_k+x)*2^n-1
 *			    1 => composite, 0 => unknown
 *	setname.pos	array of bits bit x represents (low_k+x)*2^n+1
 *			    1 => composite, 0 => unknown
 *	setname.ckpt	q checkpoint value
 *	map1		binary to complete sieve for 2^31 < q < 2^39
 *
 * These files were prepaired by init_sieve and mksieve.
 *
 * Given a set of 't' values:
 *
 *	    ** NOTE: we will use ^ to mean exponentiation **
 *
 *	{low_k*2^n-1, (low_k+1)*2^n-1, ... (low_k+k_len-1)*2^n-1}
 *	{low_k*2^n+1, (low_k+1)*2^n+1, ... (low_k+k_len-1)*2^n+1}
 *
 * We will determine the values that are multiples of odd 2^31 < q < 2^39
 */
/*
 *
 * Copyright (C) 1997 Landon Curt Noll, all rights reserved.
 *
 * Permission to use, copy, modify, and distribute this software and
 * its documentation for any purpose is hereby granted, provided that
 * the above copyright, this permission notice, and the disclaimer
 * below appear in all of the following:
 *
 *         * supporting documentation
 *         * source copies
 *         * source works derived from this source
 *
 * THE COPYRIGHT HOLDERS DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO
 * EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY SPECIAL, INDIRECT OR
 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
 * USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
 * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 * PERFORMANCE OF THIS SOFTWARE.
 *
 * Landon Curt Noll <chongo was here> /\oo/\
 *
 * chongo@{toad,sgi}.com		Share and Enjoy!
 */


#include <stdio.h>
#include <sys/types.h>
#include <sys/file.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <signal.h>
#include <fcntl.h>
#include <errno.h>
#define NO_EXTERNS		/* we don't need the prime_diff[] array */
#include "primetbl.h"
#include "jump.h"
#include "nbits.h"


/* major types */
typedef char w8;
typedef unsigned char uw8;
typedef short w16;
typedef unsigned short uw16;
#if !defined(_MIPS_SZLONG) || _MIPS_SZLONG == 32
typedef long w32;
typedef unsigned long uw32;
typedef long long w64;
typedef unsigned long long uw64;
#else
typedef int w32;
typedef unsigned int uw32;
typedef long w64;
typedef unsigned long uw64;
#endif
typedef union {
    uw64 b64;
    uw32 b32[2];	/* b32[0] is most significant */
    uw16 b16[4];	/* b16[0] is most significant */
    uw8 b8[8];		/*  b8[0] is most significant */
} large;
typedef union {
    uw64 b64[2];
    uw32 b32[4];	/* b32[0] is most significant */
    uw16 b16[8];	/* b16[0] is most significant */
    uw8 b8[16];		/*  b8[0] is most significant */
} vlarge;

/*
 * B64_OFF32 - 64 bits offset 32 bits within a 128 bit quad
 *
 * B64_OFF16 - 64 bits offset multiple of 16 bits within a 128 bit quad
 * B32_OFF16 - 32 bits offset multiple of 16 bits within a 128 bit quad
 *
 * B64_OFF8 - 64 bits offset multiple of 8 bits within a 128 bit quad
 * B32_OFF8 - 32 bits offset multiple of 8 bits within a 128 bit quad
 * B16_OFF8 - 16 bits offset multiple of 8 bits within a 128 bit quad
 */
#define B64_OFF32(x)   (*((uw64 *)&((x).b32[1])))

#define B64_OFF16(x,y) (*((uw64 *)&((x).b16[y])))
#define B32_OFF16(x,y) (*((uw32 *)&((x).b16[y])))

#define B64_OFF8(x,y) (*((uw64 *)&((x).b8[y])))
#define B32_OFF8(x,y) (*((uw32 *)&((x).b8[y])))
#define B16_OFF8(x,y) (*((uw16 *)&((x).b8[y])))

/*
 * load a vlarge with a unsigned long long
 *
 * LOAD(qd,val) performs:
 *
 *	vlarge vlrg;
 *	w64 val;
 *
 *	vlrg = val;
 */
#define LOAD(vlrg,val) ((vlrg).b64[0]=0LL, (vlrg).b64[1]=(uw64)(val))

/*
 * compute x % q in higher precision
 *
 * MOD(res,x,q) performs:
 *
 *	large res;
 *	vlarge x;
 *	w64 q;
 *
 *	res = x % q
 *
 * In in higher precision (up to 2^48), based on the sqrmod.c code:
 *
 *	(*
 *	 * mode 1:    2^31 <= q < 2^39  x2 < 2^79
 *	 *
 *	 * x % q = (((x >> 16) % q)<<16 | (x % 2^16)) % q
 *	 *)
 *	res.b64 = ((B64_OFF16(x,3) % q)<<16 | (uw64)(x.b16[7])) % q;
 */
#define MOD(res,x,q) \
     ((res).b64 = ((B64_OFF16((x),3) % (q))<<16 | (uw64)(x).b16[7]) % (q))

/*
 * compute x^2 in higher precision
 *
 * SQUARE(x,x2) performs:
 *
 *	large x;
 *	vlarge x2;
 *
 *	x2 = x * x;
 *
 * In in higher precision (up to 2^48), based on the sqrmod.c code:
 *
 *	(*
 *	 * form x^2
 *	 *)
 *	(* form the upper partial product *)
 *	x2.b64[1] = (uw64)(x.b32[1]) * (uw64)(x.b32[1]);
 *	(* form the lower partial product *)
 *	x2.b32[1] = (uw32)(x.b16[1]) * (uw32)(x.b16[1]);
 *	(* add the middle partial product *)
 *	B64_OFF32(x2) += (((uw64)x.b16[1] * (uw64)x.b32[1]) << 1);
 */
#define SQUARE(x2,x) \
    (((x2).b64[1] = (uw64)((x).b32[1]) * (uw64)((x).b32[1])), \
     ((x2).b32[1] = (uw32)((x).b16[1]) * (uw32)((x).b16[1])),\
     (B64_OFF32(x2) += (((uw64)(x).b16[1] * (uw64)(x).b32[1]) << 1)))

/*
 * compute x<<=1 in higher precision
 *
 * DOUBLE(x2) performs:
 *
 *	vlarge x2;
 *
 *	x2 <<= 1;
 *
 * in in higher precision (up to 2^63), based on the sqrmod.c code:
 *
 *	(*
 *	 * form 2*x
 *	 *)
 *	x2.b64[0] = ((x2.b64[0] << 1) | (uw64)(x2.b32[2] >> 31));
 *	x2.b64[1] <<= 1;
 */
#define DOUBLE(x2) \
    (((x2).b64[0] = (((x2).b64[0] << 1) | (uw64)((x2).b32[2] >> 31))), \
     ((x2).b64[1] <<= 1))

/*
 * store a vlarge into a unsigned long long
 *
 * STORE(val,qd) performs:
 *
 *	large lrg;
 *	w64 val;
 *
 *	val = lrg.b64[1];
 */
#define STORE(val,lrg) ((val) = (lrg).b64)

/*
 * NBIT_IS_{0,1} - inline processing to compute n^2 mod q
 *
 * These two macros perform the steps needed to compute n^2 mod q.
 * We use the fact that:
 *
 *	x = y mod q  ==>  (x^2 mod q) = y^2 mod q
 *	x = y mod q  ==>  (2*x^2 mod q) = 2*(y^2) mod q
 *
 * If we have the bits of n, then we can use the above two rules
 * to either square mod q (if the n bit is 0) or 2*square mod q
 * (if the n bit is 1).
 *
 * When this source was built by mksieve, these two macros were placed
 * in main inline such that they correspond to the bit values of n.
 *
 * The NBIT_IS_0 macro was adopted from the map0.c macro:
 *
 *	#define NBIT_IS_0(mod,unused,q) ((mod)%=(q), (mod)*=(mod))
 *
 * The macro NBIT_IS_0(x,tmp,q) performs in higher precision:
 *
 *	vlarge x;
 *	large tmp;		(temp storage holder)
 *	w64 q;
 *
 *	x = x*x % q;
 *
 * The NBIT_IS_1 macro was adopted from the map0.c macro:
 *
 *	#define NBIT_IS_1(mod,unused,q) ((mod)%=(q), (mod)*=(mod), (mod)<<=1)
 *
 * The macro NBIT_IS_1(x,tmp,q) performs in higher precision:
 *
 *	vlarge x;
 *	large tmp;		(temp storage holder)
 *	w64 q;
 *
 *	x = 2*x*x % q;
 */
#define NBIT_IS_0(x,tmp,q) (MOD(tmp,x,q), SQUARE(x,tmp))
#define NBIT_IS_1(x,tmp,q) (MOD(tmp,x,q), SQUARE(x,tmp), DOUBLE(x))


/*
 * SET_BIT(bytes,x) - set the x-th bit of a byte array
 * CHK_BIT(bytes,x) - return non-zero if x-th bit of a byte array is set
 */
#define SET_BIT(map, x) (map)[(x)>>3] |= (1<<((x) & 0x7))
#define CHK_BIT(map, x) ((map)[(x)>>3] & (1<<((x) & 0x7)))


#define VERSION "1.1"		/* format and filenames version number */

#define CHECKPT (60*60)		/* seconds between checkpoints */

/*
 * Q_LIMIT - maximum factor of k*2^n+/-1 for we will test
 */
#define Q_LIMIT ((uw64)((1ULL<<39)-1ULL))


/*
 * bit map of k*2^n-1 and k*2^n+1 values
 *
 *	kmap_neg[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n-1 is composite
 *	kmap_neg[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n-1 is unknown
 *
 *	kmap_pos[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n+1 is composite
 *	kmap_pos[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n+1 is unknown
 */
char *kmap_neg;
char *kmap_pos;

long k_len;		/* number of k to test */
w64 external_q;		/* globally acessable value of q that was just tested */
char *setname;		/* name of state set */
FILE *ckpt = NULL;	/* setname.ckpt open file stream */

char *program;		/* our name */


/*
 * forward declarations
 */
int lock_ckpt(char*);
void checkpt(void);
int write_kmap(char*, char*, int);
int parse_kmap(char*, char**, int);
int parse_parm(char*, w64*, long*);
int last_q(char*, w64*);


main(argc, argv)
    int argc;		/* arg count */
    char *argv[];	/* the args */
{
    w64 low_k;		/* initial multiplier as in low_k*2^n+/-1 */
    long n;		/* exponent n as in k*2^n+/-1 */
    w64 q;		/* current test factor */
    uw64 mval;		/* modular inverse of 2^n % q or 2^n % q or partial */
    vlarge x;		/* temp strage for x^2%q or 2*x^2%q partial */
    large tmp;		/* temp strage for x^2%q or 2*x^2%q partial */
    w64 k;		/* multiplier as in low_k*2^n+/-1 */
    unsigned char *jp;	/* jump pointer */
    int ckpt_fd;	/* setname.ckpt file descriptor */
    int i;

    /*
     * parse args
     */
    program = argv[0];
    if (argc != 2) {
	fprintf(stderr, "usage: %s setname\n", program);
	exit(1);
    }
    setname = argv[1];

    /*
     * lock the checkpoint file
     */
    ckpt_fd = lock_ckpt(setname);
    if (ckpt_fd < 0) {
	fprintf(stderr, "%s: %s.ckpt already locked\n", program, setname);
	exit(0);
    }
    ckpt = fdopen(ckpt_fd, "r+");

    /*
     * parse the name.parm file
     */
    i = parse_parm(setname, &low_k, &n);
    if (i < 0) {
	fprintf(stderr, "%s: bad parse of %s.parm, return: %d\n",
	  program, setname, i);
	exit(2);
    }

    /*
     * read name.neg kmap file
     */
    i = parse_kmap(setname, &kmap_neg, -1);
    if (i < 0) {
	fprintf(stderr, "%s: bad parse of %s.neg, return: %d\n",
	  program, setname, i);
	exit(3);
    }

    /*
     * read name.pos kmap file
     */
    i = parse_kmap(setname, &kmap_pos, +1);
    if (i < 0) {
	fprintf(stderr, "%s: bad parse of %s.pos, return: %d\n",
	  program, setname, i);
	exit(4);
    }

    /*
     * determine the next q value to test
     */
    i = last_q(setname, &external_q);
    if (i < 0) {
	fprintf(stderr, "%s: bad parse of %s.ckpt, return: %d\n",
	  program, setname, i);
	exit(5);
    }
    if (external_q % 0x2 == 0) {
	++external_q;
    }
    external_q += (w64)2;
    /* find the next test factor not divisable by a small prime */
    q = (w64)firstjmp(external_q, i); 
    external_q = q;
    /* verify that we have something to do */
    if (external_q > Q_LIMIT) {
	fprintf(stderr, "%s: q:%lld > %lld, nothing to do\n",
	  program, q, Q_LIMIT);
	exit(0);
    }

    /*
     * enable checkpoint alarms
     */
    if (sigset(SIGALRM, checkpt) == SIG_ERR) {
	fprintf(stderr, "%s: setset(SIGALRM, checkpt) failed\n", program);
	perror("sigset");
	exit(6);
    }
    {  
	struct itimerval checktime = {{CHECKPT, 0}, {CHECKPT, 0}};

	if (setitimer(ITIMER_REAL, &checktime, NULL) < 0) {
	    fprintf(stderr, "%s: setitimer failed\n", program);
	    perror("setitimer");
	    exit(7);
	}
    }

    /*
     * eliminate values disivble by primes NEXT_PRIME_DIFF <= q < 2^31
     */
    for(jp = jmp + jmpptr(q); q <= Q_LIMIT; q += nxtjmp(jp)) {

	/*
	 * compute 2^n mod q 
	 *
	 * We use the fact that:
	 *	x = y mod q  ==>  (x^2 mod q) = y^2 mod q
	 *	x = y mod q  ==>  (2*x^2 mod q) = 2*(y^2) mod q
	 *
	 * If we have the bits of n, then we can use the above two rules
	 * to either square mod q (if the n bit is 0) or 2*square mod q
	 * (if the n bit is 1).  We skip the first 5 steps because we
	 * have the top 5 bits of n processed as a power of 2 in top5.
	 *
	 * The loop below works because we know that q < 2^32 (in fact
	 * q is < NXT_MAP_PRIME).  If this were not so, we would need to 
	 * perform double precision mods and multiplications.
	 */
	LOAD(x, 1ULL<<(TOP_5_N_BITS));  /* bits 0 to 4 */
	N_BIT_5(x,tmp,q);
	N_BIT_6(x,tmp,q);
	N_BIT_7(x,tmp,q);
	N_BIT_8(x,tmp,q);
	N_BIT_9(x,tmp,q);
	N_BIT_10(x,tmp,q);
	N_BIT_11(x,tmp,q);
	N_BIT_12(x,tmp,q);
	N_BIT_13(x,tmp,q);
	N_BIT_14(x,tmp,q);
	N_BIT_15(x,tmp,q);
	N_BIT_16(x,tmp,q);
	N_BIT_17(x,tmp,q);
	N_BIT_18(x,tmp,q);
	N_BIT_19(x,tmp,q);
	N_BIT_20(x,tmp,q);
	N_BIT_21(x,tmp,q);
	N_BIT_22(x,tmp,q);
	N_BIT_23(x,tmp,q);
	N_BIT_24(x,tmp,q);
	N_BIT_25(x,tmp,q);
	N_BIT_26(x,tmp,q);
	N_BIT_27(x,tmp,q);
	N_BIT_28(x,tmp,q);
	N_BIT_29(x,tmp,q);
	N_BIT_30(x,tmp,q);
	N_BIT_31(x,tmp,q);
	MOD(tmp,x,q);
	STORE(mval,tmp);

	/*
	 * compute the modular inverse of 2^n mod q
	 *
	 * Calculate modular inverse suppressing unnecessary divisions.
	 * This is based on the Euclidian algorithm for large numbers.
	 * (Algorithm X from Knuth Vol 2, section 4.5.2. and exercise 17)
	 * Returns TRUE if there is no solution because the numbers
	 * are not relatively prime.
	 *
	 * Code cloned from calc v2.10.0t1
	 *
	 * Below is an inline equivalent:
	 *
	 *	mval = minv(mval, q);
	 *
	 * where minv() is the function found in init_sieve.
	 */
	{
	    w64 u2,v2,u3,v3;	/* Eclidian terms */
	    w64 q1, q2;		/* intermediate quotents */
	    w64 tmp1, tmp2;	/* temp values */

	    /* prep work */
	    v3 = (w64)mval;
	    u3 = q;
	    u2 = 0;
	    v2 = 1;

	    /* 
	     * We know that the easy no-solution case of (v3 == 0 && u3 != 1)
	     * cannot happen because u3 == q and q <= NEXT_PRIME_DIFF.
	     */

	    /* Eclidian process */
	    while (v3 != 0) {
		q1 = u3 / v3;
		tmp1 = v2 * q1;
		tmp2 = u2 - tmp1;
		u2 = v2;
		v2 = tmp2;
		q2 = u3 - (q1*v3);
		u3 = v3;
		v3 = q2;
	    }

	    /*
	     * We know that there will always be a modular inverse
	     * so we can skip the test of (u3 != 1).
	     */

	    /*
	     * Set mval to the multiplicitive inverse
	     */
	    mval = ((u2 < 0) ? (uw64)(q+u2) : (uw64)u2);
	}

	/*
	 * Sieve k*2^n-1
	 *
	 * Number Theory magic:
	 *
	 *	Let m be the modular inverse of 2^n mod q (m*2^n % q == 1).
	 *	Then for odd q > 1, q divides k*2^n-1 if and only if m = k%q.
	 *
	 * We will sieve out all k that are == mval % q.
	 */
	k = (w64)((w64)mval - (low_k % q));
	if (k < 0) {
	    if (k+q < k_len) {
#if defined(DEBUG)
		if (CHK_BIT(kmap_neg, (w32)(k+q)) == 0) {
		    printf("%lld %ld %lld -1\n", q, n, low_k+k+q);
		}
#endif
		SET_BIT(kmap_neg, (long)(k+q));
	    }
	} else {
	    if (k < k_len) {
#if defined(DEBUG)
		if (CHK_BIT(kmap_neg, (w32)k) == 0) {
		    printf("%lld %ld %lld  -1\n", q, n, low_k+k);
		}
#endif
		SET_BIT(kmap_neg, (long)k);
	    }
	}

	/*
	 * Sieve k*2^n+1
	 *
	 * Number Theory magic:
	 *
	 *	Let m be the modular inverse of 2^n mod q (m*2^n % q == 1).
	 *	Then for odd q > 1, q divides k*2^n+1 if and only if q-m = k%q.
	 *
	 * We will sieve out all k that are == (q-mval) % q.
	 */
	k = (w64)(q - (w64)mval - (low_k % q));
	if (k < 0) {
	    if (k+q < k_len) {
#if defined(DEBUG)
		if (CHK_BIT(kmap_pos, (w32)(k+q)) == 0) {
		    printf("%lld %ld %lld +1\n", q, n, low_k+k+q);
		}
#endif
		SET_BIT(kmap_pos, (long)(k+q));
	    }
	} else {
	    if (k < k_len) {
#if defined(DEBUG)
		if (CHK_BIT(kmap_pos, (w32)k) == 0) {
		    printf("%lld %ld %lld  +1\n", q, n, low_k+k);
		}
#endif
		SET_BIT(kmap_pos, (long)k);
	    }
	}

	/*
	 * update external q
	 */
	external_q = q;
    }

    /* 
     * disable checkpoint alarm 
     */
    (void) sigignore(SIGALRM);

    /*
     * perform a final checkpoint
     */
    checkpt();

    /*
     * all done
     */
    exit(0);
}


/*
 * lock_ckpt - lock the setname.ckpt file
 *
 * given:
 *	name	name of our set
 *
 * returns:
 *	-1  ==>  ckpt file already locked
 *	fd>=0 ==>  ckpt file was not locked, but is now
 *	does not return on error
 */
int
lock_ckpt(name)
    char *name;		/* name of state set */
{
    char filename[BUFSIZ+1];	/* filename holder */
    int fd;			/* locking file descriptor */

    /*
     * open the filename
     */
    sprintf(filename, "%s.ckpt", name);
    fd = open(filename, O_RDWR);
    if (fd < 0) {
	fprintf(stderr, "%s: error in opening %s.ckpt\n", program, setname);
	exit(8);
    }

    /*
     * lock the filename
     */
    errno = 0;
    if (flock(fd, LOCK_EX|LOCK_NB) < 0) {
	/* EWOULDBLOCK means already locked */
	if (errno == EWOULDBLOCK) {
	    /* already locked */
	    return -1;
	}
    }
    return fd;
}


/*
 * checkpt - checkpoint our kmap files and update the comp file
 */
void
checkpt()
{
    FILE *stream;		/* file stream */
    char filename[BUFSIZ+1];	/* old filename holder */
    int i;

    /*
     * update the -1 kmap
     */
    i = write_kmap(setname, kmap_neg, -1);
    if (i < 0) {
	fprintf(stderr, "%s: bad write of %s.neg, return: %d\n",
	  program, setname, i);
	exit(9);
    }

    /*
     * update the +1 kmap
     */
    i = write_kmap(setname, kmap_pos, +1);
    if (i < 0) {
	fprintf(stderr, "%s: bad write of %s.pos, return: %d\n",
	  program, setname, i);
	exit(10);
    }

    /*
     * append the current q to the setname.ckpt file
     */
    /* we tested the entire range, no need to double check */
    fprintf(ckpt, "%lld\n", external_q);
    /* flush to disk */
    if (fflush(ckpt) != 0) {
	fprintf(stderr, "%s: flush error on %s\n", program, filename);
	perror("setname.ckpt fflush");
	exit(11);
    }
    if (fsync(fileno(ckpt)) < 0) {
	fprintf(stderr, "%s: fsync error on %s\n", program, filename);
	perror("setname.ckpt fsync");
	exit(12);
    }

    /*
     * sync for good measure
     */
    if (sync() < 0) {
	fprintf(stderr, "%s: sync error in checkpt\n", program);
	perror("sync");
	exit(13);
    }
    return;
}


/*
 * write_kmap - write the name/name.neg or name.pos kmap file
 *
 * These files contain the bit maps of k values.
 *
 *	kmap_neg[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n-1 is composite
 *	kmap_neg[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n-1 is unknown
 *
 *	kmap_pos[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n+1 is composite
 *	kmap_pos[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n+1 is unknown
 */
int
write_kmap(name, kmap, posneg)
    char *name;		/* name of state set */
    char *kmap;		/* malloced kmap buffer */
    int posneg;		/* 1 ==> read name.pos, -1 ==> read name.neg */
{
    FILE *stream;		/* file stream */
    char filename[BUFSIZ+1];	/* old filename holder */
    char newfilename[BUFSIZ+1];	/* new filename holder */

    /*
     * firewall
     */
    if (posneg != 1 && posneg != -1 && (k_len&0x7) != 0) {
	fprintf(stderr, "%s: bad args to parse_kmap\n", program);
	return -1;
    }

    /*
     * open the file for truncation/writing
     */
    if (posneg == 1) {
	sprintf(filename, "%s.pos", name);
	sprintf(newfilename, "new.%s.pos", name);
    } else {
	sprintf(filename, "%s.neg", name);
	sprintf(newfilename, "new.%s.neg", name);
    }
    stream = fopen(newfilename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot write %s\n", program, filename);
	if (posneg == 1) {
	    perror("newname.pos open");
	} else {
	    perror("newname.neg open");
	}
	return -2;
    }

    /*
     * write the kmap
     */
    if (fwrite(kmap, sizeof(char), (size_t)(k_len/8), stream) !=
	(int)(k_len/8)) {
	fprintf(stderr, "%s: cannot write %s\n", program, filename);
	if (posneg == 1) {
	    perror("newname.pos fwrite");
	} else {
	    perror("newname.neg fwrite");
	}
	return -3;
    }

    /*
     * flush map to disk
     */
    if (fflush(stream) != 0) {
	fprintf(stderr, "%s: flush error on %s\n", program, filename);
	if (posneg == 1) {
	    perror("newname.pos fflush");
	} else {
	    perror("newname.neg fflush");
	}
	return -5;
    }
    if (fsync(fileno(stream)) < 0) {
	fprintf(stderr, "%s: fsync error on %s\n", program, filename);
	if (posneg == 1) {
	    perror("newname.pos fsync");
	} else {
	    perror("newname.neg fsync");
	}
	return -6;
    }

    /*
     * close the kmap
     */
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	if (posneg == 1) {
	    perror("newname.pos fclose");
	} else {
	    perror("newname.neg fclose");
	}
	return -7;
    }

    /*
     * move the kew kmap to the old kmap
     */
    if (rename(newfilename, filename) < 0) {
	fprintf(stderr, "%s: rename of %s to %s failed\n",
	  program, newfilename, filename);
	perror("rename");
	return -8;
    }
    return 0;
}


/*
 * parse_kmap - parse the name/name.neg or name.pos kmap file
 *
 * These files contain the bit maps of k values.
 *
 *	kmap_neg[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n-1 is composite
 *	kmap_neg[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n-1 is unknown
 *
 *	kmap_pos[y]&(1<<(x-1)) == 1  ==>  (low_k + 8*y + x)*2^n+1 is composite
 *	kmap_pos[y]&(1<<(x-1)) == 0  ==>  (low_k + 8*y + x)*2^n+1 is unknown
 */
int
parse_kmap(name, kmap, posneg)
    char *name;		/* name of state set */
    char **kmap;	/* malloced kmap buffer */
    int posneg;		/* 1 ==> read name.pos, -1 ==> read name.neg */
{
    FILE *stream;		/* file stream */
    struct stat statbuf;	/* status of the open file */
    char filename[BUFSIZ+1];	/* filename holder */

    /*
     * firewall
     */
    if (posneg != 1 && posneg != -1 && (k_len&0x7) != 0) {
	fprintf(stderr, "%s: bad args to parse_kmap\n", program);
	return -1;
    }

    /*
     * open the file
     */
    if (posneg == 1) {
	sprintf(filename, "%s.pos", name);
    } else {
	sprintf(filename, "%s.neg", name);
    }
    stream = fopen(filename, "r");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot read %s\n", program, filename);
	if (posneg == 1) {
	    perror("setname.pos open");
	} else {
	    perror("setname.neg open");
	}
	return -2;
    }

    /*
     * verify the size of the kmap
     */
    if (fstat(fileno(stream), &statbuf) < 0) {
	fprintf(stderr, "%s: cannot stat %s\n", program, filename);
	if (posneg == 1) {
	    perror("setname.pos stat");
	} else {
	    perror("setname.neg stat");
	}
	return -3;
    }
    if (statbuf.st_size != (off_t)(k_len/8)) {
	fprintf(stderr, "%s: expected %d bytes, found %d in %s\n",
	  program, (int)statbuf.st_size, (int)k_len, filename);
	return -4;
    }

    /*
     * malloc the kmap storage
     */
    *kmap = (char *)malloc((long)(k_len/8+1) * sizeof(char));
    if (*kmap == NULL) {
	fprintf(stderr, "%s: malloc of %d bytes for %s kmap failed\n",
	  program, (int)(k_len/8+1), filename);
	if (posneg == 1) {
	    perror("setname.pos malloc");
	} else {
	    perror("setname.neg malloc");
	}
	return -5;
    }

    /*
     * read the kmap
     */
    if (fread(*kmap, sizeof(char), (size_t)(k_len/8), stream) !=
	(int)(k_len/8)) {
	fprintf(stderr, "%s: cannot read %s\n", program, filename);
	if (posneg == 1) {
	    perror("setname.pos fread");
	} else {
	    perror("setname.neg fread");
	}
	return -6;
    }

    /*
     * all done
     */
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	if (posneg == 1) {
	    perror("setname.pos fclose");
	} else {
	    perror("setname.neg fclose");
	}
	return -7;
    }
    return 0;
}


/*
 * parse_parm - parse name/name.parm  
 *
 * This file contains version, low_k, high_k and n.
 *
 * Returns 0 if parsed ok, <0 otherwise.  Sets the low_k_p, k_len and
 * n_p to the low_k, k_len and n values found, respectively.
 */
int
parse_parm(name, low_k_p, n_p)
    char *name;		/* name of state set */
    w64 *low_k_p;	/* initial multiplier as in low_k*2^n+/-1 */
    long *n_p;		/* exponent in (k+i)*2^n+/-1 */
{
    FILE *stream;		/* file stream */
    char filename[BUFSIZ+1];	/* filename holder */
    char buf[BUFSIZ+1];		/* read buffer */
    w64 k_len_val;		/* k_len value read */

    /*
     * open file
     */
    sprintf(filename, "%s.parm", name);
    stream = fopen(filename, "r");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot read %s\n", program, filename);
	perror("setname.parm open");
	return -1;
    }

    /*
     * read version
     */
    if (fscanf(stream, "version = %s\n", buf) != 1) {
	fprintf(stderr, "%s: unable to read version from %s\n",
	  program, filename);
	return -2;
    }
    if (strcmp(buf, VERSION) != 0) {
	fprintf(stderr, "%s: expected version %s, found %s\n",
	  program, buf, VERSION);
	return -3;
    }

    /*
     * read low_k
     */
    if (fscanf(stream, "low_k = %lld\n", low_k_p) != 1) {
	fprintf(stderr, "%s: unable to read low_k from %s\n",
	  program, filename);
	return -4;
    }
    if (*low_k_p < 0) {
	fprintf(stderr, "%s: low_k must be >= 0\n", program);
	return -5;
    }
    if ((*low_k_p & 0x7) != 0) {
	fprintf(stderr, "%s: low_k must be a multiple of 8\n", program);
	return -6;
    }

    /*
     * read k_len
     */
    if (fscanf(stream, "k_len = %lld\n", &k_len_val) != 1) {
	fprintf(stderr, "%s: unable to read k_len from %s\n",
	  program, filename);
	return -7;
    }
    if (k_len_val < 0) {
	fprintf(stderr, "%s: k_len must be >= 0\n", program);
	return -8;
    }
    if ((k_len_val & 0x7LL) != 0) {
	fprintf(stderr, "%s: k_len must be a multiple of 8\n", program);
	return -9;
    }
    if (k_len_val > (1LL<<24)) {
	fprintf(stderr, "%s: k_len %lld must currently be <= 2^24 (%d)\n",
	  program, k_len_val, (1<<24));
	return -10;
    }
    k_len = (long)k_len_val;

    /*
     * read n
     */
    if (fscanf(stream, "n = %ld\n", n_p) != 1) {
	fprintf(stderr, "%s: cannot read n from %s\n", program, filename);
	return -11;
    }
    if (*n_p <= 2) {
	fprintf(stderr, "%s: n must be > 2\n", program);
	return -12;
    }

    /*
     * all done
     */
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("setname.parm close");
	return -13;
    }
    return 0;
}


/*
 * last_q - find the last q written into the setname.ckpt file
 *
 * The setname.ckpt file contains the most recent q value tested (q checkpoint
 * value).  The number of the last line represents the last q value tested.
 *
 * Returns 0 if a q value was found and sets *last_q, otherwise -1 is returned
 * to indicate error.
 */
int
last_q(name, prev_q)
    char *name;		/* name of state set */
    w64 *prev_q;	/* last q tested as found in the setname.ckpt file */
{
    FILE *stream;		/* file stream */
    char filename[BUFSIZ+1];	/* filename holder */
    char buf[BUFSIZ+1];		/* read buffer */
    struct stat statbuf;	/* status of the open file */
    char *p;

    /*
     * read the end of the file
     */
    if (fstat(fileno(ckpt), &statbuf) < 0) {
	fprintf(stderr, "%s: cannot stat %s\n", program, filename);
	return -2;
    }
    if (statbuf.st_size <= 0) {
	fprintf(stderr, "%s: %s is an empty file\n", program, filename);
	return -3;
    }
    buf[BUFSIZ] = '\0';
    if (statbuf.st_size <= BUFSIZ) {

	/* file is smaller than BUFSIZ, just read it all */
	if (fread(buf, 1, statbuf.st_size, ckpt) != statbuf.st_size) {
	    fprintf(stderr, "%s: cannot read %s\n", program, filename);
	    return -4;
	}
	p = buf+statbuf.st_size-1; /* final byte in the buffer */

    } else {
	/* file is larger than BUFSIZ, read the final BUFSIZ bytes */
	if (fseek(ckpt, (long)-BUFSIZ, SEEK_END) != 0) {
	    fprintf(stderr, "%s: cannot seek on %s\n", program, filename);
	    return -5;
	}
	if (fread(buf, 1, BUFSIZ, ckpt) != BUFSIZ) {
	    fprintf(stderr, "%s: cannot end read %s\n", program, filename);
	    return -6;
	}
	p = buf+BUFSIZ-1; /* final byte in the buffer */
    }
    if (*p != '\n') {
	fprintf(stderr,
	   "%s: %s does not end in a newline\n", program, filename);
	return -7;
    }

    /*
     * find the beginning of the last line
     */
    for (--p; p > buf && *p != '\n'; --p) {
    }
    if (p > buf) {
	++p;
    }

    /*
     * pull out the q value
     */
    if (!isdigit((int)(*p))) {
	fprintf(stderr,
	   "%s: non-numeric 1st char of last line of %s\n", program, filename);
	return -8;
    }
    if (sscanf(p, "%lld", prev_q) != 1) {
	fprintf(stderr,
	   "%s: cannot scan last q from last line of %s\n", program, filename);
	return -9;
    }

    /*
     * seek to the end of file
     */
    if (fseek(ckpt, 0, SEEK_END) != 0) {
	fprintf(stderr, "%s: cannot seek to end on %s\n", program, filename);
	return -10;
    }

    /*
     * all done
     */
    return 0;
}