usr/src/lib/libmp/common/gcd.c

   1 /*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
   2 /*        All Rights Reserved   */
   3
   4
   5 /*
   6  * Copyright (c) 1980 Regents of the University of California.
   7  * All rights reserved.  The Berkeley software License Agreement
   8  * specifies the terms and conditions for redistribution.
   9  */
  10 /*      Portions Copyright(c) 1988, Sun Microsystems Inc.       */
  11 /*      All Rights Reserved                                     */
  12
  13 /*
  14  * Copyright (c) 1997, by Sun Microsystems, Inc.
  15  * All rights reserved.
  16  */
  17
  18 #ident  "%Z%%M% %I%     %E% SMI"        /* SVr4.0 1.1   */
  19
  20 /* LINTLIBRARY */
  21
  22 #include <mp.h>
  23 #include "libmp.h"
  24 #include <sys/types.h>
  25
  26 void
  27 mp_gcd(MINT *a, MINT *b, MINT *c)
  28 {
  29         MINT x, y, z, w;
  30
  31         x.len = y.len = z.len = w.len = 0;
  32         _mp_move(a, &x);
  33         _mp_move(b, &y);
  34         while (y.len != 0) {
  35                 mp_mdiv(&x, &y, &w, &z);
  36                 _mp_move(&y, &x);
  37                 _mp_move(&z, &y);
  38         }
  39         _mp_move(&x, c);
  40         _mp_xfree(&x);
  41         _mp_xfree(&y);
  42         _mp_xfree(&z);
  43         _mp_xfree(&w);
  44 }
  45
  46 void
  47 mp_invert(MINT *x1, MINT *x0, MINT *c)
  48 {
  49         MINT u2, u3;
  50         MINT v2, v3;
  51         MINT zero;
  52         MINT q, r;
  53         MINT t;
  54         MINT x0_prime;
  55         static MINT *one = NULL;
  56
  57         /*
  58          * Minimize calls to allocators.  Don't use pointers for local
  59          * variables, for the one "initialized" multiple precision
  60          * variable, do it just once.
  61          */
  62         if (one == NULL)
  63                 one = mp_itom(1);
  64
  65         zero.len = q.len = r.len = t.len = 0;
  66
  67         x0_prime.len = u2.len = u3.len = 0;
  68         _mp_move(x0, &u3);
  69         _mp_move(x0, &x0_prime);
  70
  71         v2.len = v3.len = 0;
  72         _mp_move(one, &v2);
  73         _mp_move(x1, &v3);
  74
  75         while (mp_mcmp(&v3, &zero) != 0) {
  76                 /* invariant: x0*u1 + x1*u2 = u3 */
  77                 /* invariant: x0*v1 + x2*v2 = v3 */
  78                 /* invariant: x(n+1) = x(n-1) % x(n) */
  79                 mp_mdiv(&u3, &v3, &q, &r);
  80                 _mp_move(&v3, &u3);
  81                 _mp_move(&r, &v3);
  82
  83                 mp_mult(&q, &v2, &t);
  84                 mp_msub(&u2, &t, &t);
  85                 _mp_move(&v2, &u2);
  86                 _mp_move(&t, &v2);
  87         }
  88         /* now x0*u1 + x1*u2 == 1, therefore,  (u2*x1) % x0  == 1 */
  89         _mp_move(&u2, c);
  90         if (mp_mcmp(c, &zero) < 0) {
  91                 mp_madd(&x0_prime, c, c);
  92         }
  93         _mp_xfree(&zero);
  94         _mp_xfree(&v2);
  95         _mp_xfree(&v3);
  96         _mp_xfree(&u2);
  97         _mp_xfree(&u3);
  98         _mp_xfree(&q);
  99         _mp_xfree(&r);
 100         _mp_xfree(&t);
 101 }