text:
MS Paint is underrated
MS Paint was always a great time waster, but it was shunned by most graphic designers. Perhaps people lost interest because of the jarring color palette, or because of limited undo levels. Regardless, it's still possible to generate beautiful images with just the standard brush and default color palette.
Challenge
Using only the default brush (a 4x4 square with no corners) and the default color palette (the 28 colors below), attempt to replicate a source image using a technique based on stochastic hill climbing.

Algorithm
Every answer must follow the same basic algorithm (stochastic hillclimb). Details can be tweaked within each step. A movement is considered to be a stroke of the brush (ie. clicking in paint).

Guess next movement(s). Make a guess (of coordinates and colors) for the next movement(s) however you'd like. However the guess must not reference the source image.
Apply the guess. Apply the brush to the painting to make the movement(s).
Measure the benefit of movement(s). By referencing the source image, determine whether the movement(s) benefited the painting (ie. the painting more closely resembles the source image). If it is beneficial, keep the movement(s), otherwise discard the movement(s).
Repeat until convergence. Go to Step 1 and try the next guess until the algorithm has converged sufficiently. The painting should heavily resemble the source image at this point.

If your program does not match these four steps, it is probably not a stochastic hillclimb. I have tagged this as a popularity contest because the goal is to produce interesting painting algorithms based on the limited color palette and brush.
Contraints

The algorithm should be stochastic in some way.
The next guess should not be influenced by the source image. You are guessing each new movement, and then checking to see if it helped or not. For example, you are not allowed to determine where to place the brush based on the source image colors (that is similar to dithering the source image, which is not the goal).
You are allowed to influence the placement by tweaking the steps of the algorithm however you'd like. For example, you could start your guesses at the edges and move inwards, drag the brush to create lines for each guess, or decide to paint dark colors first. You are allowed to reference previous iteration images (but not the source image) in order to calculate the next desired movement. These can be restrictive as you'd like (ie. only make guesses within the top-left quadrant for the current iteration).
The measure of "difference" between the source image and the current iteration can be measured however you'd like, as long as it does not calculate other potential movements to determine whether this movement is considered to be the "best". It shouldn't know whether the current movement is the "best", only whether it fits within the tolerance of the acceptance criteria. For example, it can be as simple as abs(src.R - current.R) + abs(src.G - current.G) + abs(src.B - current.B) for each affected pixel, or any of the well-known color difference techniques.

Palette
You can download the palette as a 28x1 image or create it directly in the code.
Brush
The brush is a 4x4 square without corners. This is a scaled version of it:

(Your code must use the 4x4 version)
Example
Input:

Output:

You can see how the basic algorithm progresses in a short video I made (each frame is 500 iterations): The Starry Night. The initial stages are interesting to watch:

@vihan1086: The source image does not have transparency. The current image can depend on the previous iteration (like my example, in which new points are added on top of the previous) if that's what you mean.

I don't see what the stochastic hill climbing adds at all... considering that you can make guesses however you want and discard them if they're not good, that's effectively the same as just going through a bunch of guesses while checking and picking the best one.

@Sp3000: The point is that you don't know the "best" movement until it a potential movement has been made, at which point you can choose to accept it if it fits inside your own acceptance criteria (ie. it's "close enough"). The acceptance criteria shouldn't have knowledge of all possible moves (I may need to clarify this further). Essentially you shouldn't be able to determine the "best" movement ahead of time, instead you should be incrementally improving the image.

With this brush shape, is it impossible to paint the corner pixels (from your example I get this impression), or can the brush be placed partially outside the image boundaries?

Sp3000's concern perhaps includes the following deterministic algorithm: for each pixel in turn, try every colour. No stochasticity and very much like dithering, yet it seems to fit the rules.

@Oliphaunt: You can place it partially outside the image boundary to include the corners if desired. Technically it's possible to do this in MS Paint by dragging the brush.

@trichoplax: I don't see an obvious reason to disallow zero benefit movements if changes are being made, but I do want to ensure that it's not abused (ie. to determine the "best" movement at that step).

@Oliphaunt: Agreed. The acceptance criteria in that case would have to know the "best" color at that moment, which is not allowed. I've tried to clarify this in the constraints, but if you have any suggestions I can improve it.

A few quick questions: how show input be taken? How show output be given? What format is the image? What do you mean "create directly in the code" (as, can we store the colors in "variables")?

@SirPython: Input is an argument to the program or simply contained in a variable. Format of the image can be any common format (jpeg, png, etc.). You can create the palette directly in the code by using the integers of each channel or similar.

@grovesNL Awesome! Thanks for answering my questions. I just have one more: does "common format" include NetPBM format?

@SirPython: Sure, I guess you would need to use PPM for color

Here's the palette in a usable format: 000000
808080
800000
808000
008000
008080
000080
800080
808040
004040
0080FF
004080
8000FF
804000
FFFFFF
C0C0C0
FF0000
FFFF00
00FF00
00FFFF
0000FF
FF00FF
FFFF80
00FF80
80FFFF
8080FF
FF0080
FF8040

JavaScript
This solution uses the HTML5 canvas element to extract the image data, but without the need to use HTML, that means it can be run in your console. It access the color palette image as an array; I stored all the colors from the palette image in an array). It outputs to the console (after it finishes) and also stores the result in a variable.
The most updated version of the code is in the fiddle. The fiddle also uses a better algorithm to reduce noise in the pictures. The improvement in the algorithm is mostly fixing a function (max to min) which caused the inverse color to be chosen.
Code in the shape of the MS Paint Icon! (formatted code in fiddle or Stack Snippet)
eval(`                                                   function
                                                        Paint(t){fun
                                                         ction n(t){va
                                                         r n=t.toString(
                                                         16);return 1==n.
                                                         length?"0"+n:n}fu
                                                         nction e(t){return
                                                         "#"+n(t[0])+n(t[1]
                                                          )+n(t[2])}var a=ne
                                                          w Image,i=document.
                                                          createElement("canv
                                                          as"),h=null,o=docum
                                                          ent.createElement(
                                    "canvas"),r=          o.getContext("2d
                               ")     ,l=[],u=this,c      =[[0,0,0],[255
                            ,2       55,255],[ 192,192,   192],[128,12
                          8     ,128],[126,3,8],[252,13,27] ,[255,25
                       3,    56],[128,127,23],[15,127,18],[ 41,253
                      ,    46],[45,255,254],[17,128,127],[2 ,12,1
                    2    6],[ 11,36,2 51],[252,40,252],[12 7,15,1
                  2    6],[  128,127 ,68],[255,253,136],[4 2,253,
                 1   33],   [4,64,64],[23 ,131,251],[133,255,254],
               [    129   ,132,252],[6,6 6,126],[127,37,2 51],[127,
              6   4,1    3],[253,128,73],[252,22,129]];a.crossOrigin
             =   "",   a.src=t,this.done=this.done||function(){},a.o
            n   load=function(){function t(t){var n=0,e=0,a=0;return
           t  .forEach(function(t){n+=t[0],e+=t[1],a+=t[2]}),[n/t.leng
          t  h,e /t.length,a/t.length]}function n(t){for(var n=[],e=0;e
         <  t.l  ength;e+=1)n.push(t[e]);return n}function g(t,n){retur
        n  (Ma  th.abs(t[0]-n[0])/255+Math.abs(t[1]-n[1])/255+Math.abs(t
       [  2]- n[2])/255)/3}function f(t,n){for(var e=Math.floor(Math.ran
          do m()*n.length),a=n[e],i=(g(t,a),1-.8),h=56,o=[];o.length<=h&
         &g (t,a)>i;)if(o.push(a),a=n[Math.floor(Math.random()*n.length)]
     ,  o.length==h){var r=o.map(function(n){return g(t,n)});a=o[r.indexO
       f(Math.max.apply(Math,r))],o.push(a)}return a}function s(t,n){for(
    v  ar e=[];t.length>0;)e.push(t.splice(0,n).slice(0,-1));return e}i.w
   i  dth=a.width,i.height=2*a.height,h=i.getContext("2d"),h.drawImage(a,0
   ,0,a.width,a.height);for(var d=(function(t){reduce=t.map(function(t){re
  turn(t[ 0]+t[1]+t[2])/3})}(c),0),m=0,v=i.width*i.height/4,p=0;v>p;p+=1)d
  >2*Mat h.ceil(a.width/2)&&(d=0,m+=1),l.push(f(t(s(n(h.getImageData(2*d,2
  *m,4,4).data),4)),c)),d+=1;o.width=i.width,o.height=i.height;for(var d=0
 ,m=0,v=i.width*i.height/4,p=0;v>p;p+=1)d>2*Math.ceil(a.width/2)&&(d=0,m+=
 1),console.log("Filling point ("+d+", "+m+") : "+e(l[p])),r.fillStyle=e(l
 [p]),r.fillRect(2*d+1,2*m,2,1)  ,r.fillRect(2*d,2*m+1,4,2),r.fillRect(2*d
+1,2*m+3,2,1),d+=1;u.result=o      .toDataURL("image/png"),u.resultCanvas
=o,u.imageCanvas=i,u.image=a       ,u.done(),console.log(u.result)},a.one
rror=function(t){console.log       ("The image failed to load. "+t)}}/*..
............................       ......................................
. ..........................       .....................................
............................      ......................................
.............................    .......................................
.......................................................................
.......................................................................
..................  ..................................................
................     .................................................
..............       ................................................
.............       ................................................
...........        .................................................
 .........         ................................................
 .......          ................................................
  ....           ................................................
                ................................................
                ...............................................
               ...............................................
              ..............................................
              .............................................
             ............................................
             ..........................................
              .......................................
              .....................................
               .................................
                .............................
                  ......................
                                   .....
                                  .....
                                  .....
                                  ....
                                   */`
.replace(/\n/g,''))

Usage:
Paint('DATA URI');

Fiddle.
The fiddle uses crossorigin.me so you don't need to worry about cross-origin-resource-sharing.
I've also updated the fiddle so you can adjust some values to produce the best-looking painting. Some pictures' colors might be off, to avoid this, adjust the accept_rate to adjust the algorithm. A lower number means better gradients, a higher number will result in sharper colors.

Here's the fiddle as a Stack-Snippet (NOT updated, in case the fiddle doesn't work):

/* Options */

var accept_rate = 82,  // 0 (low) - 100 (high)
    attempts    = 16, // Attemps before giving up
    edge_multi  = 2;   // Contrast, 2-4

function Paint(image_url) {
    var image = new Image(), canvas = document.createElement('canvas'), context = null, result = document.createElement('canvas'), resultContext = result.getContext('2d'), final_colors = [], self = this, color_choices = [
        [0,0,0],
        [255,255,255],
        [192,192,192],
        [128,128,128],
        [126,3,8],
        [252,13,27],
        [255,253,56],
        [128,127,23],
        [15,127,18],
        [41,253,46],
        [45,255,254],
        [17,128,127],
        [2,12,126],
        [11,36,251],
        [252,40,252],
        [127,15,126],
        [128,127,68],
        [255,253,136],
        [42,253,133],
        [4,64,64],
        [23,131,251],
        [133,255,254],
        [129,132,252],
        [6,66,126],
        [127,37,251],
        [127,64,13],
        [253,128,73],
        [252,22,129]
      ];
  image.crossOrigin = "";
  image.src = image_url;

  this.done = this.done || function () {};

  function hex(c) {
    var res = c.toString(16);
    return res.length == 1 ? "0" + res : res;
  }

  function colorHex(r) {
    return '#' + hex(r[0]) + hex(r[1]) + hex(r[2]);
  }

  image.onload = function () {
    canvas.width = image.width; canvas.height = image.height * 2;
    context = canvas.getContext('2d');
    context.drawImage(image, 0, 0, image.width, image.height);

    function averageColors(colors_ar) {
      var av_r = 0,
          av_g = 0,
          av_b = 0;

      colors_ar.forEach(function (color) {
        av_r += color[0];
        av_g += color[1];
        av_b += color[2];
      });

      return [av_r / colors_ar.length,
              av_g / colors_ar.length,
              av_b / colors_ar.length];
    }

    function arrayFrom(ar) {
      var newar = [];
      for (var i = 0; i < ar.length; i += 1) {
        newar.push(ar[i]);
      }
      return newar;
    }

    function colorDif(c1,c2) {
      // Get's distance between two colors 0.0 - 1.0
      return (Math.abs(c1[0] - c2[0]) / 255 +
              Math.abs(c1[1] - c2[1]) / 255 +
              Math.abs(c1[2] - c2[2]) / 255) / 3;
    }

    var furthest = (function (cc) {
      // Determines furthest color

      // Reduces RGB into a "single value"
      reduce = cc.map(function(color) {
        return ( color[0] + color [1] + color [2] ) / 3;
      });

    }(color_choices));

    function intDif(i1,i2,t) {
        return Math.abs(i1 - i2) / t
    }

    function arrayIs(ar, int,d) {
        return intDif(ar[0],int,255) <= d &&
               intDif(ar[1],int,255) <= d &&
               intDif(ar[2],int,255) <= d
    }

    function colorLoop(c1,c2) {
        var edgeCap =  edge_multi * ((accept_rate / 100) / 50), values = c2.map(function (i) {
            return colorDif(c1,i);
        });

      return arrayIs(c1,255,edgeCap)?[255,255,255]:
             arrayIs(c1,0,edgeCap)  ?[0,0,0]:
             c2[values.indexOf(Math.min.apply(Math, values))];
    }

    function colorFilter(c1, c2) {
      // Does the color stuff
      var rand  = Math.floor( Math.random() * c2.length ), // Random number
          color = c2[rand], // Random color
          randdif = colorDif(c1, color),
          threshhold = 1 - accept_rate / 100, // If the color passes a threshhold
          maxTries   = attempts, // To avoid infinite looping, 56 is the maximum tries to reach the threshold
          tries      = [];

      // Repeat until max iterations have been reached or color is close enough

      while ( tries.length <= maxTries && colorDif( c1, color ) > threshhold ) {
        tries.push(color);
        color = c2[Math.floor(Math.random() * c2.length)]; // Tries again

        if (tries.length == maxTries) {
          // Used to hold color and location
          var refLayer = tries.map(function(guess) {
            return colorDif(c1, guess);
          });

          color = tries[refLayer.indexOf(Math.min.apply(Math, refLayer))];
            tries.push(color);
        }

      }

      var edgeCap =  edge_multi * ((accept_rate / 100) / 50), loop = colorLoop(c1, c2);

      return arrayIs(c1,255,edgeCap)?[255,255,255]:
             arrayIs(c1,0,edgeCap)  ?[0,0,0]:
             colorDif(c1,color)<accept_rate?color:
             loop;
    }

    function chunk(ar, len) {
      var arrays = [];

      while (ar.length > 0)
        arrays.push(ar.splice(0, len).slice(0, -1));

      return arrays;
    }

    var x = 0, y = 0, total = (canvas.width * canvas.height) / 4;

    for (var i = 0; i < total; i += 1) {
      if (x > (Math.ceil(image.width / 2) * 2)) {
        x = 0;
        y += 1;
      }

      final_colors.push( colorFilter( averageColors( chunk( arrayFrom(context.getImageData(x * 2, y * 2, 4, 4).data), 4 ) ), color_choices) );

      x += 1;
    }

    // Paint Image
    result.width = canvas.width;
    result.height = canvas.height;
    var x = 0, y = 0, total = (canvas.width * canvas.height) / 4;
    for (var i = 0; i < total; i += 1) {
      if (x > (Math.ceil(image.width / 2) * 2)) {
        x = 0;
        y += 1;
      }

        console.log("Filling point (" + x + ", " + y + ") : " + colorHex(final_colors[i]));

      resultContext.fillStyle = colorHex(final_colors[i]);
      resultContext.fillRect(x*2 + 1, y * 2, 2 , 1); // Top
      resultContext.fillRect(x * 2, y * 2 + 1, 4, 2); // Middle
      resultContext.fillRect(x * 2 + 1, y * 2 + 3, 2, 1); // Bottom

      x += 1;
    }

    self.result = result.toDataURL("image/png");
    self.resultCanvas = result;
    self.imageCanvas = canvas;
    self.image = image;
    self.done();

    console.log(self.result);

  };

  image.onerror = function(error) {
    console.log("The image failed to load. " + error);
  }

}

// Demo

document.getElementById('go').onclick = function () {
    var url = document.getElementById('image').value;
    if (!url.indexOf('data:') == 0) {
        url = 'http://crossorigin.me/' + url;
    }
    var example = new Paint(url);

    example.done = function () {
        document.getElementById('result').src = example.result;
        document.getElementById('result').width = example.resultCanvas.width;
        document.getElementById('result').height = example.resultCanvas.height;
        window.paint = example;
    };
};
<!--This might take a while-->
Enter the image data URI or a URL, I've used crossorigin.me so it can  perform CORS requests to the image. If you're inputting a URL, be sure to include the http(s)
<input id="image" placeholder="Image URI or URL"><button id="go">Go</button>
    <hr/>
    You can get the image URI from a website like <a href="http://jpillora.com/base64-encoder/">this one</a>
    <hr/>

    Result:
    <img id="result">
        <span id="error"></span><hr/>
        Check your console for any errors. After a second, you should see the colors that are being generated / printed getting outputted to the console.

To commemorate New Horizon's flyby of Pluto, I've inputted an image of Pluto:

For the following I've set it to make them resemble the original as close as possible:
I ran this with OS X Yosemite's default wallpaper. After leaving it run for a bit, the results are absolutely stunning. The original file was huge (26 MB) so I resized and compressed it:

The starry night (I've used a higher resolution image for better results)

A picture I found on google:

JavaScript + HTML
Random:
Random Point
Random Aligned:
Subdivides the canvas into 4x4 squares, and chooses a point randomly inside one of the squares. Offsets will move the grid, so you can fill in the little gaps.
Loop:
Creates a grid and Loops through all the points.
Offsets moves the grid.
Spacing determine the size of each cell. (They will start to overlap)
Color difference:

RGB
HSL
HSV

var draw = document.getElementById("canvas").getContext("2d");
var data = document.getElementById("data").getContext("2d");
colors = [
    [0, 0, 0],
    [255, 255, 255],
    [192, 192, 192],
    [128, 128, 128],
    [126, 3, 8],
    [252, 13, 27],
    [255, 253, 56],
    [128, 127, 23],
    [15, 127, 18],
    [41, 253, 46],
    [45, 255, 254],
    [17, 128, 127],
    [2, 12, 126],
    [11, 36, 251],
    [252, 40, 252],
    [127, 15, 126],
    [128, 127, 68],
    [255, 253, 136],
    [42, 253, 133],
    [4, 64, 64],
    [23, 131, 251],
    [133, 255, 254],
    [129, 132, 252],
    [6, 66, 126],
    [127, 37, 251],
    [127, 64, 13],
    [253, 128, 73],
    [252, 22, 129]
];
iteration = 0;
fails = 0;
success = 0;
x = 0;
y = 0;
//Init when the Go! button is pressed
document.getElementById("file").onchange = function (event) {
    document.getElementById("img").src = URL.createObjectURL(event.target.files[0]);
    filename = document.getElementById("file").value;
    /*if (getCookie("orginal") == filename) {
        console.log("Loading from Cookie");
        reload = true;
        document.getElementById("reload").src = getCookie("picture");
    }*/
};

/*function getCookie(cname) {
    var name = cname + "=";
    var ca = document.cookie.split(';');
    for (var i = 0; i < ca.length; i++) {
        var c = ca[i];
        while (c.charAt(0) == ' ') c = c.substring(1);
        if (c.indexOf(name) == 0) return c.substring(name.length, c.length);
    }
    return "";
}*/

//Run when the image has been loaded into memory
document.getElementById("img").onload = function () {
    document.getElementById("file").disable = "true";
    document.getElementById("canvas").hidden = "";
    document.getElementById("canvas").height = document.getElementById("img").height;
    document.getElementById("data").height = document.getElementById("img").height;
    document.getElementById("canvas").width = document.getElementById("img").width;
    document.getElementById("data").width = document.getElementById("img").width;

    var imgData = draw.createImageData(document.getElementById("img").width, document.getElementById("img").height);
    for (var i = 0; i < imgData.data.length; i += 4) {
        imgData.data[i + 0] = 0;
        imgData.data[i + 1] = 0;
        imgData.data[i + 2] = 0;
        imgData.data[i + 3] = 255;
    }
    draw.putImageData(imgData, 0, 0);
    data.putImageData(imgData, 0, 0);
    if (reload == true) {
        draw.drawImage(document.getElementById("reload"), 0, 0);
    }
    data.drawImage(document.getElementById("img"), 0, 0);
    setInterval(function () {
        for (var u = 0; u < document.getElementById("selectColor").value; u++) {
            doThing();
        }
    }, 0);
};

//The core function. Every time this function is called, is checks/adds a dot.
function doThing() {
    getCoords();
    paintBucket();
    console.count("Iteration");
    if (compare(x, y)) {
        draw.putImageData(imgData, x, y);
    }
}

function getCoords() {
    switch (document.getElementById("selectCord").value) {
        case "1":
            x = Math.floor(Math.random() * (document.getElementById("img").width + 4));
            y = Math.floor(Math.random() * (document.getElementById("img").height + 4));
            break;
        case "2":
            x = Math.floor(Math.random() * ((document.getElementById("img").width + 4) / 4)) * 4;
            console.log(x);
            x += parseInt(document.getElementById("allignX").value);
            console.log(x);
            y = Math.floor(Math.random() * ((document.getElementById("img").height + 4) / 4)) * 4;
            y += parseInt(document.getElementById("allignY").value);
            break;
        case "3":
            x += parseInt(document.getElementById("loopX").value);
            if (x > document.getElementById("img").width + 5) {
                x = parseInt(document.getElementById("allignX").value);
                y += parseInt(document.getElementById("loopY").value);
            }
            if (y > document.getElementById("img").height + 5) {
                y = parseInt(document.getElementById("allignY").value);
            }
    }
}

function compare(arg1, arg2) {
    var arg3 = arg1 + 4;
    var arg4 = arg2 + 4;
    imgData2 = data.getImageData(arg1, arg2, 4, 4);
    imgData3 = draw.getImageData(arg1, arg2, 4, 4);
    N = 0;
    O = 0;
    i = 4;
    addCompare();
    addCompare();
    i += 4;
    for (l = 0; l < 8; l++) {
        addCompare();
    }
    i += 4;
    addCompare();
    addCompare();
    i += 4;
    //console.log("New Score: " + N + " Old Score: " + O);
    iteration++;
    /*if(iteration>=1000){
        document.cookie="orginal="+filename;
        document.cookie="picture length="+document.getElementById("canvas").toDataURL().length;
        document.cookie="picture="+document.getElementById("canvas").toDataURL();

    }*/
    if (N < O) {
        return true;
    } else {
        return false;
    }
}

function addCompare() {
    if (document.getElementById("colorDif").value == "HSL") {
        HSLCompare();
        i += 4;
        return;
    }
    if (document.getElementById("colorDif").value == "HSV") {
        HSVCompare();
        i += 4;
        return;
    }
    N += Math.abs(imgData.data[i] - imgData2.data[i]);
    N += Math.abs(imgData.data[i + 1] - imgData2.data[i + 1]);
    N += Math.abs(imgData.data[i + 2] - imgData2.data[i + 2]);
    O += Math.abs(imgData3.data[i] - imgData2.data[i]);
    O += Math.abs(imgData3.data[i + 1] - imgData2.data[i + 1]);
    O += Math.abs(imgData3.data[i + 2] - imgData2.data[i + 2]);
    i += 4;
}

function HSVCompare() {
    var NewHue = rgbToHsv(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[0];
    var PicHue = rgbToHsv(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[0];
    var OldHue = rgbToHsv(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[0];

    var NScore = [Math.abs(NewHue - PicHue), ((NewHue < PicHue) ? NewHue + (1 - PicHue) : PicHue + (1 - NewHue))];
    var OScore = [Math.abs(OldHue - PicHue), ((OldHue < PicHue) ? OldHue + (1 - PicHue) : PicHue + (1 - OldHue))];

    NScore = Math.min(NScore[0], NScore[1]);
    OScore = Math.min(OScore[0], OScore[1]);

    NewHue = rgbToHsv(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[1];
    PicHue = rgbToHsv(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[1];
    OldHue = rgbToHsv(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[1];

    NScore += Math.abs(NewHue-PicHue);
    OScore += Math.abs(OldHue-PicHue);

    NewHue = rgbToHsv(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[2];
    PicHue = rgbToHsv(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[2];
    OldHue = rgbToHsv(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[2];

    N += Math.abs(NewHue-PicHue) + NScore;
    O += Math.abs(OldHue-PicHue) + OScore;
}

function rgbToHsv(r, g, b){
    r = r/255, g = g/255, b = b/255;
    var max = Math.max(r, g, b), min = Math.min(r, g, b);
    var h, s, v = max;

    var d = max - min;
    s = max == 0 ? 0 : d / max;

    if(max == min){
        h = 0; // achromatic
    }else{
        switch(max){
            case r: h = (g - b) / d + (g < b ? 6 : 0); break;
            case g: h = (b - r) / d + 2; break;
            case b: h = (r - g) / d + 4; break;
        }
        h /= 6;
    }

    return [h, s, v];
}

function HSLCompare() {
    var result = 0;
    rgb = false;

    var NewHue = rgbToHue(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[0];
    var PicHue = rgbToHue(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[0];
    var OldHue = rgbToHue(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[0];
    if (rgb == true) {
        N += Math.abs(imgData.data[i] - imgData2.data[i]);
        N += Math.abs(imgData.data[i + 1] - imgData2.data[i + 1]);
        N += Math.abs(imgData.data[i + 2] - imgData2.data[i + 2]);
        O += Math.abs(imgData3.data[i] - imgData2.data[i]);
        O += Math.abs(imgData3.data[i + 1] - imgData2.data[i + 1]);
        O += Math.abs(imgData3.data[i + 2] - imgData2.data[i + 2]);
        return;
    }
    var NScore = [Math.abs(NewHue - PicHue), ((NewHue < PicHue) ? NewHue + (1 - PicHue) : PicHue + (1 - NewHue))];
    var OScore = [Math.abs(OldHue - PicHue), ((OldHue < PicHue) ? OldHue + (1 - PicHue) : PicHue + (1 - OldHue))];

    NScore = Math.min(NScore[0], NScore[1]);
    OScore = Math.min(OScore[0], OScore[1]);

    NewHue = rgbToHue(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[1];
    PicHue = rgbToHue(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[1];
    OldHue = rgbToHue(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[1];

    NScore += Math.abs(NewHue-PicHue);
    OScore += Math.abs(OldHue-PicHue);

    NewHue = rgbToHue(imgData.data[i], imgData.data[i + 1], imgData.data[i + 2])[2];
    PicHue = rgbToHue(imgData2.data[i], imgData2.data[i + 1], imgData2.data[i + 2])[2];
    OldHue = rgbToHue(imgData3.data[i], imgData3.data[i + 1], imgData3.data[i + 2])[2];

    N += Math.abs(NewHue-PicHue) + NScore;
    O += Math.abs(OldHue-PicHue) + OScore;
}

function rgbToHue(r, g, b) {
    if (Math.max(r, g, b) - Math.min(r, g, b) < 50) {
        rgb = true
    }
    r /= 255, g /= 255, b /= 255;
    var max = Math.max(r, g, b),
        min = Math.min(r, g, b);
    var h, s, l = (max + min) / 2;

    if (max == min) {
        h = s = 0; // achromatic
    } else {
        var d = max - min;
        s = l > 0.5 ? d / (2 - max - min) : d / (max + min);
        switch (max) {
            case r:
                h = (g - b) / d + (g < b ? 6 : 0);
                break;
            case g:
                h = (b - r) / d + 2;
                break;
            case b:
                h = (r - g) / d + 4;
                break;
        }
        h /= 6;
    }
    return [h,s,l];
}

//Create a 4x4 ImageData object, random color selected from the colors var, transparent corners.
function paintBucket() {
    color = Math.floor(Math.random() * 28);
    imgData = draw.createImageData(4, 4);
    imgData2 = draw.getImageData(x, y, 4, 4);
    i = 0;
    createCorn();
    createColor();
    createColor();
    createCorn();
    for (l = 0; l < 8; l++) {
        createColor();
    }
    createCorn();
    createColor();
    createColor();
    createCorn();
}

function createCorn() {
    imgData.data[i] = imgData2.data[i];
    imgData.data[i + 1] = imgData2.data[i + 1];
    imgData.data[i + 2] = imgData2.data[i + 2];
    imgData.data[i + 3] = 255;
    i += 4;
}

function createColor() {
    imgData.data[i] = colors[color][0];
    imgData.data[i + 1] = colors[color][1];
    imgData.data[i + 2] = colors[color][2];
    imgData.data[i + 3] = 255;
    i += 4;
}
<canvas id="canvas" hidden></canvas>
<br>
<canvas id="data" hidden></canvas>
<br>
<input type="file" id="file"></input>
<br>
<img id="img">
<img id="reload" hidden>
<p>Algorithms:</p>
<select id="selectCord">
    <option value="1">Random</option>
    <option value="2">Random Alligned</option>
    <option value="3" selected>Loop</option>
</select>
<select id="selectColor">
    <option value="2000">Super Speedy</option>
    <option value="1000">Very Speedy</option>
    <option value="500" selected>Speedy</option>
    <option value="1">Visual</option>
</select>
<select id="colorDif">
    <option value="RGB" selected>RGB</option>
    <option value="HSL">HSL</option>
    <option value="HSV">HSV</option>
</select>
<p>Algorithm Options:
    <br>
</p>
<p>X Offset:
    <input id="allignX" type="range" min="0" max="3" value="0"></input>
</p>
<p>Y Offset:
    <input id="allignY" type="range" min="0" max="3" value="0"></input>
</p>
<p>Spacing X:
    <input id="loopX" type="range" min="1" max="4" value="2"></input>
</p>
<p>Spacing Y:
    <input id="loopY" type="range" min="1" max="4" value="2"></input>
</p>

RGB:

HSL:

HSV:

Very cool. The "Run code snippet" breaks for me when it tries to set document.cookie (after 1000 iterations) because the document is sandboxed. Is the cookie necessary?

No, once upon a time I ran the program for a few hours, but then my browser crashed. So I baked the cookie as a backup thing. But I'll remove it because it seems stack exchange dislikes cookies.

Looking at your code, I think it might benefit from the same speedup I suggested on wolfhammer's answer, except applied to doThing instead of loop. You might find the speed increase worth the extra line...

@trichoplax Thanks alot, Not only did your fix increase the speed of my program, while fixing I found and fixed a math error I made, and my program no longer generates those tiny black dots.

That's great! The new output image looks much better.

C# (reference implementation)
This is the code used to generate the images in the question. I thought it would be useful to give some people a reference for organizing their algorithm. A completely random coordinate and color are selected each movement. It performs surprisingly well considering the limitations imposed by the brush size/acceptance criteria.
I use the CIEDE2000 algorithm for measuring color differences, from the open source library ColorMine. This should give closer color matches (from a human perspective) but it doesn't seem to be a noticeable difference when used with this palette.
using System;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;
using ColorMine.ColorSpaces;
using ColorMine.ColorSpaces.Comparisons;

namespace Painter
{
    public class Painter
    {
        private readonly Bitmap _source;
        private readonly Bitmap _painting;
        private readonly int _width;
        private readonly int _height;
        private readonly CieDe2000Comparison _comparison = new CieDe2000Comparison();
        private const int BRUSHSIZE = 4;
        private readonly Random _random = new Random();
        private readonly ColorPalette _palette;

        private static readonly int[][] BRUSH = {
            new[] {1, 0}, new[] {2, 0},
            new[] {0, 1}, new[] {1, 1}, new[] {2, 1}, new[] {3, 1},
            new[] {0, 2}, new[] {1, 2}, new[] {2, 2}, new[] {3, 2},
            new[] {1, 3}, new[] {2, 3}
        };

        public Painter(string sourceFilename, string paletteFilename)
        {
            _source = (Bitmap)Image.FromFile(sourceFilename);
            _width = _source.Width;
            _height = _source.Height;

            _palette = Image.FromFile(paletteFilename).Palette;
            _painting = new Bitmap(_width, _height, PixelFormat.Format8bppIndexed) {Palette = _palette};

            // search for black in the color palette
            for (int i = 0; i < _painting.Palette.Entries.Length; i++)
            {
                Color color = _painting.Palette.Entries[i];
                if (color.R != 0 || color.G != 0 || color.B != 0) continue;
                SetBackground((byte)i);
            }
        }

        public void Paint()
        {
            // pick a color from the palette
            int brushIndex = _random.Next(0, _palette.Entries.Length);
            Color brushColor = _palette.Entries[brushIndex];

            // choose coordinate
            int x = _random.Next(0, _width - BRUSHSIZE + 1);
            int y = _random.Next(0, _height - BRUSHSIZE + 1);

            // determine whether to accept/reject brush
            if (GetBrushAcceptance(brushColor, x, y))
            {
                BitmapData data = _painting.LockBits(new Rectangle(0, y, _width, BRUSHSIZE), ImageLockMode.ReadWrite, PixelFormat.Format8bppIndexed);
                byte[] bytes = new byte[data.Height * data.Stride];
                Marshal.Copy(data.Scan0, bytes, 0, bytes.Length);

                // apply 4x4 brush without corners
                foreach (int[] offset in BRUSH)
                {
                    bytes[offset[1] * data.Stride + offset[0] + x] = (byte)brushIndex;
                }
                Marshal.Copy(bytes, 0, data.Scan0, bytes.Length);
                _painting.UnlockBits(data);
            }
        }

        public void Save(string filename)
        {
            _painting.Save(filename, ImageFormat.Png);
        }

        private void SetBackground(byte index)
        {
            BitmapData data = _painting.LockBits(new Rectangle(0, 0, _width, _height), ImageLockMode.WriteOnly, PixelFormat.Format8bppIndexed);
            byte[] bytes = new byte[data.Height * data.Stride];
            for (int i = 0; i < data.Height; i++)
            {
                for (int j = 0; j < data.Stride; j++)
                {
                    bytes[i*data.Stride + j] = index;
                }
            }
            Marshal.Copy(bytes, 0, data.Scan0, bytes.Length);
            _painting.UnlockBits(data);
        }

        private bool GetBrushAcceptance(Color brushColor, int x, int y)
        {
            double currentDifference = 0.0;
            double brushDifference = 0.0;
            foreach (int[] offset in BRUSH)
            {
                Color sourceColor = _source.GetPixel(x + offset[0], y + offset[1]);
                Rgb sourceRgb = new Rgb {R = sourceColor.R, G = sourceColor.G, B = sourceColor.B};
                Color currentColor = _painting.GetPixel(x + offset[0], y + offset[1]);

                currentDifference += sourceRgb.Compare(new Rgb {R = currentColor.R, G = currentColor.G, B = currentColor.B}, _comparison);
                brushDifference += sourceRgb.Compare(new Rgb {R = brushColor.R, G = brushColor.G, B = brushColor.B}, _comparison);
            }
            return brushDifference < currentDifference;
        }
    }
}

Then you can generate series of images (like my video) by calling an instance in a way similar to the code below (tweak based on number of iterations/frames/name desired). The first argument is the file path to the source image, the second argument is the file path to the palette (linked in the question), and the third argument is the file path for output images.
namespace Painter
{
    class Program
    {
        private static void Main(string[] args)
        {
            int i = 0;
            int counter = 1;
            Painter painter = new Painter(args[0], args[1]);
            while (true)
            {
                painter.Paint();
                if (i%500000 == 0)
                {
                    counter++;
                    painter.Save(string.Format("{0}{1:D7}.png", args[2], counter));
                }
                i++;
            }
        }
    }
}

I searched for some colorful canvas paintings online and came across the images below, which seem to be great (complicated) test images. All copyrights belong to their respective owners.

Source

Source

Source

This is how I know I don't know anything about anything. Nice solution

JavaScript Canvas
Update
Excellent suggestions in the comments. It's now faster and doesn't slow down the ui!

function previewFile() {
  var srcImage = document.getElementById('srcImage');
  var file = document.querySelector('input[type=file]').files[0];
  var reader = new FileReader();

  reader.onloadend = function() {
    srcImage.src = reader.result;
  }

  if (file) {
    reader.readAsDataURL(file);
  } else {
    srcImage.src = "";
  }
}

var buckets = []; // buckets / iterations
var iter_per_focus = 5000;

var pal = "00FFFF,0000FF,00FF00,FFFF00,\
C0C0C0,FF0000,FF00FF,FFFF78,\
FF0078,FF7848,7878FF,78FFFF,\
00FF78,784800,007800,\
007878,787800,780000,787878,\
000078,780078,004878,7800FF,\
0078FF,004848,787848,000000,FFFFFF".split(",");
var pLen = pal.length;
var p = 0;
var _R = 0;
var _G = 1;
var _B = 2;
var _CAN = 3;

// Create fast access palette with r,g,b values and
// brush image for color.
function initPal() {

  for (var i = 0; i < pal.length; i++) {
    var r = parseInt(pal[i].substr(0, 2), 16);
    var g = parseInt(pal[i].substr(2, 2), 16);
    var b = parseInt(pal[i].substr(4, 2), 16);
    var pcan = document.createElement('canvas');
    pcan.width = 4;
    pcan.height = 4;
    var pctx = pcan.getContext('2d');
    pctx.fillStyle = '#' + pal[i];
    pctx.beginPath();
    pctx.rect(1, 0, 2, 4);
    pctx.rect(0, 1, 4, 2);
    pctx.fill();

    pal[i] = [r,g,b,pcan];

  }
}
initPal();

var score = [];
var can = document.getElementById("canB");
var ctx = can.getContext('2d');
var mainDiv = document.getElementById("main");
var bCan = document.createElement('canvas');
bCan.width = can.width;
bCan.height = can.height;
var bCtx = bCan.getContext('2d');

var canA = document.getElementById("canA");
can.width = can.height = canA.width = canA.height = 200;
var ctxA = canA.getContext('2d');
var imageData;
var data;

function getSrcImage() {
  var img = document.getElementById('srcImage');
  can.width = canA.width = img.width;
  can.height = canA.height = img.height;
  ctxA.drawImage(img, 0, 0);
  imageData = ctxA.getImageData(0, 0, img.width, img.height);
  data = imageData.data;

  // adjust for brush offset
  var w = can.width - 2;
  var h = can.height - 2;

  var n = Math.floor((w * h) / iter_per_focus);
  buckets = [];
  for (var i = 0; i < n; i++) {
    var bucket = [];
    bucket.r = Math.floor(Math.random() * pLen);
    buckets.push(bucket);
  }
  var b = 0;
  var pt = 0;
  for (var y = 0; y < h; y++) {
    for (var x = 0; x < w; x++, pt+=4) {
      var r = Math.floor((Math.random() * n));
      buckets[r].push([x,y,pt,256 * 12,Math.floor(Math.random()*pLen)]);
      b %= n;
    }
    pt += 8; // move past brush offset.
  }

}

var loopTimeout = null;

function loopInit() {
  var r, y, x, pt, c, s;
  var row = can.width * 4;

  var b = 0;

  function loop() {
    clearTimeout(loopTimeout);
    var bucket = buckets[b++];
    var len = bucket.length;
    // Stepping color
    //c = pal[p];
    // Pulsing color;
    //c = pal[Math.floor(Math.random()*pLen)]
    // Pulsting step
    c = pal[bucket.r++];
    bucket.r%=pLen;
    b %= buckets.length;
    if (b === 0) {
      p++;
      p%=pLen;
    }

    for (var i = 0; i < len; i++) {

      var x = bucket[i][0]
      var y = bucket[i][1];
      var pt = bucket[i][2];
      // Random color
      //c = pal[bucket[i][4]++];
      //bucket[i][4]%=pLen;

      s = Math.abs(data[pt] - c[_R]) +
        Math.abs(data[pt + 1] - c[_G]) +
        Math.abs(data[pt + 2] - c[_B]) +
        Math.abs(data[pt + 4] - c[_R]) +
        Math.abs(data[pt + 5] - c[_G]) +
        Math.abs(data[pt + 6] - c[_B]) +
        Math.abs(data[pt + row] - c[_R]) +
        Math.abs(data[pt + row + 1] - c[_G]) +
        Math.abs(data[pt + row + 2] - c[_B]) +
        Math.abs(data[pt + row + 4] - c[_R]) +
        Math.abs(data[pt + row + 5] - c[_G]) +
        Math.abs(data[pt + row + 6] - c[_B]);
      if (bucket[i][3] > s) {
        bucket[i][3] = s;
        bCtx.drawImage(c[_CAN], x - 1, y - 1);
      }

    }
    loopTimeout = setTimeout(loop, 0);
  }

  loop();
}

// Draw function is separate from rendering. We render
// to a backing canvas first.
function draw() {
  ctx.drawImage(bCan, 0, 0);
  setTimeout(draw, 100);
}

function start() {

  getSrcImage();
  imageData = ctxA.getImageData(0, 0, can.width, can.height);
  data = imageData.data;
  bCan.width = can.width;
  bCan.height = can.height;
  bCtx.fillStyle = "black";
  bCtx.fillRect(0, 0, can.width, can.height);
  loopInit();

  draw();
}
body {
  background-color: #444444;
  color: #DDDDEE;
}
#srcImage {
  display: none;
}
#testBtn {
  display: none;
}
#canA {
  display:none;
}
<input type="file" onchange="previewFile()">
<br>
<img src="" height="200" alt="Upload Image for MS Painting">

<button onclick="genImage()" id="testBtn">Generate Image</button>

<div id="main">
  <img id="srcImage" src="" onload="start()">
  <canvas id="canA"></canvas>
  <canvas id="canB"></canvas>
</div>

@trichoplax I was having some cross site issues with loading the image. I'll see if I can figure something out.

@trichoplax No the dark wasn't intentional. It was a bug with the transparency in the generated image. The comparing code thought transparent should be black.

@trichoplax I've changed it to only compare a random color.

All works perfectly now :)

I'd love to see this converge faster. Is the s /= 4 needed? Does it make a difference since s is only compared with its own previous values?

Some example output images would help for voters who don't have time to wait for convergence.

I copied your code into a jsfiddle and tried an experiment. It made the convergence somewhat faster. You might like to try it... All I did was surround the contents of the loop function with a for loop to repeat the contents 1000 times. This means mouse and keyboard events are only checked for every 1000 iterations instead of every iteration. Your loop is fast enough that every 1000 iterations still leaves the mouse and keyboard responsive, and it saves waiting hours for convergence :)

@tricholplax Wow those suggestions made things a lot faster.  I think s /= 4 is needed, I don't get the cool color animation with it.

Mathematica
It's not really all that fast though, but at least it makes vaguely recognisable images, so I'm happy.
img = Import["https://i.stack.imgur.com/P7X6g.jpg"]
squigglesize = 20;
squiggleterb = 35;
colors = Import["https://i.stack.imgur.com/u9JAD.png"];
colist = Table[RGBColor[PixelValue[colors, {x, 1}]], {x, 28}];
imgdim0 = ImageDimensions[img];
curimg = Image[ConstantArray[0, Reverse[imgdim0]]];

rp := RandomInteger[squigglesize, 2] - squigglesize/2;
i = 0; j = 0;
Module[{imgdim = imgdim0, randimg, points, randcol, squigmid, st,
  randist, curdist = curdist0, i = 0, j = 0},

 While[j < 10,
  squigmid = {RandomInteger[imgdim[[1]]], RandomInteger[imgdim[[1]]]};
  While[i < 20,
   randcol = RandomChoice[colist];
   st = RandomInteger[squiggleterb, 2] - squiggleterb/2;
   points = {rp + squigmid + st, rp + squigmid + st, rp + squigmid + st, rp + squigmid + st};

   randimg =
    Rasterize[
     Style[Graphics[{Inset[curimg, Center, Center, imgdim],
        {randcol, BezierCurve[Table[{-1, 0}, {4}] + points]},
        {randcol, BezierCurve[Table[{-1, 1}, {4}] + points]},
        {randcol, BezierCurve[Table[{0, -1}, {4}] + points]},
        {randcol, BezierCurve[points]},
        {randcol, BezierCurve[Table[{0, 1}, {4}] + points]},
        {randcol, BezierCurve[Table[{0, 2}, {4}] + points]},
        {randcol, BezierCurve[Table[{1, -1}, {4}] + points]},
        {randcol, BezierCurve[Table[{1, 0}, {4}] + points]},
        {randcol, BezierCurve[Table[{1, 1}, {4}] + points]},
        {randcol, BezierCurve[Table[{1, 2}, {4}] + points]},
        {randcol, BezierCurve[Table[{2, 0}, {4}] + points]},
        {randcol, BezierCurve[Table[{2, 1}, {4}] + points]}
       }, ImageSize -> imgdim, PlotRange -> {{0, imgdim[[1]]}, {0, imgdim[[2]]}}],
      Antialiasing -> False], RasterSize -> imgdim];
   randist = ImageDistance[img, randimg];
   If[randist < curdist, curimg = randimg; curdist = randist; i = 0;
    j = 0;];
   i += 1;
   ]; j += 1; i = 0;];
 Print[curimg]]

Output:

Output could probably be a bit better with more iterations, and there's still a lot that I can try to speed it up/improve convergence, but for now this seems good enough.

SmileBASIC
OPTION STRICT
OPTION DEFINT

DEF MSPAINT(IMAGE,WIDTH,HEIGHT,STEPS)
 'read color data
 DIM COLORS[28]
 COPY COLORS%,@COLORS
 @COLORS
 DATA &H000000,&H808080,&H800000
 DATA &H808000,&H008000,&H008080
 DATA &H000080,&H800080,&H808040
 DATA &H004040,&H0080FF,&H004080
 DATA &H8000FF,&H804000,&HFFFFFF
 DATA &HC0C0C0,&HFF0000,&HFFFF00
 DATA &H00FF00,&H00FFFF,&H0000FF
 DATA &HFF00FF,&HFFFF80,&H00FF80
 DATA &H80FFFF,&H8080FF,&HFF0080
 DATA &HFF8040

 'create output array and fill with white
 DIM OUTPUT[WIDTH,HEIGHT]
 FILL OUTPUT,&HFFFFFFFF

 VAR K
 FOR K=1 TO STEPS
  'Pick random position/color
  VAR X=RND(WIDTH -3)
  VAR Y=RND(HEIGHT-3)
  VAR COLOR=COLORS[RND(28)]

  'Extract average (really the sum) color in a 4x4 area.
  'this is less detailed than checking the difference of every pixel
  'but it's better in some ways...
  'corners are included so it will do SOME dithering
  'R1/G1/B1 = average color in original image
  'R2/G2/B2 = average color in current drawing
  'R3/G3/B3 = average color if brush is used
  VAR R1=0,G1=0,B1=0,R2=0,G2=0,B2=0,R3=0,G3=0,B3=0
  VAR R,G,B
  VAR I,J
  FOR J=0 TO 3
   FOR I=0 TO 3
    'original image average
    RGBREAD IMAGE[Y+J,X+I] OUT R,G,B
    INC R1,R
    INC G1,G
    INC B1,B
    'current drawing average
    RGBREAD OUTPUT[Y+J,X+I] OUT R,G,B
    INC R2,R
    INC G2,G
    INC B2,B
    'add the old color to the brush average if we're in a corner
    IF (J==0||J==3)&&(I==0||I==3) THEN
     INC R3,R
     INC G3,G
     INC B3,B
    ENDIF
   NEXT
  NEXT
  'brush covers 12 pixels
  RGBREAD COLOR OUT R,G,B
  INC R3,R*12
  INC G3,G*12
  INC B3,B*12

  'Compare
  BEFORE=ABS(R1-R2)+ABS(G1-G2)+ABS(B1-B2)
  AFTER =ABS(R1-R3)+ABS(G1-G3)+ABS(B1-B3)

  'Draw if better
  IF AFTER<BEFORE THEN
   FILL OUTPUT,COLOR, Y   *WIDTH+X+1,2 ' ##
   FILL OUTPUT,COLOR,(Y+1)*WIDTH+X  ,4 '####
   FILL OUTPUT,COLOR,(Y+2)*WIDTH+X  ,4 '####
   FILL OUTPUT,COLOR,(Y+3)*WIDTH+X+1,2 ' ##
  ENDIF
 NEXT

 RETURN OUTPUT
END

MSPAINT image%[] , width% , height% , steps% OUT output%[]

image% - 2D [y,x] integer array with the image data (32 bit ARGB format (alpha is ignored))
width% - image width
height% - image height
steps% - number of iterations
output% - output array, same as image%.

can you add some examples?

Yeah, I'll add some soon (It's a lot of work to transfer images though, so I'll have to just take pictures of the screen for now)

Java + Processing
I saw this challenge and immediately had to try it.
My algorithm isn't too well optimized, so the output is a bit choppy.
color[] palette={#000000,#FFFFFF,#C0C0C0,#808080,#7E0308,#FC0D1B,#FFFD38,#807F17,#0F7F12,#29FD2E,#2DFFFE,#11807F,#020C7E,#0B24FB,#FC28FC,#7F0F7E,#807F44,#FFFD88,#2AFD85,#044040,#1783FB,#85FFFE,#8184FC,#06427E,#7F25FB,#7F400D,#FD8049,#FC1681};
PImage p;
int currTime;
int sTime=0;
int yPix;
void settings() {
  p=loadImage("arun-clarke-V-bkeYsTflY-unsplash.jpg");
  yPix = 0;

  size(p.width,p.height);
}
void setup() {
  frameRate(p.width);
    noStroke();
  rectMode(CENTER);
  //image(p,0,0);
  //for(int i=0; i<palette.length;i++) {
  //  fill(palette[i]);
  //  square(i*20,0, 20);
  //}
}

void draw() {

    int x =  int(random(p.width));
    int y = int(random(p.height));
    color randPix = p.get(x,y);
    float[] d = new float[palette.length];
    for(int k=0;k<palette.length;k++) {
    d[k]=abs(red(randPix)-red(palette[k]))+abs(blue(randPix)-blue(palette[k]))+abs(green(randPix)-green(palette[k]));
    }int b=0;
    for(int l=0;l<palette.length;l++){
      if(d[l]<=d[b]){
        b=l;
      }
    }
    fill(palette[b]);
    rect(x,y,2,4);
    rect(x,y,4,2);

}
void mouseClicked() {
  noLoop();
  save("c.png");
}

Examples


text:
Sum of smallest prime factors
SF(n) is a function which computes the smallest prime factor for a given number n.
We'll call T(N) the sum of every SF(n) with 2 <= n <= N.
T(1) = 0 (the sum is over 0 summands)
T(2) = 2 (2 is the first prime)
T(3) = 5 = 2 + 3
T(4) = 7 = 2 + 3 + 2
T(5) = 12 = 2 + 3 + 2 + 5
...
T(10000) = 5786451
The winner will be the one who manages to compute the largest T(N) in 60 seconds on my own laptop (Toshiba Satellite L845, Intel Core i5, 8GB RAM).

Current top score: Nicolás Siplis - 3.6e13 points - Nim

Pf(2) = 2, Pf(3) = 3, So, T(3) = 2 + 3 = 5. Am I right? I do programmed to find prime factors, but could you elaborate about the current requirement. Thank you

That is correct, I'll clarify it in the post.

I think this could be done relatively efficiently, using a smallest to largest sorting algorithm and then computing the sum of numbers up to length N of the sort. Doing it any other way, seems like a Greedy Method if you are working out the sum, while calculating each prime.

How is "60 seconds" to be measured? CPU seconds? Or single-core seconds? (A multi-threaded implementation might use 480 CPU seconds yet complete in 60 wall-clock seconds.) Run this on our own personal computer? Fastest CPU we can get our hands on?

@ToddLehman I'm running each code in my own laptop (Sony Vaio SVF14A16CLB), so if it takes less than 60 seconds I'll increase the number and decrease it when it takes longer.

So if our program forks multiple threads that run in parallel and still produces the correct answer in <= 60 seconds, that is acceptable?

Yes, as long as it runs on my own machine and outputs the correct answer in 60 seconds or less, it is acceptable.

How many CPU cores is the Sony Vaio SVF14A16CLB?

It has 4 threads.

Does third party Libraries are allowed? It is ok if the program is creating threads?

Third party libraries are allowed, but make sure to clarify how to install them. Creating threads is allowed as well.

@NicolásSiplis It might be helpful to order the solutions in your post by n.

Agreed, just finished doing so.

I KNOW there's a formulaic answer to this of the form ${T(N) = \sum_{i \in {primes \le N}} i * COUNT(N, i)}, where ${COUNT(N, i) \approx. \frac{N - COUNT(N, i - 1)}{i}}.  I kept getting thrown off by off-by-one errors that skewed the results  ~_~

Nim, 3.6e13
Simply sieving is not the best answer when trying to calculate the highest N possible since the memory requirements become too high. Here's a different approach (started with Nim a couple days ago and fell in love with the speed and syntax, any suggestions to make it faster or more readable are welcome!).
import math
import sequtils
import nimlongint # https://bitbucket.org/behrends/nimlongint/

proc s(n : int) : int128 =
    var x = toInt128(n)
    (x * x + x) div 2 - 1

proc sum_pfactor(N : int) : int128 =
    var
        root = int(sqrt(float(N)))
        u = newSeqWith(root+1,false)
        cntA,cntB,sumA,sumB = newSeq[int128](root+1)
        pcnt,psum,ret : int128
        interval,finish,d,q,t : int

    for i in 0..root:
        cntA[i] = i-1
        sumA[i] = s(i)

    for i in 1..root:
        cntB[i] = N div i - 1
        sumB[i] = s(N div i)

    for p in 2..root:
        if cntA[p] == cntA[p-1]:
            continue

        pcnt = cntA[p - 1]
        psum = sumA[p - 1]
        q = p * p
        ret = ret + p * (cntB[p] - pcnt)
        cntB[1] = cntB[1] - cntB[p] + pcnt
        sumB[1] = sumB[1] - (sumB[p] - psum) * p
        interval = (p and 1) + 1
        finish = min(root,N div q)

        for i in countup(p+interval,finish,interval):

            if u[i]:
                continue

            d = i * p

            if d <= root:
                cntB[i] = cntB[i] - cntB[d] + pcnt
                sumB[i] = sumB[i] - (sumB[d] - psum) * p
            else:
                t = N div d
                cntB[i] = cntB[i] - cntA[t] + pcnt
                sumB[i] = sumB[i] - (sumA[t] - psum) * p

        if q <= root:
            for i in countup(q,finish-1,p*interval):
                u[i] = true

        for i in countdown(root,q-1):
            t = i div p
            cntA[i] = cntA[i] - cntA[t] + pcnt
            sumA[i] = sumA[i] - (sumA[t] - psum) * p

    sumB[1] + ret

var time = cpuTime()
echo(sum_pfactor(int(3.6e13))," - ",cpuTime() - time)

I tried implementing the GMP wrapper for Nim into my code but just couldn't make it work (never used GMP before so that certainly didn't help).

You also don't need the return in f's definition. Single-expression procs automatically return.

This isn't the first [tag:fastest-code] that Nim has won by a noticeable margin. Might be worth investigating.

I'm curious to see how it performs when using GMP, but couldn't implement it correctly despite my efforts.

Nim is definitely going on my to-learn list!

We definitely need more answers in nim!

C, Prime Sieve: 5e9
Results:
$ time ./sieve
Finding sum of lowest divisors of n = 2..5000000000
572843021990627911

real    0m57.144s
user    0m56.732s
sys 0m0.456s

Program:
While it's a rather staightforward program, it took me a while to figure out how to get the memory management right - I only have enough ram for 1 byte per number in the range, so I had to be careful. It's a standard sieve of Erasthones.
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<assert.h>

#define LIMIT ((unsigned long long)5e9 +1)
#define ROOT_LIMIT floor(sqrt(LIMIT))

int main()
{
    printf("Finding sum of lowest divisors of n = 2..%llu\n", LIMIT - 1);
    char * found_divisor;
    found_divisor = malloc(LIMIT * sizeof(char));
    if (found_divisor == NULL) {
        printf("Error on malloc");
        return -1;
    }
    unsigned long long i;
    unsigned long long trial_div;
    unsigned long long multiple;
    unsigned long long sum = 0;

    for (i = 0; i < LIMIT; ++i) {
        found_divisor[i] = 0;
    }

    for (trial_div = 2; trial_div <= ROOT_LIMIT; ++trial_div) {
        if (found_divisor[trial_div] == 0) {
            for (multiple = trial_div * trial_div; multiple < LIMIT; multiple += trial_div) {
                if (found_divisor[multiple] == 0) {
                    found_divisor[multiple] = 1;
                    sum += trial_div;
                }
            }
        }
    }

    for (i = 2; i < LIMIT; ++i) {
        if (found_divisor[i] == 0) {
            sum += i;
        }
    }

    free(found_divisor);
    printf("%lld\n", sum);
    return 0;
}

If memory is a concern, one bit per number should be enough. You can use a bitmask to store the flags.

@RetoKoradi Unfortunately, that would probably slow the program down enough to put it over the 1 minute mark.

What do you need assert.h for?

@MaxRied It was left over from an earlie version.

Perl, brute force factoring
use ntheory ":all";
sub T {
  my $sum=0;
  for (1..$_[0]) {
    $sum += !($_%2) ? 2 : !($_%3) ? 3 : !($_%5) ? 5 : (factor($_))[0];
  }
  $sum
}
T(90_000_000);

I can get to about 9e7 in 25 seconds on my Linux machine.  It could be faster by digging into the C code, as it is saying after a check for 2/3/5, completely factor the number.
There are much more clever ways to do this using sieving. I thought a simple brute force way would be a start.  This is basically Project Euler problem 521, by the way.

If it's at all useful to know, in Python with a sieve I can only manage T(47000). I'm going to try something similar to what you're doing to see if it's faster.

Looks like not using a sieve is faster.. I was able to calculate T(493900) with a similar method to yours.

Never used Perl before but I managed to verify your answer, I'll add you to the list!

To be fair, this uses my module that does the factoring in C (you can force it to use pure Perl for everything, but of course it isn't as fast).

The answer can be computed using any combination of languages, so that's OK.

Go, 21e9
Does a sieve to find the minimum factor of each number <= N.  Spawns goroutines to count sections of the number space.
Run with "go run prime.go -P 4 -N 21000000000".
package main

import (
    "flag"
    "fmt"
    "runtime"
)

const S = 1 << 16

func main() {
    var N, P int
    flag.IntVar(&N, "N", 10000, "N")
    flag.IntVar(&P, "P", 4, "number of goroutines to use")
    flag.Parse()
    fmt.Printf("N = %d\n", N)
    fmt.Printf("P = %d\n", P)
    runtime.GOMAXPROCS(P)

    // Spawn goroutines to check sections of the number range.
    c := make(chan uint64, P)
    for i := 0; i < P; i++ {
        a := 2 + (N-1)*i/P
        b := 2 + (N-1)*(i+1)/P
        go process(a, b, c)
    }
    var sum uint64
    for i := 0; i < P; i++ {
        sum += <-c
    }
    fmt.Printf("T(%d) = %d\n", N, sum)
}

func process(a, b int, res chan uint64) {
    // Find primes up to sqrt(b).  Compute starting offsets.
    var primes []int
    var offsets []int
    for p := 2; p*p < b; p++ {
        if !prime(p) {
            continue
        }
        primes = append(primes, p)
        off := a % p
        if off != 0 {
            off = p - off
        }
        offsets = append(offsets, off)
    }

    // Allocate sieve array.
    composite := make([]bool, S)

    // Check factors of numbers up to b, a block of S at a time.
    var sum uint64
    for ; a < b; a += S {
        runtime.Gosched()
        // Check divisibility of [a,a+S) by our set of primes.
        for i, p := range primes {
            off := offsets[i]
            for ; off < S; off += p {
                if composite[off] {
                    continue // Divisible by a smaller prime.
                }
                composite[off] = true
                if a+off < b {
                    sum += uint64(p)
                }
            }
            // Remember offset for next block.
            offsets[i] = off - S
        }
        // Any remaining numbers are prime.
        for i := 0; i < S; i++ {
            if composite[i] {
                composite[i] = false // Reset for next block.
                continue
            }
            if a+i < b {
                sum += uint64(a + i)
            }
        }
    }
    res <- sum
}

func prime(n int) bool {
    for i := 2; i*i <= n; i++ {
        if n%i == 0 {
            return false
        }
    }
    return true
}

Note that the answer for N=21e9 is between 2^63 and 2^64, so I had to use unsigned 64-bit ints to count correctly...

I had to modify it to run on my machine (decreased N to 1e9) but the runtime itself is pretty fast, good work!

@NicolásSiplis: memory usage has been fixed.

Runtime was 80 seconds but 1.6e10 was calculated in almost exactly 60!

C++, 1<<34 ~ 1.7e10
Intel(R) Core(TM) i5-2410M CPU @ 2.30GHz

$ g++ -O2 test3.cpp
$ time ./a.out
6400765038917999291

real    0m49.640s
user    0m49.610s
sys 0m0.000s

#include <iostream>
#include <vector>

using namespace std;

const long long root = 1 << 17; // must be a power of two to simplify modulo operation
const long long rootd2 = root >> 1;
const long long rootd2m1 = rootd2 - 1;
const long long mult = root; // must be less than or equal to root
const long long n = root * mult; // unused constant (function argument)

int main() {
  vector < int > sieve(rootd2, 0);
  vector < int > primes;
  vector < long long > nexts;
  primes.reserve(root);
  nexts.reserve(root);
  // initialize sum with result for even numbers
  long long sum = n / 2 * 2;
  // sieve of Erathosthenes for numbers less than root
  // all even numbers are skipped
  for(long long i = 1; i < rootd2; ++i){
    if(sieve[i]){
      sieve[i] = 0;
      continue;
    }
    const long long val = i * 2 + 1;
    primes.push_back(val);
    sum += val;
    long long j;
    for(j = (val + 1) * i; j < rootd2; j += val){
      sum += val * (1 - sieve[j]); // conditionals replaced by multiplication
      sieve[j] = 1;
    }
    nexts.push_back(j);
  }
  int k = primes.size();
  long long last = rootd2;
  // segmented sieve of Erathosthenes
  // all even numbers are skipped
  for(int segment = 2; segment <= mult; ++segment){
    last += rootd2;
    for(int i = 0; i < k; ++i){
      long long next = nexts[i];
      long long prime = primes[i];
      if(next < last){
        long long ptr = next & rootd2m1; // modulo replaced by bitmasking
        while(ptr < rootd2){
          sum += prime * (1 - sieve[ptr]); // conditionals replaced by multiplication
          sieve[ptr] = 1;
          ptr += prime;
        }
        nexts[i] = (next & ~rootd2m1) + ptr;
      }
    }
    for(int i = 0; i < rootd2; ++i){
      sum += ((segment - 1) * root + i * 2 + 1) * (1 - sieve[i]);
      sieve[i] = 0;
    }
  }
  cout << sum << endl;
}

Java 8: 1.8e8 2.4e8
This entry does not compare to several of the other ones already up, but I wanted to post my answer since I had fun working on this.
The main optimizations of my approach are as follow:

Every even number has a smallest factor of 2, so these can be added for free after every odd number is processed.  Basically, if you have done the work to calculate T(N) when N % 2 == 1, you know that T(N + 1) == T(N) + 2.  This allows me to start my counting at three and to increment by iteration by twos.
I store my prime numbers in an array as opposed to a Collection type.  This more than doubled the N I can reach.
I use the prime numbers to factor a number as opposed to performing the Sieve of Eratosthenes.  This means that my memory storage is restricted almost completely to my primes array.
I store the square root of the number for which I am trying to find the smallest factor.  I tried @user1354678's approach of squaring a prime factor each time, but this cost me about than 1e7 from my score.

That's about all there is to it.  My code iterates from 3 on by twos until it detects that it has hit or exceeded the time limit, at which point it spits out the answer.
package sum_of_smallest_factors;

public final class SumOfSmallestFactors {
    private static class Result {
        private final int number;
        int getNumber() {
            return number;
        }

        private final long sum;
        long getSum() {
            return sum;
        }

        Result(int number, long sum) {
            this.number = number;
            this.sum = sum;
        }
    }

    private static final long TIME_LIMIT = 60_000_000_000L; // 60 seconds x 1e9 nanoseconds / second

    public static void main(String[] args) {
        SumOfSmallestFactors main = new SumOfSmallestFactors();
        Result result = main.run();
        int number = result.getNumber();
        long sum = result.getSum();
        System.out.format("T(%,d) = %,d\n", number, sum);
    }

    private int[] primes = new int[16_777_216];
    private int primeCount = 0;
    private long startTime;

    private SumOfSmallestFactors() {}

    private Result run() {
        startClock();
        int number;
        long sumOfSmallestFactors = 2;
        for (number = 3; mayContinue(); number += 2) {
            int smallestFactor = getSmallestFactor(number);
            if (smallestFactor == number) {
                addPrime(number);
            }
            sumOfSmallestFactors += smallestFactor + 2;
        }
        --number;

        Result result = new Result(number, sumOfSmallestFactors);
        return result;
    }

    private void startClock() {
        startTime = System.nanoTime();
    }

    private boolean mayContinue() {
        long currentTime = System.nanoTime();
        long elapsedTime = currentTime - startTime;
        boolean result = (elapsedTime < TIME_LIMIT);
        return result;
    }

    private int getSmallestFactor(int number) {
        int smallestFactor = number;
        int squareRoot = (int) Math.ceil(Math.sqrt(number));

        int index;
        int prime = 3;
        for (index = 0; index < primeCount; ++index) {
            prime = primes[index];

            if (prime > squareRoot) {
                break;
            }

            int remainder = number % prime;
            if (remainder == 0) {
                smallestFactor = prime;
                break;
            }
        }

        return smallestFactor;
    }

    private void addPrime(int prime) {
        primes[primeCount] = prime;
        ++primeCount;
    }
}

Running on a different system (Windows 8.1, Intel core i7 @ 2.5 GHz, 8 GB RAM) with the latest version of Java 8 has markedly better results with no code changes:
T(240,308,208) = 1,537,216,753,010,879

If you could replace the mayContinue() in for loop condition with just a simple condition, you could achieve higher result. And I like your way of precalculating even sum, then incrementing by two.

@user1354678, Thanks for the recommendation.  Strangely, it did not work.  I tried variations of this code on a different computer and found that the posted version is the fastest one.  Removing the clock calls from the code and using a simple threshhold number cost me a little over a second.  I even tried switching my startTime to an endTime to eliminate the ~2e7 subtractions, but that cost me 3e7 from my score!

Did you tried it with System.nanoTime() - startTime < TIME_LIMIT, coz it fastens your code a little for me. It's not blazingly fast, considering the fact, this condition is checked millions of times, it'll be a little fast. One thing I learned from your code is, don't put for inside a for.. After moving for to another method in my code, my code speed gets increased by 40%, Thanks.. One thing I'm still figuring out is, Did arrays are much efficient than ArrayList when considering the fact that it's fetched millions of times..

You could achieve x2 result if you implement MultiThreading. But it would need to precalculate the whole array before running the Prime calculation.

@user1354678, moving the check from the mayContinue() method into the for loop costs me 8e6 from my score.  This may be an issue of local optimizations.  I experimented with several data types for storing the primes when I was developing this solution.  I was only able to reach 8.8e7 with ArrayList, but I hit 1.8e8 (now 2.4e8) using an array.  There may be some performance boosts involved with the lookup, but there are definite boosts for the memory allocation.  I have thought about multi-threading the algorithm, but I ran into problems.

R, 2.5e7
Simple minded sieve of Eratosthenes, vectorised as much as possible. R's not really designed for this sort of problem and I'm pretty sure it can be made faster.
MAX <- 2.5e7
Tot <- 0
vec <- 2:MAX
while(TRUE) {
    if (vec[1]*vec[1] > vec[length(vec)]) {
        Tot <- Tot + sum(as.numeric(vec))
        break
    }

    fact <- which(vec %% vec[1] == 0)
    Tot <- Tot + vec[1]*length(vec[fact])
    vec <- vec[-fact]
}
Tot

Fair point about T.
2:MAX is a vector of integers so for large values of MAX, sum(vec) leads to an integer overflow and returns NA.
sum(as.numeric(vec)) is summing a vector of doubles which doesn't overflow (although it might not give the right answer)

Python, ~7e8
Using an incremental Sieve of Erathostenes. Some care does need to be taken that a marked value is marked with its lowest divisor, but the implementation is otherwise fairly straight forward.
Timing was taken with PyPy 2.6.0, input is accepted as a command line argument.
from sys import argv
from math import sqrt

n = int(argv[1])
sieve = {}
imax = int(sqrt(n))

t = n & -2
i = 3
while i <= n:
  divs = sieve.pop(i, [])
  if divs:
    t += divs[-1]
    for v in divs:
      sieve.setdefault(i+v+v, []).append(v)
  else:
    t += i
    if i <= imax: sieve[i*i] = [i]
  i += 2

print t

Sample Usage
$ pypy sum-lpf.py 10000
5786451

$ pypy sum-lpf.py 100000000
279218813374515

Common Lisp, 1e7
(defvar input 10000000)
(defvar numbers (loop for i from 2 to input collect i))
(defvar counter)
(defvar primes)

(setf primes (loop for i from 2 to (floor (sqrt input))
    when (loop for j in primes
        do (if (eq (mod i j) 0) (return nil))
        finally (return t))
    collect i into primes
    finally (return primes)))

(format t "~A~%"
    (loop for i in primes
        do (setf counter 0)
        summing (progn (setf numbers (remove-if #'(lambda (x) (if (eq (mod x i) 0) (progn (incf counter) t))) numbers))
                (* i counter)) into total
        finally (return (+ total (reduce #'+ numbers)))))

I've opted to first generate a list of prime numbers from 2 to (sqrt input), then test every value with the primes, whereas previously I would test against every number up to (sqrt input), which would be pointless (e.g., if a number is divisible by 4, it's also divisible by 2, so it's already accounted for.)
Thank God for side-effects while I'm at it. The remove-if both lowers the size of the list and counts how many elements were removed, so I just have to multiply that by whatever value the loop is on and add that into the running total.
(Fun fact: delete is the destructive equivalent of remove, but for whatever reason, delete is all sorts of slower than remove in this case.)

Never used Lisp before, I'm getting a compiler error when trying to run your code:
(defvar total 0)
(defvar counter 0)
(defvar input 10000)
(defvar numbers (loop for i from 2 to input collect i))
(loop for i from 2 to (floor (sqrt input))
(setf counter 0)
summing (prog2 (nsubstitute-if 0 #'(lambda (x) (if (eq (mod x i) 0) (progn (incf counter) t))) numbers)
(* i counter)
(setf numbers (remove 0 numbers))) into total
finally (return (+ total (reduce #'+ numbers))))

I'm using SBCL 1.0.38, but when I get home I'll update to the latest version and see how it goes. If you save it into a file, you can run it with "sbcl --script ".

I tried but still no luck, just in case I tried compiling online with Ideone but that didn't work either.

Oh, sorry, I forgot the "do" keyword on line 6. It should run now though, give it another shot.

Great, it computes 6e6 in 60 seconds on my machine! By the way, if I decide to enter my own code, do you know if I should submit it as an answer? I'm not sure if that would allow new submissions.

You can submit your own answer and others can still add theirs, although you could also just put it in your original question.

Perfect, thanks for the clarification!

Java 2.14e9
Pure Sieve of Eratosthenes with advantage of BitSet
I have calculated the Smallest Prime factor sum upto Integer.MAX_VALUE - 1 just in 33.89 s. But I can't proceed any larger because any further will lead to Integer Overflow on the Size of Bitset. So I'm working on to create another Bitset for the next set of Ranges. Until then, this is the fastest I'm able to generate.

T(214,74,83,646) = 109931450137817286 in 33.89 s
aka
T(2,147,483,646) = 109931450137817286 in 33.89 s

import java.util.BitSet;

public class SmallPrimeFactorSum {

    static int    limit     = Integer.MAX_VALUE - 1;

    // BitSet is highly efficient against boolean[] when Billion numbers were involved
    // BitSet uses only 1 bit for each number
    // boolean[] uses 8 bits aka 1 byte for each number which will produce memory issues for large numbers
    static BitSet primes    = new BitSet(limit + 1);
    static int    limitSqrt = (int) Math.ceil(Math.sqrt(limit));

    static long   start     = System.nanoTime();

    static long   sum       = 0;

    public static void main(String[] args) {
        genPrimes();
    }

    // Generate Primes by Sieve of Eratosthenes
    // Sieve of Eratosthenes is much efficient than Sieve of Atkins as
    // Sieve of Atkins involes Division, Modulus, Multiplication, Subtraction, Addition but
    // Sieve of Eratosthenes involves only addition
    static void genPrimes() {

        // Inverse the Bit values
        primes.flip(0, limit + 1);

        // Now all Values in primes will now be true,
        // True  represents     prime number
        // False represents not prime number

        // Set 0 and 1 as not Prime number
        primes.clear(0, 2);

        // Set all multiples of each Prime as not Prime;
        for ( int prime = 2; prime > 0 && prime <= limit && prime > 0; prime = primes.nextSetBit(prime + 1) ) {
            // Add Current Prime as its Prime factor
            sum += prime;
            // Skip marking if Current Prime > SQRT(limit)
            if ( prime > limitSqrt ) {
                continue;
            }
            // Mark Every multiple of current Prime as not Prime
            for ( int multiple = prime + prime; multiple <= limit && multiple > 0; multiple += prime ) {
                // Mark as not Prime only if it's true already
                if ( primes.get(multiple) ) {
                    // Add Current Prime as multiple's Prime factor
                    sum += prime;
                    primes.clear(multiple);
                }
            }
        }

        System.out.printf("T(%d) = %d in %.2f s", limit, sum, (System.nanoTime() - start) / 1000000000.0);
        //System.out.printf("Total Primes upto %d : %d\n", limit, primes.cardinality());
    }

}

Julia, 5e7
Surely Julia can do better but this is what I have for now. This does 5e7 in about 60 seconds on JuliaBox but I can't test locally yet. Hopefully by then I'll have thought of a more clever approach.
const PRIMES = primes(2^16)

function lpf(n::Int64)
    isprime(n) && return n
    for p in PRIMES
        n % p == 0 && return p
    end
end

function T(N::Int64)
    local x::Int64
    x = @parallel (+) for i = 2:N
        lpf(i)
    end
    x
end

Here we're creating a function lpf that iterates through sequential primes and checks the input for divisibility by each. The function returns the first divisor encountered, thereby obtaining the least prime factor.
The main function computes lpf on the integers from 2 to the input in parallel and reduces the result by summing.

Rust 1.5e9
A very naive Eratosthene sieve, but I felt Rust did not receive any love here !
// Expected (approximate) number of primes
fn hint(n:usize) -> usize {
    if n < 2 {
        1
    } else {
        n / ((n as f64).ln() as usize) + 1
    }
}

fn main() {
    let n:usize = match std::env::args().nth(1) {
        Some(s) => s.parse().ok().expect("Please enter a number !"),
        None => 10000,
    };
    let mut primes = Vec::with_capacity(hint(n));
    let mut sqrt = 2;
    let s = (2..).map(|n:u32| -> u32 {
        if (sqrt * sqrt) < n {
            sqrt += 1;
        }
        let (div, unseen) = match primes.iter().take_while(|&p| *p <= sqrt).filter(|&p| n % p == 0).next() {
            Some(p) => (*p, false),
            None => (n, true),
        };
        if unseen {
            primes.push(div);
        }
        div
    }).take(n-1).fold(0, |acc, p| acc + p);
    println!("{}", s);
}

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Fill in the Number Blanks
You are given a list of 2-tuples of positive integers, of the form:
a1 b1, a2 b2, ..., an bn
In other words, the 2-tuples are comma-separated, and the integers within each tuple are space-separated. All of the a numbers are in strictly increasing order, as well as the bs. Your goal is to "fill in the gaps." What is meant by that is the following:

If ai+1 == ai+1, don't do anything.
Else if ai+1 < ai+1, you will insert tuples into the sequence as follows: ai bi, ai+1 bi+1, ai+2 bi+1, ..., ai+1-1 bi+1, ai+1 bi+1.

In other words, if there is a "gap" in the ai, then you will insert more into the list with the inserted bs to be bi+1 (the end of the gap).
Example:
Input: 6 125, 7 170, 10 185
Output: 6 125, 7 170, 8 185, 9 185, 10 185
Scoring: the score of the submission is the number of bytes needed to construct a program to, when given a string as described above, output in the same format the "filled" list of tuples.
Lowest score wins!

Can we assume that there will be at least two tuples?

@Dennis yes you can.

So in general for input of three tuples a1 b1, a2 b2, a3 b3 the only pair that triggers an insert is a1 b1, a3 b3 requiring an insert of a2 b3, but where should it be inserted?

My interpretation is that a counts from a(min) to a(max) and outputs b for the next highest-or-equal value of a.  So because there's no a=8 and a=9, those both use the b value for a=10.

does the output have to be in that same format with space separation and no brackets?

and does it have to be comma-space separated or is comma separated okay

@PeterTaylor It should be inserted to keep the sorted ordering of the as.

@Hand-E-Food That's correct.

@Maltysen It has to be in the same format as given above.

I tried to edit a bit to make it clear. Rollback if that is not what you meant.

@jimmy23013 Thank you! I didn't know how to do it without LaTeX (why is it not supported?).

@Ryan It was disabled because of bugs.

@jimmy23013 That's unfortunate!

@ophact These seem totally unrelated to me

@RadvylfPrograms this was discussed in TNB, it was agreed that this question would be closed as a duplicate of the new one because the new one is "more usable and better specified". They are related because they both involve filling in gaps with the next seen number.

Pyth, 31 30 bytes
Jsrz7Pjds,R+ehfghTdcJ2\,}hJePJ

1 byte thanks to Jakube.
At the high level, for each number in the range of the a's, the list is filtered for a's greater than that number, then the b from the first such tuple is used as the b.
Demonstration.

CJam, 42 41 39 38 bytes
q~,]:T(\-2=),>{IST2/z~\{I<!}#=,',S}fI&

Try it online in the CJam interpreter.

Python 2, 192 180 171 Bytes
Unfortunately, it takes a ton of bytes to parse the input and format output. Also, thanks to Sp3000 for kicking my brain into gear :)
n=[map(int,x.split())for x in input().split(',')]
for i in range(len(n)-1):j=0;exec'print n[i][0]+j,`n[i+(j>0)][1]`+",",;j+=1;'*(n[i+1][0]-n[i][0])
print n[-1][0],n[-1][1]

Python 3.5.0b1+, 123 bytes
L=map(eval,(input()+",").split())
o=""
while L:
 *L,a,b=L;x=a-1
 while[x,b,*L][-2]<a:o=", %d %d"%(a,*b)+o;a-=1
print(o[2:])

Parsing the input was a pain, so I just ignored it and parsed it as is. This means that
6 125, 7 170, 10 185

gets parsed as
[6, (125,), 7, (170,), 10, (185,)]

hence the splat in the string formatting.

Why so specific as to the version?

@BetaDecay Probably relies on features introduced in the 3.5.0 beta 1 that were absent in previous releases.

@BetaDecay What Alex said, specifically PEP448 - Additional Unpacking Generalisations. Otherwise, instead of (a,*b) I'd have to do ((a,)+b) or (a,b[0]). Older versions of Python would give a SyntaxError.

JavaScript (ES6), 97 103
Finding the question wording hard to understand, I based my work on the example.
First try, there could be a shorter way with just 1 split Scanning the input string in a single pass using replace.

f=l=>
  l.replace(/\d+/g,b=>{
    for(i=p||a-1;a&&i++<a;p=a)o+=`, ${i} `+b;
    a=a?0:b
  },p=a=o='')
&&o.slice(2)

// TEST

out=x=>O.innerHTML+=x+'\n'

test=v=>out(v+'\n->'+f(v))

test('6 125, 7 170, 10 185')

<input id=I><button onclick='test(I.value)'>-></button>

.

Haskell, 152 bytes
p(x,y)=", "++show x++' ':show y
(a:b:c)%x|a<x=c%x|0<1=(x,b)
g a=p.(a%)=<<[0+head a..last$init a]
main=interact$drop 2.g.(map(read.fst).lex=<<).words

How it works
p is a simple string formatting helper function:
>>> p (6, 125)
", 6 125"

(%) is a "lookup" function for filling in blanks:
>>> let nums = [6, 125, 7, 135, 10, 185]
>>> nums % 8
(8,185)

[head a..last$init a] is the range of numbers to use: from the first integer parsed until the second-to-last integer parsed. (The 0+ helps the type-checker infer that the whole program deals with Num values; otherwise read doesn't know what to parse.) Then we concatMap (or =<<) our previous functions over the range to get one big string. That's what g does:
>>> g nums
", 6 125, 7 135, 8 185, 9 185, 10 185"

Finally, there's main. It's of the form interact$s, so all we care about is the function s :: String -> String that turns stdin into stdout.
First, to read input, we split into words:
>>> (words) "6 125, 7 170, 10 185"
["6","125,","7","170,","10","185"]

Then we concatMap the function map(read.fst).lex over this list. For one word, this does:
>>> (map(read.fst).lex) "125," :: [Int]
[125]

So for a list of words, you get concat [[6], [125], [7], ...] == [6, 125, 7...], which is the list we want to pass to g. Finally, we drop the leading ", " from the string g gives us.

CJam, 34 bytes
ri_[_q~,]:,2/::fe&:.e|+ee>Sf*", "*

I started with golfing Dennis's answer but finally it becomes nothing like the original.
Try it online.

Perl, 67 bytes
(66 characters code + 1 character command line option)
s/(\b\d+ \d+)(?=, (?!(??{$k=$1+1}) ).+( \d+))/"$1, ".$k.$2/e&&redo

Execution example:
perl -pe 's/(\b\d+ \d+)(?=, (?!(??{$k=$1+1}) ).+( \d+))/"$1, ".$k.$2/e&&redo' <<< "6 125, 7 170, 10 185"

Explanation:
General approach is to construct a regex which is capable of looking ahead to the next token to determine if it increases by one (regex is rubbish at arithmetic, so we add a perl execution to do the addition, then look for that). If it does not, we replace this match with the original text, plus the new token with the increased index. We then repeat this multiple times on the input until it does not change.
For the below explanation, a 'section' is the comma separated value, and the 'tokens' are the space separated values within this.
 s/
   (\b\d+ \d+)             # Matches the first two numbers in a section
   (?=,                 # Look ahead to the next section
     (?!(??{$k=$1+1}) ) # Make sure the first token in the next section is not increased by one, store this value in $k. Happily, although our $1 contains TWO space separated numbers, it increases the first number and ignores the second - strange behaviour in Perl.
     .+                 # Ignore the actual value of first token of the next section (i.e. the index)
     ( \d+)             # Create a backreference on the second (i.e. the value)
   )
 /                      # Start the replacement string...
   "$1, "               # Put the first section back into the result
   .$k                  # Append the increased index
   .$2                  # Append the value of the next token
 /e                     # Modifier - allows us to execute Perl inside the RegEx
 &&redo              # Repeat until it no longer changes the string

C#, 174 bytes
Indented for clarity:
class P{
    static void Main(string[]z){
        int a=0,b,i=-1;
        while(++i<z.Length){
            b=int.Parse(z[i++]);
            while(a1:l[:i+2]+=l[i]+1,;l[:i+3]+=l[i+4],;x+=1
for n in range(len(l)):print(l[n],end=','if n%2 else' ')

Bash + coreutils, 87
tr , \\n|while read a b
do for((p=p?p:a-1;p++<a;c=1)){
echo -n ${c+, }$p $b
}
p=$a
done

Test output:
$ echo "6 125, 7 170, 10 185" | ./fillblanks.sh
6 125, 7 170, 8 185, 9 185, 10 185

Java, 229 bytes
I think I've whittled this down as much as I can. I would be interested to see if there is a better approach that leads to a shorter Java answer:
String f(String s){String r="",g=" ",x[]=s.split(",\\s"),z[];int i=0,p=0,l=x.length,a,b;while(i<l){z=x[i].split(g);a=new Integer(z[0]);b=new Integer(z[1]);if(i++>0)while(a>++p)r+=p+g+b+", ";r+=a+g+b+(i==l?"":", ");p=a;}return r;}

Formatted:
String f(String s) {
    String r = "", g = " ", x[] = s.split(",\\s"), z[];
    int i = 0, p = 0, l = x.length, a, b;
    while (i < l) {
        z = x[i].split(g);
        a = new Integer(z[0]);
        b = new Integer(z[1]);
        if (i++ > 0)
            while (a > ++p)
                r += p + g + b + ", ";
        r += a + g + b + (i == l ? "" : ", ");
        p = a;
    }
    return r;
}

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
all non-empty sublists partitionings of a list
I had to write a code for finding all non-empty sublists partitionings of a list:
def f(s):
    if s:
        for j in range(1,len(s)+1):
            for p in f(s[j:]):
                yield [s[:j]]+p
    else:
        yield []

I managed to shorten it (for Python 3.4+) to:
def f(s, p=[]):
    if s:
        for j in range(1,len(s)+1):
            yield from f(s[j:],p+[s[:j]])
    else:
        yield p

using yield from and an accumulator argument p (prefix).
Any suggestions on how to shorten it even more?
Example output:
>>> for p in f([1,2,3,4]): print(p)
...
[[1], [2], [3], [4]]
[[1], [2], [3, 4]]
[[1], [2, 3], [4]]
[[1], [2, 3, 4]]
[[1, 2], [3], [4]]
[[1, 2], [3, 4]]
[[1, 2, 3], [4]]
[[1, 2, 3, 4]]
>>>

For starters, you can use single char variable names and drop a bunch of whitespace

Can we answer with code in other languages, too? This would work pretty well as a language agnostic challenge.

@JanDvorak This is asking for help in golfing the python program.

@okx hence my trying to convert this into a regular challenge via my comment. The answer would carry over just fine.

66 bytes
p=lambda l:[q+[l[i:]]for i in range(len(l))for q in p(l[:i])]or[l]

If the last part is changed from [l] to [[]], then it can also be used with tuples or strings, in addition to lists.

Try it online


text:
Implement bzip2's run-length encoding
Background
After applying the BWT (as seen in Burrows, Wheeler and Back) and the MTF (as seen in Move to the printable ASCII front), the bzip2 compressor applies a rather unique form of run-length encoding.
Definition
For the purpose of this challenge, we define the transformation BRLE as follows:
Given an input string s that consists solely of ASCII characters with code points between 0x20 and 0x7A, do the following:

Replace each run of equal characters by a single occurrence of the character and store number of repetitions after the first.
Encode the number of repetitions after the first occurrence of the character, using bijective base-2 numeration and the symbols { and }.
A non-negative integer n is encoded as the string bk…b0 such that n = 2ki(bk)+…+20i(b0), where i({) = 1 and i(}) = 2.
Note that this representation is always unique. For example, the number 0 is encoded as an empty string.
Insert the string of curly brackets that encodes the number of repetitions after the single occurrence of the corresponding character.

Step-by-step example
Input:  "abbcccddddeeeeeffffffggggggghhhhhhhh"
Step 1: "abcdefgh" with repetitions 0, 1, 2, 3, 4, 5, 6, 7
Step 2: "" "{" "}" "{{" "{}" "}{" "}}" "{{{"
Step 3: "ab{c}d{{e{}f}{g}}h{{{"

Task
Implement an involutive program or function that reads a single string from STDIN or as command-line or function argument and prints or returns either the BRLE or its inverse of the input string.
If the input contains no curly brackets, apply the BRLE. If the input contains curly brackets, apply its inverse.
Examples
INPUT:  CODEGOLF
OUTPUT: CODEGOLF

INPUT:  PROGRAMMING
OUTPUT: PROGRAM{ING

INPUT:  PUZ{LES
OUTPUT: PUZZLES

INPUT:  444488888888GGGGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
OUTPUT: 4{{8{{{G{{{{W{{{{{

INPUT:  y}}}{{
OUTPUT: yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy

Additional rules

You cannot use any built-ins that compute the BRLE or its inverse of a string.
You can use built-ins that:

Compute the RLE or RLD of a string, as long as the number of repetitions is not stored in bijective base-2.
Perform base conversion of any kind.

Your code may print a trailing newline if you choose STDOUT for output.
Your code has to work for any input of 1000 or less ASCII characters in the range from 0x20 to 0x7A, plus curly brackets (0x7B and 0x7D).
If the input contains curly brackets, you may assume that it is a valid result of applying the BRLE to a string.
Standard code golf rules apply. The shortest submission in bytes wins.

Why are builtins not allowed?

CJam, 50 48 bytes
l_{}`:T&1${_T#)_(@a?)+}%{(\2b)*}%@e`{(2b1>Tf=}%?

Thanks for Dennis for saving 2 bytes.
Try it online.
Explanation
l_              e# Read and duplicate input.
{}`:T           e# T = "{}"
&               e# If the input has '{ or '}:
    1$          e# Input.
    {           e# For each character:
        _T#)    e# If it is '{ or '}:
            _(  e# Return 0 for '{ or 1 for '}.
            @a  e# Otherwise, convert the character itself to an array (string).
        ?
        )+      e# If it is a number, increment and append to the previous array.
                e# If it is a string with at least 1 character, do nothing.
    }%
    {(\         e# For each character and bijective base 2 number:
        2b)*    e# Repeat the character 1 + that many times.
    }%
                e# Else:
    @           e# Input.
    e`          e# Run-length encoding.
    {(          e# For each character and length:
        2b1>    e# Convert the length to base 2 and remove the first bit.
        Tf=     e# Map 0 to '{ and 1 to '}.
    }%
?               e# End if.

Pyth, 48 50 bytes
J`Hs?m*hdi+1xLJtd2tczf-@zTJUz@Jzsm+ed@LJtjhd2rz8

2 bytes thanks to @Maltysen.
Demonstration. Test harness.
Explanation:
J`Hs?m*hdi+1xLJtd2tczf-@zTJUz@Jzsm+ed@LJtjhd2rz8
                                                    Implicit: z = input()
                                                    H is empty dict.
J`H                                                 J = repr(H) = "{}"
   s                                                Print the concatenation of
    ?                        @Jz                    If z and J share any chars:
                     f     Uz                       Filter range(len(z))
                      -@zTJ                         On the absence of z[T] in J.
                   cz                               Chop z at these indices.
                                                    just before each non '{}'.
                  t                                 Remove empty initial piece.
     m*hd                                           Map to d[0] *
         i       2                                  the base 2 number
            xLJtd                                   index in J mapped over d[:-1]
          +1                                        with a 1 prepended.
                                             rz8    Otherwise, run len. encode z
                                 m                  map over (len, char)
                                         jhd2       Convert len to binary.
                                        t           Remove leading 1
                                     @LJ            Map to element of J.
                                  +ed               Prepend char.
                                s                   Concatenate.

instead of "{}" you can use \H`, tied with CJam :)

@Jakube Sorry for the mixup.

OCaml, 252
t is the function that does the transformation.
#load"str.cma"open Str
let rec(!)=group_beginning and
g=function|1->""|x->(g(x/2)^[|"{";"}"|].(x mod 2))and($)i s=if g i=s then i else(i+1)$s and
t s=global_substitute(regexp"\(.\)\1*\([{}]*\)")(fun s->String.make(1$matched_group 2 s)s.[!1]^g(!2- !1))s

At first I was thinking I had to check for the presence of curly braces, but it turned out to be unnecessary. The decoding clearly has no effect on strings that are already decoded, and the encoding part proved equally harmless to ones that were already encoded.

the encoding part proved equally harmless does it? Encoding 4{{8{{{G{{{{W{{{{{ don't you get 4{{8{}G{{{W{{}?

@edc65 No, I get the answer specified in the examples. How are you testing it?

"4{{8{{{G{{{{W{{{{{" as input is not one of the examples. Did you tried it?

@edc65 It is the inverse of one of the examples and I tested them both ways. Yes, I tried it, both before posting and after your comment.

Ok, good. I pointed out the quoted sentence, as a "straightforward` encoding (as mine) would not be harmless at all with the given test case. Clearly your encoding part is more clever.

JavaScript (ES6), 162

f=s=>
(t=s[R='replace'](/[^{}][{}]+/g,n=>n[0].repeat('0b'+n[R](/./g,c=>c!='{'|0))))==s
?s[R](/(.)\1*/g,(r,c)=>c+r.length.toString(2).slice(1)[R](/./g,c=>'{}'[c])):t

// TEST
out=x=>O.innerHTML += x + '\n';

test=s=>O.innerHTML = s+' -> '+f(s) +'\n\n' + O.innerHTML;

[['CODEGOLF','CODEGOLF']
,['PROGRAMMING','PROGRAM{ING']
,['PUZ{LES','PUZZLES']
,['444488888888GGGGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW','4{{8{{{G{{{{W{{{{{']
,['y}}}{{','yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy']]
.forEach(v=>{
  w=f(v[0])
  out('Test ' + (w==v[1]?'OK':'Fail')+'\nInput:    '+v[0]+'\nExpected: '+v[1]+'\nResult:   '+w)
})
Your test: <input id=I style='width:300px'><button onclick='test(I.value)'>-></button>

Some explanation
Number n to BB2 using 0 and 1: (n+1).toString(2).slice(1)
String in BB2 to number: to_number ("0b1"+string) - that is, add a leftmost 1 binary digit and convert from binary (and decrease by 1, not needed in this specific instance).
Regular expression to find any character followed by { or }: /[^{}][{}]+/g
Regular expression to find repeated characters: /(.)\1*/g
Using that regexp in replace, first param is the "repeated" char (eventually repeat just 1 time), second param is the total repeated string, whose length is the number that I need to encode in BB2 already incremented by 1
... then put all together ...