Zufallszahlen

SideWinder

Hi!

Background:
Ich kann mit einer mir vorliegenden rand()-Funktion Zufallszahlen zwischen [0, 65535) erzeugen Ich benötige aber 4-Byte-Unsigned-Zufallszahlen, also bis zum max. Intervall [0, 4294967295] (Achtung: hier inklusive des letzten Integers im Intervall!)

Ich will nun mit wiederholtem Aufrufen der rand()-Funktion die notwendigen 4-Byte-Zufallszahlen erzeugen. Da das Bit-Pattern 0xFFFF offenbar nicht erzeugt wird, kann ich nur 15 zufällige Bits mit meinem rand() kreieren. Für 32bit-Zahlen benötige ich also zumindest 3 Anläufe.

Gibt's dafür einen guten Algorithmus?

Ich habe zwei Versuche gestartet (basiernd darauf mir für jedes Byte eine Zufallszahl zu erzeugen und das über Bit-Shfiting zusammenzubauen, aber wenn das höhere Byte == höheresByteErlaubtesMax ist, darf man weiter unten allerhöchstens nur niedrigeresByteErlaubtesMax haben, falls nicht, darf man 0-BYTE_MAX haben - es ist hier baer serh schwer die Gleichverteiltheit der Zahlen beizubehalten).

Any ideas?

SideWinder

Hm, geht das mathematisch überhaupt so einfach...

Nehmen wir an ich kann bis zu 0b11 zufällig erzeugen lassen (00, 01, 10 oder 11) und ich möchte aber gerne die Zahlen 0-5 (0b000 - 0b101) zufällig erzeugen. Ich kann es nicht durch 2 HalfBytes erzeugen. Erzeuge ich zuerst das höhere HalfByte und dann abhängig von dessen Ausgang (0 -> [00-11], 1 -> [00-01]) ein niedriges HalfByte erhalte ich 1/4 und 1/8 Chancen statt gleichverteilte 1/6 Chancen. Genauso wenn ich es umgekehrt mache....

D.h. mein Multi-Wort-Ansatz wird nicht funktionieren.

Andere Ideen wie das richtig funktioniert?

SideWinder

Hmm, muss mir morgen mal https://github.com/dials/boost/blob/master/boost/random/uniform_int_distribution.hpp reinziehen...

Christoph

Ich würde auch empfehlen, was fertiges zu nehmen wie uniform_int_distribution, denn das selber zu machen klingt nur nach einer unnötigen und schwer zu findenden Fehlerquelle, vor allem wenn der Fehler subtil genug ist.

Wenn du 3 Mal rand() aufrufst und jeder Aufruf dir 15 zufällige Bits gibt (indem du jeweils auf 15 Bit runtermaskierst), kannst du die aber auch einfach alle verketten zu 45 Bits insgesamt und dann auf 32 Bit maskieren. Oder überseh ich was?

Ich würde aber egal was du machst sehr sicher sein wollen, dass deine rand()-Funktion gut ist und sogar in Frage stellen, warum du eine eigene rand()-Funktion benutzen möchtest, anstatt was fertiges. Hängt aber natürlich auch davon ab, ob es für Kryptografie oder für ein Spiel oder sowas ist (bei ersterem wär eine eigene rand()-Funktion nämlich ein Warnzeichen).

SideWinder

@christoph sagte in Zufallszahlen:

Wenn du 3 Mal rand() aufrufst und jeder Aufruf dir 15 zufällige Bits gibt (indem du jeweils auf 15 Bit runtermaskierst), kannst du die aber auch einfach alle verketten zu 45 Bits insgesamt und dann auf 32 Bit maskieren. Oder überseh ich was?

Naja, ja. Ich hab nicht dazugesagt, dass ich manchmal auch andere Intervalle benötige. Dann geht das nicht mehr. Aber auch hier hilft mir die uniform_int_distribution von Boost weiter.

Ich würde aber egal was du machst sehr sicher sein wollen, dass deine rand()-Funktion gut ist und sogar in Frage stellen, warum du eine eigene rand()-Funktion benutzen möchtest, anstatt was fertiges.

Mag ich ja nicht, aber die, die ich hab liefert zu kleine Zahlen. Aber mit uniform_int_distribution löst sich mein Problem.

Falls jemand noch meine Rückfragen zu uniform_int_distribution beantworten kann wäre mir sehr geholfen: https://github.com/boostorg/random/issues/43 Ich muss das nämlich in eine andere Sprache übersetzen.

MfG SideWinder

Christoph

@sidewinder Mir ist auch gerade aufgefallen, dass mein Vorschlag mit dem Maskieren auf 15 Bits gar nicht funktioniert, weil es auch die Gleichverteilung kaputt macht.

Dein rand() so oft aufrufen bis eine Zahl mit <= 15 Bits rauskommt, sollte gleichverteilte 15-Bit-Zahlen produzieren, ist aber vielleicht nicht optimal.

hustbaer

@SideWinder Ein Generator wie xorshift ist gut und sehr einfach selbst zu implementieren. Dazu eine Funktion ala uniform_int_distribution und ... tadaa.

Die Kernidee eines Zufallsgenerators sieht in der ur Mimik etwa so aus:

 ... Code ... 

char zufall(short bis)
{
  static unsigned char Z[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20} ;  
  static short i=2;

  if((i++) > 20)i=2;  
    
  Z[i] = (Z[i-1]+Z[i-2]) % (bis); 

       
   return( Z[i] );
}





void main()
{

  while(getch() != 32)
  {
    printf("wert = %d\n",zufall(2));  
  }

}

Wenn man es sehr schnell benötigt kann man mit Asm weiter kommen:

 ... Code ... 
short random();
#pragma aux random            = "mov esi,offset seed    "\
                                "mov     eax,1107030247 "\
                                "mul     dword ptr [esi]"\
                                "add     eax,97177      "\
                                "mov     [esi],eax      "\
                                "shr     eax,15         "\
                                "and     ax,8191        "\
                                "sub     ax,4096        "\
                                value    [AX]            \

Oder etwas ausgeprägter mit timer:

 ... Code ... 


;±±±±±±±±±±±±±±±± rand ±±±±±±±±±±±±±±±±
;returns a random value in range -4096..4095
rand    PROC NEAR
        mov     eax,1107030247
        mul     ds:seed
        add     eax,97177
        mov     ds:seed,eax
        shr     eax,15
        and     ax,8191
        sub     ax,4096
;size optimizatin, some code moved from after all rand calls
        add     bx,2
        mov     ds:[bx],ax
        ret
rand    ENDP

;±±±±±±±±±±±±±±±± timer ±±±±±±±±±±±±±±±±
inittimer PROC NEAR
        mov     eax,fs:[8*4]
        mov     ds:oldint8,eax
        mov     ax,cs
        shl     eax,16
        mov     ax,OFFSET intti8
        mov     dx,17000 ;70hz
        jmp     @@1
deinittimer:
        mov     eax,ds:oldint8
        xor     dx,dx
@@1:    cli
        mov     fs:[8*4],eax
        mov     al,036h
        out     43h,al
        mov     al,dl
        out     40h,al
        mov     al,dh
        out     40h,al
        sti
        ret
inittimer ENDP

intti8  PROC FAR ;timer interrupt
        push    ax
        mov     al,20h
        out     20h,al
        inc     cs:framecounter
        pop     ax
        iret
intti8  ENDP

Der beste Zufallsgenerator ist jedoch das weiße rauschen aus einem DVBT-Stick

indem man Stichproben aus den Daten abgreift.

Grüße
Karsten

rustyoldguy

Hallo SideWinder!

Vielleicht teilst du die ganze Sache auf. "Normale" rand-Funktionen gehen von 0 bis 65535.
Bei einem 64-Bit-System verarbeitet.
define RAND_MAX /implementation defined/
ist bei 32767

4294967295 / 65535 = 65537

Du müsstest also das Ganze in 65537 Felder unterteilen. Da aber eine Null auch vorkommen kann(sollte)
sind das genau 65536^2.
also zuerst mit rand eine Zahl zwischen 0 und 65535 dann noch eine und beide Multiplizieren.
Praktisch wie früher in DOS, also segment und Offset, wenn du weißt was ich meine.
Auf einigen Seiten
http://man7.org/linux/man-pages/man3/rand.3.html
ist das aber anders zu lesen. Einige Systeme arbeiten bis etwa 199506.
Da ich erst kürzlich meine CD-Sammlung entsorgt habe, hier mal ein Code:

/************************************************************************
 This random number generator originally appeared in "Toward a Universal
 Random Number Generator" by George Marsaglia and Arif Zaman.
 Florida State University Report: FSU-SCRI-87-50 (1987)

 It was later modified by F. James and published in "A Review of Pseudo-
 random Number Generators"

 Converted from FORTRAN to C by Phil Linttell, James F. Hickling
 Management Consultants Ltd, Aug. 14, 1989.

 THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
       (However, a newly discovered technique can yield
         a period of 10^600. But that is still in the development stage.)

 It passes ALL of the tests for random number generators and has a period
   of 2^144, is completely portable (gives bit identical results on all
   machines with at least 24-bit mantissas in the floating point
   representation).

 The algorithm is a combination of a Fibonacci sequence (with lags of 97
   and 33, and operation "subtraction plus one, modulo one") and an
   "arithmetic sequence" (using subtraction).

 On a Vax 11/780, this random number generator can produce a number in
    13 microseconds.
************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define TRUE    1
#define FALSE   0

float u[97], c, cd, cm;
int i97, j97, test;

int rmarin(int ij, int kl);
int ranmar(float rvec[], int len);
void sow( int *seed1, int *seed2 );


int main()
{

        float temp[100];
        int i;
        int ij, kl, len;

        /*These are the seeds needed to produce the test case results*/
        sow( &ij, &kl);
        //ij = 1802;
        //kl = 9373;

        /*Do the initialization*/

        if (1 == rmarin(ij,kl))
                return 1;

        /*Generate 20000 random numbers*/

        len = 100;
        for ( i=0; i<=199 ; i++)
                if (1 == ranmar(temp, len))
                        return 1;

        /*If the random number generator is working properly,
          the next six random numbers should be:

            6533892.0  14220222.0   7275067.0
            6172232.0   8354498.0  10633180.0
        */

        len = 6;
        if (1 == ranmar(temp, len))
                return 1;

        for ( i=0; i<=5; i++)
                printf("%12.1f\n",4096.0*4096.0*temp[i]);

  return 0;
}


/************************************************************************
 This is the initialization routine for the random number generator RANMAR()
 NOTE: The seed variables can have values between:    0 <= IJ <= 31328
                                                      0 <= KL <= 30081
 The random number sequences created by these two seeds are of sufficient
 length to complete an entire calculation with. For example, if several
 different groups are working on different parts of the same calculation,
 each group could be assigned its own IJ seed. This would leave each group
 with 30000 choices for the second seed. That is to say, this random
 number generator can create 900 million different subsequences -- with
 each subsequence having a length of approximately 10^30.

 Use IJ = 1802 & KL = 9373 to test the random number generator. The
 subroutine RANMAR should be used to generate 20000 random numbers.
 Then display the next six random numbers generated multiplied by 4096*4096
 If the random number generator is working properly, the random numbers
 should be:
           6533892.0  14220222.0   7275067.0
           6172232.0   8354498.0  10633180.0
************************************************************************/

int rmarin(int ij, int kl)
{

        float s, t;
        int i, j, k, l, m;
        int ii, jj;

        /* Change FALSE to TRUE in the next statement to test the
           random routine.*/

        test = TRUE;

        if ( ( ij < 0 || ij > 31328 ) ||
                ( kl < 0 || kl > 30081 ) )
        {
                printf ("RMARIN: The first random number seed must have a "
                        "value between 0 and 31328\n");
                printf ("        The second random number seed must have a "
                        "value between 0 and 30081");
                return 1;
        }

        i = fmod(ij/177.0, 177.0) + 2;
        j = fmod(ij      , 177.0) + 2;
        k = fmod(kl/169.0, 178.0) + 1;
        l = fmod(kl      , 169.0);

        for ( ii=0; ii<=96; ii++ )
        {
                s = 0.0;
                t = 0.5;
                for ( jj=0; jj<=23; jj++ )
                {
                        m = fmod( fmod(i*j,179.0)*k , 179.0 );
                        i = j;
                        j = k;
                        k = m;
                        l = fmod( 53.0*l+1.0 , 169.0 );
                        if ( fmod(l*m,64.0) >= 32)
                                s = s + t;
                        t = 0.5 * t;
                }
                u[ii] = s;
        }

        c  =   362436.0 / 16777216.0;
        cd =  7654321.0 / 16777216.0;
        cm = 16777213.0 / 16777216.0;

        i97 = 96;
        j97 = 32;

        test = TRUE;

        return 0;
}

int ranmar(float rvec[], int len)
{
        float uni;
        int ivec;

        if ( !test )
        {
                printf ("RANMAR: Call the initialization routine (RMARIN) "
                        "before calling RANMAR.\n");
                return 1;
        }

        for ( ivec=0; ivec < len; ivec++)
        {
                uni = u[i97] - u[j97];
                if ( uni < 0.0 )
                        uni = uni + 1.0;
                u[i97] = uni;
                i97--;
                if ( i97 < 0 )
                        i97 = 96;
                j97--;
                if ( j97 < 0 )
                        j97 = 96;
                c = c - cd;
                if ( c < 0.0 )
                        c = c + cm;
                uni = uni - c;
                if ( uni < 0.0 )
                        uni = uni + 1.0;
                rvec[ivec] = uni;
        }
        return 0;
}

/* I use the following procedure in TC to generate seeds:

  The sow() procedure calculates two seeds for use with the random number
  generator from the system clock.  I decided how to do this myself, and
  I am sure that there must be better ways to select seeds; hopefully,
  however, this is good enough.  The first seed is calculated from the values
  for second, minute, hour, and year-day; weighted with the second most
  significant and year-day least significant.  The second seed weights the
  values in reverse.
*/

void sow( int *seed1, int *seed2 )
{
        struct tm *tm_now;
        float s_sig, s_insig, maxs_sig, maxs_insig;
        long secs_now;
        int s, m, h, d, s1, s2;

        time(&secs_now);
        tm_now = localtime(&secs_now);

        s = tm_now->tm_sec + 1;
        m = tm_now->tm_min + 1;
        h = tm_now->tm_hour + 1;
        d = tm_now->tm_yday + 1;

        maxs_sig   = 60.0 + 60.0/60.0 + 24.0/60.0/60.0 + 366.0/24.0/60.0/60.0;
        maxs_insig = 60.0 + 60.0*60.0 + 24.0*60.0*60.0 + 366.0*24.0*60.0*60.0;

        s_sig      = s + m/60.0 + h/60.0/60.0 + d/24.0/60.0/60.0;
        s_insig    = s + m*60.0 + h*60.0*60.0 + d*24.0*60.0*60.0;

        s1 = s_sig   / maxs_sig   * 31328.0;
        s2 = s_insig / maxs_insig * 30081.0;

        *seed1 = s1;
        *seed2 = s2;
}

Vielleicht kannst du den Code ja umarbeiten.

rustyoldguy

Genau genommen meine ich so etwas

int randinitcount = 0;

void RANDOMIZE(void)
{
long bwer;
srand( (unsigned)(time(&bwer) % 20));
if (randinitcount < 2)
 randinitcount++;
}

int getrandom( int min, int max)
{
int zuza, zuza2, zuza3, rewer;
zuza = rand();
zuza2 = (max - min + 1);
zuza3 = zuza % zuza2;
rewer = zuza3 + min;
return rewer;
}

long getlongrand(long min, long max)
{
long rewer = 0, rest, bereichsfeld;
int bereich, limit = 32767;

bereichsfeld = max / 32767;

bereich = getrandom((int)min, (int)bereichsfeld);
if (bereich > 1) rewer = 32767 * (long)bereich;
rest = max - rewer;
if (rest < 32767) limit = (int)rest;

rewer += (long)getrandom(0, limit);

return(rewer);
}

Christoph

Nur um das nicht so stehenzulassen: Was gutes aus dem Standard wie Mersenne twister (mt19937) und uniform_int_distribution nehmen ist schon sinnvoll, außer man hat sehr spezielle Anforderungen und kennt die Mathematik dahinter.

Die ganzen Funktionen, die hier in den letzten 3 Postings gepostet wurden, haben Probleme mit der dem Algorithmus und der Gleichverteilung. Modulo funktioniert nicht, wenn man eine Gleichverteilung haben möchte.

Swordfish

@KahnSoft sagte in Zufallszahlen:

Die Kernidee eines Zufallsgenerators sieht in der ur Mimik etwa so aus:

Mhm. Das ist der erste rng den ich sehe der ein statisches Array braucht um zu funktionieren. Du könntest wiklich mal aufhören geistigen Dünnpfiff zu verbreiten. Ein paar Iterationen und dein "Zuffalsgenerator" erledigt sich zu

für zufall(100).

Edit: Die Zitierfunktion scheint jetzt noch kaputter als zu Anfang. But what shells, death is near.

uint32_t num;
uint32_t i;

uint8_t *ptr;

ptr = (uint8_t *) &num;

for(i = 0; i < sizeof(num); i++)
{
     *ptr = rand();
     ptr++;
}

unter unix macht man das übrigens eigentlich mit

uint32_t num;

int file = open("/dev/urandom", O_RDONLY);
read(file, &num, sizeof(num));
close(file);

unter windows wirds etwas komplizierter, geht aber auch.