/* @(#) $Revision: 1.5 $ */
/* @(#) RCS control in //prime.corp/usr/local/src/cmd/prime/src/minv.c */
/*
 * minv - multiplicitive inverse
 */
/*
 *
 * 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!
 */

#if defined(STANDALONE)
#include <stdio.h>
#endif

typedef long long large;


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

    /* prep work */
    if (u < 0 || v < 0 || 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 (large)0;
    }

    /* Eclidian process */
    while (v3 != 0) {
        q1 = u3 / v3;
#if defined(STANDALONE)
	if (((v2 > 0 && v2 >= (1LL<<32)) || (v2 < 0 && v2 <= -(1LL<<32))) &&
	    ((q1 > 0 && q1 >= (1LL<<32)) || (q1 < 0 && q1 <= -(1LL<<32)))) {
	    fprintf(stderr,
	      "multiply overflow v2:0x%llx q1:0x%llx\n", v2, q1);
	}
#endif
        tmp1 = v2 * q1;
        tmp2 = u2 - tmp1;
        u2 = v2;
        v2 = tmp2;
#if defined(STANDALONE)
	if (((v3 > 0 && v3 >= (1LL<<32)) || (v3 < 0 && v3 <= -(1LL<<32))) &&
	    ((q1 > 0 && q1 >= (1LL<<32)) || (q1 < 0 && q1 <= -(1LL<<32)))) {
	    fprintf(stderr,
	      "multiply overflow v3:0x%llx q1:0x%llx\n", v3, q1);
	}
#endif
        q2 = u3 - (q1*v3);
        u3 = v3;
        v3 = q2;
    }

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


#if defined(STANDALONE)

char *program;		/* our name */

/*
 * minv - calculate modular inverse
 *
 * usage:
 *	minv number mod
 */
main(argc, argv)
    int argc;		/* arg count */
    char *argv[];	/* the args */
{
    large u;		/* compute the modular inverse of this value */
    large v;		/* compute the inverse with this modulus */
    large res;		/* modular inverse of u mod v */

    /*
     * parse args
     */
    program = argv[0];
    if (argc != 3) {
	fprintf(stderr, "usage: %s number mod\n", program);
	exit(1);
    }
    sscanf( argv[1], "%lld", &u );
    sscanf( argv[2], "%lld", &v );

    /*
     * print modular inverse of u mod v
     */
    res = minv(u,v);
    printf("%lld\n", res);
    exit(0);
}
#endif /* STANDALONE */