text:
Write a FAST word equation solver
This question is similar to Write a word equation solver by David Frank but different in three important ways: (1) David asked for just the first solution found, and I am asking for all the solutions.  (2) I am asking to output the time required to solve the puzzle.  And (3), it is a completely different undertaking to write a program designed to be fast compared with one designed to be short.  A brute force method is probably best for code golf, but will not beat my ?amazing? time of 0.15 seconds.  -I say "amazing" only to encourage other to try it out.  I'm sure mine won't be the fastest.
(I'm going to copy a lot of the text from that question.)
Introduction
Consider the following example:
  CODE
+ GOLF
——————
 GREAT

This is an equation where each letter represents a decimal digit and the words represent natural numbers (similar letters represent similar digits and different letters represent different digits). The task is to match each letter with its digit value so that the equation is correct. One of the four solutions for the equation above is:
  9265
+ 1278
——————
 10543

Your task
Your task is to write a program or function to find ALL of the answers to such a word equation.  Your program should output all the answers and the time it took to find them in seconds.
Input
The input will be something like this, two words added together to make a third word.
Example:

CODE+GOLF=GREAT
AA+BB=CC

Spaces are omitted and only letters between capital A and Z will be used (no special or small letters).  This string can be read from the standard input, from a file, or as a function parameter.
Output
The output should look something like this.  To take luck out of the competition, your program must list all the solutions.
time = XX.XXX seconds
(1) 9265+1278=10543
(2) 9275+1268=10543
(3) 9428+1437=10865
(4) 9438+1427=10865

Rules

Similar letters represent similar digits and different letters represent different digits
To make things easier, you can assume that numbers can or don't start with a leading zero, it's up to you.  The four answers I got for the sample problem assume numbers do not start with leading zeros.
You can use any language and the standard library of the chosen language (no external libs)
You cannot connect to any resources on the internet (why would you anyways?)
So the goal here is speed, but of coarse that depends on what computer you use.  In order to account for the different speed of different computers, you could run my example on your computer and then report how many times faster or slower yours is than my sample answer.
Please explain how your code works.

Example in JavaScript
http://goo.gl/TdQczX

For fairness, submissions to fastest-code challenges should all be tested on the same machine - which usually means your machine. Are you prepared to install all of those submissions on your own machine? Alternatively, another option might be to have everyone run your own solution and their own, the score being the ratio between the two.

Also I think your example problem may be too small to accurately measure the timing for well-optimised solutions. You might want to provide a larger test case as the benchmark.

The same question with a different winning criterion is a duplicate. (Even if the original question has a really stupid winning criterion).

@PeterTaylor I don't see how you read that into the answers of that meta question. What they say (and what the usual consensus is) is that it depends on whether answers from one challenge are competitive and valid solutions to the other. Seeing that the other one is code golf, I suspect that no solution will be optimised for efficiency (most of them will probably brute-force it). Hence, I'd say the difference in winning criterion is quite significant.

@MartinBüttner, it's the plain meaning of the most voted answer, as further explained in Tyzoid's comment on his answer.

@PeterTaylor Right, I didn't see his comment. I personally think (and I've seen answers along those lines on meta, too - and the comment with 3 upvotes is one of them) that we should take into account whether any submission from the original challenge would be competitive against new submissions which tackle the new challenge in a way according to the new winning criterion. If three optimised answers to this are posted, who would copy over a golfed brute force solution? I'm pretty sure to win this one, you do have to solve it in a novel way.

(Furthermore, copying a golfed submission into a fastest-code challenge comes very close to justifying a delete-vote for "This does not provide an answer to the question.")

Hello Martin and Peter.  I am very new to this.  I can say that the answers for the code-golf question seemed to take many seconds or even minutes to run, and they did not find all the solutions.  They just found one.  The method I used in my example is pretty unique - I think.  But then again my code golf answers are typically 5-10 times longer than the shortest ones.  As for judging the speed of an algorithm, I would agree that it's approximate.

@MartinBüttner, even taking that view, I think this is one of the clearest cases: the scoring system on the other question is so badly exploitable that there's very little stopping an answer from scoring well on both that and speed.

@PeterTaylor That's a fair point, although in that case I'd rather view the old one a duplicate of this one, as this here is definitely the more interesting challenge (answers to which can be submitted to the other one by encoding it in a whitespace string). One more thing about Tyzoid's opinion: there is a reasonable chance that his strong opinion and the upvotes were correlated to the fact that the two winning criteria in question were much more similar than code-golf vs fastest-code.

Also, for comparison sake, my answer can be run on anyone's computer because it's JavaScript on a web page -not that this completely solves the problem.  (I'll add a link to it.)

Can multiple letters have the same digit? or are there only a maximum of 10 different letters per equation.

The "O" in GOLF and the "O" in CODE have the same value. This is what David meant by "Similar letters represent similar digits and different letters represent different digits."  Also two different letters can not have the same value.

I adjusted rule number 5 to try to standardize the speed of different computers.  You can answer in seconds if you want to be general/approximate, or you can answer as a factor of the sample answer if you want to be more precise.

Why not use most elegant solution? This is determine by popularity, also is applicable to each language.

@JeffSB If I had ABCDEFGHIJK, only 10 of 11 letters would be able to be assigned digits. Would K have one of the previously assigned digits, or can there only be 10 letters?

Adam, I'm fine with most votes.  It's about the fun of it.

Ourous, You can only have ten different letters, but letters can be repeated.  CODE+GOLF=GREAT has 13 letters, but only 10 unique letters - so it works.

@AdamSpeight Because why turn something into a popularity contest if you don't have to? (@JeffSB, please, please, stick to fastest code.)

Okay, fastest code it is.

I inserted the 2nd paragraph to address the duplicate issue.

Python3, 0.03 seconds:
def f(a, b, r, c = {}, s = [], k = 0, l = [], A = [], B = []):
   if not a and not b:
      yield f"{''.join(map(str, A))}+{''.join(map(str, B))}={''.join(map(str, l))}"
      return
   for x in c.get(a[-1], range(10)):
      if a[-1] in c or x not in s:
         _c, _s = {**c, a[-1]:[x]}, s+[x]
         for y in _c.get(b[-1], range(10)):
             if b[-1] in _c or y not in _s:
                n_c, v = {**_c, b[-1]:[y]}, k+x+y
                n_s = _s+[y]
                if a[:-1] and b[:-1]:
                    _v, _k = v%10, v//10
                    if [_v] == n_c.get(r[-1], [_v])and (r[-1] in n_c or _v not in n_s):
                        yield from f(a[:-1], b[:-1], r[:-1], {**n_c, r[-1]:[_v]}, n_s+[_v], _k, [_v]+l, [x]+A, [y]+B)
                else:
                    if len(str(v)) == len(r):
                        F = 1
                        for J, K in zip(r[::-1], str(v)[::-1]):
                            if n_c.get(J, [int(K)]) == [int(K)] and (J in n_c or int(K) not in n_s):
                               n_c[J] = [int(K)]
                               n_s = n_s + [int(K)]
                            else:
                               F = 0
                               break
                        if F: yield from f(a[:-1], b[:-1], '', n_c, n_s, 0, [v]+l, [x]+A, [y]+B)

t = time.time()
result = [*f('CODE', 'GOLF', 'GREAT')]
print(time.time() - t)
print(result)

Try it online!

I think this answer qualifies for 400 rep from this bounty and 500 rep if you add an explanation.

does not work on
print([*f('FIFTY', 'STATES', 'AMERICA')])

C++, 10000× faster (~20 µs)
#include <stdio.h>
#include <time.h>
#include <algorithm>

// digits
typedef int digit_t;
const digit_t BASE = 10;
typedef unsigned digit_mask_t;
inline char format_digit(digit_t d) { return '0' + d; }

// letters
typedef int letter_t;
const letter_t ALPHABETSIZE = 26;
typedef unsigned letter_mask_t;
inline bool is_letter(char c) { return unsigned(c - 'A') < 26u; }
inline letter_t parse_letter(char c) { return letter_t(c - 'A'); }

// terms
typedef int term_t;
inline term_t term_abs(term_t x) { return x < 0 ? -x : x; }

bool parse_formula(const char* formula, term_t coeffs[ALPHABETSIZE], letter_mask_t* ppresent, letter_mask_t* pleading) {
    // Convert BA+BB=CCC into A + 21*B - 111*C = 0.
    letter_mask_t present = 0, leading = 0;
    const char* p = formula;
    term_t side = +1, place = +1;
    for (;;) {
        const char* r;
        for (r = p; is_letter(*r); ++r) {}
        if (p < r)
            leading |= letter_mask_t(1) << parse_letter(*p);
        for (const char* q = r; q != p; --q) {
            letter_t l = parse_letter(q[-1]);
            coeffs[l] += place;
            present |= letter_mask_t(1) << l;
            place *= BASE;
        }
        p = r;
        switch (*p++) {
            case 0:
                *ppresent = present; *pleading = leading; return true;
            case '+':
                place = +side; break;
            case '-':
                place = -side; break;
            case '=':
                if (side < 0)
                    return false;
                place = side = -1;
                break;
            default:
                return false;
        }
    }
}

struct Var {
    digit_t value;
    term_t coeff;
    digit_t minval;  // 0 or 1
    letter_t name;
};

class Solver {
    const char* formula;
    letter_t num_vars;
    Var vars[BASE];
    void solve2(Var*, digit_mask_t, term_t, term_t);
public:
    bool solve(const char* formula);
};

bool Solver::solve(const char* formula) {
    this->formula = formula;
    term_t coeffs[ALPHABETSIZE] = {0};
    letter_mask_t present, leading;
    if (!parse_formula(formula, coeffs, &present, &leading))
        return false;
    // sort the used letters by magnitude of their coefficient
    letter_t letters[BASE];
    letter_t num_vars = 0;
    for (letter_t i = 0; i < ALPHABETSIZE; ++i) {
        if (present & (letter_mask_t(1) << i)) {
            if (num_vars == BASE)
                return false;
            letters[num_vars++] = i;
        }
    }
    if (num_vars == 0)
        return false;
    this->num_vars = num_vars;
    std::sort(letters, letters + num_vars, [&](letter_t a, letter_t b) { return term_abs(coeffs[a]) < term_abs(coeffs[b]); });
    // initialize tables
    term_t bounds[2] = {0};
    for (letter_t i = 0; i < num_vars; ++i) {
        vars[i].name = letters[i];
        vars[i].minval = digit_t((leading >> letters[i]) & 1);
        term_t coeff = vars[i].coeff = coeffs[letters[i]];
        bounds[coeff > 0] += coeff * (BASE - 1);
    }
    // call the recursive solver
    solve2(vars + num_vars - 1, 0, bounds[0], bounds[1]);
    return true;
}

void Solver::solve2(Var* pos, digit_mask_t digits_used, term_t lower_bound, term_t upper_bound) {
    term_t coeff = pos->coeff;
    digit_t d = BASE - 1, e = pos->minval;
    // adjust bounds to reflect a digit equal to BASE-1, instead of anything in [0,BASE)
    if (coeff < 0)
        upper_bound += coeff * d;
    else
        lower_bound += coeff * d;
    do {
        if (lower_bound <= 0 && upper_bound >= 0 && !(digits_used & (digit_mask_t(1) << d))) {
            pos->value = d;
            if (pos == vars) {
                // print the solution
                char digit_chars[ALPHABETSIZE];
                for (letter_t i = 0; i < num_vars; ++i)
                    digit_chars[vars[i].name] = format_digit(vars[i].value);
                for (const char* p = formula; *p; ++p)
                    putchar(is_letter(*p) ? digit_chars[parse_letter(*p)] : *p);
                putchar('\n');
            } else {
                solve2(pos - 1, digits_used | (digit_mask_t(1) << d), lower_bound, upper_bound);
            }
        }
        lower_bound -= coeff;
        upper_bound -= coeff;
    } while (--d >= e);
}

int main(int argc, char **argv) {
    int code = 0;
    for (int i = 1; i < argc; ++i) {
        puts(argv[i]);
        clock_t t = clock();
        code |= !Solver().solve(argv[i]);
        t = clock() - t;
        printf("%f\n", double(t) / CLOCKS_PER_SEC);
    }
    return code;
}

Attempt This Online!
This is roughly 10,000 times faster than JeffSB's solver and 2,500 times faster than Ajax1234's solver on my machine on the problem CODE+GOLF=GREAT. The speed ratio varies widely with the problem. For FIFTY+STATES=AMERICA, the JeffSB/me ratio is 25,000, while for CORRECTS+REJECTED=MATTRESS it's about 25.
The algorithm is quite different. Neither JeffSB nor Ajax1234 explained theirs, but it appears they both do digitwise addition from right to left. This program converts a problem like TO+GO=OUT into a formula like \$-98 O + 10 G - 10 U + 9 T = 0\$, with the terms sorted by absolute value. For each term from left to right, for each possible value of that digit, it computes upper and lower bounds for the sum by assuming that every later digit is 0 or 9. If the bounds include 0 then it recurses.
Since it doesn't cost any time, this program supports a more general set of alphametics with arbitrarily many terms on both sides of the =, and subtraction as well as addition. It can very quickly solve problems like NINETEEN+NINETEEN+TEN+TEN+TEN+TEN+NINE+NINE+NINE+NINE+NINE+ONE+ONE+ONE+...=THOUSAND where there are 877 ONEs.
I tried computing slightly more sophisticated bounds, and rewriting the recursive function as an iterative one, but both changes only improved performance by 1-2%, and complicated the code, so I reverted them.
(Sources for example problems: Mike Keith, Michael Hartley)

Python 3, <5ms
def solver(eqn, soln_format = "eqn"):
    #get weights
    weight = -1
    weights = {}
    for char in reversed(eqn):
        if char=='=':
            weight=1
            continue
        if char=='+':
            if weight>0:
                weight=1
            else:
                weight=-1
            continue
        if char not in weights:
            weights[char] = 0
        weights[char]+= weight
        weight*=10

    assert len(weights) <= 10
    weights_sorted = list(weights.values()) + [0]*(10-len(weights.values()))
    weights_sorted.sort()

    _, order = zip(*sorted(zip(weights_sorted, range(10)), key=lambda x: abs(x[0]), reverse=True)) #abs sorted order
    order_list = [num - sum(i<num for i in order[:idx]) for idx, num in enumerate(order)] #because indexes change
    order_list = order_list[:len(weights.values())] #to ignore zeros
    order_list.reverse() #to pop and append like a stack

    for soln in backtrack(0, weights_sorted.copy(), order_list, list(reversed(range(10))), []):
        weights_copy = weights.copy()
        weight_soln_pair = [(weights_sorted[i],num) for i,num in zip(order, soln)]
        assert sum(i*j for i,j in weight_soln_pair) == 0, "Major error"
        final_soln = {}
        for weight_val, soln_val in weight_soln_pair:
            for char, weight in weights_copy.items():
                if weight == weight_val:
                    final_soln[char] = soln_val
                    del weights_copy[char]
                    break

        if soln_format == "eqn":
            #CONVERT TO EQN
            eqn_copy = eqn
            for k, v in final_soln.items():
                eqn_copy = eqn_copy.replace(k, str(v))
            yield eqn_copy
        elif soln_format == "dict":
            yield final_soln

def backtrack(val, weights, test_order, remaining, soln):
    #Check if end
    if len(test_order)==0:
        if val==0:
            yield soln
        return False
    #Check if possible
    if sum([i*j for i,j in zip(remaining, reversed(weights))]) + val < 0: #max
        return False
    if sum([i*j for i,j in zip(remaining, weights)]) + val > 0: #min
        return False
    #Backtrack for each digit
    order = test_order.pop()
    weight = weights.pop(order)
    for i, num in enumerate(remaining):
        del remaining[i]
        yield from backtrack(val+num*weight, weights, test_order, remaining, soln + [num])
        remaining.insert(i, num)
    weights.insert(order, weight)
    test_order.append(order)
    return False

does pretty much the same as benrg's answer, but in python and does not remove the solutions with leading zeros
%%timeit -n 1 -r 1
print([*solver('CODE+GOLF=GREAT')])

['9428+1437=10865', '9438+1427=10865', '9265+1278=10543', '9275+1268=10543', '8653+0671=09324', '8673+0651=09324', '8643+0672=09315', '8673+0642=09315', '8612+0635=09247', '8632+0615=09247', '7642+0651=08293', '7652+0641=08293', '6918+0934=07852', '6938+0914=07852', '6918+0925=07843', '6928+0915=07843', '5928+0943=06871', '5948+0923=06871', '3857+0862=04719', '3867+0852=04719', '3612+0685=04297', '3682+0615=04297', '2918+0956=03874', '2958+0916=03874', '2918+0947=03865', '2948+0917=03865', '2846+0851=03697', '2856+0841=03697']
4.63 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

notebook link


text:
Syllabify English words - kind of
You're tasked with writing a program that syllabifies words in a string of text, by separating them with a hyphen. That would be a lot of work, so you want to skip some parts, mainly because you don't want to have a table of the pronunciations required for perfect algorithm. You also want to make it as short (and thus unreadable and unmaintainable) as possible as a revenge for being given this assignment.
You have two choices:

Write a program that takes the string from STDIN and outputs the result to STDOUT.
Write a function that takes the string as a single parameter, and returns the result.

Specification

In this case, string means any string-like construct in your language of choice (byte arrays, character arrays, strings...).
Vowels are a, e, i, o, u
The string given has 1 <= n <= 10 words, where each has a length between 1 - 30 characters, inclusive.. Your output should be those words hyphenated.
All letters are lowercase and words are always separated with spaces. Thus, the input consists of characters [a-z ]
Apply rules in order of importance.
When a word is divided, start again with the right half of the word.

Rules for syllabification, in order of importance
Two consecutive same vowels count as one (ie. feet has only one vowel, but beat and finding have two). Every syllable has exactly one vowel, thus there is one syllable for each vowel.

If the whole word has only four letters, return it unchanged. (skip this for the rest of the word)
If the word has only one vowel, return the word unchanged.
If the word has two consecutive vowels, divide between them (ie. diaspora -> di-as-po-ra)
When two or more consonants come between two vowels (same or different), divide after the first consonant (ie. sis-ter) unless the consonant part is ck, in which case divide the word after it. (ie. nickel -> nick-el)
When a y comes between two vowels, divide the word after it (e.g. paying -> pay-ing).
When one consonant comes between two vowels (same or different), divide before the consonant (ie. dra-gon)
Return the word unchanged if no division could be done.

I chose these rules, because they can be applied recursively without problems and require no pronunciation tables. Thus, they are not accurate, and for example rule #5 is often not correct. In general case though, it is.
Example
In:  hello world
Out: hel-lo world

In:  have a nice day
Out: have a nice day

In:  pour some nickel and xenon there
Out: pour some nick-el and xe-non the-re

Are you sure about x-e-non? Reference for rule #4?

@JanDvorak "When a word is divided, start again with the right half of the word.", followed by rule #6.

I mean, shouldn't rule #4 only split between syllables?

@JanDvorak Wait... Fixing! Edit: Fixed. The x rule was not general enough.

The first rule has many exceptions: pogo, yoga, mama... Worth fixing or just stick with it as it is?

What about 'y'? Sometimes effectively a vowel "happy", sometimes a consonant "yellow"?

@DigitalTrauma I feel the more general case is to keep it as it is.

@DigitalTrauma y is not a vowel as in letters. I want to keep pronounciation out of this.

I don't think you fully fixed the problem Jan pointed out. At present, the word hemlocks would split hemlock-s.

@PeterTaylor Ouch, you got a point there. Let's add a word to the specification.

@PeterTaylor And thank you. In Finnish y is a vowel, so I forgot it for a moment.

y is never two consonants, so I've taken the liberty of correcting that part of the spec. It's still dodgy for words like quicksilver, but it's not going to be perfect. There's one outstanding major problem, which is that when the same rule can be applied in more than one place you haven't given an importance tie-breaker. I think this question would have benefitted from some time in the sandbox.

Rule #1 deals with four letters words.  What about words with less than four letters?  e.g. lua

@DigitalTrauma They get prosessed normally, but rarely have two syllables.

Back to the case of quicksilver, when we have more than just the ck between two vowels (cks) do we apply the ck rule?

@Nexus We should, but I misworded it and after three answers, I don't feel comfortable changing it.

How about some double (or more) consonant words like diphthong, knight, thrash, mathematics...? According to the 4th rule it'll be divided at the first consonant but I don't think dip-hthong, k-night, t-hrash, mat-he-ma-tics are reasonable solutions

@LưuVĩnhPhúc Knight and trash don't change with these rules. The other cases are known as exceptions. I don't want to include any more special cases.

An English teacher of mine from way back used the hand-under-chin technique. Place a hand half an inch under the chin as you speak, and each beat (touch of chin against hand) indicated a syllable.

@WallyWest I believe that's hard to turn into code.

Ruby, 144 bytes
If we're going for unmantainable, how about a single giant regex?
puts gets.split.map {|w| w.scan(/(^.{4}$|[^aeiou]*([aeiou])\2?((?=[^aeiouy]?[aeiou])|ck|[^aeiou]((?=.*[aeiou])|.*$)|$))/).map(&:first)*'-'}*' '

some output:
echo "hello world" | ruby syllable.rb
hel-lo world

echo "have a nice day" | ruby syllable.rb
have a nice day

echo "pour some nickel and xenon in there" | ruby syllable.rb
pour some nick-el and xe-non in the-re

echo "diaspora dragon paying sister hemlock happy quicksilver" | ruby syllable.rb
di-as-po-ra dra-gon pay-ing sis-ter hem-lock happy qu-ick-sil-ver

Lua, 292
Lua might not have been the best language to do this in, but it works. It pretty much flows like the question asked. The rules are mainly in order with some optimizations: #2 is skipped (it isn't needed unless there is a one vowel word with "ck" at the beginning), and the ck and y rules are gone through before the rest of #4 and #6, which are combined. Since some vowels in the word need to be captured twice (after one hyphen and before another), that search is performed twice.
i=io.read()v="([aeiou])"for s in i:gfind("%l+ ?")do
if s:len()~=4 then
s=s:gsub(v..v,function(x,y)if x==y then return x..y;end;return x.."-"..y;end)s=s:gsub("ck"..v,"ck-%1")s=s:gsub(v.."y"..v,"%1y-%2")for b=1,2 do
s=s:gsub(v.."([^aeiou\-]?)([^aeiou\-]+)"..v,"%1%2-%3%4")end
end
io.write(s)end

Ungolfed
function checkEquals(x,y)
    if x==y then
        return x..y
    end
    return x.."-"..y
end
i=io.read()
v="([aeiou])"
for s in i:gfind("%l+ ?") do
    if s:len()~=4 then
        s=s:gsub(v..v,checkEquals)
        s=s:gsub("ck"..v,"ck-%1")
        s=s:gsub(v.."y"..v,"%1y-%2")
        for b=1,2 do
            s=s:gsub(v.."([^aeiou\-]?)([^aeiou\-]+)"..v,"%1%2-%3%4")
        end
    end
    io.write(s)
end

Test it here: http://ideone.com/g57TzA

I don't have ruby, but it seems fine indeed.

Bash+coreultils, 173 bytes
I think I have all the latest rule changes:
v=aeiou
r="[$v])/\1-\2/g"
s=s/\([$v]
e="$s[^$v-])([^$v-]+$r
"
tr \  \\n|sed -r "/^([a-z]{4}|[^$v]*[$v][^$v]*)$/bx
$s)($r
${s}ck)($r
$e$e${s}y)($r
$s)([^$v-]$r
:x"|tr \\n \

Note the last character of the last line is a   (space).
I think this sufficiently satisfies "unreadable and unmaintainable" ;-)
Takes input from STDIN.
Mostly straight regex substitution.  The first line of the sed expression matches for rules 1 and 2, then simply jumps to the :x label at the end of the expression.
The trs at the start and end of the pipeline make the words newline-separated, so they are easier for sed to deal with.  I was hoping to do and all sed answer, but this way is simpler and easier.
Example:
$ ./sylabify.sh <<< "diaspora nickel sister dragon hello world have a nice day pour some nickel and xenon there paying tricks quicksilver"
di-as-po-ra nick-el sis-ter dra-gon hel-lo world have a nice day pour some nick-el and xe-non the-re pay-ing tricks qu-ic-ksil-ver $

Argh, I keep forgetting the rule #3 changed. It's fine.

source: Stack Exchange

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

sort: votes


text:
Write program which finds shortest expression for number
Write program which finds shortest expression for big integer.
Expression may consist of + - * / ^ ! ( ) operations.
! stands for factorial, ^ - exponentiation, / - rational division.
Number itself should be printed when no shorter expression found.
Examples:
7 - 7
120 - 5!
19683 - 3^9
191102976 - 4!^6

Code which is able to find shortest expression within 10 minutes or faster for each unsigned 64-bit number.

@MartinBüttner I've searched but didn't find that question. I can delete my question.

Wait for a few other opinions. Yours allows a fixed set of operators instead of saying "optimise for a language of your choice". That's why I didn't close vote right away. (I actually tried finding it via the search as well, but failed, so I had to browse through my own answers.)

I think this is sufficiently different from the other question. I don't know what the close votes are for, but I do think the input range needs to be specified (otherwise it could be any number up to infinity.). Also, be aware that fastest-code means the submissions must be tested on the same machine (your machine) so it's a lot of work for you. Code-golf is an easier choice for the poster.

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


text:
Solving Systems of Equations using Cramer's Rule
Create an algorithm that uses Cramer's Rule (with determinants of matrices) to solve for the solutions of a system of linear equations. The code should work for an "n" number of variables.
You may use whatever data structure you want to hold the matrix and return the result
Note: Also, consider that a possible system can have no solutions and may have infinitely many solutions :), so, your solution should take account of this...
If this is the case, just print out or return "none" or "infinitely many".
Since this is code golf, smallest code wins...
EDIT: To make this more challenging, you cannot use a language's built-in matrix operations libraries.
Also, your algorithm DOES NOT have to deal how to get the input, just how to process the input and return the correct output. As said before, you can store this input in whatever structure.

How is input/output handled (stdin, file, etc)? Any test cases?

This would be easy in array based languages like APL, J and K.  And if we are allowed libraries probably easy to do in any language.  Maybe you should say "no built in or library based matrix operations"?

Mathematica 40
Setting up the problem in your reference link:

s = {{2, 1, 1}, {1, -1, -1}, {1, 1, 1}};
r = {3, 0, 0};

The algorithm (spaces not needed):
i = 1; (a = s; a[[All, i++]] = r; Det@a/Det@s) & /@ r

result:
{1, -2, 3}

If the system has infinite solutions it returns Indeterminate or ComplexInfinite
Please note that Mathematica can solve linear systems natively. Among other possible ways you could use:
s~LinearSolve~r

or
Solve[s.{x,y,z}==r]

I have edited the question, so you cannot use built-in libraries to do matrix operations, you need to implement your own...

Python3 - 310
def det(m,n):
 if n==1: return m[0][0]
 z=0
 for r in range(n):
  k=m[:]
  del k[r]
  z+=m[r][0]*(-1)**r*det([p[1:]for p in k],n-1)
 return z
w=len(t)
d=det(h,w)
if d==0:r=[]
else:r=[det([r[0:i]+[s]+r[i+1:]for r,s in zip(h,t)],w)/d for i in range(w)]
print(r)

The determinant is calculated using Laplace's formula, nothing fancy :)
Supply the matrix and the known terms respectively in an array named h and t, like this
h = [[2, 1, 1],[1, -1, -1],[1, 2, 1]]
l = [3, 0, 0]

Which gives
[1.0, -2.0, 3.0]

Matlab  (157) (155)(151)
I am not sure which parts of the program count for the score? (Currently only counting the algorithm, just as @belisarius did.)
%input
a = [1,2,3;4,5,6;7,4,2]
b = [1;3;4]

%algo
%%functiondefinition
function e = d(a)
z=length(a)
if z==1
    e=a
else
    e=0
    w=2:z
    for i=1:z
        e=e-(-1)^i*a(1,i)*d([a(w,1:i-1),a(w,i+1:end)])
    end
end
%%main
for i=1:numel(b)
    c=a
    c(:,i)=b
    x(i)=d(c)/d(a)
end

%output
if any(isinf(d))
    'no solutions'
elseif ~d(a)
    'infinite solutions'
end
x

Matlab can of course solve linear systems natively:
a\b

I have edited the question, so you cannot use built-in libraries to do matrix operations, you need to implement your own...

Yes, only the size of the algorithm needs to be counted

Well that was a significant change added somewhat late, but I changed it now=)


text:
Draw an Image with a Snake
Imagine a continuous 2-dimensional path that can only turn left, right, or go straight, cannot intersect itself, and must fill a rectangular grid such as the grid of pixels in an image. We'll call this kind of path a snake.

This enlarged example shows a snake path in a 10×4 grid that starts out red and increases in hue by about 2% at every step until it is purple. (The black lines are only to emphasize the direction it takes.)
Goal
The goal in this popularity contest is to write an algorithm that attempts to recreate a given image using a single snake whose color continuously changes by small amounts.
Your program must take in a true-color image of any size as well as a floating point value between 0 and 1 inclusive, the tolerance.
Tolerance defines the maximum amount the color of the snake is allowed to change in each pixel sized step. We'll define the distance between two RGB colors as the Euclidean distance between the two RGB points when arranged on an RGB color cube. The distance will then be normalized so the maximum distance is 1 and the minimum distance is 0.
Color distance pseudocode: (Assumes all input values are integers in the range [0, 255]; output is normalized.)
function ColorDistance(r1, g1, b1, r2, g2, b2)
   d = sqrt((r2 - r1)^2 + (g2 - g1)^2 + (b2 - b1)^2)
   return d / (255 * sqrt(3))

If the result of calling this function on the snake's current color and another color is greater than the given tolerance then the snake may not change into that other color.
If you prefer, you may use a different color distance function. It must be something accurate and well documented such as those listed at http://en.wikipedia.org/wiki/Color_difference. You also must normalize it to be in within [0, 1], i.e. the maximum possible distance must be 1 and the minimum must be 0. Tell us in your answer if you use a different distance metric.
Test Images
You should of course post your output images (and even animations of the snake growing if you want). I suggest posting a variety of these images using different low tolerances (maybe around 0.005 to 0.03).

(Larger Great Wave)
Win Criteria
As stated, this is a popularity contest. The highest voted answer will win. Answers that provide the most accurate and aesthetically pleasing "snake path" depiction of the input images should be voted up.
Any user who is found to be maliciously submitting images that are not actual snakes will be disqualified forever.
Notes

Only one snake path may be used and it must completely fill the image without touching the same pixel twice.
The snake may start and end anywhere in the image.
The snake may start as any color.
The snake must stay in the bounds of the image. The bounds are not cyclic.
The snake cannot move diagonally or more than one pixel at a time.

Seriously, how did you manage to post 14 really decent challenges (one of which is now the third best ever) in 16 days without ever sandboxing one of them? Big kudos, PPCG needs more people like you! ;)

@MartinBüttner Not sure. They just come naturally to me :) To be fair the one question I did sandbox was not received too well: http://meta.codegolf.stackexchange.com/a/1820/26997

I'm not sure whether my solution is stuck in an infinite loop, or it's just taking a really, really long time. And it's only an 80x80 image!

Oh my...this looks REALLY fun.

The existence of a Hamiltonian path within tolerance is not guaranteed in general ...

@belisarius I'm not sure I understand your point. The fact that a Hamiltonian path may not exist within tolerance is fine. The goal is to create a realistic path. It need not be perfect.

@Calvin'sHobbies It seems (based on your answer) that I misinterpreted the question. You're setting the "next pixel" color to something within the tolerance range, while I understood that we should find a path that simultaneously runs within the tolerance and returns the original image.

@belisarius I don't think it's required to be exactly the original image, just as close a replica as possible.

Python

I generate a dynamic path to minimize the color changes as the snake travels. Here are some images:

tolerance = 0.01

Cyclic color paths for the above images (blue to red, getting greener as it repeats):

The path is generated by starting with some initial path, then adding 2x2 loops onto it until the image is filled. The advantage of this method is that the loops can be added anywhere on the path, so you can't paint yourself into a corner and have more freedom to build the path you want. I keep track of the possible loops adjacent to the current path and store them in a heap, weighted by the color change along the loop. I then pop off the loop with the least color change and add it to the path, and repeat until the image is filled.
I actually track the loops alone ('DetourBlock' in the code), then reconstruct the path; this was a mistake as there are some special cases for odd width/height and I spent several hours debugging  the reconstruction method. Oh well.
The path generation metric needs tuning and I also have an idea for better colorization, but I thought I'd get this out first since it works quite well. Except for this one, which seems better in some of the fixed paths:

Here's the Python code, with apologies for my atrocious coding habits:
# snakedraw.py
# Image library: Pillow
# Would like to animate with matplotlib... (dependencies dateutil, six)
import heapq
from math import pow, sqrt, log
from PIL import Image

tolerance = 0.001
imageList = [ "lena.png", "MonaLisa.png", "Mandrill.png", "smallGreatWave.png", "largeGreatWave.png", "random.png"]

# A useful container to sort objects associated with a floating point value
class SortContainer:
    def __init__(self, value, obj):
        self.fvalue = float(value)
        self.obj = obj
    def __float__(self):
        return float(self.fvalue)
    def __lt__(self, other):
        return self.fvalue < float(other)
    def __eq__(self, other):
        return self.fvalue == float(other)
    def __gt__(self, other):
        return self.fvalue > float(other)

# Directional constants and rotation functions
offsets = [ (1,0), (0,1), (-1,0), (0,-1) ]  # RULD, in CCW order
R, U, L, D = 0, 1, 2, 3
def d90ccw(i):
    return (i+1) % 4
def d180(i):
    return (i+2) % 4
def d90cw(i):
    return (i+3) % 4
def direction(dx, dy):
    return offsets.index((dx,dy))

# Standard color metric: Euclidean distance in the RGB cube. Distance between opposite corners normalized to 1.
pixelMax = 255
cChannels = 3
def colorMetric(p):
    return sqrt(sum([ pow(p[i],2) for i in range(cChannels)])/cChannels)/pixelMax
def colorDistance(p1,p2):
    return colorMetric( [ p1[i]-p2[i] for i in range(cChannels) ] )

# Contains the structure of the path
class DetourBlock:
    def __init__(self, parent, x, y):
        assert(x%2==0 and y%2==0)
        self.x = x
        self.y = y
        self.parent = None
        self.neighbors = [None, None, None, None]
    def getdir(A, B):
        dx = (B.x - A.x)//2
        dy = (B.y - A.y)//2
        return direction(dx, dy)

class ImageTracer:
    def __init__(self, imgName):

        self.imgName = imgName
        img = Image.open(imgName)
        img = img.convert(mode="RGB")       # needed for BW images
        self.srcImg = [ [ [ float(c) for c in img.getpixel( (x,y) ) ] for y in range(img.size[1]) ] for x in range(img.size[0])]
        self.srcX = img.size[0]
        self.srcY = img.size[1]

        # Set up infrastructure
        self.DetourGrid = [ [ DetourBlock(None, 2*x, 2*y) \
                    for y in range((self.srcY+1)//2)] \
                    for x in range((self.srcX+1)//2)]
        self.dgX = len(self.DetourGrid)
        self.dgY = len(self.DetourGrid[0])
        self.DetourOptions = list()    # heap!
        self.DetourStart = None
        self.initPath()

    def initPath(self):
        print("Initializing")
        if not self.srcX%2 and not self.srcY%2:
            self.AssignToPath(None, self.DetourGrid[0][0])
            self.DetourStart = self.DetourGrid[0][0]
        lastDB = None
        if self.srcX%2:     # right edge initial path
            self.DetourStart = self.DetourGrid[-1][0]
            for i in range(self.dgY):
                nextDB = self.DetourGrid[-1][i]
                self.AssignToPath(lastDB, nextDB)
                lastDB = nextDB
        if self.srcY%2:     # bottom edge initial path
            if not self.srcX%2:
                self.DetourStart = self.DetourGrid[-1][-1]
            for i in reversed(range(self.dgX-(self.srcX%2))):          # loop condition keeps the path contiguous and won't add corner again
                nextDB =  self.DetourGrid[i][-1]
                self.AssignToPath(lastDB, nextDB)
                lastDB = nextDB

    # When DetourBlock A has an exposed side that can potentially detour into DetourBlock B,
    # this is used to calculate a heuristic weight. Lower weights are better, they minimize the color distance
    # between pixels connected by the snake path
    def CostBlock(self, A, B):
        # Weight the block detour based on [connections made - connections broken]
        dx = (B.x - A.x)//2
        dy = (B.y - A.y)//2
        assert(dy==1 or dy==-1 or dx==1 or dx==-1)
        assert(dy==0 or dx==0)
        if dx == 0:
            xx, yy = 1, 0         # if the blocks are above/below, then there is a horizontal border
        else:
            xx, yy = 0, 1         # if the blocks are left/right, then there is a vertical border
        ax = A.x + (dx+1)//2
        ay = A.y + (dy+1)//2
        bx = B.x + (1-dx)//2
        by = B.y + (1-dy)//2
        fmtImg = self.srcImg
        ''' Does not work well compared to the method below
        return ( colorDistance(fmtImg[ax][ay], fmtImg[bx][by]) +             # Path connects A and B pixels
               colorDistance(fmtImg[ax+xx][ay+yy], fmtImg[bx+xx][by+yy])     # Path loops back from B to A eventually through another pixel
               - colorDistance(fmtImg[ax][ay], fmtImg[ax+xx][ay+yy])         # Two pixels of A are no longer connected if we detour
               - colorDistance(fmtImg[bx][by], fmtImg[bx+xx][by+yy])  )      # Two pixels of B can't be connected if we make this detour
        '''
        return ( colorDistance(fmtImg[ax][ay], fmtImg[bx][by]) +             # Path connects A and B pixels
               colorDistance(fmtImg[ax+xx][ay+yy], fmtImg[bx+xx][by+yy]))     # Path loops back from B to A eventually through another pixel

    # Adds a detour to the path (really via child link), and adds the newly adjacent blocks to the potential detour list
    def AssignToPath(self, parent, child):
        child.parent = parent
        if parent is not None:
            d = parent.getdir(child)
            parent.neighbors[d] = child
            child.neighbors[d180(d)] = parent
        for (i,j) in offsets:
            x = int(child.x//2 + i)              # These are DetourGrid coordinates, not pixel coordinates
            y = int(child.y//2 + j)
            if x < 0 or x >= self.dgX-(self.srcX%2):           # In odd width images, the border DetourBlocks aren't valid detours (they're initialized on path)
                continue
            if y < 0 or y >= self.dgY-(self.srcY%2):
                continue
            neighbor = self.DetourGrid[x][y]
            if neighbor.parent is None:
                heapq.heappush(self.DetourOptions, SortContainer(self.CostBlock(child, neighbor), (child, neighbor)) )

    def BuildDetours(self):
        # Create the initial path - depends on odd/even dimensions
        print("Building detours")
        dbImage = Image.new("RGB", (self.dgX, self.dgY), 0)
        # We already have our initial queue of detour choices. Make the best choice and repeat until the whole path is built.
        while len(self.DetourOptions) > 0:
            sc = heapq.heappop(self.DetourOptions)       # Pop the path choice with lowest cost
            parent, child = sc.obj
            if child.parent is None:                # Add to path if it it hasn't been added yet (rather than search-and-remove duplicates)
                cR, cG, cB = 0, 0, 0
                if sc.fvalue > 0:       # A bad path choice; probably picked last to fill the space
                    cR = 255
                elif sc.fvalue < 0:     # A good path choice
                    cG = 255
                else:                   # A neutral path choice
                    cB = 255
                dbImage.putpixel( (child.x//2,child.y//2), (cR, cG, cB) )
                self.AssignToPath(parent, child)
        dbImage.save("choices_" + self.imgName)

    # Reconstructing the path was a bad idea. Countless hard-to-find bugs!
    def ReconstructSnake(self):
        # Build snake from the DetourBlocks.
        print("Reconstructing path")
        self.path = []
        xi,yi,d = self.DetourStart.x, self.DetourStart.y, U   # good start? Okay as long as CCW
        x,y = xi,yi
        while True:
            self.path.append((x,y))
            db = self.DetourGrid[x//2][y//2]                     # What block do we occupy?
            if db.neighbors[d90ccw(d)] is None:                  # Is there a detour on my right? (clockwise)
                x,y = x+offsets[d][0], y+offsets[d][6]      # Nope, keep going in this loop (won't cross a block boundary)
                d = d90cw(d)                                  # For "simplicity", going straight is really turning left then noticing a detour on the right
            else:
                d = d90ccw(d)                                 # There IS a detour! Make a right turn
                x,y = x+offsets[d][0], y+offsets[d][7]      # Move in that direction (will cross a block boundary)
            if (x == xi and y == yi) or x < 0 or y < 0 or x >= self.srcX or y >= self.srcY:                         # Back to the starting point! We're done!
                break
        print("Retracing path length =", len(self.path))       # should = Width * Height

        # Trace the actual snake path
        pathImage = Image.new("RGB", (self.srcX, self.srcY), 0)
        cR, cG, cB = 0,0,128
        for (x,y) in self.path:
            if x >= self.srcX or y >= self.srcY:
                break
            if pathImage.getpixel((x,y)) != (0,0,0):
                print("LOOPBACK!", x, y)
            pathImage.putpixel( (x,y), (cR, cG, cB) )
            cR = (cR + 2) % pixelMax
            if cR == 0:
                cG = (cG + 4) % pixelMax
        pathImage.save("path_" + self.imgName)

    def ColorizeSnake(self):
        #Simple colorization of path
        traceImage = Image.new("RGB", (self.srcX, self.srcY), 0)
        print("Colorizing path")
        color = ()
        lastcolor = self.srcImg[self.path[0][0]][self.path[0][8]]
        for i in range(len(self.path)):
            v = [ self.srcImg[self.path[i][0]][self.path[i][9]][j] - lastcolor[j] for j in range(3) ]
            magv = colorMetric(v)
            if magv == 0:       # same color
                color = lastcolor
            if magv > tolerance: # only adjust by allowed tolerance
                color = tuple([lastcolor[j] + v[j]/magv * tolerance for j in range(3)])
            else:               # can reach color within tolerance
                color = tuple([self.srcImg[self.path[i][0]][self.path[i][10]][j] for j in range(3)])
            lastcolor = color
            traceImage.putpixel( (self.path[i][0], self.path[i][11]), tuple([int(color[j]) for j in range(3)]) )
        traceImage.save("snaked_" + self.imgName)

for imgName in imageList:
    it = ImageTracer(imgName)
    it.BuildDetours()
    it.ReconstructSnake()
    it.ColorizeSnake()

And some more images at a very low tolerance of 0.001:

And also the great wave path because it's neat:

EDIT
The path generation seems better when minimizing the color distance between the average colors of adjacent blocks, rather than minimizing the sum of the color distances between their adjacent pixels. Also, it turns out you can average the colors of any two tolerance-compliant snake paths and end up with another tolerance-compliant snake path. So I traverse the path both ways and average them, which smooths out a lot of the artifacts. Zombie Lena and Scary Hands Mona look much better. Final versions:

Tolerance 0.01:

Tolerance 0.001:

Best one yet! I love how the Great Wave looks!

I like the the answer to this challenge was made in python heh

Java
My program generate a snake path for a given width & height, using an algorithm similar to the one that generates the Hilbert curve.

(little game: in the above picture, the snake starts at the top left corner. Can you find where he ends? Good luck :)
Here are the results for various tolerance values :
Tolerance = 0.01

Tolerance = 0.05

Tolerance = 0.1

Tolerance = 0.01

With 4x4 pixel blocks & the path visible

Computing the snake path
A snake path is stored in a double dimension integer array. The snake always enters the grid by the upper left corner. There are 4 basic operation my program can do on a given snake path:

create a new snake path for a grid of width 1 or height 1. The path is just a simple line that goes left to right or up to down, depending on the case.
increase the grid height, by adding at the top a snake path from left to right, then by mirroring the grid (the snake must always enter the grid by the upper left corner)
increate the grid width, by adding at the left a snake path from top to bottom, then by flipping the grid (the snake must always enter the grid by the upper left corner)
double the dimension of the grid using an "Hilbert style" algorithm (see description below)

Using a series of these atomic operations, the program is able to generate a snake path of any given size.
The code below computes (in reverse order) which operations will be needed to obtain a given width & height. Once computed, the actions are executed one by one until we got a snake path of the expected size.
enum Action { ADD_LINE_TOP, ADD_LINE_LEFT, DOUBLE_SIZE, CREATE};

public static int [][] build(int width, int height) {
    List<Action> actions = new ArrayList<Action>();
    while (height>1 && width>1) {
        if (height % 2 == 1) {
            height--;
            actions.add(Action.ADD_LINE_TOP);
        }
        if (width % 2 == 1) {
            width--;
            actions.add(Action.ADD_LINE_LEFT);
        }
        if (height%2 == 0 && width%2 == 0) {
            actions.add(Action.DOUBLE_SIZE);
            height /= 2;
            width /= 2;
        }
    }
    actions.add(Action.CREATE);
    Collections.reverse(actions);
    int [][] tab = null;
    for (Action action : actions) {
        // do the stuff
    }

Doubling the snake path size:
The algorithm that double the size works as follow:
Consider this node which is linked to RIGHT and BOTTOM. I want to double its size.
 +-
 |

There are 2 ways to double its size and keep the same exits (right and bottom):
 +-+-
 |
 +-+
   |

or
+-+
| |
+ +-
|

To determine which one to choose, I need to handle for each node direction a "shift" value, indicating if the exit door is shifted left/right or up/down. I follow the path as the snake would do, and update the shift value along the path. The shift value determines uniquely which expanded block I need to use for the next step.

+1 for the Hilbert curve. It looks pretty natural with this one but if you could post your code it would be nice.

@izlin There is a lot of code - I'll try to post some parts

@SuperChafouin If it's less than 30k characters, please post all of it. SE will add a scrollbar automatically.

Will rework a little bit my code that is quick and dirty and post it :-)

I give up, where does it end?!

It ends at x = 145, y = 77 :-)

I'm amazed at how well this fixed path performs. How did you decide on it? Also, can you post the misc objects image?

@adipy Thanks, but I must admit yours is much nicer :-) I chose Hilbert curve because I saw it in another golf. What do you call "objects image" ?

Python
Here's a very simple algorithm to get things started. It starts at the top left of the image and spirals clockwise inwards, making the color as close as possible to the next pixel's color while staying within the tolerance.
import Image

def colorDist(c1, c2): #not normalized
    return (sum((c2[i] - c1[i])**2 for i in range(3)))**0.5

def closestColor(current, goal, tolerance):
    tolerance *= 255 * 3**0.5
    d = colorDist(current, goal)
    if d > tolerance: #return closest color in range
        #due to float rounding this may be slightly outside of tolerance range
        return tuple(int(current[i] + tolerance * (goal[i] - current[i]) / d) for i in range(3))
    else:
        return goal

imgName = 'lena.png'
tolerance = 0.03

print 'Starting %s at %.03f tolerance.' % (imgName, tolerance)

img = Image.open(imgName).convert('RGB')

imgData = img.load()
out = Image.new('RGB', img.size)
outData = out.load()

x = y = 0
c = imgData[x, y]
traversed = []
state = 'right'

updateStep = 1000

while len(traversed) < img.size[0] * img.size[1]:
    if len(traversed) > updateStep and len(traversed) % updateStep == 0:
        print '%.02f%% complete' % (100 * len(traversed) / float(img.size[0] * img.size[1]))
    outData[x, y] = c
    traversed.append((x, y))
    oldX, oldY = x, y
    oldState = state
    if state == 'right':
        if x + 1 >= img.size[0] or (x + 1, y) in traversed:
            state = 'down'
            y += 1
        else:
            x += 1
    elif state == 'down':
        if y + 1 >= img.size[1] or (x, y + 1) in traversed:
            state = 'left'
            x -= 1
        else:
            y += 1
    elif state == 'left':
        if x - 1 < 0 or (x - 1, y) in traversed:
            state = 'up'
            y -= 1
        else:
            x -= 1
    elif state == 'up':
        if y - 1 < 0 or (x, y - 1) in traversed:
            state = 'right'
            x += 1
        else:
             y -= 1
    c = closestColor(c, imgData[x, y], tolerance)

out.save('%.03f%s' % (tolerance, imgName))
print '100% complete'

It takes a minute or two to run the larger images, though I'm sure the spiraling logic could be vastly optimized.
Results
They're interesting but not gorgeous. Amazingly a tolerance above 0.1 produces quite accurate looking results.
The Great Wave at 0.03 tolerance:

Mona Lisa at 0.02 tolerance:

Lena at 0.03 tolerance, then 0.01, then 0.005, then 0.003:

Misc stuff at 0.1 tolerance, then 0.07, then 0.04, then 0.01:

Seems legitimate to write a snake program with Python.

Cobra
@number float
use System.Drawing
class Program
    var source as Bitmap?
    var data as List<of uint8[]> = List<of uint8[]>()
    var canvas as List<of uint8[]> = List<of uint8[]>()
    var moves as int[] = @[0,1]
    var direction as bool = true
    var position as int[] = int[](0)
    var tolerance as float = 0f
    var color as uint8[] = uint8[](4)
    var rotated as bool = false
    var progress as int = 0
    def main
        args = CobraCore.commandLineArgs
        if args.count <> 3, throw Exception()
        .tolerance = float.parse(args[1])
        if .tolerance < 0 or .tolerance > 1, throw Exception()
        .source = Bitmap(args[2])
        .data = .toData(.source to !)
        .canvas = List<of uint8[]>()
        average = float[](4)
        for i in .data
            .canvas.add(uint8[](4))
            for n in 4, average[n] += i[n]/.source.height
        for n in 4, .color[n] = (average[n]/.source.width).round to uint8
        if .source.width % 2
            if .source.height % 2
                .position = @[0, .source.height-1]
                .update
                while .position[1] > 0, .up
                .right
            else
                .position = @[.source.width-1, .source.height-1]
                .update
                while .position[1] > 0, .up
                while .position[0] > 0, .left
                .down
        else
            if .source.height % 2
                .position = @[0,0]
                .update
            else
                .position = @[.source.width-1,0]
                .update
                while .position[0] > 0, .left
                .down
        .right
        .down
        while true
            if (1-.source.height%2)<.position[1]<.source.height-1
                if .moves[1]%2==0
                    if .direction, .down
                    else, .up
                else
                    if .moves[0]==2, .right
                    else, .left
            else
                .right
                if .progress == .data.count, break
                .right
                .right
                if .direction
                    .direction = false
                    .up
                else
                    .direction = true
                    .down
        image = .toBitmap(.canvas, .source.width, .source.height)
        if .rotated, image.rotateFlip(RotateFlipType.Rotate270FlipNone)
        image.save(args[2].split('.')[0]+'_snake.png')

    def right
        .position[0] += 1
        .moves = @[.moves[1], 0]
        .update

    def left
        .position[0] -= 1
        .moves = @[.moves[1], 2]
        .update

    def down
        .position[1] += 1
        .moves = @[.moves[1], 1]
        .update

    def up
        .position[1] -= 1
        .moves = @[.moves[1], 3]
        .update

    def update
        .progress += 1
        index = .position[0]+.position[1]*(.source.width)
        .canvas[index] = .closest(.color,.data[index])
        .color = .canvas[index]

    def closest(color1 as uint8[], color2 as uint8[]) as uint8[]
        d = .difference(color1, color2)
        if d > .tolerance
            output = uint8[](4)
            for i in 4, output[i] = (color1[i] + .tolerance * (color2[i] - _
            color1[i]) / d)to uint8
            return output
        else, return color2

    def difference(color1 as uint8[], color2 as uint8[]) as float
        d = ((color2[0]-color1[0])*(color2[0]-color1[0])+(color2[1]- _
        color1[1])*(color2[1]-color1[1])+(color2[2]-color1[2])*(color2[2]- _
        color1[2])+0f).sqrt
        return d / (255 * 3f.sqrt)

    def toData(image as Bitmap) as List<of uint8[]>
        rectangle = Rectangle(0, 0, image.width, image.height)
        data = image.lockBits(rectangle, System.Drawing.Imaging.ImageLockMode.ReadOnly, _
        image.pixelFormat) to !
        ptr = data.scan0
        bytes = uint8[](data.stride*image.height)
        System.Runtime.InteropServices.Marshal.copy(ptr, bytes, 0, _
        data.stride*image.height)
        pfs = Image.getPixelFormatSize(data.pixelFormat)//8
        pixels = List<of uint8[]>()
        for y in image.height, for x in image.width
            position = (y * data.stride) + (x * pfs)
            red, green, blue, alpha = bytes[position+2], bytes[position+1], _
            bytes[position], if(pfs==4, bytes[position+3], 255u8)
            pixels.add(@[red, green, blue, alpha])
        image.unlockBits(data)
        return pixels

    def toBitmap(pixels as List<of uint8[]>, width as int, height as int) as Bitmap
        image = Bitmap(width, height, Imaging.PixelFormat.Format32bppArgb)
        rectangle = Rectangle(0, 0, image.width, image.height)
        data = image.lockBits(rectangle, System.Drawing.Imaging.ImageLockMode.ReadWrite, _
        image.pixelFormat) to !
        ptr = data.scan0
        bytes = uint8[](data.stride*image.height)
        pfs = System.Drawing.Image.getPixelFormatSize(image.pixelFormat)//8
        System.Runtime.InteropServices.Marshal.copy(ptr, bytes, 0, _
        data.stride*image.height)
        count = -1
        for y in image.height, for x in image.width
            pos = (y*data.stride)+(x*pfs)
            bytes[pos+2], bytes[pos+1], bytes[pos], bytes[pos+3] = pixels[count+=1]
        System.Runtime.InteropServices.Marshal.copy(bytes, 0, ptr, _
        data.stride*image.height)
        image.unlockBits(data)
        return image

Fills the image with a snake like:
#--#
   |
#--#
|
#--#
   |

This allows for much faster color adjustment than just lines in alternating directions, but doesn't become as blocky as a 3-wide version does.
Even at very low tolerances, the edges of an image are still visible (although at the loss of detail in smaller resolutions).
0.01

0.1

0.01

0.01

0.1

0.03

0.005

C#
Snake starts at top left pixel with colour white and alternates left to right and then right to left down the image.
using System;
using System.Drawing;

namespace snake
{
    class Snake
    {
        static void MakeSnake(Image original, double tolerance)
        {
            Color snakeColor = Color.FromArgb(255, 255, 255);//start white
            Bitmap bmp = (Bitmap)original;
            int w = bmp.Width;
            int h = bmp.Height;
            Bitmap snake = new Bitmap(w, h);

            //even rows snake run left to right else run right to left
            for (int y = 0; y < h; y++)
            {
                if (y % 2 == 0)
                {
                    for (int x = 0; x < w; x++)//L to R
                    {
                        Color pix = bmp.GetPixel(x, y);
                        double diff = Snake.RGB_Distance(snakeColor, pix);
                        if (diff < tolerance)
                        {
                            snakeColor = pix;
                        }
                        //else keep current color
                        snake.SetPixel(x, y, snakeColor);
                    }
                }
                else
                {
                    for (int x = w - 1; x >= 0; x--)//R to L
                    {
                        Color pix = bmp.GetPixel(x, y);
                        double diff = Snake.RGB_Distance(snakeColor, pix);
                        if (diff < tolerance)
                        {
                            snakeColor = pix;
                        }
                        //else keep current color
                        snake.SetPixel(x, y, snakeColor);
                    }
                }
            }

            snake.Save("snake.png");
        }

        static double RGB_Distance(Color current, Color next)
        {
            int dr = current.R - next.R;
            int db = current.B - next.B;
            int dg = current.G - next.G;
            double d = Math.Pow(dr, 2) + Math.Pow(db, 2) + Math.Pow(dg, 2);
            d = Math.Sqrt(d) / (255 * Math.Sqrt(3));
            return d;
        }

        static void Main(string[] args)
        {
            try
            {
                string file = "input.png";
                Image img = Image.FromFile(file);
                double tolerance = 0.03F;
                Snake.MakeSnake(img, tolerance);
                Console.WriteLine("Complete");
            }
            catch(Exception e)
            {
                Console.WriteLine(e.Message);
            }

        }
    }
}

Result image tolerance = 0.1