> restart;
> with (numtheory,ifactors):
> with (combinat,fibonacci):
> N:=100;
> 
> for i from 1 to 1000 do a[i]:=0;b[i]:=0; od: 
> 
> # Phase 1 : Anlegen der Primzahlliste: 
> a[1]:=2; Nmax:=1;
> 

for i from 2 to N do           # Teste alle natuerlichen Zahlen
#if( isprime(i)=false) then    # Schliesse Primzahlen aus 
if (true) then
mm:=fibonacci(i);

s:=ifactors(mm);
anzahl:=nops(s[2]):          # Anzahl Primfaktoren   
   for j from 1 to anzahl do # Fuer gefundene Primfaktoren
   pzahl:=(s[2][j][1]):
   wertig:=(s[2][j][2]):
   lbl:=0; 
   for k from 1 to Nmax do   # Primzahlliste durchsuchen
    if (a[k]=pzahl) then
     lbl:=1;                 # Bereits eingetragen,Abbruch 
     b[k]:=b[k]+wertig;
     #printf(`   n=%g a=%g b=%g\n`,i,a[k],wertig);
    k:=Nmax;
    fi;           
    od;                      # Ende k Liste durchsuchen     
    if (lbl=0) then          # Primzahl nicht eingetragen
     Nmax:=Nmax+1;
     a[Nmax]:=pzahl;
     b[Nmax]:=b[Nmax]+wertig;
     #printf(`NEU n=%g a=%g b=%g\n`,i,a[Nmax],wertig);
   fi;    
   od;                       # end j..Anzahl                    
fi;                          # end ist nicht prim
 
#printf(`%g `,i);

od:                          # end natuerliche Zahlen
> 
> # ifactors liefert die Primzahlzerlegung der Zahl i  
  # dabei ist das Ergebnis ein Array dass die Primzahlen 
  # sowie deren Wertigkeit, Exponent enthaelt 
> #nops(s[2]) lefert die anzahl gefundener Primzahlen
> #s[2][1..anzahl][1]liefert die Primzahl
> #s[2][1..anzahl][2]liefert deren Wertigkeit
> #In Nmax stellt die Anzahl der gefundenen Primzahlen
> 
> 
> #Normieren der Primzahlmesswerte
> s:=0;
> for i from 1 to Nmax do s:=s+b[i]; od: 
> for i from 1 to Nmax do b[i]:=b[i]/s; od: 
> 
> #b sortieren nach Wert(Bubblesort) 
> for i from 1 to Nmax-1 do  
> for j from i to Nmax do 
> if (b[i]<b[j]) then 
> mem:=b[i];b[i]:=b[j];b[j]:=mem;
> mem:=a[i];a[i]:=a[j];a[j]:=mem;
> 
> fi;
> od;od;    
> 
> 
> # Erstellen der Zipf Verteilung 1/i 
> czipf:=1/evalf(sum(1/kk,kk=1..Nmax));
> for i from 1 to Nmax do zipf[i]:=czipf*1/i; od:
> 
> druck1:=seq([log(i),log(b[i])],i=1..Nmax):
> druck2:=seq([log(i),log(zipf[i])],i=1..Nmax):
> 
> plot([[druck1],[druck2]]);
> 
> druck1:=seq([(i),(b[i])],i=1..Nmax):
> druck2:=seq([(i),(zipf[i])],i=1..Nmax):
> 
> 
> plot([[druck1],[druck2]]);