mpz_multi_invert.cc

Go to the documentation of this file.

// compute a[i]^(-1) mod n, i=0..k-1 
// written by Thorsten Reinecke, 1999-11-21
// last change: 2003-11-10

/*

  Invertierung nach folgendem Prinzip:
   (1/x)*(1/y)=1/(x*y)  ["Additivität": f(x)*f(y)=f(x*y)]
   (1/x)=y*(1/(x*y))    ["???": f(x)=y*f(x*y)]
   also: anstatt 1/x und 1/y direkt zu berechnen,
   berechnen wir stattdessen z=1/(x*y) und erhalten 1/x=y*z, 1/y=x*z.
   Mehrere Kehrwerte können also zusammen schneller berechnet werden:
  
  a,b sind Felder der Größe k (a[0]..a[k-1], b[0]..b[k-1]),
  in a stehen die Argumente, in b anschließend die Kehrwerte,
  a und b müssen disjunkten Speicherplatz belegen;
  das Resultat ist undefiniert, falls eines der zu invertierenden Elemente
  nicht invertierbar ist (mithin also das Produkt der zu invertierenden
  Elemente nicht invertierbar ist).
 
  Algorithmus:

  multi_invert(var array of nat b, const var array of nat a, int k, nat n)
    {
      b[0]:=a[0];
      for i from 1 to k-1 do b[i]:=b[i-1]*a[i] mod n; od;
      x:=1/b[k-1] mod n; // x = 1/(a[0]*a[1]*...*a[k-1]) mod n
      for i from k-1 downto 1 do b[i]:=b[i-1]*x mod n; x:=x*a[i]; od;
      b[0]:=x;
    }

  Kosten: 1 Invertierung, 3*(k-1) modulare Multiplikationen
          anstelle von k Invertierungen.
  
  mögliche Anwendungsprobleme:
     Für die Anwendung müssen die k zu invertierenden Zahlen
     bereits bekannt sein. In vielen Algorithmen gilt leider
     a[i] = f(1/a[j]) für i>j, f ist irgendeine nichttriviale Funktion,
     d.h. die zu invertierenden Elemente sind voneinander abhängig und
     die Abhängigkeit kann ohne Kehrwertberechnung nicht beseitigt werden.

  Referenz:
    Peter L. Montgomery, "Speeding the Pollard and Elliptic Curve Methods
     of Factorization", Mathematics Of Computation, Vol. 48 (177), January 1987,
     pages 243-264, (hier: insbesondere Seite 260)  


 Bemerkung am Rande:
  Wie bei so vielen Algorithmen, auf die man sicher auch selbst gekommen wäre,
  wenn man sich "nur" das richtige Problem gestellt und konzentriert
  nachgedacht hätte, so findet sich auch dieser Algorithmus in der Literatur.
  Was mal wieder beweist, dass die eigentliche Leistung oftmals nicht in einer
  neuen algorithmischen Idee an sich besteht, sondern vielmehr in der Anwendung
  hinlänglich bekannter Eigenschaften auf eine neue, alternative
  Problemstellung gesucht werden kann.
  -> Wer will schon viele Zahlen auf einmal invertieren? Aber wenn man
  nun auf den Gedanken kommt, dass dies (tatsächlich) effizient möglich ist?
  Das müßte sich doch irgendwie anwenden lassen! -> Und wenn das der Fall ist,
  dann will man viele Zahlen auf einmal invertieren! -> Diesen Zusammenhang,
  der eine Problemtransformation ermöglicht, erkannt (und bekannt gemacht) zu
  haben, darin besteht die eigentliche wissenschaftliche Leistung...
  In dem oben angegebenen, sehr lesenswerten Paper hagelt es nur so von
  solchen (teilweise zitierten) "Problemtransformationen"...

*/



// Implementierung:

#include <gmp.h>

int mpz_multi_invert(mpz_t *b, const mpz_t *a, int k, const mpz_t n)
{
 /* 
    a and b MUST be disjoint arrays!
    returns 1, if inverse of (a[0]*...*a[k-1]) (mod n) exist
               and sets b[i]=a[i]^(-1) mod n, for i=0...k-1
    returns 0, otherwise
 */
  mpz_set(b[0],a[0]);
  for (int i=1; i<k; ++i)
   { 
     mpz_mul(b[i],b[i-1],a[i]); mpz_mod(b[i],b[i],n);
   }
  mpz_t x;
  mpz_init(x);
  if (mpz_invert(x,b[k-1],n)==0) { mpz_clear(x); return 0; } // inverse does not exist
  for (int i=k-1; i>0; --i)
   {
     mpz_mul(b[i],b[i-1],x); mpz_mod(b[i],b[i],n);
     mpz_mul(x,x,a[i]); mpz_mod(x,x,n);
   }
  mpz_set(b[0],x);
  mpz_clear(x);
  return 1; // inverse exist
}


#if 0

#include <iostream>

int main()
{
  using std::cout;
  using std::endl;
  const unsigned int p = 23;
  mpz_t n, x, a[p], b[p];
  mpz_init(x);
  mpz_init(n); mpz_set_ui(n,p);
  for (unsigned int i=0; i<p; ++i)
   {
     mpz_init(a[i]);
     mpz_set_ui(a[i],i);
     mpz_init(b[i]);
   }
  if (mpz_multi_invert(&b[1],&a[1],p-1,n)==0) cout << "Inverse doesn't exist!" << endl;
  for (unsigned int i=1; i<p; ++i)
   {
     cout << i << "^(-1) mod " << p << " = ";
     mpz_out_str(stdout,10,b[i]);
     cout << " --check--> ";
     mpz_mul_ui(x,b[i],i); mpz_mod(x,x,n);
     mpz_out_str(stdout,10,x);
     cout << endl;
   }
}
#endif

---