Teil 1 - Gleichverteilung
Eins vorweg, das Thema Zufallszahlen ist ein sehr umfangreiches Thema mit teilweisen sehr komplexen Algorithmen. Ich werde daher nur einen Ansatz aufzeigen können. Allerdings muß ich eingestehen, daß dieser Artikel doch etwas groß ausgefallen ist als eigentlich gewollt war. Wer mehr Interesse an diesem Thema haben sollte, den verweise ich als erste Anlaufstelle an http://de.wikipedia.org/wiki/Zufallszahl. Dort bekommt man auch weiterführende Links.
So, los geht's...
Überall dort, wo Simulationen erzeugt werden, findet man ein scheinbar simples und harmloses mathematisches Werkzeug welches Zufallszahlen berechnet. Die Aufgabe dieses Werkzeugs ist die Ausgabe von Millionen oder Milliarden Zufallszahlen am Computer, die als Eingabewerte für Simulationsmodelle dienen. Diese Werkzeuge oder auch Zufallsgeneratoren genannt arbeiten mit deterministischen Algorithmen (definierte und reproduzierbare Zustände), weshalb man sie auch Pseudo-Zufallszahlengeneratoren nennt. Bei den meisten Computersprachen gibt es bereits ein Funktion zur Ermittlung von Zufallszahlen. In PHP z.B. heißt diese Funktion rand bzw. mt_rand (mt = "Mersenne Twister"). Da es in AutoLISP keine fertige Funktion zur Ermittlung von Zufallszahlen gibt, werden wir uns eben diese selbst programmieren.
Hier im 1. Tutorial zeige ich wie man mit AutoLISP sogenannte gleichverteilte Zufallszahlen erzeugt. Dazu verwenden wir den "Linear Kongruenz-Generator" (LCG = Linear Congruence Generator), der wohl am meisten eingesetzte Generator in der Stochastik (die Lehre der Häufigkeit und Wahrscheinlichkeit). Dieser Algorithmus erzeugt aus einer endlichen Menge "zuf&¨llige" Zahlen.
Achtung! Zufall auf dem Computer gibt es nicht, alles wird "nur" berechnet!
Formel für die lineare Kongruenzmethode: x = (a * x + c) % m
% m steht für Modulo und ist eine in der Informatik verbreitete Funktion, die den Rest aus der Division zweier Ganzzahlen angibt.
Hier nun die Formel in AutoLISP ausgedrückt : *x* = (rem (+ (* a *x*) c) m))
Die Formel sieht einfach aus und ist es auch, die "Kunst" ist aber geeignete Konstanten zu finden. Diese Konstanten sollten ganzzahlig (Integer) sein und gewissen mathematischen Regeln unterliegen, auf die ich weiter unten kurz eingehen werde. Fangen wir also ganz von vorne an und erzeugen eine Lisp Funktion. Dabei setzen wir folgende mathematische Regeln ein:
m ist eine Primzahl
a-1 ist ein vielfaches von m
c ist KEIN vielfaches von m
(defun rand-tmp (#len / a c lst)
(if (null *x*)
(setq *x* 0)
)
(repeat #len
(setq lst (cons (setq *x* (rem (+ (* 50 *x*) 9) 7)) lst))
)
lst
)
Rufen wir diese Funktion auf (dabei ist *x* noch nicht definiert) mit einer Länge von 16.
- Aufruf... (rand-tmp 16) => (4 2 0 5 3 1 6 4 2 0 5 3 1 6 4 2)
- Aufruf... (rand-tmp 16) => (1 6 4 2 0 5 3 1 6 4 2 0 5 3 1 6)
Bevor Sie weiterlesen, sollten Sie ein wenig mit den Zahlen spielen, ersetzen Sie in der Funktion die Zahlen 50, 9 und 7 mit anderen Ganzzahlen und ändern zusätzlich die Länge (Empfehlung für L&¨nge: ≥ Modulo-Konstante = obiges Beispiel 7). Beobachten Sie dabei welche Ergebnisse zurückgegeben werden.
Was können wir erkennen?
Wenn Sie ausreichend getestet haben, werden Sie feststellen, daß die Wahl der Konstanten SEHR entscheidend ist um eine brauchbare Liste von Zufallszahlen zu bekommen. Sicherlich haben Sie auch festgestellt, daß durch Änderung der Modulo-Konstante höhere Zahlen herauskommen. Daraus können wir nun schliessen, daß durch die Modulo-Konstante die "Obere Grenze" der Zahlen bestimmt wird, in unserem Fall also 7-1. Das heißt weiterhin, daß sich damit eine zyklische Wiederholung generiert. Im ersten Beispiel (*x* = nil) beginnt der Zyklus bei der Zahl 4 (+7) und im zweiten Beispiel bei der Zahl 1 (+7). Wir wissen nun, daß wir für eine lange Periode eine große Modulo-Konstante benötigen.
Um eine volle Periodenlänge (wir erinnern uns: zyklische Wiederholung) mit gleichverteilten Zufallszahlen zu erreichen, sollte folgendes gelten:
0 < m
a < m
c < m
Weiterhin ist ein auf der linearen Kongruenzmethode basierender Zufallszahlengenerator genau dann ein "full-period" Generator, wenn gilt:
- Der größte gemeinsame Teiler von m und c ist 1 => ggT(c,m) = 1
- Wenn q (seed) Primteiler von m ist, dann muss a-1 auch durch q teilbar sein
- Wenn 4 Teiler von m ist, dann muss 4 auch Teiler von a-1 sein
Da es bereits viele Lösungen hierzu gibt, habe ich mich für folgenden Konstanten enschieden, die von Doug Cooper und Mike Clancy stammen (Buch Oh! Pascal!;W.W. Norton & Co. 1982, 1987)
Modulus (m) = 2^16 = 65536
multiplier (a) = 25173
increment (c) = 13849
Seed (x) = <TimeStamp> -> (getvar "date")
Damit nicht immer dieselbe Zufallszahlenfolge erzeugt wird, verwendet man für den initialen Wert x einen Timestamp (auch Seed genannt).
Hierzu schreiben wir uns eine Hilfsfunktion.
;Erzeugung eines Seed Wertes als initialen Wert x
;für die Generierung von Zufallszahlen
(defun :T-Seed ()
(cond (*seed*)
((setq *seed* (getvar "date")))
)
)
Diese Funktion prüft mittels der Bedingungsfunktion cond, ob *seed* bereits gesetzt ist. Ist *seed* bereits gesetzt gibt die Funktion *seed* zurück,
andernfalls wird *seed* mittels der Systemvariablen "DATE" neu initialisiert. Die beiden Sternchen vor und nach seed kennzeichnen diese Variable als GLOBAL, dass heisst, dass die Variable auch außerhalb der Funktion :T-Seed verfügbar ist.
Nun unsere rand() Funktion:
(defun :m-rand (#len / retlst)
(repeat #len
(setq retlst (cons (setq *seed*
(fix (rem (+ (* 25173 (:T-Seed)) 13849) 65536))
)
retlst
)
)
)
)
Erzeugen wir nun damit eine Liste mit Zahlen:
(:M-RAND 16) => (25822 27313 59960 47267 26210 4645 61404 46007 20774 48601 43200 21771 2346 49613 6372 37023)
Hmmm, Sie sind immer noch da und wollen mehr wissen? Freut mich, daß ich Ihr Interesse geweckt habe. Unsere Funktion erzeugt nun gleichverteilte Zufallszahlen und zwar eine sogenannte "Volle Periode". Prüfen wir diesen Algorithmus:
(setq *TmpLst* (:m-rand (* 65536 2)))
Sehen wir uns die ersten Zahlen an:
=> (33139 64498 33141 4588 64903 25526 ...)
Prüfen wir nun die Periode indem wir an die 65536'ste Stelle unserer Liste springen und wir erkennen, daß exakt an dieser Stelle sich die Zahlen wiederholen:
(length(member (car *TmpLst*) (cdr *TmpLst*))) => 65536 => (33139 64498 33141 4588 64903 25526 ...)
Für unsere spätere Bibliotheksfunktion ist diese Funktion jedoch noch nicht optimal, denn wir wollen mit dieser Funktion nicht nur Ganzzahlen von 0 bis 65536, sondern auch Zufallszahlen in einem bestimmten Bereich erzeugen und diese Funktion soll später im 2. Tutorial dazu verwendet werden um Standardnormalverteilte Zufallszahlen zu erzeugen.
Für diese Aufgabe legen wir folgende Randbedingungen fest:
Funktion :m-rand sollte gleichverteilte Zufallszahlen zwischen 0 und 1 erzeugen. Hierzu werden wir nicht den Modulusoperator verändern, da das Ergebnis eine nicht mehr gesicherte gleichverteilte Zufallszahlenfolge ist, sondern vielmehr auf eine Divison und Rundung zurückgreifen.
Benötigen wir eine weitere Funktion die Zufallszahlen in einem bestimmten Bereich generiert.
Los geht's.´Änderung der Funktion :m-rand
;Erzeugung einer Liste mit gleichverteilten Zufallszahlen zwischen 0 und 1
;Argumente:
;#len = Länge der zu generierenden Liste
(defun :m-rand (#len / retlst)
(repeat #len
(setq retlst
(cons
(/ (setq *seed* (rem (+ (* 25173 (:t-seed)) 13849) 65536))
65536.0
)
retlst
)
)
)
)
Durch diese Änderung erhalten wir nun Zahlen zwischen 0 und 1...
Dieses Ergebniss erhalten wir, indem wir den errechneten Realwert durch den Modulusoperator (!) dividieren.
(:M-RAND 10) => (0.651596 0.499481 0.100357 0.864136 0.873581 0.942902 0.441971 0.823669 0.0418549 0.297455)
Und jetzt schreiben wir uns eine Funktion, die aus diesen Zahlen Ganzzahlen mit einer Ober- und Untergrenze erzeugt. Die Formel darfür lautet als AutoLISP-Ausdruck inklusive Auf- oder Abrundung:
(runde (+ (* (- Obergrenze Untergrenze) Zahl_zwischen_0_und_1)Untergrenze))
Dazu schreiben wir drei neue Funktionen.
Ober- und Untergrenzen Festlegung durch Übergabe einer Liste mit Zufallszahlen.
Um kleinere statistische Fehler beim kaufmännischen Runden von Zahlen in Ganzzahlen zu vermeiden (Aufrunden bei 0.5, jedoch nie ein Abrunden bei 0.5) schreiben wir uns eine Funktion, die ein unverzerrtes Runden ermöglicht. Dieses Verfahren rundet bei X+0.5 auf X ab, wenn X eine gerade Zahl ist. Ist X ungerade wird bis zur nächsten geraden Zahl aufgerundet! Beispiel: 2.5 => 2; -2.5=>-2; 5.5=>6; -5.5=-6; 4.51=>5; -4.51=>5
Zwecks einfacherem Handling wird für die Funktion :m-RoundToEven die Zahlen als Liste übergeben.
;Agrumente:
;#lrnd = Liste mit Gleichverteilten Zufallszahlen zwischen 0 und 1
;#LowerLimit = Untere Begrenzung
;#UpperLimit = Obere Begrenzung
(defun :m-RndLimit (#lrnd #LowerLimit #UpperLimit)
(mapcar '(lambda (rnd)
(+ (* (- #UpperLimit #LowerLimit) rnd) #LowerLimit)
)
#lrnd
)
)
;Kaufmännisches Runden von Realzahlen in Ganzzahlen
;Argumente: Realzahl
(defun :m-Round (#real)
(cond ((= 'real (type #real))
(if (<= 0.5 (abs (rem #real 1)))
(if (> #real 0)
(1+ (fix #real))
(1- (fix #real))
)
(fix #real)
)
)
(#real)
)
)
;Unverzerrtes Runden von Realzahlen in Ganzzahlen
;Argumente: Liste mit Realzahlen
(defun :m-RoundToEven (#numlst)
(mapcar '(lambda (n)
(if (/= (rem n 1) 0.5)
(:m-Round n)
(progn
(if (> n 0)
(* 2 (:m-Round (/ n 2)))
(* 2 (:m-Round (/ n 2)) -1)
)
)
)
)
#numlst
)
)
OK, dann wollen wir uns ein paar Zufallszahlen generieren lassen...let's go
(:m-RoundToEven(:m-RndLimit (:m-rand 10) 34 99)) => (75 52 77 68 94 85 40 65 59 66)
(:m-RoundToEven(:m-RndLimit (:m-rand 10) 1 15)) => (2 1 1 10 3 15 6 10 9 8)
Wollen wir Realwerte haben, dann lassen wir einfach die Funktion :m-RoundToEven weg:
(:m-RndLimit (:m-rand 10) 65 122) => (96.3285 122.425 74.2687 70.3384 121.325 82.7798 87.1721 117.722 68.8932 86.9182)
Yeeahh, nun haben wir alle nötigen Werkzeuge um gleichverteile Zufallszahlen inklusive einer Ober- und Untergrenze zu erzeugen. Was nun noch weiterhin eingebaut werden kann, ist eine Liste von Ganzzahlen ohne doppelte Zahlen.
Das überlasse ich aber nun Ihnen ;-)
Zum Schluß: Sammeln der hier vorgestellten Funktionen
Beispiele
Sammeln der Funktionen...
;Erzeugung eines Seed Wertes als initialen Wert x
;für die Generierung von Zufallszahlen
(defun :T-Seed ()
(cond (*seed*)
((setq *seed* (getvar "date")))
)
)
;Erzeugung einer Liste mit gleichverteilten Zufallszahlen zwischen 0 und 1
;Argumente:
;#len = Länge der zu generierenden Liste
;sexy coded by Rolf Wischnewski
(defun :m-rand (#len / retlst)
(repeat #len
(setq retlst
(cons
(/ (setq *seed* (rem (+ (* 25173 (:t-seed)) 13849) 65536))
65536.0
)
retlst
)
)
)
)
;Argumente:
;#lrnd = Liste mit Zahlen zwischen 0 und 1
;#LowerLimit = Untere Begrenzung
;#UpperLimit = Obere Begrenzung
(defun :m-RndLimit (#lrnd #LowerLimit #UpperLimit)
(mapcar '(lambda (rnd)
(+ (* (1+ (- #UpperLimit #LowerLimit)) rnd) #LowerLimit)
)
#lrnd
)
)
;Kaufmännisches Runden von Realzahlen in Ganzzahlen
;Argumente: Realzahl
;sexy coded by Rolf Wischnewski
(defun :m-Round (#real)
(cond ((= 'real (type #real))
(if (<= 0.5 (abs (rem #real 1)))
(if (> #real 0)
(1+ (fix #real))
(1- (fix #real))
)
(fix #real)
)
)
(#real)
)
)
;Unverzerrtes Runden von Realzahlen in Ganzzahlen
;Argumente: Liste mit Realzahlen
;sexy coded by Rolf Wischnewski
(defun :m-RoundToEven (#numlst)
(mapcar '(lambda (n)
(if (/= (rem n 1) 0.5)
(:m-Round n)
(progn
(if (> n 0)
(* 2 (:m-Round (/ n 2)))
(* 2 (:m-Round (/ n 2)) -1)
)
)
)
)
#numlst
)
)
Beispiele
Damit wir uns ein Bild machen können, was eigentlich gleichverteilte Zufallszahlen sind, generieren wir aus einer Zufalls-Punkteliste Punktobjekte. Dazu schreiben wir einfach ein paar zusätzlich neue Funktionen, die ich kurz beschreibe (natürlich werden die obigen Funktionen auch benötigt!)
:M-3d2d = Wandelt einen 3D-Punkt (getpoint) in einen 2D-Punkt um
:PO-addPoint = Erzeugt einen grafischen Punkt in Autocad mit den aktuellen Attributen
PointList = Erzeugt aus den übergebenen Argumenten eine Punkteliste
DrawPoints = Zeichnet Punke
;Argumente
;#PointSet = Punktmenge (Anzahl der zu zeichnenden Punkte)
;#Point = Punkt als Liste (getpoint)
;#Limit = Punktausbreitung
;sexy coded by Rolf Wischnewski
(defun DrawPoints (#PointSet #Point #Limit / :M-3d2d :PO-addPoint PointList)
;; local functions
(defun :M-3d2d (3dpt)
(list (float (car 3dpt)) (float (cadr 3dpt)))
)
(defun :PO-addPoint (#p)
(entmake (list (cons 0 "POINT") (cons 10 #p)))
)
(defun PointList (#PointSet #Point #Limit / retval)
(repeat #PointSet
(setq retval (cons (mapcar '+
(:M-3d2d #point)
(:m-rndlimit
(:m-rand 2)
0
(fix #Limit)
)
)
retval
)
)
)
)
;;body
(foreach p (PointList #PointSet #point #Limit)
(:PO-addPoint p)
)
(prin1)
)
Aufruf:
(DrawPoints
300
(setq *point* (getpoint " Bitte einen Punkt klicken"))
(getdist *point* " Ausbreitung festlegen")
)
Wie wir auf dem Bild erkennen, werden die Zufallszahlen offensichtlich gleichverteilt erzeugt. Somit haben wir einen erneuten Beleg für unsere gleichverteilte Zufallszahlen-Funktion.
Als Letztes Beispiel erzeugen wir mit Hilfe einer Funktion ein Art Paßwortliste. Diese Liste sollte nur aus Buchstaben und Ziffern bestehen.
Die Funktion besitzt eine lokale Funktion zur Erzeugung eines Paßwortes mit einer bestimmten Länge. Die Rückgabe ist eine Liste mit der gew&¨nschten Anzahl von Paßwörtern.
;Argumente
; #len-of-passwort = Ganzzahl für die Länge des Paßwortes
; #how-often = wieviele Paßwörter sollen erzeugt werden
;sexy coded by Rolf Wischnewski
(defun MakePasswordList (#len-of-passwort #how-often / ret-lst table password)
;; local Function
(defun password (#len-of-passwort / table)
(setq table
"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
)
(apply 'strcat
(mapcar '(lambda (a) (substr table a 1))
(:m-RoundToEven
(:m-RndLimit (:m-rand #len-of-passwort) 1 62)
)
)
)
)
;; body
(repeat #how-often
(setq ret-lst (cons (password #len-of-passwort) ret-lst))
)
ret-lst
)
Aufruf
(MAKEPASSWORDLIST 6 10) => ("MiYK4i" "D9JqBK" "Pky7Vu" "zrx9GB" "dhUBbr" "TS3xyX" "21i6z5" "YvSwtE" "cJeOM5" "HPKSOc")
(MAKEPASSWORDLIST 15 5) => ("UqZ2WdWD2yun3VV" "TZLuktwQ4mfyEvQ" "ABnI7jbqVQwb2Q6" "Zgsh1c9YKtifU0Y" "f4chH8kIXB1Cafl")
Mit diesem letzten Beispiel schließe ich dieses 1. Tutorial ab, hoffe bei Ihnen das Interesse für die Welt der Zufallszahlen geweckt zu haben und verabschiede mich mit einem Zitat aus der Chaostheorie von Edward N. Lorenz:
"Der Flügelschlag eines Schmetterlings im Amazonas-Urwald kann einen Orkan in Europa auslösen."
Weiter geht's hier =>Zufallszahlen erzeugen 2. Teil
Mit freundlicher Genehmigung von Rolf Wischnewski. Originalbeitrag im Juli 2006, CADmaro.de