User:Anatoly E. Voevudko/MyMixedPariScripts.gp

My PARI Scripts Published on Rosettacode Wiki

Sequences, etc.

Ludic numbers

• [1] Ludic numbers. Version #1.

 \\ Creating Vlf - Vector of ludic numbers' flags,
 \\ where the index of each flag=1 is the ludic number.
 \\ 2/28/16 aev
 ludic(maxn)={my(Vlf=vector(maxn,z,1),n2=maxn\2,k,j1);
 for(i=2,n2, 
     if(Vlf[i], k=0; j1=i+1;
        for(j=j1,maxn, if(Vlf[j], k++); if(k==i, Vlf[j]=0; k=0))
       ); 
    );
 return(Vlf);
 }
 {
 \\ Required tests:
 my(Vr,L=List(),k=0,maxn=25000);
 Vr=ludic(maxn);
 print("The first 25 Ludic numbers: ");
 for(i=1,maxn, if(Vr[i]==1, k++; print1(i," "); if(k==25, break)));
 print("");print("");
 k=0;
 for(i=1,999, if(Vr[i]==1, k++));
 print("Ludic numbers below 1000: ",k);
 print("");
 k=0;
 print("Ludic numbers 2000 to 2005: ");
 for(i=1,maxn, if(Vr[i]==1, k++; if(k>=2000&&k<=2005, listput(L,i)); if(k>2005, break)));
 for(i=1,6, print1(L[i]," "));
 print(""); print("");
 print("Ludic Triplets below 250: ");
 for(i=1,250, if(Vr[i]&&Vr[i+2]&&Vr[i+6], print1("(",i," ",i+2," ",i+6,") ")));
 }
  
 Output:
 The first 25 Ludic numbers:
 1 2 3 5 7 11 13 17 23 25 29 37 41 43 47 53 61 67 71 77 83 89 91 97 107
 
 Ludic numbers below 1000: 142
 
 Ludic numbers 2000 to 2005:
 21475 21481 21487 21493 21503 21511
 
 Ludic Triplets below 250:
 (1 3 7) (5 7 11) (11 13 17) (23 25 29) (41 43 47) (173 175 179) (221 223 227) (233 235 239)

• [2] Ludic numbers. Version #2.

 \\Upgraded script from A003309 to meet task requirements.
 \\ Creating Vl - Vector of ludic numbers.
 \\ 2/28/16 aev
 ludic2(maxn)={my(Vw=vector(maxn, x, x+1),Vl=Vec([1]),vwn=#Vw,i);
 while(vwn>0, i=Vw[1]; Vl=concat(Vl,[i]);
 Vw=vector((vwn*(i-1))\i,x,Vw[(x*i+i-2)\(i-1)]); vwn=#Vw
 ); return(Vl);
 }
 {
 \\ Required tests:
 my(Vr,L=List(),k=0,maxn=22000,vrs,vi);
 Vr=ludic2(maxn); vrs=#Vr;
 print("The first 25 Ludic numbers: ");
 for(i=1,25, print1(Vr[i]," "));
 print("");print("");
 k=0;
 for(i=1,vrs, if(Vr[i]<1000, k++, break));
 print("Ludic numbers below 1000: ",k);
 print("");
 k=0;
 print("Ludic numbers 2000 to 2005: ");
 for(i=2000,2005, print1(Vr[i]," "));
 print("");print("");
 print("Ludic Triplets below 250: ");
 for(i=1,vrs, vi=Vr[i]; if(i==1,print1("(",vi," ",vi+2," ",vi+6,") "); next); if(vi+6<250,if(Vr[i+1]==vi+2&&Vr[i+2]==vi+6, print1("(",vi," ",vi+2," ",vi+6,") "))));
 }
  
 Output:
 The first 25 Ludic numbers:
 1 2 3 5 7 11 13 17 23 25 29 37 41 43 47 53 61 67 71 77 83 89 91 97 107
 
 Ludic numbers below 1000: 142
 
 Ludic numbers 2000 to 2005:
 21475 21481 21487 21493 21503 21511
 
 Ludic Triplets below 250:
 (1 3 7) (5 7 11) (11 13 17) (23 25 29) (41 43 47) (173 175 179) (221 223 227) (233 235 239)

Hailstone sequence

• [3] Hailstone sequence. Version #2.

Works with PARI/GP 2.7.4 and above

Different kind of PARI scripts for Collatz sequences you can find in OEIS, e.g.: A070165

 \\ Get vector with Collatz sequence for the specified starting number.
 \\ Limit vector to the lim length, or less, if 1 (one) term is reached (when lim=0).
 \\ 3/26/2016 aev
 Collatz(n,lim=0)={
 my(c=n,e=0,L=List(n)); if(lim==0, e=1; lim=n*10^6); 
 for(i=1,lim, if(c%2==0, c=c/2, c=3*c+1); listput(L,c); if(e&&c==1, break));
 return(Vec(L)); } 
 Collatzmax(ns,nf)={
 my(V,vn,mxn=1,mx,im=1);
 print("Search range: ",ns,"..",nf);
 for(i=ns,nf, V=Collatz(i); vn=#V; if(vn>mxn, mxn=vn; im=i); kill(V)); 
 print("Hailstone/Collatz(",im,") has the longest length = ",mxn);
 } 

 {
 \\ Required tests:
 print("Required tests:");
 my(Vr,vrn);
 Vr=Collatz(27); vrn=#Vr;
 print("Hailstone/Collatz(27): ",Vr[1..4]," ... ",Vr[vrn-3..vrn],"; length = ",vrn);
 Collatzmax(1,100000);
 }

 Output:
 Required tests:
 Hailstone/Collatz(27): [27, 82, 41, 124] ... [8, 4, 2, 1]; length = 112
 Search range: 1..100000
 Hailstone/Collatz(77031) has the longest length = 351

 (15:32) gp > ##
   ***   last result computed in 15,735 ms.

Stern-Brocot sequence

• [4] Stern-Brocot sequence

Works with PARI/GP 2.7.4 and above

 \\ Stern-Brocot sequence
 \\ 5/27/16 aev
 SternBrocot(n)={
 my(L=List([1,1]),k=2);
 if(n<3,return(L));
 for(i=2,n, listput(L,L[i]+L[i-1]); if(k++>=n, break); listput(L,L[i]);if(k++>=n, break));
 return(Vec(L));
 \\ Find the first item in any list starting with sind index (return 0 or index).
 \\ 9/11/2015 aev
 findinlist(list, item, sind=1)={
 my(idx=0, ln=#list); if(ln==0 || sind<1 || sind>ln, return(0)); 
 for(i=sind, ln, if(list[i]==item, idx=i; break;)); return(idx);
 }
 }
 {
 \\ Required tests:
 my(v,j);
 v=SternBrocot(15); 
 print1("The first 15: "); print(v);
 v=SternBrocot(1200); 
 print1("The first i@n: "); \\print(v);
 for(i=1,10, if(j=findinlist(v,i), print1(i,"@",j,", ")));
 if(j=findinlist(v,100), print(100,"@",j));
 v=SternBrocot(10000); 
 print1("All GCDs=1?: ");
 j=1; for(i=2,10000, j*=gcd(v[i-1],v[i]));
 if(j==1, print("Yes"), print("No"));
 }

 Output:

 The first 15: [1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 5, 2, 5, 3, 4]
 The first i@n: 1@1, 2@3, 3@5, 4@9, 5@11, 6@33, 7@19, 8@21, 9@35, 10@39, 100@1179
 All GCDs=1?: Yes

Fractran

• [5] Fractran

In this version ideas were borrowed from C++, Java and JavaScript.

Works with PARI/GP 2.7.4 and above

 \\ FRACTRAN 
 \\ 4/27/16 aev
 fractran(val,ft,lim)={
 my(ftn=#ft,fti,di,L=List(),j=0);
 while(val&&j<lim, listput(L,val);
   for(i=1,ftn, fti=ft[i]; di=denominator(fti);
       if(val%di==0, break));\\fend i
   val= numerator(fti)*val/di; j++);\\wend j
 return(Vec(L));
 }

 {\\ Executing:
 my(v=[17/91,78/85,19/51,23/38,29/33,77/29,95/23,77/19,1/17,11/13,13/11,15/14,15/2,55/1]);
 print(fractran(2,v,15));
 }

 Output:

 [2, 15, 825, 725, 1925, 2275, 425, 390, 330, 290, 770, 910, 170, 156, 132]

Convert seconds to compound duration

• [6] Convert seconds to compound duration

Works with PARI/GP 2.7.4 and above

 \\ Convert seconds to compound duration
 \\ 4/11/16 aev
 secs2compdur(secs)={
 my(us=[604800,86400,3600,60,1],ut=[" wk, "," d, "," hr, "," min, "," sec"],
    cd=[0,0,0,0,0],u,cdt="");
 for(i=1,5, u=secs\us[i]; if(u==0, next, cd[i]=u; secs-=us[i]*cd[i]));
 for(i=1,5, if(cd[i]==0, next, cdt=Str(cdt,cd[i],ut[i])));
 if(ssubstr(cdt,#cdt-1,1)==",", cdt=ssubstr(cdt,1,#cdt-2));
 return(cdt);
 }
 {\\ Required tests:
 print(secs2compdur(7259));
 print(secs2compdur(86400));
 print(secs2compdur(6000000));
 }

 Output:

 2 hr, 59 sec
 1 d
 9 wk, 6 d, 10 hr, 40 min

String functions

Repeat a string

• [7] Repeat a string. Version #2.

 \\ Repeat a string str the specified number of times ntimes and return composed string.
 \\ 3/3/2016 aev
 srepeat(str,ntimes)={
 my(srez=str,nt=ntimes-1);
 if(ntimes<1||#str==0,return("")); 
 if(ntimes==1,return(str)); 
 for(i=1,nt, srez=concat(srez,str));
 return(srez);
 }
  
 {
 \\ TESTS
 print(" *** Testing srepeat:");
 print("1.",srepeat("a",5));
 print("2.",srepeat("ab",5));
 print("3.",srepeat("c",1));
 print("4.|",srepeat("d",0),"|");
 print("5.|",srepeat("",5),"|");
 print1("6."); for(i=1,10000000, srepeat("e",10));
 }
  
 Output:
  *** Testing srepeat:
 1.aaaaa
 2.ababababab
 3.c
 4.||
 5.||
 6.
 (16:00) gp > ##
   ***   last result computed in 1min, 2,939 ms.

Reverse a string

• [8] Reverse a string. Version #2.

 \\ Return reversed string str.
 \\ 3/3/2016 aev
 sreverse(str)={return(Strchr(Vecrev(Vecsmall(str))))}
  
 {
 \\ TEST1 
 print(" *** Testing sreverse from Version #2:");
 print(sreverse("ABCDEF"));
 my(s,sr,n=10000000);
 s="ABCDEFGHIJKL";
 for(i=1,n, sr=sreverse(s));
 }

 Output:
  *** Testing sreverse from Version #2:
 FEDCBA
 (17:28) gp > ##
   ***   last result computed in 8,642 ms.

 \\ Version #1 upgraded to complete function. Practically the same.
 reverse(str)={return(concat(Vecrev(str)))}
  
 {
 \\ TEST2 
 print(" *** Testing reverse from Version #1:");
 print(reverse("ABCDEF"));
 my(s,sr,n=10000000);
 s="ABCDEFGHIJKL";
 for(i=1,n, sr=reverse(s));
 }
  
 Output:
  *** Testing reverse from Version #1:
 FEDCBA
 (17:31) gp > ##
   ***   last result computed in 11,814 ms.
 NOTE: Version #2 is faster then Version #1.

Substring of a string

• [9] Substring

Works with PARI/GP 2.7.4 and above

 \\ Returns the substring of string str specified by the start position s and length n.
 \\ If n=0 then to the end of str.
 \\ ssubstr() 3/5/16 aev
 ssubstr(str,s=1,n=0)={
 my(vt=Vecsmall(str),ve,vr,vtn=#str,n1);
 if(vtn==0,return(""));
 if(s<1||s>vtn,return(str));
 n1=vtn-s+1; if(n==0,n=n1); if(n>n1,n=n1);
 ve=vector(n,z,z-1+s); vr=vecextract(vt,ve); return(Strchr(vr));
 }
 {\\ TEST
 my(s="ABCDEFG",ns=#s);
 print(" *** Testing ssubstr():");
 print("1.",ssubstr(s,2,3));
 print("2.",ssubstr(s));
 print("3.",ssubstr(s,,ns-1));
 print("4.",ssubstr(s,2));
 print("5.",ssubstr(s,,4));
 print("6.",ssubstr(s,0,4));
 print("7.",ssubstr(s,3,7));
 print("8.|",ssubstr("",1,4),"|");
 }

 Output:
  *** Testing ssubstr():
 1.BCD
 2.ABCDEFG
 3.ABCDEF
 4.BCDEFG
 5.ABCD
 6.ABCDEFG
 7.CDEFG
 8.||

Tokenize a string

• [10] Tokenize a string.

Version #1. Simple version, like the most custom ones here (for this task). This version has 1 character delimiter, which is not allowed in the beginning and at the end of string, in addition, double, triple, etc., delimiters are not allowed too.

Works with PARI/GP 2.7.4 and above

 \\ Tokenize a string str according to 1 character delimiter d. Return a list of tokens.
 \\ Using ssubstr() from http://rosettacode.org/wiki/Substring#PARI.2FGP
 \\ tokenize() 3/5/16 aev
 tokenize(str,d)={
 my(str=Str(str,d),vt=Vecsmall(str),d1=sasc(d),Lr=List(),sn=#str,v1,p1=1);
 for(i=p1,sn, v1=vt[i]; if(v1==d1, listput(Lr,ssubstr(str,p1,i-p1)); p1=i+1)); 
 return(Lr);
 }
 {
 \\ TEST 
 print(" *** Testing tokenize from Version #1:");
 print("1.", tokenize("Hello,How,Are,You,Today",","));
 \\ BOTH 2 & 3 are NOT OK!!
 print("2.",tokenize("Hello,How,Are,You,Today,",","));
 print("3.",tokenize(",Hello,,How,Are,You,Today",","));
 }

 Output:
  *** Testing tokenize from Version #1:
 1.List(["Hello", "How", "Are", "You", "Today"])
 2.List(["Hello", "How", "Are", "You", "Today", ","])
 3.List([",Hello,,How,Are,You,Today,", "Hello", ",How,Are,You,Today,", "How", "Ar
 e", "You", "Today"])

Version #2. Advanced version. Delimiter is allowed in any place. In addition, multiple delimiters are allowed too. This is really useful for considering omitted data. This version can be used for positional parameters processing, or for processing data from tables with string rows.

Works with PARI/GP 2.7.4 and above

 \\ Tokenize a string str according to 1 character delimiter d. Return a list of tokens.
 \\ Using ssubstr() from http://rosettacode.org/wiki/Substring#PARI.2FGP
 \\ stok() 3/5/16 aev
 stok(str,d)={
 my(d1c=ssubstr(d,1,1),str=Str(str,d1c),vt=Vecsmall(str),d1=sasc(d1c),
    Lr=List(),sn=#str,v1,p1=1,vo=32);
 if(sn==1, return(List(""))); if(vt[sn-1]==d1,sn--);
 for(i=1,sn, v1=vt[i];
     if(v1!=d1, vo=v1; next);
     if(vo==d1||i==1, listput(Lr,""); p1=i+1; vo=v1; next);
     if(i-p1>0, listput(Lr,ssubstr(str,p1,i-p1)); p1=i+1);
     vo=v1;
    ); 
 return(Lr);
 }
 {
 \\ TEST 
 print(" *** Testing stok from Version #2:");
 \\ pp - positional parameter(s)
 print("1. 5 pp: ", stok("Hello,How,Are,You,Today",","));
 print("2. 5 pp: ", stok("Hello,How,Are,You,Today,",","));
 print("3. 9 pp: ", stok(",,Hello,,,How,Are,You,Today",","));
 print("4. 6 pp: ", stok(",,,,,,",","));
 print("5. 1 pp: ", stok(",",","));
 print("6. 1 pp: ", stok("Hello-o-o??",","));
 print("7. 0 pp: ", stok("",","));
 }

 Output:
  *** Testing stok from Version #2:
 1. 5 pp: List(["Hello", "How", "Are", "You", "Today"])
 2. 5 pp: List(["Hello", "How", "Are", "You", "Today"])
 3. 9 pp: List(["", "", "Hello", "", "", "How", "Are", "You", "Today"])
 4. 6 pp: List(["", "", "", "", "", ""])
 5. 1 pp: List([""])
 6. 1 pp: List(["Hello-o-o??"])
 7. 0 pp: List([""])

Jaro distance

• [11] Jaro distance

Works with PARI/GP 2.7.4 and above

 \\ Jaro distance between 2 strings s1 and s2.
 \\ This version was translated from Java and Perl on RC.
 \\ 5/12/16 aev
 jaroDist(s1,s2)={
 my(vt1=Vecsmall(s1),vt2=Vecsmall(s2),n1=#s1,n2=#s2,d,
    md=max(n1,n2)\2-1,cs,ce,mc=0,tr=0,k=1,ds,
    s1m=vector(n1,z,0),s2m=vector(n2,z,0));
 if(!n1||!n2, return(0));
 for(i=1,n1,
   cs=max(1,i-md);
   ce=min(i+md+1,n2);
   for(j=cs,ce,
     if(s2m[j],next);
     if(vt1[i]!=vt2[j], next);
     mc++; s1m[i]=1; s2m[j]=1; break;
   );\\fend j
 );\\fend i
 if(!mc, return(0));
 for(i=1,n1,
   if(!s1m[i], next);
   while(!s2m[k], k++);
   if(vt1[i]!=vt2[k], tr++);
   k++
 );\\fend i
 d=(mc/n1+mc/n2+(mc-tr/2)/mc)/3.0;
 ds=Strprintf("%.5f",d);
 print(" *** Jaro distance is: ",ds," for strings: ",s1,", ",s2);
 return(d);
 }
 { \\ Testing:
 jaroDist("MARTHA","MARHTA"); 
 jaroDist("DIXON","DICKSONX");
 jaroDist("JELLYFISH","SMELLYFISH");
 jaroDist("DWAYNE","DUANE");
 }

 Output:

  *** Jaro distance is: 0.94444 for strings: MARTHA, MARHTA
  *** Jaro distance is: 0.76667 for strings: DIXON, DICKSONX
  *** Jaro distance is: 0.89630 for strings: JELLYFISH, SMELLYFISH
  *** Jaro distance is: 0.82222 for strings: DWAYNE, DUANE

Levenshtein distance

• [12] Levenshtein distance between two words

Translated from JavaScript.
Works with PARI/GP 2.7.4 and above

 \\ Levenshtein distance between two words
 \\ 6/21/16 aev
 levensDist(s1,s2)={
 my(n1=#s1,n2=#s2,v1=Vecsmall(s1),v2=Vecsmall(s2),c,
    n11=n1+1,n21=n2+1,t=vector(n21,z,z-1),u0=vector(n21),u=u0);
 if(s1==s2, return(0)); if(!n1, return(n2)); if(!n2, return(n1));
 for(i=2,n11, u=u0; u[1]=i-1;
   for(j=2,n21,
     if(v1[i-1]==v2[j-1], c=t[j-1], c=vecmin([t[j-1],t[j],u[j-1]])+1);
     u[j]=c;
   );\\fend j
   t=u;
 );\\fend i 
 print(" *** Levenshtein distance = ",t[n21]," for strings: ",s1,", ",s2);
 return(t[n21]);
 }
 { \\ Testing:
 levensDist("kitten","sitting"); 
 levensDist("rosettacode","raisethysword");
 levensDist("Saturday","Sunday");
 levensDist("oX","X");
 levensDist("X","oX");
 }

 Output:
  *** Levenshtein distance = 3 for strings: kitten, sitting
  *** Levenshtein distance = 8 for strings: rosettacode, raisethysword
  *** Levenshtein distance = 3 for strings: Saturday, Sunday
  *** Levenshtein distance = 1 for strings: oX, X
  *** Levenshtein distance = 1 for strings: X, oX

Trees and other plotting

Draw a cuboid

• [13] Draw a cuboid

Plotting lines and scaling in PARI/GP is not designed for "cuboids". But you are welcome to play with parameters of this Cuboid() function.

Output Cuboid1.png

Output Cuboid2.png

Works with PARI/GP 2.7.4 and above

 \\ Simple "cuboid". Try different parameters of this Cuboid() function.
 \\ 4/11/16 aev
 Cuboid(a,b,c,u=10)={
 my(dx,dy,ttl="Cuboid AxBxC: ",size=200,da=a*u,db=b*u,dc=c*u);
 print(" *** ",ttl,a,"x",b,"x",c,"; u=",u);
 plotinit(0);
 plotscale(0, 0,size, 0,size); 
 plotcolor(0,7); \\grey
 plotmove(0, 0,0);
 plotrline(0,dc,da\2); plotrline(0,db,0); plotrline(0,-db,0);
 plotrline(0,0,da); 
 plotcolor(0,2); \\black
 plotmove(0, db,da);
 plotrline(0,0,-da); plotrline(0,-db,0);
 plotrline(0,0,da); plotrline(0,db,0);
 plotrline(0,dc,da\2); plotrline(0,-db,0); plotrline(0,-dc,-da\2);
 plotmove(0, db,0);
 plotrline(0,dc,da\2); plotrline(0,0,da);
 plotdraw([0,size,size]);
 }
 
 {\\ Executing:
 Cuboid(2,3,4,20); \\Cuboid1.png
 Cuboid(5,3,1,20); \\Cuboid2.png
 }
 
 Output:
 
 > Cuboid(2,3,4,20); \\Cuboid1.png
 *** Cuboid AxBxC: 2x3x4; u=20
 > Cuboid(5,3,1,20); \\Cuboid2.png
 *** Cuboid AxBxC: 5x3x1; u=20

Barnsley fern fractal

• [14] Barnsley fern fractal

Translated from zkl.
Works with PARI/GP 2.7.4 and above

Output BarnsleyFern.png

 \\ Barnsley fern fractal
 \\ 6/17/16 aev
 pBarnsleyFern(size,lim)={
 my(X=List(),Y=X,x=y=xw=yw=0.0,r,c=0.625);
 print(" *** Barnsley Fern, size=",size," lim=",lim);
 plotinit(0); plotcolor(0,6); \\green
 plotscale(0, -3,3, 0,10); plotmove(0, 0,0);
 for(i=1, lim,
   r=random(100);
   if(r<=1, xw=0;yw=0.16*y,
     if(r<=8, xw=0.2*x-0.26*y;yw=0.23*x+0.22*y+1.6,
       if(r<=15, xw=-0.15*x+0.28*y;yw=0.26*x+0.24*y+0.44,
         xw=0.85*x+0.04*y;yw=-0.04*x+0.85*y+1.6)));
   x=xw;y=yw; listput(X,x*c); listput(Y,y);
 );\\fend i
 plotpoints(0,Vec(X),Vec(Y));
 plotdraw([0,-3,-0]);
 }
 {\\ Executing:
 pBarnsleyFern(530,100000);  \\ BarnsleyFern.png
 }

 Output:
 > pBarnsleyFern(530,100000);  \\ BarnsleyFern.png
  *** Barnsley Fern, size=530 lim=100000

Barnsley fern fractal in JavaScript

• [15] Barnsley fern fractal in JavaScript

Translated from PARI/GP.

Output BarnsleyFernjs.png

 // Barnsley fern fractal
 //6/17/16 aev
 function pBarnsleyFern(canvasId,lim) {
   // DCLs
   var canvas = document.getElementById(canvasId);
   var ctx = canvas.getContext("2d");
   var w = canvas.width;
   var h = canvas.height;
   var x=0.,y=0.,xw=0.,yw=0.,r;
   // Like in PARI/GP: return random number 0..max-1
   function randgp(max) {return Math.floor(Math.random()*max)}
   // Cleaning canvas and setting plotting color 
   ctx.fillStyle="white"; ctx.fillRect(0,0,w,h);
   ctx.fillStyle="green";
   // MAIN LOOP
   for(var i=0; i<lim; i++) {
     r=randgp(100);
     if (r<=1) {xw=0;yw=0.16*y;}
     else if (r<=8) {xw=0.2*x-0.26*y;yw=0.23*x+0.22*y+1.6;}
     else if (r<=15) {xw=-0.15*x+0.28*y;yw=0.26*x+0.24*y+0.44;}
     else {xw=0.85*x+0.04*y;yw=-0.04*x+0.85*y+1.6;}
     x=xw;y=yw; ctx.fillRect(x*50+260,-y*50+540,1,1);
   }//fend i
 }

Executing:

 <html>
  <head><script src="BarnsleyFern.js"></script></head>
  <body onload="pBarnsleyFern('canvas', 100000)">
    <xbr /><xh3>Barnsley fern fractal</xh3>
    <canvas id="canvas" width="540" height="540" style="border: 2px inset;"></canvas>
  </body>
 </html>

Output:

 Page with BarnsleyFernjs.png

Barnsley fern fractal in R

• [16] Barnsley fern fractal in R

Translated from PARI/GP.

Output BarnsleyFernR.png

 ## pBarnsleyFern(fn, n, clr, ttl, psz=600): Plot Barnsley fern fractal.
 ## Where: fn - file name; n - number of dots; clr - color; ttl - plot title;
 ## psz - picture size.
 ## 7/27/16 aev
 pBarnsleyFern <- function(fn, n, clr, ttl, psz=600) {
   cat(" *** START:", date(), "n=", n, "clr=", clr, "psz=", psz, "\n");
   cat(" *** File name -", fn, "\n");
   pf = paste0(fn,".png"); # pf - plot file name
   A1 <- matrix(c(0,0,0,0.16,0.85,-0.04,0.04,0.85,0.2,0.23,-0.26,0.22,-0.15,0.26,0.28,0.24), ncol=4, nrow=4,   byrow=TRUE);
   A2 <- matrix(c(0,0,0,1.6,0,1.6,0,0.44), ncol=2, nrow=4, byrow=TRUE);
   P <- c(.01,.85,.07,.07);
   # Creating matrices M1 and M2.
   M1=vector("list", 4); M2 = vector("list", 4);
   for (i in 1:4) {
     M1i <- matrix(c(A1[i,1:4]), nrow=2);
     M2i <- matrix(c(A2[i, 1:2]), nrow=2);
   }
   x <- numeric(n); y <- numeric(n);
   x[1] <- y[1] <- 0;
   for (i in 1:(n-1)) {
     k <- sample(1:4, prob=P, size=1);
     M <- as.matrix(M1k);
     z <- M%*%c(x[i],y[i]) + M2k; 
     x[i+1] <- z[1]; y[i+1] <- z[2];
   }
   plot(x, y, main=ttl, axes=FALSE, xlab="", ylab="", col=clr, cex=0.1);
   # Writing png-file
   dev.copy(png, filename=pf,width=psz,height=psz);
   # Cleaning 
   dev.off(); graphics.off();
   cat(" *** END:",date(),"\n");
 }
 ## Executing:
 pBarnsleyFern("BarnsleyFernR", 100000, "dark green", "Barnsley Fern Fractal", psz=600)

 Output:

  > pBarnsleyFern("BarnsleyFernR", 100000, "dark green", "Barnsley Fern Fractal", psz=600)
  *** START: Wed Jul 27 13:50:49 2016 n= 1e+05 clr= dark green psz= 600 
  *** File name - BarnsleyFernR 
  *** END: Wed Jul 27 13:50:56 2016 
  + BarnsleyFernR.png file

Sierpinski carpet fractal

• [17] Sierpinski carpet fractal

Works with PARI/GP 2.7.4 and above

Translation of: Python

Output SierpC5.png

 \\ Sierpinski carpet fractal
 \\ Note: plotmat() can be found here on my VoeLib.gp page or on 
 \\ http://rosettacode.org/wiki/Brownian_tree#PARI.2FGP page.
 \\ 6/10/16 aev
 inSC(x,y)={
 while(1, if(!x||!y,return(1)); 
          if(x%3==1&&y%3==1, return(0));
          x\=3; y\=3;);\\wend
 }
 pSierpinskiC(n,pflg=0)={
 my(n3=3^n-1,M);
 if(pflg<0||pflg>1, pflg=0); if(pflg, M=matrix(n3+1,n3+1));
 for(i=0,n3, for(j=0,n3,
     if(inSC(i,j), 
     if(pflg, M[i+1,j+1]=1, print1("* ")), if(!pflg, print1("  ")));
    ); if(!pflg, print(""));
    );\\fend i
 if(pflg, plotmat(M));
 }
 {\\ Test:
 pSierpinskiC(3);
 pSierpinskiC(5,1); \\ SierpC5.png
 }

 Output:

 > pSierpinskiC(3);
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *   * *   * *   * *   * *   * *   * *   * *   * *   *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * * *       * * * * * *       * * * * * *       * * *
 *   *       *   * *   *       *   * *   *       *   *
 * * *       * * * * * *       * * * * * *       * * *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *   * *   * *   * *   * *   * *   * *   * *   * *   *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * * * * * * * * *                   * * * * * * * * *
 *   * *   * *   *                   *   * *   * *   *
 * * * * * * * * *                   * * * * * * * * *
 * * *       * * *                   * * *       * * *
 *   *       *   *                   *   *       *   *
 * * *       * * *                   * * *       * * *
 * * * * * * * * *                   * * * * * * * * *
 *   * *   * *   *                   *   * *   * *   *
 * * * * * * * * *                   * * * * * * * * *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *   * *   * *   * *   * *   * *   * *   * *   * *   *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * * *       * * * * * *       * * * * * *       * * *
 *   *       *   * *   *       *   * *   *       *   *
 * * *       * * * * * *       * * * * * *       * * *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 *   * *   * *   * *   * *   * *   * *   * *   * *   *
 * * * * * * * * * * * * * * * * * * * * * * * * * * *
 > pSierpinskiC(5,1); \\ SierpC5.png
   *** matrix(243x243) 32768 DOTS

Sierpinski triangle fractal

• [18] Sierpinski triangle fractal

Output SierpT9.png

Works with PARI/GP 2.7.4 and above

 \\ Sierpinski triangle fractal
 \\ Note: plotmat() can be found here on my VoeLib.gp page or on 
 \\ http://rosettacode.org/wiki/Brownian_tree#PARI.2FGP page.
 \\ 6/3/16 aev
 pSierpinskiT(n)={
 my(sz=2^n,M=matrix(sz,sz),x,y);
 for(y=1,sz, for(x=1,sz, if(!bitand(x,y),M[x,y]=1);));\\fends
 plotmat(M);
 }
 {\\ Test:
 pSierpinskiT(9); \\ SierpT9.png
 }

   Output:
 > pSierpinskiT(9); \\ SierpT9.png
  *** matrix(512x512) 19682 DOTS

Dragon curve

• [19] Dragon curve (with pictures on Rosettacode Wiki)

Version #3. This is actualy Version #1 upgraded to the reusable function.

Works with PARI/GP 2.7.4 and above

 \\ Dragon curve
 \\ 4/8/16 aev
 Dragon(level)={my(p=[0,1],end);
 print(" *** Dragon curve, level ",level);
 for(i=1,level, end=(1+I)*p[#p];
     p=concat(p,apply(z->(end-I*z),Vecrev(p[^-1]))) );
 plothraw(apply(real,p),apply(imag,p), 1);
 }

 {\\ Executing/Testing:
 Dragon(13); \\ Dragon13.png
 Dragon(17); \\ Dragon17.png
 Dragon(21); \\ Dragon21.png
 Dragon(23); \\ No result
 }

 Output: 
 *** Dragon curve, level 13
   ***   last result computed in 282 ms.
   *** Dragon curve, level 17
 ***   last result computed in 453 ms.
   *** Dragon curve, level 21
   ***   last result computed in 7,266 ms.
 *** Dragon curve, level 23
   *** concat: the PARI stack overflows !
   ***   last result computed in 0 ms.

The Archimedean spiral

• [20] The Archimedean spiral

Works with PARI/GP 2.7.4 and above

Output ArchiSpiral1.png

Output ArchiSpiral2.png

 \\ The Archimedean spiral  
 \\ ArchiSpiral() - Where: lps is a number of loops, c is a direction 0/1
 \\ (counter-clockwise/clockwise). 6/6/16 aev
 \\ Note: cartes2() can be found here on 
 \\ http://rosettacode.org/wiki/Polyspiral#PARI.2FGP page.
 ArchiSpiral(size,lps,c=0)={
 my(a=.0,ai=.1,r=.0,ri=.1,as=lps*2*Pi,n=as/ai,x,y,vc,vx=List(.0),vy=vx);
 if(c<0||c>1, c=0); if(c, ai*=-1);
 print(" *** The Archimedean spiral: size=",size," loops=",lps," c=",c);
 for(i=1, n, vc=cartes2(r,a); x=vc[1]; y=vc[2];
     listput(vx,x); listput(vy,y);
     r+=ri; a+=ai;
 );\\fend i
 plothraw(Vec(vx),Vec(vy));
 }
 {\\ Executing:
 ArchiSpiral(640,5);   \\ArchiSpiral1.png
 ArchiSpiral(640,5,1); \\ArchiSpiral2.png
 }

 Output:

 > ArchiSpiral(640,5);  \\ArchiSpiral1.png
  *** The Archimedean spiral: size=640 loops=5 c=0
 > ArchiSpiral(640,5,1);  \\ArchiSpiral2.png
  *** The Archimedean spiral: size=640 loops=5 c=1

Polyspiral

• [21] Polyspiral: a spiral made of multiple line segments

Works with PARI/GP 2.7.4 and above

Plotting helper functions

Both versions #1 and #2 are based on using my own small plotting helper functions. You can find a few others on OEIS Wiki.

 \\ Plot the line from x1,y1 to x2,y2.
 plotline(x1,y1,x2,y2,w=0)={plotmove(w, x1,y1);plotrline(w,x2-x1,y2-y1);}
 \\ Convert degrees to radians.
 rad2(degs)={return(degs*Pi/180.0)}
 \\ Convert Polar coordinates to Cartesian.
 cartes2(r,a,rndf=0)={my(v,x,y); x=r*cos(a); y=r*sin(a);
 if(rndf==0, return([x,y]), return(round([x,y])))}

Version #1. Polyspiral (a spiral made of multiple line segments).

In this version function plotpspiral() was translated from Java and J, plus my own plotline() function was used. Some tweaks and options were added to make it reusable and outputting differently looking polyspirals. There are no animation features in PARI/GP.

Output Polyspiral1.png

Output Polyspiral2.png

Output Polyspiral3.png

Output Polyspiral3b.png

Output Polyspiral4.png

 \\ Polyspiral() - Where: ai is an angle increment (in radians), d is a distance/length,
 \\ c is a direction 0/1 (clockwise/counter-clockwise); other parameters are self explanative. 
 \\ 4/15/16 aev Last updated: 4/18/16
 plotpspiral(size,lim,ai,d,di,c=0)={
 my(x1,y1,x2,y2,air=ai*Pi,a,sai=Strprintf("%.3f",ai));
 if(c<0||c>1, c=0);
 print(" *** Polyspiral, size=",size," lim=",lim," ai=",sai," d=",d," di=",di," c=",c);
 x1=0; y1=0;\\from J OK!! the best
 a=air;
 for(i=0, lim,
     if(c==0, x2=x1+cos(a)*d; y2=y1-sin(a)*d,
              x2=x1-sin(a)*d; y2=y1+cos(a)*d;);
     plotline(x1,y1,x2,y2);
     x1=x2; y1=y2; d+=di; a+=air;
    );\\fend i
 }
 
 polyspiral(size,lim,ai,d,di,c=0)={
 plotinit(0);
 plotcolor(0,3); \\blue
 plotscale(0, -size,size, -size,size); 
 plotmove(0, 0,0);
 plotpspiral(size,lim,ai,d,di,c);
 plotdraw([0,size,size]);
 }
 
 {\\ Executing:
 polyspiral(1500,1500,0.25,9,5);  \\Polyspiral1.png
 polyspiral(1500,1500,0.25,3,2);  \\Polyspiral2.png
 polyspiral(10000,10000,0.03,3,2);  \\Polyspiral3.png
 polyspiral(10000,10000,0.03,3,2,1);\\Polyspiral3b.png
 polyspiral(100000,100000,0.03,3,2);\\Polyspiral4.png
 }
  
 Output:
 
 > polyspiral(1500,1500,0.25,9,5);  \\Polyspiral1.png
 *** Polyspiral, size=1500 lim=1500 ai=0.250 d=9 di=5 c=0
  
 > polyspiral(1500,1500,0.25,3,2);  \\Polyspiral2.png
 *** Polyspiral, size=1500 lim=1500 ai=0.250 d=3 di=2 c=0
  
 > polyspiral(10000,10000,0.03,3,2);  \\Polyspiral3.png
 *** Polyspiral, size=100000 lim=100000 ai=0.030 d=3 di=2 c=0
 
 > polyspiral(10000,10000,0.03,3,2,1);  \\Polyspiral3b.png
 *** Polyspiral, size=100000 lim=100000 ai=0.030 d=3 di=2 c=1

 > polyspiral(100000,100000,0.03,3,2);  \\Polyspiral4.png
 *** Polyspiral, size=100000 lim=100000 ai=0.030 d=3 di=2 c=0

Version #2. Multi-spiral figure translated from zkl.

This is definitely not a polyspiral, but a very nice "multi-spiral" figure similar to shown in zkl and in a few other languages. Also, there is a very nice and impressive animation created in zkl, but not possible in PARI/GP.

Output Spiralz.png

 \\ plotpspiralz() Multi-spiral figure translated from zkl using my own ploting functions.
 \\ 4/15/16 aev
 plotpspiralz(size,lim,ai,di,lim2)={
 my(x1,y1,u1,v1,air=rad2(ai),a,sai=Strprintf("%.3f",ai),sdi=Strprintf("%.3f",di),
    sz2=size\2,aj,inc,ao,x,y,u,v,vc,r2i=rad2(130.0),d=0.0);
 print(" *** Spiralz: size=",size," lim=",lim," ai=",sai," di=",sdi," lim2=",lim2);
 x1=0; y1=0; u1=0; v1=0;
 for(i=1, lim,
   r=0.0; a=0.0;ao=0.0;
   if(i>1, inc=air+r2i, inc=air);
   for(j=1, lim2,
     d=r+di; aj=a+inc;
     vc=cartes2(r,a); x=vc[1]; y=vc[2];
     vc=cartes2(r,aj); u=vc[1]; v=vc[2];
     plotline(ao+x,ao+y,ao+u,ao+v);
     r=d; a=aj;
   );\\fend j
   air+=0.05;
 );\\fend i
 }

 \\ Spiralz() - Where: ai is an angle increment (in radians), di is a distance/length
 \\ increment, other parameters are self explanative. 
 \\ 4/15/16 aev
 Spiralz(size,lim,ai,di,lim2)={
 plotinit(0); plotcolor(0,3); \\blue
 plotscale(0, -size,size, -size,size); 
 \\plotscale(0, 0,size, 0,size); 
 plotmove(0, 0,0);
 plotpspiralz(size,lim,ai,di,lim2);
 plotdraw([0,size,size]);
 }

 {\\ Executing:
 Spiralz(640,2,3.0,3.0,128);  \\Spiralz1.png
 }

 Output:

 > Spiralz(640,2,3.0,3.0,128);  \\Spiralz1.png
  *** Spiralz: size=640 lim=2 ai=3.000 di=3.000 lim2=128

Ulam spiral

• [22] Ulam spiral: from small presentation to very big

In this version function plotulamspir() was translated from VB, plus upgraded to plot/print different kind of Ulam spirals. My own plotting helper functions and string functions were used and made it possible. See Voelib.gp from my OEIS user page.

Output ULAMspiral1.png

Output ULAMspiral2.png

Works with PARI/GP 2.7.4 and above

 \\ Ulam spiral (plotting/printing)
 \\ plotulamspir(n,pflg=0): Creating and plotting/printing matrix nxn according to the p-flag
 \\ (0-real plot,1-"*",2-primes)
 \\ 4/19/16 aev
 plotulamspir(n,pflg=0)={
 my(n=if(n%2==0,n++,n),M=matrix(n,n),x,y,xmx,ymx,cnt,dir,n2=n*n,pch,sz=#Str(n2),pch2=srepeat(" ",sz));
 if(pflg<0||pflg>2,pflg=0);
 print(" *** Ulam spiral: ",n,"x",n," matrix, p-flag=",pflg);
 x=y=n\2+1; xmx=ymx=cnt=1; dir="R";
 for(i=1,n2,
     if(isprime(i), if(!insm(M,x,y), break); if(pflg==2, M[y,x]=i, M[y,x]=1));
     if(dir=="R", if(xmx>0, x++;xmx--, dir="U";ymx=cnt;y--;ymx--); next); 
     if(dir=="U", if(ymx>0, y--;ymx--, dir="L";cnt++;xmx=cnt;x--;xmx--); next); 
     if(dir=="L", if(xmx>0, x--;xmx--, dir="D";ymx=cnt;y++;ymx--); next); 
     if(dir=="D", if(ymx>0, y++;ymx--, dir="R";cnt++;xmx=cnt;x++;xmx--); next); 
    );\\fend
 \\Plot/Print according to the p-flag (0-real plot,1-"*",2-primes)
 if(pflg==0, plotmat(M));
 if(pflg==1, for(i=1,n, 
             for(j=1,n, if(M[i,j]==1, pch="*", pch=" ");
                        print1(" ",pch)); print(" ")));
 if(pflg==2, for(i=1,n, 
             for(j=1,n, if(M[i,j]==0, pch=pch2, pch=spad(Str(M[i,j]),sz,,1));
                        print1(" ",pch)); print(" ")));
 }
 
 {\\ Executing:
 plotulamspir(9,1); \\ (see output)
 plotulamspir(9,2); \\ (see output)
 plotulamspir(100); \\ ULAMspiral1.png
 plotulamspir(200); \\ ULAMspiral2.png
 }
 
 Output:
 
 > plotulamspir(9,1);
   *** Ulam spiral: 9x9 matrix, p-flag=1
   
           *   *
     *           *
   *   *       *
         *   *   *
       *     * *   *
     *   *
   *       *
     *       *
   *           *
  
 > plotulamspir(9,2);
   *** Ulam spiral: 9x9 matrix, p-flag=2
  
              61    59
     37                31
  67    17          13
            5     3    29
        19        2 11    53
     41     7
  71          23
     43          47
  73                79
  
 > plotulamspir(100); \\ ULAMspiral1.png
   *** Ulam spiral: 101x101 matrix, p-flag=0
   *** matrix(101x101) 1252 DOTS

 > plotulamspir(200); \\ ULAMspiral2.png
   *** Ulam spiral: 201x201 matrix, p-flag=0
   *** matrix(201x201) 4236 DOTS

Fibonacci word/fractal

• [23] Fibonacci word/fractal

Version #1.

In this version only function plotfibofract() was translated from C++, plus upgraded to plot different kind/size of Fibonacci word/fractals.

Output Fibofrac1.png

Output Fibofrac2.png

Works with PARI/GP 2.7.4 and above

 \\ Fibonacci word/fractals
 \\ 4/25/16 aev
 fibword(n)={
 my(f1="1",f2="0",fw,fwn,n2);
 if(n<=4, n=5);n2=n-2;
 for(i=1,n2, fw=Str(f2,f1); f1=f2;f2=fw;); fwn=#fw;
 fw=Vecsmall(fw);
 for(i=1,fwn,fw[i]-=48);
 return(fw);
 }
 
 nextdir(n,d)={
 my(dir=-1);
 if(d==0, if(n%2==0, dir=0,dir=1)); \\0-left,1-right
 return(dir);
 }
 
 plotfibofract(n,sz,len)={
 my(fwv,fn,dr,px=10,py=420,x=0,y=-len,g2=0,
    ttl="Fibonacci fractals: n=");
 plotinit(0); plotcolor(0,6); \\green 
 plotscale(0, -sz,sz, -sz,sz);
 plotmove(0, px,py);
 fwv=fibword(n); fn=#fwv;
 for(i=1,fn, 
     plotrline(0,x,y);
     dr=nextdir(i,fwv[i]);
     if(dr==-1, next);
     \\up
     if(g2==0, y=0; if(dr, x=len;g2=1, x=-len;g2=3); next);
     \\right
     if(g2==1, x=0; if(dr, y=len;g2=2, y=-len;g2=0); next);
     \\down
     if(g2==2, y=0; if(dr, x=-len;g2=3, x=len;g2=1); next);
     \\left
     if(g2==3, x=0; if(dr, y=-len;g2=0, y=len;g2=2); next);
    );\\fend i 
 plotdraw([0,-sz,-sz]);
 print(" *** ",ttl,n," sz=",sz," len=",len," fw-len=",fn);
 }
 
 {\\ Executing:
 plotfibofract(11,430,20); \\ Fibofrac1.png
 plotfibofract(21,430,2);  \\ Fibofrac2.png
 }
 
 Output:
 
 > plotfibofract(11,430,20); \\ Fibofrac1.png
  *** Fibonacci fractals: n=11 sz=430 len=20 fw-len=89
 
 > plotfibofract(21,430,2);  \\ Fibofrac2.png
  *** Fibonacci fractals: n=21 sz=430 len=2 fw-len=10946

Version #2.

In this version only function plotfibofract1() was translated from Java, plus upgraded to plot different kind/size of Fibonacci word/fractals.

Output Fibofrac3.png

Output Fibofrac4.png

Works with PARI/GP 2.7.4 and above

 \\ Fibonacci word/fractals
 \\ 4/26/16 aev
 fibword(n)={
 my(f1="1",f2="0",fw,fwn,n2); \\check n2 in v2 ADD it!!
 if(n<=4, n=5); n2=n-2;
 for(i=1,n2, fw=Str(f2,f1); f1=f2;f2=fw;); fwn=#fw;
 fw=Vecsmall(fw);
 for(i=1,fwn,fw[i]-=48);
 return(fw);
 }

 plotfibofract1(n,sz,len)={
 my(fwv,fn,dx=len,dy=0,nr,ttl="Fibonacci word/fractal, n=");
 plotinit(0); plotcolor(0,5); \\red 
 plotscale(0, -sz,sz, -sz,sz); plotmove(0, 0,0);
 fwv=fibword(n); fn=#fwv;
 for(i=1,fn, plotrline(0,dx,dy);
     if(fwv[i]==0, tx=dx; nr=i%2; if(!nr,dx=-dy;dy=tx, dx=dy;dy=-tx));
    );\\fend i 
 plotdraw([0,0,0]); 
 print(" *** ",ttl,n," sz=",sz," len=",len," fw-len=",fn);
 }

 {\\ Executing:
 plotfibofract1(17,500,6); \\ Fibofrac3.png
 plotfibofract1(21,600,1); \\ Fibofrac4.png
 }

 Output:
 > plotfibofract1(17,500,6); \\ Fibofrac3.png
  *** Fibonacci word/fractal: n=17 sz=500 len=6 fw-len=1597
 > plotfibofract1(21,600,1); \\ Fibofrac4.png
  *** Fibonacci word/fractal: n=21 sz=600 len=1 fw-len=10946

Brownian tree

• [24] Brownian tree (with pictures on Rosettacode Wiki)

All versions #1 - #4 are based on using 2 small plotting helper functions, which are allowing to unify all BrownianTreeX() functions and make them shorter.

Works with PARI/GP 2.7.4 and above

Plotting helper functions

 \\ 2 plotting helper functions 3/2/16 aev
 \\ insm(): x,y are inside matrix mat (+/- p deep).
 insm(mat,x,y,p=0)={my(xz=#mat[1,],yz=#mat[,1]);
   return(x+p>0 && x+p<=xz && y+p>0 && y+p<=yz && x-p>0 && x-p<=xz && y-p>0 && y-p<=yz)}
 \\ plotmat(): Simple plotting using matrix mat (filled with 0/1).
 plotmat(mat)={
 my(xz=#mat[1,],yz=#mat[,1],vx=List(),vy=vx,x,y);
 for(i=1,yz, for(j=1,xz, if(mat[i,j]==0, next, listput(vx,i); listput(vy,j))));
 plothraw(Vec(vx),Vec(vy));
 print(" *** matrix(",xz,"x",yz,") ",#vy, " DOTS");
 }
 

Version #1. Translated from AutoHotkey.

 \\ 3/8/2016
 BrownianTree1(size,lim)={
 my(Myx=matrix(size,size),sz=size-1,sz2=sz\2,x,y,ox,oy);
 x=sz2; y=sz2; Myx[y,x]=1;  \\ seed in center
 print(" *** START: ",x,"/",y);
 for(i=1,lim,
   x=random(sz)+1; y=random(sz)+1;
   while(1,
     ox=x; oy=y;
     x+=random(3)-1; y+=random(3)-1;
     if(insm(Myx,x,y)&&Myx[y,x],
        if(insm(Myx,ox,oy), Myx[oy,ox]=1; break));
     if(!insm(Myx,x,y), break);
   );\\wend
 );\\ fend i
 plotmat(Myx);
 }
 {\\ Executing:
 BrownianTree1(400,15000);
 }
 
 Output:
  *** START: 199/199
  *** matrix(400x400) 3783 DOTS
  ***   last result computed in 26min, 1,295 ms.
 

Version #2. Translated from Octave.

 \\ 3/11/2016
 BrownianTree2(size,lim)={
 my(Myx=matrix(size,size),sz=size-1,dx,dy,x,y);
 x=random(sz); y=random(sz); Myx[y,x]=1; \\ random seed
 print(" *** START: ",x,"/",y);
 for(i=1,lim,
   x=random(sz)+1; y=random(sz)+1;
   while(1,
     dx=random(3)-1; dy=random(3)-1;
     if(!insm(Myx,x+dx,y+dy), x=random(sz)+1; y=random(sz)+1,
        if(Myx[y+dy,x+dx], Myx[y,x]=1; break, x+=dx; y+=dy));
   );\\wend
 );\\fend i
 plotmat(Myx);
 }
 {\\ Executing:
 BrownianTree2(1000,3000);
 }
 
 Output:
  *** START: 697/753
  *** matrix(1000x1000) 2984 DOTS
  ***   last result computed in 6h, 41min, 15,610 ms.
 

Version #3. Translated from Seed7.

 \\ 3/14/2016
 BrownianTree3(size,lim)={
 my(Myx=matrix(size,size),sz=size-2,x,y,dx,dy,b=0);
 x=random(sz); y=random(sz); Myx[y,x]=1; \\ random seed
 print("*** START: ", x,"/",y);
 for(i=1,lim,
     x=random(sz); y=random(sz);
     b=0; \\ bumped not
     while(!b,
        dx=random(3)-1; dy=random(3)-1;
        if(!insm(Myx,x+dx,y+dy), x=random(sz); y=random(sz),
           if(Myx[y+dy,x+dx]==1, Myx[y,x]=1; b=1, x+=dx; y+=dy);
          );
     );\\wend
 );\\fend i
 plotmat(Myx);
 }
 {\\ Executing:
 BrownianTree3(400,5000);
 }
 
 Output:
 *** Start: 275/183
 *** matrix(400x400) 4859 DOTS
 ***   last result computed in 18min, 30,110 ms.
 

Version #4. Translated from PureBasic.

 \\ 3/17/2016
 \\ s=1/2(random seed/seed in the center); p=0..n (level of the "deep" checking).
 BrownianTree4(size,lim,s=1,p=0)={
 my(Myx=matrix(size,size),sz=size-3,x,y);
 \\ seed s=1 for BTPB1, s=2 for BTPB2, BTPB3
 if(s==1,x=random(sz); y=random(sz), x=sz\2; y=sz\2); Myx[y,x]=1; 
 print(" *** START: ",x,"/",y);
 for(i=1,lim,
   if(!(i==1&&s==2), x=random(sz)+1; y=random(sz)+1);
   while(insm(Myx,x,y,1)&&
         (Myx[y+1,x+1]+Myx[y+1,x]+Myx[y+1,x-1]+Myx[y,x+1]+
          Myx[y-1,x-1]+Myx[y,x-1]+Myx[y-1,x]+Myx[y-1,x+1])==0,
     x+=random(3)-1; y+=random(3)-1;
     \\ p=0 for BTPB1, BTPB2; p=5 for BTPB3
     if(!insm(Myx,x,y,p), x=random(sz)+1; y=random(sz)+1;);
   );\\wend
   Myx[y,x]=1;  
 );\\fend i
 plotmat(Myx);
 }
 
 {
 \\ Executing:
 BrownianTree4(200,4000);      \\BTPB1.png
 BrownianTree4(200,4000,2);    \\BTPB2.png
 BrownianTree4(200,4000,2,5);  \\BTPB3.png
 }
 
 Output:
 > BrownianTree4(200,4000)
 *** START: 4/4
 *** matrix(200x200) 3824 DOTS
 ***   last result computed in 1min, 17,891 ms.
 
 > BrownianTree4(200,4000,2)
 *** START: 98/98
 *** matrix(200x200) 3806 DOTS
 ***   last result computed in 38,580 ms.

 > BrownianTree4(200,4000,2,5)
 *** START: 98/98
 *** matrix(200x200) 3645 DOTS
 ***   last result computed in 1min, 12,109 ms.

Fractal tree

• [25] Fractal tree

Output FracTree1.gif

Output FracTree2.gif

Output FracTree3.gif

This version with recursion, in general, is a translation of JavaScript version. Some tweaks and options were added to make it reusable and outputting different size of a tree.

Works with PARI/GP 2.7.4 and above

 \\ Fractal tree (w/recursion)
 \\ 4/10/16 aev
 plotline(x1,y1,x2,y2)={plotmove(0, x1,y1);plotrline(0,x2-x1,y2-y1);}
 
 plottree(x,y,a,d)={
 my(x2,y2,d2r=Pi/180.0,a1=a*d2r,d1);
 if(d<=0, return(););
 if(d>0, d1=d*10.0;
    x2=x+cos(a1)*d1;
    y2=y+sin(a1)*d1;
    plotline(x,y,x2,y2);
    plottree(x2,y2,a-20,d-1);
    plottree(x2,y2,a+20,d-1),
    return();
   );
 }
 
 FractalTree(depth,size)={
 my(dx=1,dy=0,ttlb="Fractal Tree, depth ",ttl=Str(ttlb,depth));
 print1(" *** ",ttl); print(", size ",size);
 plotinit(0);
 plotcolor(0,6); \\green
 plotscale(0, -size,size, 0,size); 
 plotmove(0, 0,0);
 plottree(0,0,90,depth);
 plotdraw([0,size,size]);
 }
 
 {\\ Executing:
 FractalTree(9,500);     \\FracTree1.gif
 FractalTree(12,1100);   \\FracTree2.gif
 FractalTree(15,1500);   \\FracTree3.gif
 }

 Output:
  *** Fractal Tree, depth 9, size 500
  ***   last result computed in 140 ms.
 
  *** Fractal Tree, depth 12, size 1100
  ***   last result computed in 236 ms.

  *** Fractal Tree, depth 15, size 1500
  ***   last result computed in 1,095 ms

Pythagoras tree

• [26] Pythagoras_tree

Output PythTree1.png

This version with recursion, in general, is a translation of zkl version. Almost "as is", so, outputting upside-down tree. This was fixed manually.

Works with PARI/GP 2.7.4 and above

 \\ Pythagoras Tree (w/recursion)
 \\ 4/11/16 aev
 plotline(x1,y1,x2,y2)={plotmove(0, x1,y1);plotrline(0,x2-x1,y2-y1);}
 
 pythtree(ax,ay,bx,by,d=0)={
 my(dx,dy,x3,y3,x4,y4,x5,y5);
 if(d>10, return());
 dx=bx-ax; dy=ay-by;
 x3=bx-dy; y3=by-dx;
 x4=ax-dy; y4=ay-dx;
 x5=x4+(dx-dy)\2; y5=y4-(dx+dy)\2;
 plotline(ax,ay,bx,by);
 plotline(bx,by,x3,y3);
 plotline(x3,y3,x4,y4);
 plotline(x4,y4,ax,ay);
 pythtree(x4,y4,x5,y5,d+1);
 pythtree(x5,y5,x3,y3,d+1);
 }
 
 PythagorTree(x1,y1,x2,y2,depth=9,size)={
 my(dx=1,dy=0,ttlb="Pythagoras Tree, depth ",ttl=Str(ttlb,depth));
 print1(" *** ",ttl); print(", size ",size);
 print(" *** Start: ",x1,",",y1,",",x2,",",y2);
 plotinit(0);
 plotcolor(0,6); \\green
 plotscale(0, -size,size, 0,size); 
 plotmove(0, 0,0);
 pythtree(x1,y1, x2,y2);
 plotdraw([0,size,size]);
 }
 
 {\\ Executing:
 PythagorTree(275,500,375,500,9,640);    \\PythTree1.png
 }
 
 Output:
 
  *** Pythagoras Tree, depth 9, size 640
  *** Start: 275,500,375,500