API tools faq
paste
A-KouZ1

Mandelbrot Perturbation

A-KouZ1
May 11th, 2020
141
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Java 16.15 KB | None | 0 0
raw download clone embed print report
  1. // Mandelbrot Perturbation Demonstration written by Kouzeru,
  2. // credits goes to pianomiles, and CSA in partial contribution.
  3. // sukop, implementation of double-double arithmetics reference.
  4.  
  5. // Instructions:
  6. // Use mouse, hover to the spot and scroll to zoom, click to drag.
  7.  
  8. // Increase roughness as far as you zoomed to see more in details,
  9. // as it costs more computation work, keep the roughness low as
  10. // possible while keeping details.
  11.  
  12. // When the image goes pixelated or glitched, right-click to it.
  13. // Or just switch precision mode into appropriate choice.
  14. // Increase the color speed when the spots hard to distinguish.
  15.  
  16. float zoom,scl,maxIT,cmag,bailOut,logBOut,log2,deg,fpfAvg=0;
  17. int lMillis,mHover,frames,maxITl,zooml,refIT,precMode;
  18. double[] xPos,yPos,xRef,yRef,xRefO,yRefO;
  19. boolean infoShow;String precName[]={
  20.   "Single-Precision delta","Double-Precision delta",
  21.   "Single-Precision main","Double-Precision main",
  22.   "Quadruple-Precision main"};
  23. PFont font;
  24.  
  25. //Double-double precision arithmetic methods.
  26. //I'm going to call them quad here to avoid confusion.
  27. //A reference which is actually rewritten from Python to Java:
  28. //https://github.com/sukop/doubledouble/blob/master/doubledouble.py
  29. static class dd{private static double r,e,s,t,u,v,f,g;
  30.     private static double[] qs(double x,double y){
  31.     r = x+y;e = y-(r-x);double[] z={r,e}; return z;}
  32.     private static void ts(double x,double y){
  33.     r = x+y;t = r-x;e = x-(r-t);e+= y-t;}
  34.     private static void pr(double x,double y){
  35.     u = x*134217729;v = y*134217729;
  36.         u-= u-x; v-= v-y; f = x-u; g = y-v;
  37.     r = x*y; e = u*v-r; e+= u*g+f*v; e+= f*g; }
  38.     private static void sq(double x){
  39.     u = x*134217729; u-= u-x; f = x-u;
  40.     r = x*x; e = u*u-r; e+= 2*u*f; e+= f*f; }
  41.     static double[] add(double[] x,double[] y){
  42.     ts(x[0], y[0]); e+=x[1]+y[1];return qs(r,e);}
  43.     static double[] sub(double[] x,double[] y){
  44.     ts(x[0],-y[0]); e+=x[1]-y[1];return qs(r,e);}
  45.     static double[] mul(double[] x,double[] y){
  46.     pr(x[0],y[0]); e+=x[0]*y[1]+x[1]*y[0]; return qs(r,e);}
  47.     static double[] div(double[] x,double[] y){
  48.     t = x[0]/y[0]; pr(t,y[0]); s = x[0]-r-e+x[1];
  49.     s = (s-t*y[1])/y[0]; return qs(t,s);}
  50.     static double[] squ(double[] x){
  51.         sq(x[0]); e+=2*x[0]*x[1]; return qs(r,e);}
  52.     static double[] sqr(double[] x){
  53.         if(x[0]<=0){ if(x[0]<0||x[1]<0){
  54.         double[] z={Double.NaN,Double.NaN}; return z;}}
  55.         t = java.lang.Math.sqrt(x[0]); sq(x[0]);
  56.         s = x[0]-r-e+x[1]; s/= t+t; return qs(t,s);  };
  57.     static double[] neg(double[] x){
  58.         double[] z={-x[0],-x[1]}; return z;}
  59.     static double[] abs(double[] x){double[] z={0,0};
  60.         if(x[0]<0){z[0]=-x[0]; z[1]=-x[1];
  61.         }else{z[0]=x[0]; z[1]=x[1];} return z;}
  62.     static double[] doub(double[] x){
  63.         double[] z = { x[0]*2, x[1]*2};return z;}
  64.     static double[] half(double[] x){
  65.         double[] z = { x[0]/2, x[1]/2};return z;}
  66.     static double[] copy(double[] x){
  67.         double[] z = { x[0], x[1]};return z;}
  68.     static double[] parse(double x){
  69.         double[] z = { x,0 }; return z;}
  70.     //AAAAA BASIQUE DOUBLE DOUBLE TO STRING VICEVERSA CONVERSION
  71.     static boolean dSgn; static int dExp; static int[] dDgt;
  72.     private static int p(char x){ //just to make the code short
  73.         return Integer.parseInt(String.valueOf(x));}
  74.     private static void digitsStr(String x){
  75.         int[] d=new int[40]; int b,e=d.length,f,g=0,h=0,
  76.         i=x.charAt(0)=='-'||x.charAt(0)=='+'?1:0,j=x.lastIndexOf('e');
  77.         if(j<0){j=x.lastIndexOf('E');if(j<0){j=x.length();}}
  78.         f=x.lastIndexOf('.');f=f<0?j:f;
  79.         while(i<f){if(p(x.charAt(i))>0){break;}i++;}
  80.         while(i<f){d[h++]=p(x.charAt(i++));g++;}
  81.         if(g>0){g--;}if(f<j){i++;}
  82.         if(h==0){while(i<j){g--;if(p(x.charAt(i))>0){
  83.             break;}i++;}if(i==j){g=0;}}
  84.         while(i<j&&h<e){d[h++]=p(x.charAt(i++));}
  85.         while(h<e){d[h++]=0;}i=j+1;j=x.length();b=0;if(i<j){
  86.         if(x.charAt(i++)=='-'){while(i<j){b=b*10-p(x.charAt(i++));}
  87.           }else{while(i<j){b=b*10+p(x.charAt(i++));}}}
  88.         dSgn=x.charAt(0)=='-';dExp=b+g;dDgt=d;
  89.     }
  90.     private static String digitsFixed(int s,int e,int[] d){
  91.         int i=d.length-1,j,k;d[--i]|=d[i+1]<5?0:1;
  92.         while(d[i]==0&&i>0){i--;};
  93.         String a=(s!=0)?"-":"";int b=0,c=i,f=0;
  94.         if(e<0){a+="0.";f++;while(++e!=0&&b<c){a+="0";b++;f++;}b=0;
  95.         }else{e++;while(b<e&b<c){a+=d[b++];f++;}
  96.         while(b++<e&b<c){a+="0";f++;}a+=".";}
  97.         while(b<c){a+=d[b++];f++;}return a;
  98.     }
  99.     static double[] parse(String x){
  100.         int g,h,i,j; double a=Double.parseDouble(x),b=0,k;
  101.         digitsStr(String.format("%1.40e",a<0?-a:a));
  102.         int e=dExp; int[] c=dDgt; digitsStr(x); int[] d=dDgt;
  103.         for(g=0,h=i=c.length-1;i>=0;i--){
  104.         j=d[i]-c[i]-g+10;g=j<10?1:0;c[i]=j%10;}
  105.         if(g!=0){k=-1;while(++i<h){c[i]=9-c[i];}
  106.         c[i]=10-c[i];}else{k=1;}
  107.         for(j=0;j<h;j++){if(c[j]>0){break;}}
  108.         b=0;while(j<=h){b+=c[h--]*k;k*=10;}
  109.         k/=10;e-=j;i=0;if(e<0){while(i-->e){
  110.         k*=10;}}else{while(i++<e){k/=10;}}
  111.         b=(a<0?-b:b)/(k<0?-k:k);
  112.         double[] z={a,b}; return z;
  113.     }
  114.     static String string(double[] x){
  115.         digitsStr(String.format("%1.40g",x[0]));
  116.         boolean a=dd.dSgn;int c=dd.dExp; int[] e=dDgt;
  117.         digitsStr(String.format("%1.40g",x[1]));
  118.         boolean b=dd.dSgn;int d=dd.dExp; int[] f=dDgt;
  119.         int g=0,h=e.length,i,j;
  120.         if(f[0]>0){d=c-d;i=h;if(a!=b){
  121.             while(--i>=0){j=i>=d?f[i-d]:0;
  122.             j=e[i]-j-g+10;g=j<10?1:0;f[i]=j%10;}
  123.             if(g!=0){while(++i<h){f[i]=9-f[i];}
  124.                 f[i-1]++;}i=j=0;
  125.             while(i<h){if(f[i]>0){break;}i++;}c-=i;
  126.             }else{while(--i>=0){j=i>=d?f[i-d]:0;
  127.             j=e[i]+j-g+1 ;g=j<10?1:0;f[i]=j%10;}
  128.             i=j=0;if(g==0){f[0]=1;j++;c++;}g=0;}
  129.             while(i<h){e[j++]=f[i++];}
  130.             while(j<h){e[j++]=0;}
  131.         }return digitsFixed((a?1:0)^g,c,e);
  132.     };
  133. }
  134.  
  135. int[] HSV2RGB(float h,float s,float v){
  136.     float r,g,b,k;s=s<0?0:s<1?s:1;v=v<0?0:v<1?v:1;
  137.         k=0.5-sin(deg*(((((h+=3)%1)*60)-30)));
  138.         r=(h+5)%6<3?1:0;r+=(((h+4)%6<3?1:0)-r)*k;
  139.         g=(h+3)%6<3?1:0;g+=(((h+2)%6<3?1:0)-g)*k;
  140.         b=(h+1)%6<3?1:0;b+=(((h+0)%6<3?1:0)-b)*k;
  141.         r=1+(r-1)*s;g=1+(g-1)*s;b=1+(b-1)*s;
  142.         r*=255*v;g*=255*v;b*=255*v;
  143.     int[] col={(int)r,(int)g,(int)b};
  144.     return col;
  145. }
  146.  
  147. int setColour(float IT){if(IT<maxIT){
  148.         IT = IT>0?sqrt(IT*cmag):0;
  149.         float col = (sin(deg*IT*10)+1)/2;
  150.         float br = (cos(deg*IT*40) + 1)/2;
  151.         float h = (((col+0.05))%1)*6;
  152.         float s = min(br*2,1);
  153.         float v = (1-(br*br));
  154.         int[] c = HSV2RGB(h,s,v);
  155.         return color(c[0],c[1],c[2]);
  156.     }else{return color(0);}
  157. }
  158.  
  159. double[] refX,refY;
  160. float[] refXf,refYf;
  161.  
  162. //precalculate single point with quad precision
  163. float precalc(double[] x,double[] y){
  164.     double[] a=dd.copy(x),b=dd.copy(y),
  165.     a2=dd.squ(a),b2=dd.squ(b);int IT=0;
  166.     if(precMode==0){
  167.         refXf= new float[(int) maxIT+4];
  168.         refYf= new float[(int) maxIT+4];
  169.         refXf[0]=(float)a[0];refYf[0]=(float)b[0];
  170.     }else{
  171.         refX = new double[(int) maxIT+4];
  172.         refY = new double[(int) maxIT+4];
  173.         refX[0]=a[0];refY[0]=b[0];
  174.     }
  175.     while(IT++<maxIT){b=dd.add(dd.doub(dd.mul(a,b)),y);
  176.         a=dd.add(dd.sub(a2,b2),x);a2=dd.squ(a);b2=dd.squ(b);
  177.         if(precMode!=0){refX[IT]=a[0];refY[IT]=b[0];}else{
  178.             refXf[IT]=(float)a[0];refYf[IT]=(float)b[0];}
  179.         if(a2[0]+b2[0]>=4){break;}
  180.     }return (refIT=IT);
  181. }
  182. //precalculate single point with double precision
  183. float precalcd(double x,double y){
  184.     double a=x,b=y,a2=a*a,b2=b*b;int IT=0;
  185.     if(precMode==0){
  186.         refXf= new float[(int) maxIT+4];
  187.         refYf= new float[(int) maxIT+4];
  188.         refXf[0]=(float)a;refYf[0]=(float)b;
  189.     }else{
  190.         refX = new double[(int) maxIT+4];
  191.         refY = new double[(int) maxIT+4];
  192.         refX[0]=a; refY[0]=b;
  193.     }
  194.     while(IT++<maxIT){b=2*a*b+y;a=a2-b2+x;
  195.         if(precMode!=0){refX[IT]=a;refY[IT]=b;}else{
  196.             refXf[IT]=(float)a;refYf[IT]=(float)b;}
  197.         if((a2=a*a)+(b2=b*b)>=4){break;}
  198.     }return (refIT=IT);
  199. }
  200. //iterate delta with double precision
  201. float iteratedd(double x,double y){
  202.     int IT=0;x+=xRefO[0];y+=yRefO[0];
  203.     double a=x,b=y,e,f,c=refX[0],d=refY[0],
  204.     g=c+x,h=d+y,i,j;while(IT++<refIT){
  205.         i =a*c;i-=b*d;j =a+b;j*=a-b;
  206.         b*=a+c;b+=a*d;b+=b+y;a =i+i+j+x;
  207.         c=refX[IT];d=refY[IT];
  208.         e=a+c;f=b+d;i=e*e;j=f*f;
  209.         if(i+j>=bailOut){
  210.             i=log((float)(i+j));e-=g;f-=h;
  211.             j=log((float)(e*e+f*f))/2;
  212.             i=(logBOut-j)/(i-j);
  213.             return (float)IT+log((float)i+1)/log2;
  214.     }}return IT;
  215. }
  216. //iterate delta with single precision
  217. //wait wait, even with less type conversion,
  218. //it still running slow than double, what?
  219. float iterateds(double dx,double dy){
  220.     int IT=0;dx+=xRefO[0];dy+=yRefO[0];
  221.     float x=(float)dx,y=(float)dy,a=x,b=y,
  222.     e,f,c=refXf[0],d=refYf[0],
  223.     g=c+x,h=d+y,i,j;while(IT++<refIT){
  224.         i =a*c;i-=b*d;j =a+b;j*=a-b;
  225.         b*=a+c;b+=a*d;b+=b+y;a =i+i+j+x;
  226.         c=refXf[IT];d=refYf[IT];
  227.         e=a+c;f=b+d;i=e*e;j=f*f;
  228.         if(i+j>=bailOut){
  229.             i=log(i+j);e-=g;f-=h;
  230.             j=log(e*e+f*f)/2;
  231.             i=(logBOut-j)/(i-j);
  232.             return (float)IT+log(i+1)/log2;
  233.     }}return IT;
  234. }
  235. //iterate the set with quad precision
  236. float iterateq(double dx,double dy){
  237.     double[] x = dd.add(dd.parse(dx),xPos),
  238.     y = dd.add(dd.parse(dy),yPos),
  239.     a=dd.copy(x),b=dd.copy(y),
  240.     a2=dd.squ(a),b2=dd.squ(b);
  241.     float m,n,o,p,q;int IT=0;
  242.     while(IT++<maxIT){
  243.         b=dd.add(dd.doub(dd.mul(a,b)),y);
  244.         a=dd.add(dd.sub(a2,b2),x);
  245.         a2=dd.squ(a);b2=dd.squ(b);
  246.         if(a2[0]+b2[0]>=bailOut){
  247.             m=(float)(a[0]-x[0]);
  248.             n=(float)(b[0]-y[0]);
  249.             p=log((float)(a2[0]+b2[0]));
  250.             q=log(m*m+n*n)/2;
  251.             o=(logBOut-q)/(p-q);
  252.             return (float)IT+log(o+1)/log2;
  253.         }
  254.     }
  255.     return IT;
  256. }
  257. //iterate the set with double precision
  258. float iterated(double x,double y){
  259.     int IT=0; x+=xPos[0];y+=yPos[0];
  260.     double a=x,b=y,a2=a*a,b2=b*b;
  261.     while(IT++<maxIT){
  262.         b=2*a*b+y;a=a2-b2+x;
  263.         if((a2=a*a)+(b2=b*b)>bailOut){
  264.             a2=log((float)(a2+b2));a-=x;b-=y;
  265.             b2=log((float)(a*a+b*b))/2;
  266.             a2=((double)logBOut-b2)/(a2-b2);
  267.             return (float)IT+log((float)a2+1)/log2;
  268.     }}return IT;
  269. }
  270. //iterate the set with single precision
  271. float iterates(double xd,double yd){
  272.     int IT=0; xd+=xPos[0];yd+=yPos[0];
  273.     float x=(float)xd,y=(float)yd;
  274.     float a=x,b=y,a2=a*a,b2=b*b;
  275.     while(IT++<maxIT){b=2*a*b+y;a=a2-b2+x;
  276.         if((a2=a*a)+(b2=b*b)>bailOut){
  277.             a2=log(a2+b2);a-=x;b-=y;
  278.             b2=log(a*a+b*b)/2;
  279.             a2=(logBOut-b2)/(a2-b2);
  280.             return (float)IT+log(a2+1)/log2;
  281.     }}return IT;
  282. }
  283.  
  284. void setup(){
  285.     size(640,480);background(0);
  286.     font = createFont("Monospaced.bold",18);
  287.     deg = atan(1)*4/180;
  288.     lMillis=millis();
  289.     scl = height;
  290.     logBOut = log(4);
  291.     log2 = log(2);
  292.     cmag = 1;
  293.     mHover = 0;
  294.     precMode = 1;
  295.     infoShow = true;
  296.     xPos = dd.parse( "0.26544602163041845079010564709539");
  297.     yPos = dd.parse("-0.003176395242688470419416969446745");
  298.     xRef = dd.parse( "0.26544602163041845079010785861635");
  299.     yRef = dd.parse("-0.0031763952426884704194213468422653");
  300.     zooml= (int)(16*log(1.4432772E-23)/log2);
  301.     maxITl= (int)(16*log(6400)/log2);
  302.     zoom = exp((float)(zooml*log2/16));
  303.     loadPixels();//frameRate(60);
  304.     textFont(font,18);
  305. }
  306.  
  307. void mouseWheel(MouseEvent e){
  308.     zooml+=e.getCount();
  309. };
  310.  
  311. void draw(){
  312.     if(mousePressed){switch(mouseButton){
  313.         case LEFT:{
  314.             xPos =dd.add(xPos,dd.parse((double)(pmouseX-mouseX)*zoom/scl));
  315.             yPos =dd.sub(yPos,dd.parse((double)(pmouseY-mouseY)*zoom/scl));
  316.             break;
  317.         }
  318.         case RIGHT:{
  319.             if(precMode<2){
  320.                 xRef=dd.add(xPos,dd.parse((double)((mouseX-width/2)/scl)*zoom));
  321.                 yRef=dd.sub(yPos,dd.parse((double)((mouseY-height/2)/scl)*zoom));
  322.             }
  323.             break;
  324.         }
  325.     }}
  326.     if(keyPressed){switch(key){
  327.         case 'z':{maxITl-=2;break;}
  328.         case 'x':{maxITl+=2;break;}
  329.         case 'a':{logBOut-=0.125;break;}
  330.         case 's':{logBOut+=0.125;break;}
  331.         case 'q':{cmag/=1.125;break;}
  332.         case 'w':{cmag*=1.125;break;}
  333.         case 'k':{precMode=0;break;} //use single delta
  334.         case 'l':{precMode=1;break;} //use double delta
  335.         case 'i':{precMode=2;break;} //use single iter
  336.         case 'o':{precMode=3;break;} //use double iter
  337.         case 'p':{precMode=4;break;} //use quad iter
  338.     }}
  339.     maxIT=exp((float)(maxITl*log2/16));bailOut = exp(logBOut);
  340.     float lzoom=zoom;zoom=exp((float)(zooml*log2/16));
  341.     double[] xTar = dd.add(xPos,dd.parse((double)(lzoom*((mouseX-width/2)/scl))));
  342.     double[] yTar = dd.add(yPos,dd.parse((double)(lzoom*((height/2-mouseY)/scl))));
  343.     xPos = dd.add(dd.mul(dd.sub(xPos,xTar),dd.parse(zoom/lzoom)),xTar);
  344.     yPos = dd.add(dd.mul(dd.sub(yPos,yTar),dd.parse(zoom/lzoom)),yTar);
  345.     xRefO= dd.sub(xPos,xRef);yRefO= dd.sub(yPos,yRef);
  346.     float rx = scl*(float)dd.sub(xRef,xPos)[0]/zoom+width/2;
  347.     float ry = height/2-scl*(float)dd.sub(yRef,yPos)[0]/zoom;
  348.    
  349.     int framesPerFrame = -frames;
  350.     if(precMode<2){precalc(xRef,yRef);}
  351.     for(lMillis=millis()+50;millis()<lMillis;){
  352.         int fr = frames++&255;
  353.         int sx = (fr>>0&1)<<3|(fr>>2&1)<<2|(fr>>4&1)<<1|(fr>>6&1)<<0;
  354.         int sy =((fr>>1&1)<<3|(fr>>3&1)<<2|(fr>>5&1)<<1|(fr>>7&1)<<0)^sx;
  355.         for(int j=sy;j<height;j+=16){
  356.             for(int i=sx;i<width;i+=16){
  357.                 float x=zoom*(i-width/2)/scl;
  358.                 float y=zoom*(height/2-j)/scl;
  359.                 float IT=0;
  360.                 switch(precMode){
  361.                     case 0: IT = iterateds(x,y); break;
  362.                     case 1: IT = iteratedd(x,y); break;
  363.                     case 2: IT = iterates(x,y); break;
  364.                     case 3: IT = iterated(x,y); break;
  365.                     case 4: IT = iterateq(x,y); break;
  366.                 }
  367.                 pixels[i+j*width]=setColour(IT);
  368.             }
  369.         }
  370.     }
  371.     framesPerFrame += frames;
  372.     fpfAvg+=((float)framesPerFrame-fpfAvg)/4;
  373.     updatePixels();
  374.     for(int i=0;i<2;i++){fill(240*i);
  375.         if(infoShow){
  376.             textAlign(LEFT,TOP);
  377.             text("Real: "+(xPos[0]<0?"":" ")+dd.string(xPos), 4-i,  4-i);
  378.             text("Imag: "+(yPos[0]<0?"":" ")+dd.string(yPos), 4-i, 20-i);
  379.             text("RefR: "+(xRef[0]<0?"":" ")+dd.string(xRef), 4-i, 36-i);
  380.             text("RefI: "+(yRef[0]<0?"":" ")+dd.string(yRef), 4-i, 52-i);
  381.             text("Magn:  "+Float.toString(zoom), 4-i, 68-i);
  382.             text("Iter:  "+Float.toString(maxIT), 4-i, 84-i);
  383.             text("BOut:  "+Float.toString(bailOut), 4-i, 100-i);
  384.             text("Prec:  "+precName[precMode],4-i, 116-i);
  385.             text("FPFA:  "+String.format("%1.3f",fpfAvg),4-i, 132-i);
  386.             if(mouseX!=pmouseX||mouseY!=pmouseY){mHover=millis()+5000;}
  387.             if(millis()<mHover){textAlign(CENTER,BOTTOM);
  388.                 text("press 'm' key to hide info", width/2-i, height-100-i);
  389.                 text("scroll up/down to zoom out/in", width/2-i, height-84-i);
  390.                 text("left-click drag to move the view", width/2-i, height-68-i);
  391.                 text("z/x to increase/decrease roughness", width/2-i, height-52-i);
  392.                 text("right-click to choose reference point", width/2-i, height-36-i);
  393.                 text("a/s for bailOut, q/w for the color speed", width/2-i, height-20-i);
  394.                 text("k/l/i/o/p to switch between precision mode", width/2-i, height-6-i);
  395.             }
  396.         }
  397.     }
  398.     if(precMode<2){fill(255);stroke(0);strokeWeight(1);
  399.         rectMode(CENTER);rect(rx,ry,3,9);rect(rx,ry,9,3);
  400.         stroke(255);rect(rx,ry,1,3);}
  401. }
  402.  
  403. void keyPressed(){
  404.     switch(key){
  405.         case 'm': infoShow = !infoShow; break;
  406.     }
  407.     switch(keyCode){
  408.         case 32 : background(0); loadPixels(); break;
  409.     }
  410. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!