/* @(#) $Revision: 1.6 $ */
/* @(#) RCS control in //prime.corp/usr/local/src/cmd/prime/src/sqrmod.c */
/*
 * quomod - compute 2*x^2 mod q
 *
 * usage:
 *	quomod x q
 *
 *	where 0 <= x < q  and  1 < q < 2^47
 *
 * Print the value 2*x^2 mod q (where ^ means exponentation)
 *
 * This test simulates the double precision mod that will be needed
 * for sieveing out w64 factors q < 2^47.  The value x may be thought
 * of as an intermediate value in the 2^n mod q calculation.   As such
 * we will assume that 0 <= x < q.
 */
/*
 *
 * 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>


/* 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 vlarge
 *
 * B64_OFF16 - 64 bits offset multiple of 16 bits within a 128 bit vlarge
 * B32_OFF16 - 32 bits offset multiple of 16 bits within a 128 bit vlarge
 *
 * B64_OFF8 - 64 bits offset multiple of 8 bits within a 128 bit vlarge
 * B32_OFF8 - 32 bits offset multiple of 8 bits within a 128 bit vlarge
 * B16_OFF8 - 16 bits offset multiple of 8 bits within a 128 bit vlarge
 */
#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])))


char *program;		/* our name */


/* forward declarations */
void printlarge(char*, large);
void printvlarge(char*, vlarge);


main(argc, argv)
    int argc;		/* arg count */
    char *argv[];	/* the args */
{
    large x;		/* the value to 2*square */
    large qlarge;	/* modulus as a large structure */
    uw64 q;		/* modulus */
    vlarge x2;		/* x^2 and 2*x^2 */
    large r;		/* 2*x^2 mod q */
    w64 xarg;		/* x input value */
    w64 qarg;		/* q input value */
    int i;

    /*
     * parse args
     */
    program = argv[0];
    if (argc != 3) {
	fprintf(stderr, "usage %s x q\n", program);
	exit(1);
    }
    sscanf(argv[1], "%lld", &xarg);
    if (xarg < 0 || xarg >= 0x7fffffffffffULL) {
	fprintf(stderr, "%s: x must be >= 0 and x < 2^47\n", program);
	exit(2);
    }
    sscanf(argv[2], "%lld", &qarg);
    if (qarg < 2 || qarg >= 0x7fffffffffffULL) {
	fprintf(stderr, "%s: q must be > 1 and q < 2^47\n", program);
	exit(3);
    }
    if (xarg >= qarg) {
	fprintf(stderr, "%s: for this test, x < q\n", program);
	exit(4);
    }
    x.b64 = (uw64)xarg;
    q = (uw64)qarg;
    qlarge.b64 = q;

    /*
     * detail values
     */
    printlarge("x", x);
    printlarge("q", qlarge);

    /*
     * 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);
    printvlarge("x^2", x2);

    /*
     * form 2*x^2
     */
    x2.b64[0] = ((x2.b64[0] << 1) | (uw64)(x2.b32[2] >> 31));
    x2.b64[1] <<= 1;
    printvlarge("2*x^2", x2);

    /*
     * form 2*x^2 mod q
     */
    if (q < (1ULL<<31)) {

	/*
	 * mode 0:      0 <= q < 2^31  x2 < 2^63 
	 *
	 * x2 % q = x2 % q
	 */
	printf("mode 0\n");
	r.b64 = x2.b64[1] % q;


    } else if (q < (1ULL<<39)) {

	/* 
	 * mode 1:    2^31 <= q < 2^39  x2 < 2^79 
	 *
	 * x2 % q = (((x2 >> 16) % q)<<16 | (x2 % 2^16)) % q
	 */
	printf("mode 1\n");
	r.b64 = ((B64_OFF16(x2,3) % q)<<16 | (uw64)(x2.b16[7])) % q;

    } else if (q < (1ULL<<40)) {

	/* 
	 * mode 1.5:  2^39 <= q < 2^40  x2 < 2^81 
	 *
	 * x2 % q = (((x2 >> 24) % q)<<24 | (x2 % 2^24)) % q
	 */
	printf("mode 1.5\n");
	r.b64 = ((B64_OFF8(x2,5) % q)<<24 | (uw64)(x2.b32[3]&0xffffff)) % q;

    } else if (q < (1ULL<<47)) {

	/* 
	 * mode 2:    2^40 <= q < 2^47  x2 < 2^95 
	 *
	 * x2 % q = ( ((((x2 >> 32) % q)<<16 | ((x2 % 2^32)<<16)) % q)<<16 |
	 *	      (x2 % 2^16)) % q
	 */
	printf("mode 2\n");
	r.b64 = ( ((((B64_OFF32(x2)%q)<<16 | (uw64)x2.b16[6]) % q)<<16) |
	         (uw64)x2.b16[7]) % q;

    } else {

	/* 
	 * mode 3:    2^47 <= q 
	 *
	 * not done
	 */
	printf("mode 3\n");
	fprintf(stderr, "%s: how did e get into mode 3?\n", program);
	exit(5);

    }
    printlarge("2*x^2 % q", r);

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


/*
 * printlarge - print a 64 bit value
 *
 * given:
 *	name - what to call the value
 *	x - 64 bit large to print
 */
void
printlarge(char *str, large x)
{
    printf("%s = 0x%016llx (%lld)\n", str, x.b64, x.b64);
    printf("%s = 0x%08x %08x\n", str, x.b32[0], x.b32[1]);
    printf("%s = 0x%04x %04x %04x %04x\n", str,
      x.b16[0], x.b16[1], x.b16[2], x.b16[3]);
    printf("%s = 0x%02x %02x %02x %02x %02x %02x %02x %02x\n", str,
      x.b8[0], x.b8[1], x.b8[2], x.b8[3], x.b8[4], x.b8[5], x.b8[6], x.b8[7]);
    return;
}


/*
 * printvlarge - print a 128 bit value
 *
 * given:
 *	name - what to call the value
 *	x - 128 bit large to print
 */
void
printvlarge(char *str, vlarge x)
{
    printf("%s = 0x%016llx%016llx\n", str, x.b64[0], x.b64[1]);

    printf("%s = 0x%08x %016llx %08x\n", str,
      x.b32[0], B64_OFF32(x), x.b32[3]);

    printf("%s = 0x%08x %08x %08x %08x\n", str,
      x.b32[0], x.b32[1], x.b32[2], x.b32[3]);

    printf("%s = 0x%04x %08x %08x %08x %04x\n", str,
      x.b16[0], B32_OFF16(x,1), B32_OFF16(x,3), B32_OFF16(x,5), x.b16[7]);

    printf("%s = 0x%04x %04x %04x %04x %04x %04x %04x %04x\n", str,
      x.b16[0], x.b16[1], x.b16[2], x.b16[3],
      x.b16[4], x.b16[5], x.b16[6], x.b16[7]);

    printf("%s = 0x%02x %02x %02x %02x %02x %02x %02x %02x", str,
      x.b8[0], x.b8[1], x.b8[2], x.b8[3], x.b8[4], x.b8[5], x.b8[6], x.b8[7]);

    printf(" %02x %02x %02x %02x %02x %02x %02x %02x\n",
      x.b8[8], x.b8[9], x.b8[10], x.b8[11], x.b8[12], x.b8[13], x.b8[14],
      x.b8[15]);

    return;
}