1+ // jscottpilgrim
2+ // greyscale julia by distance estimator
3+ // distance estimation info: http://www.iquilezles.org/www/articles/distancefractals/distancefractals.htm
4+
5+ // fragment shader
6+
7+ #ifdef GL_ES
8+ precision highp float ;
9+ precision mediump int ;
10+ #endif
11+
12+ uniform vec2 resolution;
13+
14+ uniform vec2 center;
15+ uniform float zoom;
16+
17+ uniform vec2 juliaParam;
18+
19+ const int maxIterations = 250 ;
20+
21+ // double precision not natively supported in webgl/gl_es
22+ // emulated double precision from: https://gist.github.com/LMLB/4242936fe79fb9de803c20d1196db8f3
23+ // performance drops dramatically from single precision. drop number of iterations to balance performance and max zoom. maybe ~300
24+
25+ float times_frc(float a, float b) {
26+ return mix (0.0 , a * b, b != 0.0 ? 1.0 : 0.0 );
27+ }
28+ float plus_frc(float a, float b) {
29+ return mix (a, a + b, b != 0.0 ? 1.0 : 0.0 );
30+ }
31+ float minus_frc(float a, float b) {
32+ return mix (a, a - b, b != 0.0 ? 1.0 : 0.0 );
33+ }
34+ // Double emulation based on GLSL Mandelbrot Shader by Henry Thasler (www.thasler.org/blog)
35+ //
36+ // Emulation based on Fortran-90 double-single package. See http://crd.lbl.gov/~dhbailey/mpdist/
37+ // Addition: res = ds_add(a, b) => res = a + b
38+ vec2 add (vec2 dsa, vec2 dsb) {
39+ vec2 dsc;
40+ float t1, t2, e;
41+ t1 = plus_frc(dsa.x, dsb.x);
42+ e = minus_frc(t1, dsa.x);
43+ t2 = plus_frc(plus_frc(plus_frc(minus_frc(dsb.x, e), minus_frc(dsa.x, minus_frc(t1, e))), dsa.y), dsb.y);
44+ dsc.x = plus_frc(t1, t2);
45+ dsc.y = minus_frc(t2, minus_frc(dsc.x, t1));
46+ return dsc;
47+ }
48+ // Subtract: res = ds_sub(a, b) => res = a - b
49+ vec2 sub (vec2 dsa, vec2 dsb) {
50+ vec2 dsc;
51+ float e, t1, t2;
52+ t1 = minus_frc(dsa.x, dsb.x);
53+ e = minus_frc(t1, dsa.x);
54+ t2 = minus_frc(plus_frc(plus_frc(minus_frc(minus_frc(0.0 , dsb.x), e), minus_frc(dsa.x, minus_frc(t1, e))), dsa.y), dsb.y);
55+ dsc.x = plus_frc(t1, t2);
56+ dsc.y = minus_frc(t2, minus_frc(dsc.x, t1));
57+ return dsc;
58+ }
59+ // Compare: res = -1 if a < b
60+ // = 0 if a == b
61+ // = 1 if a > b
62+ float cmp(vec2 dsa, vec2 dsb) {
63+ if (dsa.x < dsb.x) {
64+ return - 1 .;
65+ }
66+ if (dsa.x > dsb.x) {
67+ return 1 .;
68+ }
69+ if (dsa.y < dsb.y) {
70+ return - 1 .;
71+ }
72+ if (dsa.y > dsb.y) {
73+ return 1 .;
74+ }
75+ return 0 .;
76+ }
77+ // Multiply: res = ds_mul(a, b) => res = a * b
78+ vec2 mul (vec2 dsa, vec2 dsb) {
79+ vec2 dsc;
80+ float c11, c21, c2, e, t1, t2;
81+ float a1, a2, b1, b2, cona, conb, split = 8193 .;
82+ cona = times_frc(dsa.x, split);
83+ conb = times_frc(dsb.x, split);
84+ a1 = minus_frc(cona, minus_frc(cona, dsa.x));
85+ b1 = minus_frc(conb, minus_frc(conb, dsb.x));
86+ a2 = minus_frc(dsa.x, a1);
87+ b2 = minus_frc(dsb.x, b1);
88+ c11 = times_frc(dsa.x, dsb.x);
89+ c21 = plus_frc(times_frc(a2, b2), plus_frc(times_frc(a2, b1), plus_frc(times_frc(a1, b2), minus_frc(times_frc(a1, b1), c11))));
90+ c2 = plus_frc(times_frc(dsa.x, dsb.y), times_frc(dsa.y, dsb.x));
91+ t1 = plus_frc(c11, c2);
92+ e = minus_frc(t1, c11);
93+ t2 = plus_frc(plus_frc(times_frc(dsa.y, dsb.y), plus_frc(minus_frc(c2, e), minus_frc(c11, minus_frc(t1, e)))), c21);
94+ dsc.x = plus_frc(t1, t2);
95+ dsc.y = minus_frc(t2, minus_frc(dsc.x, t1));
96+ return dsc;
97+ }
98+ // create double-single number from float
99+ vec2 set(float a) {
100+ return vec2 (a, 0.0 );
101+ }
102+ float rand(vec2 co) {
103+ // implementation found at: lumina.sourceforge.net/Tutorials/Noise.html
104+ return fract (sin (dot (co.xy ,vec2 (12.9898 ,78.233 ))) * 43758.5453 );
105+ }
106+ vec2 complexMul(vec2 a, vec2 b) {
107+ return vec2 (a.x* b.x - a.y* b.y,a.x* b.y + a.y * b.x);
108+ }
109+ // double complex multiplication
110+ vec4 dcMul(vec4 a, vec4 b) {
111+ return vec4 (sub(mul(a.xy,b.xy),mul(a.zw,b.zw)),add(mul(a.xy,b.zw),mul(a.zw,b.xy)));
112+ }
113+ // double complex addition
114+ vec4 dcAdd(vec4 a, vec4 b) {
115+ return vec4 (add(a.xy,b.xy),add(a.zw,b.zw));
116+ }
117+ // Length of double complex
118+ vec2 dcLength(vec4 a) {
119+ return add(mul(a.xy,a.xy),mul(a.zw,a.zw));
120+ }
121+ // create double-complex from complex float
122+ vec4 dcSet(vec2 a) {
123+ return vec4 (a.x,0 .,a.y,0 .);
124+ }
125+ // create double-complex from complex double
126+ vec4 dcSet(vec2 a, vec2 ad) {
127+ return vec4 (a.x, ad.x,a.y,ad.y);
128+ }
129+ // Multiply double-complex with double
130+ vec4 dcMul(vec4 a, vec2 b) {
131+ return vec4 (mul(a.xy,b),mul(a.wz,b));
132+ }
133+
134+
135+ void main()
136+ {
137+ vec4 uv = dcSet( gl_FragCoord .xy - resolution.xy * 0.5 );
138+ vec4 dcCenter = dcSet( center );
139+
140+ vec4 c = dcSet( juliaParam );
141+ vec4 z = dcAdd( dcMul( uv, set( zoom ) ), dcCenter );
142+ vec4 dz = dcSet( vec2 ( 1.0 , 0.0 ) );
143+ bool escape = false;
144+
145+ float dotZZ = 0.0 ;
146+
147+ for ( int i = 0 ; i < maxIterations; i++ ) {
148+ // z' = 2 * z * z'
149+ // dz = 2.0 * vec2( z.x * dz.x - z.y * dz.y, z.x * dz.y + z.y * dz.x );
150+ dz = dcMul( dcMul( z, dz ), set( 2.0 ) );
151+
152+ // mandelbrot function on z
153+ // z = c + complexSquare( z );
154+ z = dcAdd( c, dcMul( z, z ) );
155+
156+ // escape radius 32^2
157+ // if ( dot( z, z ) > 1024.0 )
158+ // lowered to compensate for emulated double computation speed
159+ dotZZ = z.x * z.x + z.z * z.z; // extract high part
160+ if ( dotZZ > 400.0 )
161+ {
162+ escape = true;
163+ break ;
164+ }
165+ }
166+
167+ vec3 color = vec3 ( 1.0 );
168+
169+ if ( escape )
170+ {
171+ // distance
172+ // d(c) = (|z|*log|z|)/|z'|
173+
174+ float dotDZDZ = dz.x * dz.x + dz.z * dz.z; // extract high part
175+
176+ float d = sqrt ( dotZZ );
177+ d *= log ( d );
178+ d /= sqrt ( dotDZDZ );
179+
180+ d = clamp ( pow ( 4.0 * d, 0.1 ), 0.0 , 1.0 );
181+
182+ color = vec3 ( d );
183+ }
184+ else
185+ {
186+ // set inside points to inside color
187+ color = vec3 ( 0.0 );
188+ }
189+
190+ gl_FragColor = vec4 ( color, 1.0 );
191+ }
0 commit comments