lib/silcmath/modinv.c

   1 /*
   2
   3   modinv.h
   4
   5   Author: Pekka Riikonen <priikone@poseidon.pspt.fi>
   6
   7   Copyright (C) 1997 - 2000 Pekka Riikonen
   8
   9   This program is free software; you can redistribute it and/or modify
  10   it under the terms of the GNU General Public License as published by
  11   the Free Software Foundation; either version 2 of the License, or
  12   (at your option) any later version.
  13
  14   This program is distributed in the hope that it will be useful,
  15   but WITHOUT ANY WARRANTY; without even the implied warranty of
  16   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17   GNU General Public License for more details.
  18
  19 */
  20 /* $Id$ */
  21
  22 #include "silcincludes.h"
  23
  24 /* Table for finding multiplicative inverse */
  25 typedef struct {
  26   SilcMPInt x;
  27 } ModInv;
  28
  29 #define plus1   (i == 2 ? 0 : i + 1)
  30 #define minus1  (i == 0 ? 2 : i - 1)
  31
  32 /* Find multiplicative inverse using Euclid's extended algorithm.
  33    Computes inverse such that a * inv mod n = 1, where 0 < a < n.
  34    Algorithm goes like this:
  35
  36    g(0) = n    v(0) = 0
  37    g(1) = a    v(1) = 1
  38
  39    y = g(i-1) / g(i)
  40    g(i+1) = g(i-1) - y * g(i) = g(i)-1 mod g(i)
  41    v(i+1) = v(i-1) - y * v(i)
  42
  43    do until g(i) = 0, then inverse = v(i-1). If inverse is negative then n,
  44    is added to inverse making it positive again. (Sometimes the algorithm
  45    has a variable u defined too and it behaves just like v, except that
  46    initalize values are swapped (i.e. u(0) = 1, u(1) = 0). However, u is
  47    not needed by the algorithm so it does not have to be included.)
  48 */
  49
  50 void silc_mp_modinv(SilcMPInt *inv, SilcMPInt *a, SilcMPInt *n)
  51 {
  52   int i;
  53   SilcMPInt y;
  54   SilcMPInt x;
  55
  56   ModInv g[3];
  57   ModInv v[3];
  58
  59   /* init MP vars */
  60   silc_mp_init(&y);
  61   silc_mp_init(&x);
  62   silc_mp_init(&v[0].x);
  63   silc_mp_init(&v[1].x);
  64   silc_mp_set_ui(&v[0].x, 0L);          /* v(0) = 0 */
  65   silc_mp_set_ui(&v[1].x, 1L);          /* v(1) = 1 */
  66   silc_mp_init(&v[2].x);
  67   silc_mp_init(&g[0].x);
  68   silc_mp_init(&g[1].x);
  69   silc_mp_set(&g[0].x, n);              /* g(0) = n */
  70   silc_mp_set(&g[1].x, a);              /* g(1) = a */
  71   silc_mp_init(&g[2].x);
  72
  73   i = 1;
  74   while(silc_mp_cmp_ui(&g[i].x, 0) != 0) {
  75     silc_mp_div(&y, &g[minus1].x, &g[i].x);     /* y = n / a */
  76     silc_mp_mod(&g[plus1].x, &g[minus1].x, &g[i].x); /* remainder */
  77     silc_mp_mul(&x, &y, &v[i].x);
  78     silc_mp_set(&v[plus1].x, &v[minus1].x);
  79     silc_mp_sub(&v[plus1].x, &v[plus1].x, &x);
  80     i = plus1;
  81   }
  82
  83   /* set the inverse */
  84   silc_mp_set(inv, &v[minus1].x);
  85
  86   /* if inverse is negative, add n to inverse */
  87   if (silc_mp_cmp_ui(inv, 0) < 0)
  88     silc_mp_add(inv, inv, n);
  89
  90   /* clear the vars */
  91   memset(&g, 0, sizeof(g));
  92   memset(&v, 0, sizeof(v));
  93   silc_mp_uninit(&y);
  94   silc_mp_uninit(&x);
  95   silc_mp_uninit(&g[0].x);
  96   silc_mp_uninit(&g[1].x);
  97   silc_mp_uninit(&g[2].x);
  98   silc_mp_uninit(&v[0].x);
  99   silc_mp_uninit(&v[1].x);
 100   silc_mp_uninit(&v[2].x);
 101 }