Diskussion:Mersenne-Twister

Letzter Kommentar: vor 2 Jahren von Megatherium in Abschnitt Bessere Initialisierung
Diese Diskussionsseite dient dazu, Verbesserungen am Artikel „Mersenne-Twister“ zu besprechen. Persönliche Betrachtungen zum Thema gehören nicht hierher. Für allgemeine Wissensfragen gibt es die Auskunft.

Füge neue Diskussionsthemen unten an:

Klicke auf Abschnitt hinzufügen, um ein neues Diskussionsthema zu beginnen.
Archiv
Wie wird ein Archiv angelegt?
Auf dieser Seite werden Abschnitte ab Überschriftenebene 2 automatisch archiviert, die seit 10 Tagen mit dem Baustein {{Erledigt|1=--~~~~}} versehen sind.

Pseudo-Code statt C? Bearbeiten

Schön wäre ein Link auf den Original-Code, und statt diesem auf der Tastatur herumgewälzten Gürteltier ein verständlicher Pseudo-Code. Oder von mir aus Pascal, das versteht jeder. -- Maxus96 00:37, 9. Jun. 2009 (CEST)Beantworten

Der englische Wiki-Artikel hat Pseudocode. (nicht signierter Beitrag von 93.237.23.57 (Diskussion) 06:42, 7. Dez. 2020 (CET))Beantworten

Du kannst den Code ja von mir aus in Pascal umschreiben, wenn du meinst dass das der Lesbarkeit gut tut... er ist IMHO schon vereinfacht, zu sauberem Stil gehört z.B. dass man nach Möglichkeit statische Variablen vermeidet, und statt den #defines kann man auch const int, oder was weiß ich, nutzen... --178.15.218.194 23:11, 27. Jul. 2013 (CEST)Beantworten
Ist das besser lesbar? Ist weiterhin in (sogar kompilierbarem!) C geschrieben, aber in Unterfunktionen zerlegt. --202.79.203.59 21:29, 14. Dez. 2013 (CET)Beantworten

MT 19937

#include <stdint.h>

#define N     624
#define M     397

/*
 * Initialisiere den Vektor mit Pseudozufallszahlen.
 * Siehe Donald E. Knuth: The Art of Computer Programming, Band 2, 3. Ausgabe, S. 106 für geeignete Werte für den Faktor mult
 *
 * 2002/01/09 modifiziert von Makoto Matsumoto:
 *   In der vorhergehenden Version die MSBs des Seeds haben auch nur die MSBs des Arrays beeinflusst.
 */

static void
mersenne_twister_vector_init (uint32_t* const p, const int len)
{
  const uint32_t  mult = 1812433253ul;
  uint32_t        seed = 5489ul;
  int             i;

  for (i = 0; i < len; i++)
  {
    p[i] = seed;
    seed = mult * (seed ^ (seed >> 30)) + (i+1);
  }
}

/*
 * Berechne neuen Zustandsvektor aus altem Zustandsvektor
 * Dies erfolgt periodisch alle N=624 Aufrufe von mersenne_twister().
 *
 * Der folgende Code stellt eine optimierte Version von
 *   ...
 *   for (i = 0; i < N; i++)
 *     Proc3 (p[i], p[(i+1) % N], p[(i+M) % N]);
 *   ...
 * mit
 *   void Proc3 (uint32_t &a0, uint32_t a1, uint32_t aM)
 *   {
 *     a0 = aM ^ (((a0 & 0x80000000) | (a1 & 0x7FFFFFFF)) >> 1) ^ A[a1 & 1];
 *   }
 * dar, in der die (zeitintensiven) Modulo N-Operationen vermieden werden.
 */

static void
mersenne_twister_vector_update (uint32_t* const p)
{
  static const uint32_t  A[2] = { 0, 0x9908B0DF };
  int                    i = 0;

  for (; i < N-M; i++)
    p[i] = p[i+(M  )] ^ (((p[i  ] & 0x80000000) | (p[i+1] & 0x7FFFFFFF)) >> 1) ^ A[p[i+1] & 1];
  for (; i < N-1; i++)
    p[i] = p[i+(M-N)] ^ (((p[i  ] & 0x80000000) | (p[i+1] & 0x7FFFFFFF)) >> 1) ^ A[p[i+1] & 1];
  p[N-1] = p[M-1]     ^ (((p[N-1] & 0x80000000) | (p[0  ] & 0x7FFFFFFF)) >> 1) ^ A[p[0  ] & 1];
}

/*
 * Die eigentliche Funktion:
 * - beim 1. Aufruf wird der Vektor mittels mersenne_twister_vector_init() initialisiert.
 * - beim 1., 625., 1249., ... Aufruf wird ein neuer Vektor mittels mersenne_twister_vector_update berechnet.
 * - bei jedem Aufruf wird eine Zufallszahl aus dem Vektor bestimmt, der noch einem Tempering unterzogen wird.
 */

uint32_t
mersenne_twister ()
{
  static uint32_t  vektor [N];  /* Zustandsvektor */
  static int       idx = N+1;   /* Auslese-Index; idx=N: neuer Vektor muß berechnet werden, idx>N: Vektor muß überhaupt erst mal initialisiert werden */
  uint32_t         e;

  if (idx >= N)
  {
    if (idx > N)
      mersenne_twister_vector_init (vektor, N);
    mersenne_twister_vector_update (vektor);
    idx = 0;
  }

  e  = vektor [idx++];
  e ^= (e >> 11);             /* Tempering */
  e ^= (e <<  7) & 0x9D2C5680;
  e ^= (e << 15) & 0xEFC60000;
  e ^= (e >> 18);

  return e;
}

#undef N
#undef M

TT 800

#include <stdint.h>

#define N  25
#define M   7

/* 
 * Initialisiere den Vektor mit Pseudozufallszahlen.
 */

static void 
TT800_vector_init (uint32_t* const p, const int len)
{
  const uint32_t  mult = 509845221ul;
  const uint32_t  add  =         3ul;
  uint32_t        seed1 =    9;
  uint32_t        seed2 = 3402;
  int       i;
  
  for (i = 0; i < len; i++) 
  {
    seed1  = seed1 * mult + add;
    seed2 *= seed2 + 1;
    p[i]   = seed2 + (seed1 >> 10);
  }
}

/* 
 * Berechne neuen Zustandsvektor aus altem Zustandsvektor
 * Dies erfolgt periodisch alle N=25 Aufrufe von TT800().
 *
 * Der folgende Code stellt eine optimierte Version von
 *   ...
 *   for (i = 0; i < N; i++) 
 *     Proc2 (p[i], p[(i+M) % N]);
 *   ...
 * mit
 *   void Proc2 (uint32_t &a0, uint32_t aM)
 *   {
 *     a0 = aM ^ (a0 >> 1) ^ A[a0 & 1];
 *   }
 * dar, in der die (zeitintensiven) Modulo N-Operationen vermieden werden.
 */

static void
TT800_vector_update (uint32_t* const p)
{
  static const uint32_t  A[2] = { 0, 0x8ebfd028 };
  int                    i = 0;

  for (; i < N-M; i++)
    p[i] = p[i+(M  )] ^ (p[i] >> 1) ^ A[p[i] & 1];
  for (; i < N  ; i++)
    p[i] = p[i+(M-N)] ^ (p[i] >> 1) ^ A[p[i] & 1];
}

/*
 * Die eigentliche Funktion:
 * - beim 1. Aufruf wird der Vektor mittels TT800_vector_init() initialisiert.
 * - beim 1., 26., 51., ... Aufruf wird ein neuer Vektor mittels TT800_vector_update berechnet.
 * - bei jedem Aufruf wird eine Zufallszahl aus dem Vektor bestimmt, der noch einem Tempering unterzogen wird.
 */
 
uint32_t 
TT800 () 
{
  static uint32_t  vektor [N];	/* Zustandsvektor */
  static int       idx = N+1;	/* Auslese-Index; idx=N: neuer Vektor muß berechnet werden, idx>N: Vektor muß überhaupt erst mal initialisiert werden */
  uint32_t         e;

   if (idx >= N) 
   {
     if (idx > N)
	    TT800_vector_init (vektor, N);
	 TT800_vector_update (vektor);
     idx = 0;
   }
 
   e  = vektor [idx++];
   e ^= (e <<  7) & 0x2b5b2500;             /* Tempering */
   e ^= (e << 15) & 0xdb8b0000;
   e ^= (e >> 16);
   return e;
}

#undef N
#undef M

Verallgemeinerung Bearbeiten

Wenn dieser Algorithmus auf Mersenne-Primzahlen aufbaut, dann sollte man ihn doch auch auf andere   - Primzahlen verallgemeinern können. Oder wenn nicht, warum? Leider erklärt der Artikel das Funktionsprinzip überhaupt nicht. -- Maxus96 19:15, 10. Jun. 2009 (CEST)Beantworten

Der Algorithmus basiert nicht auf Primzahlen. Es ist nur so, dass eine seiner Eigenschaften (die Periodenlänge) einem Primzahlwert entspricht: "Er hat die extrem lange Periode von 219937-1 (≈ 4,3·106001). Diese Periodenlänge erklärt auch den Namen des Algorithmus: Sie ist eine Mersenne-Primzahl, und einige Eigenschaften des Algorithmus resultieren daraus." (Zitat aus dem Artikel) 19937 ist übrigens ebenfalls eine Primzahl und sie entspricht 32·623+1. Natürlich könnte man versuchen, andere Primzahlen auf diese Weise zu finden, aber ob man daraus einen ähnlich guten Algorithmus konstruieren kann, lässt sich nicht aus der Definition des MT ableiten. —Wormbo 10:32, 15. Jan. 2011 (CET)Beantworten

Bessere Initialisierung Bearbeiten

static void
mersenne_twister_vector_init (uint32_t* const p, const int len)
{
  const uint32_t  mult = 1812433253;
  uint32_t        seed = 5489;
  int             i, j;

  for (i=j=0; i < len; i++)
  {
    seed *= mult;
    seed ^= seed >> 16 ^ i;
    p[j] = seed;
    j += 17; // Inkrement muss zu len teilerfremd sein
    if (j >= len) j -= len;
  }
}

--Megatherium (Diskussion) 17:51, 2. Jul. 2021 (CEST)Beantworten