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.