Sunday, 7 January 2018

One of my greatest fears is that one day I'll wake up in a gridlock I predicted ages ago

Einstein field equations (EFE)



describe the four dimensional geometry g associated with the points of spacetime and corresponding stress-energy tensor T. The differential terms are a function of the metric as described below.


Variable g with upper indices (contravariant) is the inverse of the (covariant) metric.


T describes flux of four momentum between different dimensions, including time. The first term (with upper indices tt) is simply the invariant mass-energy. For a single isolated particle the following (in natural units) applies (v is 4-dimensional velocity vector).


Contravariant stress-energy tensor can be converted into covariant form by the following.


EFE are a kind of generalization of Newtons gravitational potential (Poisson equation).


EFE describe a kind of 4-dimensional array of matrices (in naive numerical sense) corresponding to the structure of such spacetime. This sort of block universe can be seen as static. Theoretically it's possible that in this sort of block there could exist "flows" and geometries that form a closed loop in what we would normally call time (allowing a kind of time travel into the past). The naive interpretation of EFE would suggest that if such loops exist, they ought to form consistent spacetime. That is to say in order for them to be a solution to the equations, no grandfather paradoxes can form. You will (for one reason or another) only do such things in the "past" that are consistent with the block as a whole. Quantum mechanics may significantly complicate the situation if something along the lines of these loops are possible, but instead lead to alternative realities. The block of all universes must still be consistent as a whole, but as there are now several alternative realities and bridges between them, the consistency requirement can still apply and allow some particular futures to be consistent with time traveller "causing" them.

How to actually compute a solution to some practical problem with EFE is a topic for some other time.
--
The end justifies the means as long as the means are factored into the end.

If you could actually rationalize your opinion, it wouldn't be an opinion.

I don't really care about subjective changes all that much, it's the objective ones that I find interesting.

Every January I remember my mortality, but then I forget it for a while again.

Sometimes I consider will a flaw.

What to do when you reach a conclusion that is highly likely to be sound, yet highly undesirable?

...it's just not in my nature.

A true explanation must come to an end where no further explanations are required and none can exist, but can it be?
--
BTW. The Great Philosophers by Stephen Law is pretty nice and compact book about philosophy.
--
double g[4][4][N*N*N*N]; /* spacetime metric tensors */
double T[4][4][N*N*N*N]; /* spacetime stress-energy tensors */

int npow(int i) {
  int val = 1;
  for(int j=0; j<i; j++)
    val = val*N;
  return val;
}

/* computes values of individual Christoffel symbols */
double christoffel(int i, int j, int l, int loc) {
  double val = 0.0, d;
  for(int k=0; k<4; k++) {
    d = \
      (g[k][i][loc + npow(j)] - g[k][i][loc - npow(j)]) +  \
      (g[k][j][loc + npow(i)] - g[k][j][loc - npow(i)]) -  \
      (g[i][j][loc + npow(k)] - g[i][j][loc - npow(k)]);
    val = val + g[l][k][loc]*d;
  }
  return 0.5*val;
}

/* computes values of individual components of Ricci curvature tensor */
double ricci(int i, int j, int loc) {
  double val = 0.0;
  double a = 0.0, b = 0.0;
  for(int l=0; l<4; l++) val += \
    christoffel(i, j, l, loc + npow(l)) - christoffel(i, j, l, loc - npow(l));
  for(int l=0; l<4; l++) val -= \
    christoffel(i, l, l, loc + npow(j)) - christoffel(i, j, l, loc - npow(j));
  for(int m=0; m<4; m++)
    for(int l=0; l<4; l++) val += \
      christoffel(i, j, m, loc)*christoffel(m, l, l, loc);
  for(int m=0; m<4; m++)
    for(int l=0; l<4; l++) val -= \
      christoffel(i, l, m, loc)*christoffel(m, j, l, loc);
  return val;
}

/* computes scalar curvature */
double scalar_curvature(int i, int j, int loc) {
  double gi[4][4], A[4][4], val = 0.0;
  double A2[4][4], A3[4][4], trA, trA2, trA3, det = 0.0;

  /* find metric tensor and initialize gi to zero */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++) {
      A[i][j] = g[i][j][loc];
      gi[i][j] = 0.0;
    }

  /* compute g^2 and g^3 */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      for(int k=0; k<4; k++)
A2[i][j] += A[i][k]*A[k][j]; 
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      for(int k=0; k<4; k++)
A3[i][j] += A2[i][k]*A[k][j]; 

  /* use Cayley-Hamilton method to compute inverse of g */
  trA  =  A[0][0] +  A[1][1] +  A[2][2] +  A[3][3];
  trA2 = A2[0][0] + A2[1][1] + A2[2][2] + A2[3][3];
  trA3 = A3[0][0] + A3[1][1] + A3[2][2] + A3[3][3];
  for(int i=0; i<4; i++) {
    gi[i][i] = (trA*trA*trA - 3*trA*trA2 + 2*trA3)/6.0;
    for(int j=0; j<4; j++)
      gi[i][j] += A2[i][j]*trA - A3[i][j] - 0.5*A[i][j]*(trA*trA-trA2);
  }

  /* division of gi such that (g*gi)_{00} = 1 */
  for(int k=0; k<4; k++)
    det += A[0][k]*gi[k][0];
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      gi[i][j] /= det;

  /* compute Ricci scalar */
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      val += gi[i][j]*ricci(i, j, loc);

  return val;
}

Sunday, 19 November 2017

Spacetime tells matter how to move, matter tells spacetime how to curve.

"Beyond the corridor of our spacetime there are infinite number of universes, each of them governed by their own set of laws."
Not general relativity, maybe next time if my computer has enough memory and time for that.
Just Maxwell's equations with two similar "forcefully" orbiting charges.

That's all Folks!
% ffmpeg -r 30 -f image2 -i %04d.png -ss 00:00:5 -vcodec libx264 -crf 2 output.mp4
clear all; close all; map = parula(256);
N    = 256;
phi  = zeros(N, N);
dphi = zeros(N, N);
ddphi= zeros(N, N);
rho  = zeros(N, N);
[x y] = meshgrid(linspace(-1, 1, N));
f = 0;
for t = 1:3*N
    A = 0.05;
    if t>N
        f = 0.01;
    end
    x1 = cos(2*pi*f*t)*A;
    y1 = sin(2*pi*f*t)*A;
    x2 = -cos(2*pi*f*t)*A;
    y2 = -sin(2*pi*f*t)*A;
    rho = exp(-2000*((x+x1).^2+(y+y1).^2)) + exp(-2000*((x+x2).^2+(y+y2).^2));

    ddphi = rho + del2(phi);
    dphi = dphi + ddphi;
    
     phi(1, :)   =   phi(2, :);
    dphi(1, :)   = -dphi(2, :);
     phi(end, :) =   phi(end-1, :);
    dphi(end, :) = -dphi(end-1, :);
     phi(:, 1)   =   phi(:, 2);
    dphi(:, 1)   = -dphi(:, 2);
     phi(:, end) =   phi(:, end-1);
    dphi(:, end) = -dphi(:, end-1);

    phi = phi + dphi;
    
    cla;
    c = phi/max(max(phi));
    surface(c, 'edgecolor', 'none');
    colormap(map);
    caxis([0 1]);
    axis([1 N 1 N]);
    axis square;
    axis off;
    shading interp;
    drawnow;
    imwrite(c*255, map, ['png/' num2str(t, '%04.f') '.png']);
end

Saturday, 11 November 2017

50 Things I'd like to do if I had all the time in the world

 1. Solve quantum gravity.
 2. Crack RSA (factoring) or prove it impossible (polynomial time).
 3. Build a working quantum computer (capable of factoring RSA) or prove it impossible (a separate question from the possibility of classical factoring).
 4. Learn to play piano (I can a bit, but learn it well).
 5. Build my own CPU from transistor level up.
 6. Build my own compiler.
 7. Build my own operating system.
 8. Build my own AI that rivals humans in general creative tasks.
 9. Build my own house.
10. Build my own car.
11. Acquire a few airplanes.
12. Build my own spaceship.
13. Build my own research satellite.
14. Visit the moon.
15. P vs. NP.
16. Climb Aconcagua.
17. Climb Mt. Denali.
18. Climb Mt. Everest.
19. Learn to fly a plane.
20. Visit Mt. Thor (Baffin island).
21. Sail to Ball's Pyramid from Sydney.
22. Visit Svalbard.
23. Directly observe the event horizon of the Milky Way central black hole.
24. Solve the mystery of baryon asymmetry.
25. Solve the mystery of dark energy and dark matter.
26. Crack FTL/time travel or prove them impossible.
27. Find extraterrestrial life.
28. Learn creative drawing.
29. Walk across the Canadian wilderness (~5000 km).
30. If there is a mystery related to human consciousness, solve it, if not, make it obvious for everone there is none.
31. See the day when religions basically are no more.
32. See the day when aging and all disease are cured.
33. See the Milky Way-Andromeda collision.
34. Visit a planet inhabited by estraterrestrial life.
35. Develop more meaningful social relationships.
36. Acquire a nice pacific island and build my own winter retreat there.
37. Acquire an underground missile base from the US and make it nice.
38. Acquire some weapons and explosives for fun and distribute them to my retreats all around the world.
39. Stockpile the retreats with supplies and build bunkers.
40. Visit Antarctica and Cape Horn area, Greenland, Siberia, India, China, New Zealand, Jordan, Africa...
41. Build an internationally successful company that manufactures useful hardware, solves a number of practical problems and employs a number of good people.
42. Make a game that tells an interesting story, teaches a few things and presents practical problems requiring decent problem solving skills.
43. Work for a while as national park ranger, chemist, philosopher, politician and as a neuroscientist.
44. Work for a while on genetic research and related biological issues.
45. Work for a while for NASA, nvidia, AI research, weapons research, some startups...
46. Live a few years in the wilderness and a monestary.
47. Try living in different cities around the world.
48. Get a cabin from Norway from some nice remote mountain slope next to a small lake with a small airstrip where I can fly from Finland with my own small plane.
49. Get a nice place from somewhere along the Californian cost, not too far from San Francisco, but in a not too densely populated location either.
50. Write a book that will stir up people and teach them some critical thinking.

* The word "prove" used in colloquial sense here.

Though I didn't like this particular turn of events in that particular episode as in my opinion while it is not possible to avoid all mistakes given the initial conditions, I never the less believe that losing is due to mistakes that under perfect knowledge and logic should not exist. Choose your mistakes wisely and you'll be fine.
Unexamined life is not worth living is often taken to imply that rather than not being worth living some lives are simply missing something important. I would, however, argue that all lives are missing something important and on the other than that no single thing in any life is something so important that one should consider it essential. Most life has intrinsic value that is independent of its fullness.

It is often also said that the saying refers to one being able to say why things are important, right or wrong rather than simply insisting they are. However, this too is a bit arbitrary as in the end there never seems to be any why. We are what we are and want what we want, because that's the way we have come to be. Some goals are better served by some actions, but no action is right or wrong without reference to motivations which in end are arbitrary. Without such motivations we would not be humans, but rather simple machines simply waiting to do what we are told.

It is of course worthwhile to know how to influence other people by justifying actions and inactions for the betterment of humanity in general.

Wednesday, 1 November 2017

Origin of the power of quantum computing

Predicting the behavior of many particle quantum system requires solving the Schrödinger equation for a wavefunction that has complex dimensionality proportinal to the number of particles. Such a wavefunctions contains a value for all possible combinations of particle positions, i.e. it describes a kind of "multiverse" where everything conceivable happens simultaneously.


The system evolves as a whole and thus in order to predict phenomena such as entanglement, it is necessary to solve the time evolution of the "multiverse" rather than time evolution of separate wavefunctions corresponding to single particles. Single separate wavefunctions cannot yield predictions corresponding to observations in every circumstance.
The power of quantum computing is based on the complexity of the time evolution of this kind of multiparticle wavefunction. The classical computational resources required to do full simulations scale extremely poorly, but it is possible to demonstrate the idea with a simple code (below).

Full simulation of the time evolution of three electrons starting at rest relative to each other. Variable containing the wavefunction is a 6D grid of 26 complex points (3D double complex 26x26) and requires about 2.3 GBs of memory. Adding one more particle would make it 8D and require about 1556 GBs of memory. Simulating 30 particles in this configuration would require more bytes of memory than there are atoms in the observable universe. Many problems can be approximated effectively without conducting such simulations, but there are a few cases like quantum computers which cannot. Qubit are typically in a simple static grid rather than free space such as electrons in this simulation so larger number of qubits can still be simulated. The point at which quantum computers are expected to become more powerful than classical supercomputers (so called quantum supremacy) is a the moment 56 qubits.
n = 0.05;
N = 26;
x10 = -0.4; y10 =  0.4; vx1 = 0; vy1 = 0;
x20 =  0.4; y20 =  0.4; vx2 = 0; vy2 = 0;
x30 =  0.0; y30 = -0.4; vx3 = 0; vy3 = 0;
frequency_spread = 40;
psi = zeros(N, N, N, N, N, N, 'single');
V = zeros(N, N, N, N, N, N, 'single');
S = linspace(-1.0, 1.0, N);
for y1 = 1:N
    for x1 = 1:N
        for y2 = 1:N
            for x2 = 1:N
                for y3 = 1:N
                    for x3 = 1:N
                        envelope = exp(-frequency_spread*((S(x1)+x10)^2 + (S(y1)+y10)^2 + ...
                                                          (S(x2)+x20)^2 + (S(y2)+y20)^2 + ...
                                                          (S(x3)+x30)^2 + (S(y3)+y30)^2));
                        velocity = exp(-j*(vx1*S(x1) + vy1*S(y1) + ...
                                           vx2*S(x2) + vy2*S(y2) + ...
                                           vx3*S(x3) + vy3*S(y3)));
                        psi(x1, y1, x2, y2, x3, y3) =  envelope.*velocity;
                        V(x1, y1, x2, y2, x3, y3) = 1/sqrt(n + (S(x1)-S(x2))^2 + ...
                                                               (S(y1)-S(y2))^2)^2 + ...
                                                    1/sqrt(n + (S(x1)-S(x3))^2 + ...
                                                               (S(y1)-S(y3))^2)^2 + ...
                                                    1/sqrt(n + (S(x2)-S(x3))^2 + ...
                                                               (S(y2)-S(y3))^2)^2;
                    end
                end
            end
        end
    end
end
for t = 1:500
    PSI1 = squeeze(mean(mean(mean(mean(psi, 1), 2), 3), 4));
    PSI2 = squeeze(mean(mean(mean(mean(psi, 1), 2), 5), 6));
    PSI3 = squeeze(mean(mean(mean(mean(psi, 3), 4), 5), 6));
    a = abs(PSI1).^2; a = a/max(max(a)); b = abs(PSI2).^2; b = b/max(max(b));
    c = abs(PSI3).^2; c = c/max(max(c)); val = a + b + c;
    imwrite(val/max(max(val)), ['png/' num2str(t, '%04.f') '.png']);
    
    z = zeros(N, N, N, N, N, N, 'single');
    z(2:N-1, :, :, :, :, :) = z(2:N-1, :, :, :, :, :) + diff(psi, 2, 1);
    z(:, 2:N-1, :, :, :, :) = z(:, 2:N-1, :, :, :, :) + diff(psi, 2, 2);
    z(:, :, 2:N-1, :, :, :) = z(:, :, 2:N-1, :, :, :) + diff(psi, 2, 3);
    z(:, :, :, 2:N-1, :, :) = z(:, :, :, 2:N-1, :, :) + diff(psi, 2, 4);
    z(:, :, :, :, 2:N-1, :) = z(:, :, :, :, 2:N-1, :) + diff(psi, 2, 5);
    z(:, :, :, :, :, 2:N-1) = z(:, :, :, :, :, 2:N-1) + diff(psi, 2, 6);
    % that's just del2(psi)

    psi = psi + 0.005*j*(z - V.*psi);
end

--

mergesort(int a[], int n) { // O(n log n)
  int *b, *c;
  int x=0, y=0, z;
  if(n>3) {
    b = (int*)malloc(4*n/2);
    c = (int*)malloc(4*n/2);
    for(z=0; z<n/2; z++) {
      b[z] = a[z];
      c[z] = a[z+n/2];
    }
    mergesort(b, n/2);
    mergesort(c, n/2);
    for(z=0; z<n; z++)
      if(y<n/2 && x<n/2)
        if(b[x]>c[y]) a[z] = c[y++];
        else a[z] = b[x++];
      else
        if(y>x) a[z] = b[x++];
        else a[z] = c[y++];
    free(b);
    free(c);
  } else
    if(a[0]>a[1]) { // swap a[0] and a[1] (xor swap)                            
      a[0] = a[0]^a[1];
      a[1] = a[1]^a[0];
      a[0] = a[0]^a[1];
    }
}

--
Standard bicubic upscaler (4x).
A = imread('smb.png');
A = double(A)/255;

s = size(A);
S = 2;
for y = 1:S*s(2)-1
    for x = 1:S*s(1)-1
        r(x, y) = A(1+floor(x/S), 1+floor(y/S), 1);
        g(x, y) = A(1+floor(x/S), 1+floor(y/S), 2);
        b(x, y) = A(1+floor(x/S), 1+floor(y/S), 3);
    end
end

r = r + del2(r);
g = g + del2(g);
b = b + del2(b);

...repeat...

[x y] = gradient(r); r = r + abs(x) + abs(y);
[x y] = gradient(g); g = g + abs(x) + abs(y);
[x y] = gradient(b); b = b + abs(x) + abs(y);
[x y] = gradient(r); r = r - abs(x) - abs(y);
[x y] = gradient(g); g = g - abs(x) - abs(y);
[x y] = gradient(b); b = b - abs(x) - abs(y);
Nabla squared filter applied in two stages to nearest neighbourhood doubler and sharpen in post prosessing. Pixel edges in 45 degree angle are totally eliminated and regular raster becomes solid color.

Tuesday, 24 October 2017

Schrödinger sims, raytracing and Julia sets

Particle colliding with an obstacle.
N = 256;
[x y] = meshgrid(linspace(-1, 1, N), linspace(-1, 1, N));
psi = exp(-10*((x+0.5).^2 + (y+0.5).^2)).*exp(-j*(-50*x -50*y));
psi = psi/norm(psi);
V = zeros(N, N);
V(N/2-20:N/2-10, N/2:N/2+20) = 1;
V(N/2+10:N/2+20, N/2-10:N/2+10) = 1;
S = conv2(diff([0 0 1 0 0], 2), eye(N)); % del2 matrix form
for n = 1:N % H = V - del2
    Ex{n} = expm(-j*(diag(V(n, :)) - S(:, 2:end-1)));
    Ey{n} = expm(-j*(diag(V(:, n)) - S(:, 2:end-1)));
end
for t = 0:150 % time evolution psi(t+1) = psi(t)*exp(-j*H)
    for n = 1:N
        psi(:, n) = Ex{n}*psi(:, n);
    end
    for n = 1:N
        psi(n, :) = psi(n, :)*Ey{n};
    end
    A = abs(psi).^2/4.5e-4;
    cla; imagesc(A); caxis([0 1]); truesize; drawnow;
    imwrite(imresize(A, 0.5), ['png/' num2str(t+1, '%04.f') '.png']);
    % ffmpeg -r 30 -f image2 -i %04d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p test.mp4
    % convert -resize 50% -layers Optimize -delay 3 *.png test.gif
end
Three particles scattering, simplified approximate solution only.
N = 64;
[x y] = meshgrid(linspace(-1, 1, N), linspace(-1, 1, N));
A = 1.4; B = -0.2; C = 1.2;
psi{1} = exp(-50*((x-0.5*A-B).^2 + (y-0.4*A).^2)).*exp(-j*( 20*x*C + 10*y*C));
psi{2} = exp(-50*((x+0.2*A-B).^2 + (y+0.3*A).^2)).*exp(-j*(-20*x*C - 10*y*C));
psi{3} = exp(-50*((x-0.4*A-B).^2 + (y+0.3*A).^2)).*exp(-j*( 10*x*C - 20*y*C));
v = fft2(0.3./(x.^2+y.^2+0.01));
for x = 1:3
    psi{x} = psi{x}/norm(psi{x});
end
in = 1;
for t = 0:15000
    if mod(t, 150)==0
        B = 0;
        for x = 1:3
            A = abs(psi{x}).^2;
            B = B + A/max(max(A));
        end
        figure(1); cla; imagesc(imresize(B, 4)); truesize; drawnow;
        imwrite(B, ['png/' num2str(in, '%04.f') '.png']);
        in = in + 1
    end
    for x = 1:3 % fast 2D convolution with fft
        V{x} = fftshift(ifft2(v.*fft2(conj(psi{x}).*psi{x})));
    end
    for x = 1:3 % separate wave functions
        W = zeros(N, N);
        for y = 1:3
            if y~=x
                W = W + V{y};
            end
        end
        psi{x} = psi{x} + 1e-3*j*(4*del2(psi{x}) - W.*psi{x});
    end   
end
Traces ever pixel of the screen to z-direction until ray within the object and
then chooses color based on surface normal and distance.
The edges to infinity are funny due to "cutting corners" in the way the normal is computed.
/* gcc test.c -lGL -lSDL */
#include <GL/gl.h>
#include <SDL/SDL.h>

char *d =
  "float t = gl_Color.x*1000.;"
  "vec3 scene(vec2 p) {"
  "  float z = 5.;"
  "  mat3 X = mat3(1., 0., 0., 0., cos(2.*t), sin(2.*t), 0., -sin(2.*t), cos(2.*t));"
  "  mat3 Y = mat3(cos(3.*t), 0., -sin(3.*t), 0., 1., 0., sin(3.*t), 0., cos(3.*t));"
  "  vec3 c = X*Y*vec3(1.5*cos(.5*5.*t), 1.5*sin(.5*7.*t), 7.);"
  "  while(z<10.) {"
  "    vec3 r = X*Y*(vec3(sqrt(z)*p.x, sqrt(z)*p.y, z));"
  "    if(abs(r.x-c.x)<.8 && abs(r.y-c.y)<.8 && abs(r.z-c.z)<.8) return r;"
  "    if(distance(r, c)<1.) return r;"
  "    z += .05;"
  "  }"
  "  return vec3(10.);"
  "}"
  "void main() {"
  "  gl_FragColor = vec4(0.);"
  "  vec2 p = vec2(gl_FragCoord.x/256.-1., gl_FragCoord.y/256.-1.);"
  "  vec3 q = scene(p);"
  "  if(length(q)<10.) {"
  "    vec3 dx = scene(vec2(p.x + .05, p.y)) - scene(vec2(p.x - .05, p.y));"
  "    vec3 dy = scene(vec2(p.x, p.y + .05)) - scene(vec2(p.x, p.y - .05));"
  "    float a = 1.5*dot(normalize(cross(dx, dy)), normalize(q))/sqrt(length(q));"
  "    gl_FragColor = vec4(a);"
  "  }"
  "}";

main(c) {
  SDL_Event e;
  SDL_SetVideoMode(512, 512, 0, SDL_OPENGL);
  glUseProgram(glCreateShaderProgramv(GL_FRAGMENT_SHADER, 1, &d));

/*
  c = 0;
  GLubyte **pixels = malloc(3*512*512);
  int cur;
  char filename[80];
*/

  while(e.type!=SDL_QUIT) {
    glColor3us(c++, 0, 0);
    glRecti(-1, -1, 1, 1);
    SDL_GL_SwapBuffers();
    SDL_PollEvent(&e);
/*
    sprintf(filename, "%04d.ppm", c);
    FILE *f = fopen(filename, "w");
    fprintf(f, "P3\n%d %d\n%d\n", 512, 512, 255);
    *pixels = realloc(*pixels, 3*512*512);
    glReadPixels(0, 0, 512, 512, GL_RGB, GL_UNSIGNED_BYTE, *pixels);
    for(int i=0; i<512; i++) {
        for(int j=0; j<512; j++) {
            cur = 3*((512-i-1)*512+j);
            fprintf(f, "%3d %3d %3d ", (*pixels)[cur], (*pixels)[cur+1], (*pixels)[cur+2]);
        }
        fprintf(f, "\n");

    }
*/
  }
}
Zoom into the Julia set (one of them).
N = 256; B = 0.75;
X = linspace(-1, 1, N);
Y = linspace(-1, 1, N);
while true
    B = B*sqrt(sqrt(0.75))
    y0 = 0.636;
    x0 = 0.107;
    X = linspace(x0-B, x0+B, N);
    Y = linspace(y0-B, y0+B, N);
    for y = 1:N
        for x = 1:N
            % z = 0.0;                  % Mandelbrot
            % z0 = X(x) + Y(y)*j;       % Mandelbrot
            z = X(x) + Y(y)*j;          % Julia
            z0 = 0.5 + 0.5*j;           % Julia
            iteration = 0;
            while abs(z) < 2 && iteration < 128
                z = z^2 + z0;
                iteration = iteration + 1;
            end
            s(x, y) = log(iteration);
        end
    end
    imagesc(s); drawnow;
end
#include <stdio.h>
#include <math.h>
main() {
  char filename[80];
  int N = 128, x, y, c=1;
  float f[128*128];  float p[128*128];
  float u[128*128];  float v[128*128];
  float ux[128*128]; float vx[128*128];
  float uy[128*128]; float vy[128*128];
  float u2[128*128]; float v2[128*128];
  float us[128*128]; float vs[128*128];
  float px[128*128]; float py[128*128];
  for(x=0; x<128*128; x++) {
    f[x] = 0.0;    p[x] = 0.0;
    u[x] = 0.0;    v[x] = 0.0;
    ux[x] = 0.0;   vx[x] = 0.0;
    uy[x] = 0.0;   vy[x] = 0.0;
    u2[x] = 0.0;   v2[x] = 0.0;
    us[x] = 0.0;   vs[x] = 0.0;
    px[x] = 0.0;   py[x] = 0.0;  }
  for(int t=0; t<4500; t++) {
    for(y=-5; y<5; y++)
      for(x=-5; x<5; x++) {
        u[(int)(0.75*N+x)+(N/2+y)*N] = 0.15;
        v[(int)(0.75*N+x)+(N/2+y)*N] = 0.0;
      }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        ux[x+y*N] = 0.50*(u[(x+1)+y*N] -                u[(x-1)+y*N]);
        vx[x+y*N] = 0.50*(v[(x+1)+y*N] -                v[(x-1)+y*N]);
        uy[x+y*N] = 0.50*(u[x+(y+1)*N] -                u[x+(y-1)*N]);
        vy[x+y*N] = 0.50*(v[x+(y+1)*N] -                v[x+(y-1)*N]);
        u2[x+y*N] = 0.25*(u[(x+1)+y*N] - 2.0*u[x+y*N] + u[(x-1)+y*N]);
        v2[x+y*N] = 0.25*(v[(x+1)+y*N] - 2.0*v[x+y*N] + v[(x-1)+y*N]);
        u2[x+y*N]+= 0.25*(u[x+(y+1)*N] - 2.0*u[x+y*N] + u[x+(y-1)*N]);
        v2[x+y*N]+= 0.25*(v[x+(y+1)*N] - 2.0*v[x+y*N] + v[x+(y-1)*N]);
      }
    for(x=0; x<N*N; x++) {
      us[x] = u[x] + u[x]*ux[x] + v[x]*uy[x] + u2[x]/20.0 - px[x]/20.0 + 0.00000;
      vs[x] = v[x] + u[x]*vx[x] + v[x]*vy[x] + v2[x]/20.0 - py[x]/20.0 - 0.00005;
    }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        f[x+y*N] = 0.5*(us[(x+1)+y*N] - us[(x-1)+y*N] + vs[x+(y+1)*N] - vs[x+(y-1)*N]);
        p[x+y*N] = 0.25*(p[(x+1)+y*N] + p[x+(y+1)*N] + p[(x-1)+y*N] + p[x+(y-1)*N] - f[x+y*N]);
    }
    for(y=1; y<N-1; y++)
      for(x=1; x<N-1; x++) {
        px[x+y*N] = 0.5*(p[(x+1)+y*N] - p[(x-1)+y*N]);
        py[x+y*N] = 0.5*(p[x+(y+1)*N] - p[x+(y-1)*N]);
         u[x+y*N] = us[x+y*N] - px[x+y*N];
         v[x+y*N] = vs[x+y*N] - py[x+y*N];
      }
    if(t%40==0) {
      printf("%d\n", t);
      sprintf(filename, "ppm/%04d.pgm", c);
      FILE *f = fopen(filename, "w");
      fprintf(f, "P2\n%d %d\n%d\n", 128, 128, 255);
      c++;
      for(int j=0; j<128; j++) {
          for(int i=0; i<128; i++) {
             float val = 1500.0*sqrt( u[i+j*N]*u[i+j*N] + v[i+j*N]*v[i+j*N]);
             if(val>255.0) val=255.0;
             fprintf(f, "%3d ", (unsigned char)val);
          }
          fprintf(f, "\n");
      }
      close(f);
    }
  }
}