/* @(#) $Revision: 1.19 $ */
/* @(#) RCS control in //prime.corp/usr/local/src/cmd/prime/src/init_sieve.c */
/*
 * init_sieve - initialize sieve state
 *
 * usage:
 *	init_sieve low_k k_len n setname
 *
 * 	low_k	low k value in the k*2^n+/-1 sequence
 * 	k_len	high k value in the k*2^n+/-1 sequence is low_k+k_len-1
 *	n	exponent in the k*2^n+/-1 sequence
 *	setname	base state filename
 *
 * We will setup state files for sieve to work on.  The state consists
 * of several files:
 *
 *	setname/setname.parm	version, low_k, k_len, n global parameters
 *	setname/setname.neg	array of bits bit x represents (low_k+x)*2^n-1
 *				    1 => composite, 0 => unknown
 *	setname/setname.pos	array of bits bit x represents (low_k+x)*2^n+1
 *				    1 => composite, 0 => unknown
 *	setname/setname.ckpt	q checkpoint value
 *
 * After this program complete, the build.mk makefile needs be executed
 * (usually by do_init) as follows:
 *
 *	cd ../../src
 *	make -f build.mk setup all SETNAME="setname" TOPDIR=/usr/local/lib/map
 *
 * This makefile will build:
 *
 *	setname/map0		binary to complete sieve for q < 2^31
 *	setname/map1		binary to sieve 2^31 <= q < 2^39
 *
 * These files are not created by build.mk, but could be later on:
 *
 *	setname/map1a		binary to sieve 2^39 <= q < 2^40
 *	setname/map2		binary to sieve 2^40 <= q < 2^47
 */
/*
 *
 * 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/stat.h>
#include "primetbl.h"

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


/*
 * 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)))


/* 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

#define NO_MAKE		/* set do_init run the build.mk file */


/*
 * 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;

char *program;			/* our name */


/* forward declarations */
#if !defined(NO_MAKE)
void makeit(char*);
#endif
uw64 minv(uw64, uw64);


main(argc, argv)
    int argc;		/* arg count */
    char *argv[];	/* the args */
{
    struct stat statbuf;	/* inode status */
    char filename[BUFSIZ+1];	/* filename holder */
    FILE *stream;	/* file stream */
    w64 low_k;		/* initial multiplier as in low_k*2^n+/-1 */
    w64 k_len;		/* number of k to test */
    w64 end_k;		/* final multiplier is (end_k-1)*2^n+/-1 */
    char *name;		/* base state filename */
    unsigned char const *prm;	/* pointer into prime table */
    uw64 q;		/* test factor */
    w32 top5;		/* top 5 bits of n */
    uw64 modq;		/* 2^x mod q work value ending in 2^n mod q */
    uw64 minv_val;	/* modular inverse of 2^n mod q */
    w32 n;		/* exponent in (k+i)*2^n+/-1 */
    char nstep[32-5];	/* bits of n in reverse order (leading 0's ignored) */
    char step_cnt;	/* number of steps in the nstep array */
    int nlog2;		/* smallest power of 2 >= n */
    w64 k;		/* k in k*2^n+/-1 to sieve out */
    long tmp;
    int i;
    int j;

    /*
     * parse args
     */
    program = argv[0];
    if (argc != 5) {
	fprintf(stderr, "usage: %s low_k k_len n name\n", program);
	fprintf(stderr, "\tlow_k\tlow k in the k*2^n+/-1 sequence\n");
	fprintf(stderr,
	  "\thigh_k\thigh k in the k*2^n+/-1 sequence is low_k+k_len-1\n");
	fprintf(stderr, "\tn\texponent in the k*2^n+/-1 sequence\n");
	fprintf(stderr, "\tname\tbase state filename\n");
	exit(1);
    }
    sscanf(argv[1], "%lld", &low_k);
    if (low_k < 0) {
	fprintf(stderr, "%s: low_k must be >= 0\n", program);
	exit(2);
    }
    if ((low_k & 0x7) != 0) {
	fprintf(stderr, "%s: low_k must be a multiple of 8\n", program);
	exit(3);
    }
    sscanf(argv[2], "%lld", &k_len);
    if (k_len < 8) {
	fprintf(stderr, "%s: k_len must be >= 8\n", program);
	exit(4);
    }
    if ((k_len & 0x7LL) != 0) {
	fprintf(stderr, "%s: k_len must be a multiple of 8\n", program);
	exit(5);
    }
    if (k_len > (1LL<<24)) {
	fprintf(stderr, "%s: k_len must currently be <= 2^24 (%d)\n",
	  program, (1<<24));
	exit(6);
    }
    end_k = low_k+k_len;
    sscanf(argv[3], "%d", &n);
    if (n <= 2) {
	fprintf(stderr, "%s: n must be > 2\n", program);
	exit(7);
    }
    name = argv[4];

    /*
     * make the directory if needed
     */
    if (stat(name, &statbuf) < 0) {
	/* no such inode, make a directory */
	if (mkdir(name, 0755) < 0) {
	    fprintf(stderr, "%s: cannot mkdir %s\n", program, name);
	    perror("mkdir");
	    exit(8);
	}
    } else {
	/* object if not a writable directory */
	if (! S_ISDIR(statbuf.st_mode)) {
	    fprintf(stderr, "%s: %s exists and is not a directory\n",
	      program, name);
	    exit(9);
	}
    }

    /*
     * initialize the kmap
     */
    kmap_neg = (char *)malloc((long)(k_len/8+1) * sizeof(char));
    if (kmap_neg == NULL) {
	fprintf(stderr,
	 "%s: cannot malloc kmap_neg of %d bytes\n",
	 program, (long)(k_len/8+1));
	exit(10);
    }
    kmap_pos = (char *)malloc((long)(k_len/8+1) * sizeof(char));
    if (kmap_pos == NULL) {
	fprintf(stderr,
	 "%s: cannot malloc kmap_pos of %d bytes\n",
	 program, (long)(k_len/8+1));
	exit(11);
    }
    memset(kmap_neg, 0, (long)(k_len/8+1));
    memset(kmap_pos, 0, (long)(k_len/8+1));

    /*
     * If the first bit in our kmaps represent 0*2^n+/-1 then mark it
     * as composite (1 and -1 are not prime).
     */
    if (low_k <= 0) {
	SET_BIT(kmap_neg, 0);
	SET_BIT(kmap_pos, 0);
    }

    /*
     * pre-digest n
     *
     * We will form top5 the top 5 bits of n.
     * The valie 2^top5 will be the initial stage for the 2^n mod q computation.
     *
     * We will form the reverse bit values of n minus the top 5 bits.  This
     * will be the state array for other stages of the 2^n mod q computation.
     */
    /* compute the "log" of n */
    for (nlog2=0; (nlog2 < 31) && ((1<<nlog2) < n); ++nlog2) {
    }
    /* form the top5 value and nstep array */
    if (nlog2 <= 5) {
	/* empty nstep array */
	top5 = n;
	step_cnt = 0;	
    } else {
	top5 = ((n >> (nlog2-5))&0x1f);
	for (i=(nlog2-6); i >= 0; --i) {
	    nstep[(nlog2-6)-i] = ((n & (1<<i)) ? 1 : 0);
	}
	step_cnt = nlog2 - 5;
    }
#if 0
#if defined(DEBUG)
    /* print some information about n */
    printf("n=%d(%x)  nlog2=%d  top5=%lx  step_cnt=%d  nstep=",
      n, n, nlog2, top5, step_cnt); 
    for(i=0; i < step_cnt; ++i) {
	printf("%1d", (int)nstep[i]);
    }
    printf("\n\n");
#endif
#endif

    /*
     * eliminate values disivble by primes < NEXT_PRIME_DIFF
     */
    q = (uw64)3;
    prm = prime_diff;
    do {
	/*
	 * 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.
	 */
	modq = (1ULL << top5);
	for (i=0; i < step_cnt; ++i) {

	     /*
	      * perform the mod q step for either x^2 or 2*x^2
	      */
	     modq %= q;

	     /*
	      * start the next stage depending on the next bit of n
	      *
	      * nstep[i] == 0  ==>  (x^2 mod q) = y^2 mod q
	      * nstep[i] == 1  ==>  (2*x^2 mod q) = 2*(y^2) mod q
	      *
	      * We always to x^2, and only 2* if nstep[i] == 1.
	      */
	    modq *= modq;
	    if (nstep[i]) {
		modq <<= 1;	/* times 2 */
	    }
	}
	/* performt the final mod */
	modq %= q;

	/*
	 * compute the modular inverse of 2^n mod q
	 */
	minv_val = minv(modq, q);

	/*
	 * 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 == minv_val % q.
	 */
	k = (w64)((w64)minv_val - (low_k % q));
	for (k = ((k < 0) ? k+q : k); k < k_len; k += q) {
#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, (w32)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-minv_val) % q.
	 */
	k = (w64)(q - (w64)minv_val - (low_k % q));
	for (k = ((k < 0) ? k+q : k); k < k_len; k += q) {
#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, (w32)k);
	}

	/* get the next test factor */
	q += (uw64)*prm;

    } while (*prm++ != 1);

    /*
     * form name/name.parm	
     *
     * This file contains version, low_k, high_k and n.
     */
    sprintf(filename, "%s/%s.parm", name, name);
    stream = fopen(filename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot create/truncate %s\n", program, filename);
	perror("name.parm open");
	exit(12);
    }
    fprintf(stream, "version = %s\n", VERSION);
    fprintf(stream, "low_k = %lld\n", low_k);
    fprintf(stream, "k_len = %lld\n", k_len);
    fprintf(stream, "n = %ld\n", n);
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("name.parm close");
	exit(13);
    }

    /*
     * form name/name.neg
     *
     * array of bits bit x represents (low_k+x)*2^n-1
     * 1 => composite, 0 => unknown
     */
    sprintf(filename, "%s/%s.neg", name, name);
    stream = fopen(filename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot create/truncate %s\n", program, filename);
	perror("name.neg open");
	exit(14);
    }
    if (fwrite(kmap_neg, sizeof(kmap_neg[0]), (size_t)(k_len/8), stream) !=
	(int)(k_len/8)) {
	fprintf(stderr, "%s: cannot write %s\n", program, filename);
	perror("name.neg write");
	exit(15);
    }
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("name.q close");
	exit(16);
    }

    /*
     * form name/name.pos	
     *
     * array of bits bit x represents (low_k+x)*2^n+1
     * 1 => composite, 0 => unknown
     */
    sprintf(filename, "%s/%s.pos", name, name);
    stream = fopen(filename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot create/truncate %s\n", program, filename);
	perror("name.pos open");
	exit(17);
    }
    if (fwrite(kmap_pos, sizeof(kmap_pos[0]), (size_t)(k_len/8), stream) !=
	(int)(k_len/8)) {
	fprintf(stderr, "%s: cannot write %s\n", program, filename);
	perror("name.pos write");
	exit(18);
    }
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("name.q close");
	exit(19);
    }

    /*
     * form name/name.ckpt	
     *
     * This file contains the q checkpoint (last q tested) value.
     *
     * We will print just the q checkpoint value of LAST_PRIME_DIFF.
     */
    sprintf(filename, "%s/%s.ckpt", name, name);
    stream = fopen(filename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot create/truncate %s\n", program, filename);
	perror("name.ckpt open");
	exit(20);
    }
    fprintf(stream, "%d\n", LAST_PRIME_DIFF);
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("name.ckpt close");
	exit(21);
    }

    /*
     * form the name/nbits.h file
     */
    sprintf(filename, "%s/nbits.h", name);
    stream = fopen(filename, "w");
    if (stream == NULL) {
	fprintf(stderr, "%s: cannot read %s\n", program, filename);
	perror("nbits.h create");
	exit(22);
    }
    /* form start of the file */
    fputs("/* created by mksieve -- DO NOT EDIT */\n", stream);
    fputc('\n', stream);
    /* form the define values */
    fprintf(stream, "#define N_VALUE %d\t\t/* exponent in k*2^n+/-1 */\n", n);
    fputc('\n', stream);
    fprintf(stream, "#define TOP_5_N_BITS %d\t/* top 5 bits of n */\n", top5);
    fputc('\n', stream);
    /* process n (skipping bit 0 to 4) from bit 5 on to bit 31 if needed */
    for (i=5; i < 32; ++i) {
	if (i-5 < step_cnt) {
	    if (nstep[i-5]) {
		fprintf(stream,
		  "#define N_BIT_%d(mod,tmp,q) NBIT_IS_1(mod,tmp,q)\n", i);
	    } else {
		fprintf(stream,
		  "#define N_BIT_%d(mod,tmp,q) NBIT_IS_0(mod,tmp,q)\n",i);
	    }
	} else {
	    fprintf(stream, "#define N_BIT_%d(mod,tmp,q) /* undefined */\n", i);
	}
    }
    /* close up */
    if (fclose(stream)) {
	fprintf(stderr, "%s: close error on %s\n", program, filename);
	perror("nbits.h fclose");
	exit(23);
    }

    /*
     * all done
     */
#if defined(NO_MAKE)
    exit(0);
#else
    /*
     * make everything
     */
    fflush(stdout);
    fflush(stderr);
    makeit(name);
    /* how did we get here? */
    exit(24);
#endif
}


#if !defined(NO_MAKE)
/*
 * makeit - exec make to build the map0, map1, map1a and map2 binaries
 *
 * We do not return.
 */
void
makeit(setname)
    char *setname;	/* name of state set */
{
    char buf[BUFSIZ+1];		/* setname symbol */

    /*
     * form the SETNAME=arg
     */
    sprintf(buf, "SETNAME=%s", setname);

    /*
     */
    execl("/sbin/make", "make", "-f", "build.mk", buf, "setup", "all", NULL);
    fputs("exec error!\n", stderr);
    perror("exec");
    exit(26);
}
#endif


/*
 * minv - calculate modular inverse
 *
 * 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
 *
 * given:
 *	u	compute the modular inverse of this value (must be < 2^63)
 *	v	compute the inverse with this modulus (must be < 2^63)
 *
 * return:
 *	multiplicitive inverse of u mod v or 0 if no solution
 */
uw64
minv(u, v)
    uw64 u;		/* compute the modular inverse of this value */
    uw64 v;		/* compute the inverse with this modulus */
{
    w64 u2,v2,u3,v3;	/* Eclidian terms */
    w64 q1, q2;		/* intermediate quotents */
    w64 tmp1, tmp2;	/* temp values */

    /* prep work */
    if (u >= v) {
       v3 = u % v;
    } else {
       v3 = u;
    }
    u3 = v;
    u2 = 0;
    v2 = 1;

    /* catch an easy no-solution */
    if (v3 == 0 && u3 != 1) {
       /* no solution */
       return (uw64)0;
    }

    /* 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;
    }

    /* obtain our solution, if any exists */
    if (u3 != 1) {
       /* no solution */
       return (uw64)0;
    }
    return ((u2 < 0) ? (uw64)(v+u2) : (uw64)u2);
}