2N/A/* Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
2N/A/* All Rights Reserved */
2N/A
2N/A
2N/A/*
2N/A * Copyright (c) 1980 Regents of the University of California.
2N/A * All rights reserved. The Berkeley software License Agreement
2N/A * specifies the terms and conditions for redistribution.
2N/A */
2N/A/* Portions Copyright(c) 1988, Sun Microsystems Inc. */
2N/A/* All Rights Reserved */
2N/A
2N/A/*
2N/A * Copyright (c) 1997, by Sun Microsystems, Inc.
2N/A * All rights reserved.
2N/A */
2N/A
2N/A#ident "%Z%%M% %I% %E% SMI" /* SVr4.0 1.1 */
2N/A
2N/A/* LINTLIBRARY */
2N/A
2N/A#include <mp.h>
2N/A#include "libmp.h"
2N/A#include <sys/types.h>
2N/A
2N/Avoid
2N/Amp_gcd(MINT *a, MINT *b, MINT *c)
2N/A{
2N/A MINT x, y, z, w;
2N/A
2N/A x.len = y.len = z.len = w.len = 0;
2N/A _mp_move(a, &x);
2N/A _mp_move(b, &y);
2N/A while (y.len != 0) {
2N/A mp_mdiv(&x, &y, &w, &z);
2N/A _mp_move(&y, &x);
2N/A _mp_move(&z, &y);
2N/A }
2N/A _mp_move(&x, c);
2N/A _mp_xfree(&x);
2N/A _mp_xfree(&y);
2N/A _mp_xfree(&z);
2N/A _mp_xfree(&w);
2N/A}
2N/A
2N/Avoid
2N/Amp_invert(MINT *x1, MINT *x0, MINT *c)
2N/A{
2N/A MINT u2, u3;
2N/A MINT v2, v3;
2N/A MINT zero;
2N/A MINT q, r;
2N/A MINT t;
2N/A MINT x0_prime;
2N/A static MINT *one = NULL;
2N/A
2N/A /*
2N/A * Minimize calls to allocators. Don't use pointers for local
2N/A * variables, for the one "initialized" multiple precision
2N/A * variable, do it just once.
2N/A */
2N/A if (one == NULL)
2N/A one = mp_itom(1);
2N/A
2N/A zero.len = q.len = r.len = t.len = 0;
2N/A
2N/A x0_prime.len = u2.len = u3.len = 0;
2N/A _mp_move(x0, &u3);
2N/A _mp_move(x0, &x0_prime);
2N/A
2N/A v2.len = v3.len = 0;
2N/A _mp_move(one, &v2);
2N/A _mp_move(x1, &v3);
2N/A
2N/A while (mp_mcmp(&v3, &zero) != 0) {
2N/A /* invariant: x0*u1 + x1*u2 = u3 */
2N/A /* invariant: x0*v1 + x2*v2 = v3 */
2N/A /* invariant: x(n+1) = x(n-1) % x(n) */
2N/A mp_mdiv(&u3, &v3, &q, &r);
2N/A _mp_move(&v3, &u3);
2N/A _mp_move(&r, &v3);
2N/A
2N/A mp_mult(&q, &v2, &t);
2N/A mp_msub(&u2, &t, &t);
2N/A _mp_move(&v2, &u2);
2N/A _mp_move(&t, &v2);
2N/A }
2N/A /* now x0*u1 + x1*u2 == 1, therefore, (u2*x1) % x0 == 1 */
2N/A _mp_move(&u2, c);
2N/A if (mp_mcmp(c, &zero) < 0) {
2N/A mp_madd(&x0_prime, c, c);
2N/A }
2N/A _mp_xfree(&zero);
2N/A _mp_xfree(&v2);
2N/A _mp_xfree(&v3);
2N/A _mp_xfree(&u2);
2N/A _mp_xfree(&u3);
2N/A _mp_xfree(&q);
2N/A _mp_xfree(&r);
2N/A _mp_xfree(&t);
2N/A}