text:
Four Squares Together
Lagrange's four square theorem tells us any natural number can be represented as the sum of four square numbers. Your task is to write a program that does this.
Input: A natural number (below 1 billion)
Output: Four numbers whose squares sum to that number (order doesn't matter)
Note: You don't have to do a brute force search! Details here and here. If there  is a function that trivializes this problem (I will determine), it's not allowed. Automated prime functions and square root are allowed. If there is more than one representation any one is fine. If you chose to do brute force, it must run within reasonable time (3 minutes)
sample input
123456789

sample output (either is fine)
10601 3328 2 0
10601 3328 2

I may do brute force though if it makes my code shorter?

@m.buettner Yes, but it should handle large numbers

@m.buettner Read the post, any natural number below 1 billion

Ah sorry overlooked that.

Does "natural numbers" include 0?

@Dennis it certainly looks like there's no agreement on this in the maths world, so I gave myself the benefit of the doubt.  http://en.wikipedia.org/wiki/Natural_number Many programs including mine will need a patch to work with n=0. The first issue is division by 4 (you can get into an infinite loop trying to divide zero by 4 until you reach a number not divisible by 4) and more problems occur later.

@Dennis Natural numbers in this case do not include 0.

This theorem is impressive.

The natural numbers include 0. In maths.

@NicolasBarbulesco There is no agreement on whether 0 should be included, but in this case I have excluded it to make golfing easier

If my language of choice only supports numbers up to 9999 but can theoretically work for numbers larger, does it count?

CJam, 50 bytes
li:NmF0=~2/#:J(_{;)__*N\-[{_mqJ/iJ*__*@\-}3*])}g+p

My third (and final, I promise) answer. This approach is based heavily on primo's answer.
Try it online in the CJam interpreter.
Usage
$ cjam 4squares.cjam <<< 999999999
[189 31617 567 90]

Background

After seeing primo's updated algorithm, I had to see how a CJam implementation would score:
li{W):W;:N4md!}g;Nmqi)_{;(__*N\-[{_mqi__*@\-}3*])}g+2W#f*p

Only 58 bytes! This algorithm is performs in nearly constant time and doesn't exhibit much variation for different values of N. Let's change that...
Instead of starting at floor(sqrt(N)) and decrementing, we can start at 1 and increment. This saves 4 bytes.
li{W):W;:N4md!}g;0_{;)__*N\-[{_mqi__*@\-}3*])}g+2W#f*p

Instead of expressing N as 4**a * b, we can express it as p**(2a) * b – where p is the smallest prime factor of N – to save 1 more byte.
li_mF0=~2/#:J_*/:N!_{;)__*N\-[{_mqi__*@\-}3*])}g+Jf*p

The previous modification allows us to slightly change the implementation (without touching the algorithm itself): Instead of dividing N by p**(2a) and multiply the solution by p**a, we can directly restrict the possible solutions to multiples of p**a. This saves 2 more bytes.
li:NmF0=~2/#:J!_{;J+__*N\-[{_mqJ/iJ*__*@\-}3*])}g+`

Not restricting the first integer to multiples of p**a saves an additional byte.
li:NmF0=~2/#:J(_{;)__*N\-[{_mqJ/iJ*__*@\-}3*])}g+`

Final algorithm

Find a and b such that N = p**(2a) * b, where b is not a multiple of p**2 and p is the smallest prime factor of N.
Set j = p**a.
Set k = floor(sqrt(N - j**2) / A) * A.
Set l = floor(sqrt(N - j**2 - k**2) / A) * A.
Set m = floor(sqrt(N - j**2 - k**2 - l**2) / A) * A.
If N - j**2 - k**2 - l**2 - m**2 > 0, set j = j + 1 and go back to step 3.

This can be implemented as follows:
li:N          " Read an integer from STDIN and save it in “N”.                        ";
mF            " Push the factorization of “N”. Result: [ [ p1 a1 ] ... [ pn an ] ]    ";
0=~           " Push “p1” and “a1”. “p1” is the smallest prime divisor of “N”.        ";
2/#:J         " Compute p1**(a1/2) and save the result “J”.                           ";
(_            " Undo the first two instructions of the loop.                          ";
{             "                                                                       ";
  ;)_         " Pop and discard. Increment “J” and duplicate.                         ";
  _*N\-       " Compute N - J**2.                                                     ";
  [{          "                                                                       ";
    _mqJ/iJ*  " Compute K = floor(sqrt(N - J**2)/J)*J.                                ";
    __*@      " Duplicate, square and rotate. Result: K   K**2   N - J**2             ";
    \-        " Swap and subtract. Result: K   N - J**2 - K**2                        ";
  }3*]        " Do the above three times and collect in an array.                     ";
  )           " Pop the array. Result: N - J**2 - K**2 - L**2 - M**2                  ";
}g            " If the result is zero, break the loop.                                ";
+p            " Unshift “J” in [ K L M ] and print a string representation.           ";

Benchmarks
I've run all 5 versions over all positive integers up to 999,999,999 on my Intel Core i7-3770, measured execution time and counted the iterations required to find a solution.
The following table shows the average number of iterations and execution time for a single integer:
Version               |    1    |    2    |    3    |    4    |    5
----------------------+---------+---------+---------+---------+---------
Number of iterations  |  4.005  |  28.31  |  27.25  |  27.25  |  41.80
Execution time [µs]   |  6.586  |  39.69  |  55.10  |  63.99  |  88.81

At only 4 iterations and 6.6 microseconds per integer, primo's algorithm is incredibly fast.
Starting at floor(sqrt(N)) makes more sense, since this leaves us with smaller values for the sum of the remaining three squares. As expected, starting at 1 is a lot slower.
This is a classical example of a badly implemented good idea. To actually reduce the code size, we rely on mF, which factorizes the integer N. Although version 3 requires less iterations than version 2, it is a lot slower in practice.
Although the algorithm does not change, version 4 is a lot slower. This is because it performs an additional floating point division and an integer multiplication in each iteration.
For the input N = p**(2a) ** b, algorithm 5 will require (k - 1) * p**a + 1 iterations, where k is the number of iterations algorithm 4 requires. If k = 1 or a = 0, this makes no difference.
However, any input of the form 4**a * (4**c * (8 * d + 7) + 1) may perform quite badly. For the starting value j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7), so it cannot be expressed as a sum of three squares. Thus, k > 1 and at least p**a iterations are required.
Thankfully, primo's original algorithm is incredibly fast and N < 1,000,000,000. The worst case I could find by hand is 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1), which requires 6,145 iterations. Execution time is less than 300 ms on my machine. On average, this version is 13.5 times slower than the implementation of primo's algorithm.

"Instead of expressing N as 4**a * b, we can express it as p**(2a) * b." This is actually an improvement. I would have liked to have included this, but it was a lot longer (the ideal is to find the largest perfect square factor). "Starting with 1 and incrementing saves 4 bytes." This is definitely slower. The runtime for any given range is 4-5 times as long. "All positive integers up to 999,999,999 took 24.67 hours, giving an average execution time of 0.0888 milliseconds per integer." Perl only took 2.5 hours to crunch the whole range, and the Python translation is 10x faster ;)

@primo: Yes, you're right. Dividing by p**a is an improvement, but it's a small one. Dividing by the largest perfect square factor makes a big difference when starting from 1; it's still an improvement when starting from the integer part of the square root. Implementing it would cost only two more bytes. The abysmal execution time seems to be due to my unimprovements, not CJam. I'll rerun the tests for all algorithms (including the one you proposed), counting iterations rather than measuring wall time. Let's see how long that takes...

Finding the largest square factor only costs 2 additional bytes?! What kind of sorcery is this?

@primo: If the integer is on the stack, 1\ swaps it with 1 (accumulator), mF pushes its factorization and {~2/#*}/ raises every prime factor to its exponent divided by two, then multiplies it with the accumulator. For the direct implementation of your algorithm, that only adds 2 bytes. The small difference is mainly due to the awkward way I had to find the exponent of 4, since CJam doesn't (seem to) have a while loop...

Anyway, the benchmark finished. The total number of iterations required to factorize all 1,000,000 integers without finding the largest square factor is 4,004,829,417, with an execution time of 1.83 hours. Dividing by the largest square factor reduces the iteration count to 3,996,724,799, but it increases the time to 6.7 hours. Looks like factorizing takes a lot more time than finding the squares...

Complete integer factorization has always been a difficult problem. I suppose this means that since removing factors of four is sufficient - and fast - it's also the best in practice. I would have thought that the improvement would have been more than 0.2%, though... But still, 4 iterations on average, that's what I call amortized O(1) ;)

This is code golf, not fastest code!

Perl - 116 bytes 87 bytes (see update below)
#!perl -p
$.<<=1,$_>>=2until$_&3;
{$n=$_;@a=map{$n-=$a*($a-=$_%($b=1|($a=0|sqrt$n)>>1));$_/=$b;$a*$.}($j++)x4;$n&&redo}
$_="@a"

Counting the shebang as one byte, newlines added for horizontal sanity.
Something of a combination code-golf fastest-code submission.
The average (worst?) case complexity seems to be O(log n) O(n0.07). Nothing I've found runs slower than 0.001s, and I've checked the entire range from 900000000 - 999999999. If you find anything that takes significantly longer than that, ~0.1s or more, please let me know.

Sample Usage
$ echo 123456789 | timeit perl four-squares.pl
11110 157 6 2

Elapsed Time:     0:00:00.000

$ echo 1879048192 | timeit perl four-squares.pl
32768 16384 16384 16384

Elapsed Time:     0:00:00.000

$ echo 999950883 | timeit perl four-squares.pl
31621 251 15 4

Elapsed Time:     0:00:00.000

The final two of these seem to be worst case scenerios for other submissions. In both instances, the solution shown is quite literally the very first thing checked. For 123456789, it's the second.
If you want to test a range of values, you can use the following script:
use Time::HiRes qw(time);

$t0 = time();

# enter a range, or comma separated list here
for (1..1000000) {
  $t1 = time();
  $initial = $_;
  $j = 0; $i = 1;
  $i<<=1,$_>>=2until$_&3;
  {$n=$_;@a=map{$n-=$a*($a-=$_%($b=1|($a=0|sqrt$n)>>1));$_/=$b;$a*$i}($j++)x4;$n&&redo}
  printf("%d: @a, %f\n", $initial, time()-$t1)
}
printf('total time: %f', time()-$t0);

Best when piped to a file. The range 1..1000000 takes about 14s on my computer (71000 values per second), and the range 999000000..1000000000 takes about 20s (50000 values per second), consistent with O(log n) average complexity.

Update
Edit: It turns out that this algorithm is very similar to one that has been used by mental calculators for at least a century.
Since originally posting, I have checked every value on the range from 1..1000000000. The 'worst case' behavior was exhibited by the value 699731569, which tested a grand total of 190 combinations before arriving at a solution. If you consider 190 to be a small constant - and I certainly do - the worst case behavior on the required range can be considered O(1). That is, as fast as looking up the solution from a giant table, and on average, quite possibly faster.
Another thing though. After 190 iterations, anything larger than 144400 hasn't even made it beyond the first pass. The logic for the breadth-first traversal is worthless - it's not even used. The above code can be shortened quite a bit:
#!perl -p
$.*=2,$_/=4until$_&3;
@a=map{$=-=$%*($%=$=**.5-$_);$%*$.}$j++,(0)x3while$=&&=$_;
$_="@a"

Which only performs the first pass of the search. We do need to confirm that there aren't any values below 144400 that needed the second pass, though:
for (1..144400) {
  $initial = $_;

  # reset defaults
  $.=1;$j=undef;$==60;

  $.*=2,$_/=4until$_&3;
  @a=map{$=-=$%*($%=$=**.5-$_);$%*$.}$j++,(0)x3while$=&&=$_;

  # make sure the answer is correct
  $t=0; $t+=$_*$_ for @a;
  $t == $initial or die("answer for $initial invalid: @a");
}

In short, for the range 1..1000000000, a near-constant time solution exists, and you're looking at it.

Updated Update
@Dennis and I have made several improvements this algorithm. You can follow the progress in the comments below, and subsequent discussion, if that interests you. The average number of iterations for the required range has dropped from just over 4 down to 1.229, and the time needed to test all values for 1..1000000000 has been improved from 18m 54s, down to 2m 41s. The worst case previously required 190 iterations; the worst case now, 854382778, needs only 21.
The final Python code is the following:
from math import sqrt

# the following two tables can, and should be pre-computed

qqr_144 = set([  0,   1,   2,   4,   5,   8,   9,  10,  13,
                16,  17,  18,  20,  25,  26,  29,  32,  34,
                36,  37,  40,  41,  45,  49,  50,  52,  53,
                56,  58,  61,  64,  65,  68,  72,  73,  74,
                77,  80,  81,  82,  85,  88,  89,  90,  97,
                98, 100, 101, 104, 106, 109, 112, 113, 116,
               117, 121, 122, 125, 128, 130, 133, 136, 137])

# 10kb, should fit entirely in L1 cache
Db = []
for r in range(72):
  S = bytearray(144)
  for n in range(144):
    c = r

    while True:
      v = n - c * c
      if v%144 in qqr_144: break
      if r - c >= 12: c = r; break
      c -= 1

    S[n] = r - c
  Db.append(S)

qr_720 = set([  0,   1,   4,   9,  16,  25,  36,  49,  64,  81, 100, 121,
              144, 145, 160, 169, 180, 196, 225, 241, 244, 256, 265, 289,
              304, 324, 340, 361, 369, 385, 400, 409, 436, 441, 481, 484,
              496, 505, 529, 544, 576, 580, 585, 601, 625, 640, 649, 676])

# 253kb, just barely fits in L2 of most modern processors
Dc = []
for r in range(360):
  S = bytearray(720)
  for n in range(720):
    c = r

    while True:
      v = n - c * c
      if v%720 in qr_720: break
      if r - c >= 48: c = r; break
      c -= 1

    S[n] = r - c
  Dc.append(S)

def four_squares(n):
  k = 1
  while not n&3:
    n >>= 2; k <<= 1

  odd = n&1
  n <<= odd

  a = int(sqrt(n))
  n -= a * a
  while True:
    b = int(sqrt(n))
    b -= Db[b%72][n%144]
    v = n - b * b
    c = int(sqrt(v))
    c -= Dc[c%360][v%720]
    if c >= 0:
      v -= c * c
      d = int(sqrt(v))

      if v == d * d: break

    n += (a<<1) - 1
    a -= 1

  if odd:
    if (a^b)&1:
      if (a^c)&1:
        b, c, d = d, b, c
      else:
        b, c = c, b

    a, b, c, d = (a+b)>>1, (a-b)>>1, (c+d)>>1, (c-d)>>1

  a *= k; b *= k; c *= k; d *= k

  return a, b, c, d

This uses two pre-computed correction tables, one 10kb in size, the other 253kb. The code above includes the generator functions for these tables, although these should probably be computed at compile time.
A version with more modestly sized correction tables can be found here: http://codepad.org/1ebJC2OV This version requires an average of 1.620 iterations per term, with a worst case of 38, and the entire range runs in about 3m 21s. A little bit of time is made up for, by using bitwise and for b correction, rather than modulo.

Improvements
Even values are more likely to produce a solution than odd values.
The mental calculation article linked to previously notes that if, after removing all factors of four, the value to be decomposed is even, this value can be divided by two, and the solution reconstructed:

Although this might make sense for mental calculation (smaller values tend to be easier to compute), it doesn't make much sense algorithmically. If you take 256 random 4-tuples, and examine the sum of the squares modulo 8, you will find that the values 1, 3, 5, and 7 are each reached on average 32 times. However, the values 2 and 6 are each reached 48 times. Multiplying odd values by 2 will find a solution, on average, in 33% less iterations. The reconstruction is the following:

Care needs to be taken that a and b have the same parity, as well as c and d, but if a solution was found at all, a proper ordering is guaranteed to exist.
Impossible paths don't need to be checked.
After selecting the second value, b, it may already be impossible for a solution to exist, given the possible quadratic residues for any given modulo. Instead of checking anyway, or moving onto the next iteration, the value of b can be 'corrected' by decrementing it by the smallest amount that could possibly lead to a solution. The two correction tables store these values, one for b, and the other for c. Using a higher modulo (more accurately, using a modulo with relatively fewer quadratic residues) will result in a better improvement. The value a doesn't need any correction; by modifying n to be even, all values of a are valid.

This is incredible! The final algorithm is probably the simplest of all the answers, yet 190 iterations are all it takes...

@Dennis I would be very suprised if it hasn't been mentioned elsewhere. It seems too simple to have been over-looked.

I'm curious: Did any of the test values in your complexity analysis require the breadth-first traversal? 2. The Wikipedia article you linked to is a little confusing. It mentions the Rabin-Shallit algorithm, but provides an example for an entirely different one. 3. It would be interesting to see when exactly the Rabin-Shallit algorithm would outperform yours, I'd imagine the primality tests are rather expensive in practice.

Not one. 2. This is where I got my information (i.e. that this algorithm exists); I haven't seen the analysis, or even read the paper. 3. The curve becomes so steep at around 1e60, that it really wouldn't matter how 'slow' the O(log²n) is, it will still cross at about that point.

The second link in the question explains how to implement Rabin-Shallit, but it doesn't talk about the complexity. This answer on MathOverflow gives a nice summary of the paper. By the way, you rediscovered an algorithm used by Gottfried Ruckle in 1911 (link).

if N/2 = a²+b²+c²+d² then N = (a+b)²+(a-b)²+(c+d)²+(c-d)². Interesting. I wonder how big of a difference this would make? Seems like it might be fairly significant, as one quarter of the range has an odd factor of 2.

I did a short sampling (10 million integers), and the algorithm seems to perform much better for even numbers than for odd ones (30 % less iterations). By the way, (log n)² is less ambiguous/confusing, but usually log²n means just that.

It actually requires more iterations on average, which is disappointing. Factoring out 13 (3a+2b)²+(2a-3b)²+... and 17 (4a+b)²+(a-4b)²+... does give a meager improvement, but not enough to warrant use. It's interesting that even numbers perform better than odd... perhaps multiplying odd values by two could produce an improvement ;)

Multiplying by 2 is an improvement, but only if the original number was divisible by four (otherwise, you need to ensure even parity of a+b and c+d, which seems to eliminate too many combinations): http://codepad.org/yy1oBCqW

Well, the general idea still works: If (a,b,c,d) is a solution for 2N and the four integers are ordered by parity, 1/2 * (a + b, a - b, c + d, c - d) is a solution for N. Now that I think about it, it makes perfect sense: If N ≡ 2 mod 4, since a² ≡ 0 mod 4 or a² ≡ 1 mod 4, N - a² ≡ 2 mod 4 or N - a² ≡ 1 mod 4. This avoid both possible cases in which N - a² is not a sum of two squares. Another way of achieving the same might be the following: Force a to be even if N ≡ 1 mod 4 and to be odd if N ≡ 3 mod 4.

Even better ;) http://codepad.org/njnixbTd

Multiplying all odd values by 2 requires even less iterations, but sorting the array seems to take too much time. http://codepad.org/0QovLcAo There's probably a better way. I've never used Python before...

Function calls are relatively expensive in Python. But in this case, we don't need a full sort (with a specialized comparator function, no less), we can determine a proper ordering in no more than two comparisons: http://codepad.org/HklDWPGe

Much better, but still a tad slower than the other approach.

For values less than 10^9, it is a bit slower (presumably because most values would have found a solution on the first or second iteration anyway). For 10^10 and above, it's noticably faster. I'm curious, though, would it be possible to further specialize between n ≡ 1 (mod 4) and n ≡ 3 (mod 4)?

Let us continue this discussion in chat.

Mathematica  61 66 51
Three methods are shown.  Only the first approach meets the time requirement.

1-FindInstance (51 char)
This returns a single solution the equation.
FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &

Examples and timings
FindInstance[a^2 + b^2 + c^2 + d^2 == 123456789, {a, b, c, d}, Integers] // AbsoluteTiming

{0.003584, {{a -> 2600, b -> 378, c -> 10468, d -> 2641}}}

FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &[805306368]

{0.004437, {{a -> 16384, b -> 16384, c -> 16384, d -> 0}}}

2-IntegerPartitions
This works also, but is too slow to meet the speed requirement.
f@n_ := Sqrt@IntegerPartitions[n, {4}, Range[0, Floor@Sqrt@n]^2, 1][[1]]

Range[0, Floor@Sqrt@n]^2 is the set of all squares less than the square root of n (the largest possible square in the partition).
{4} requires the integer partitions of n consist of 4 elements from the above-mentioned set of squares.
1, within the function IntegerPartitions returns the first solution.
[[1]] removes the outer braces; the solution was returned as a set of one element.

f[123456]

{348, 44, 20, 4}

3-PowerRepresentations
PowerRepresentations returns all of the solutions to the 4 squares problem.  It can also solve for sums of other powers.
PowersRepresentations returns, in under 5 seconds,  the 181 ways to express 123456789 as the sum of 4 squares:
n= 123456;
PowersRepresentations[n, 4, 2] //AbsoluteTiming

However, it is far too slow for other sums.

Wow, Mathematica does the brute force fast. Is IntegerPartitions doing something much more clever than trying every combination, like DFT convolution on the sets? The specs ask for the numbers, by the way, not their squares.

I think Mathematica uses brute force, but probably has optimized IntegerPartitions.  As you can see from the timings, the speed varies greatly depending upon whether the first (largest) number is close to the square root of n. Thanks for catching the spec violation in the earlier version.

Could you benchmark f[805306368]? Without dividing by powers of 4 first, my solution takes 0.05 s for 999999999; I've aborted the benchmark for 805306368 after 5 minutes...

f[805306368] returns {16384, 16384, 16384} after 21 minutes.  I used {3} in place of {4}.  The attempt to solve it with a sum of 4 non-zero squares was unsuccessful after several hours of running.

I don't have access to Mathematica, but from what I've read in the documentation center, IntegerPartitions[n,4,Range[Floor@Sqrt@n]^2 should work as well. However, I don't think you should use method 1 for your score, since it doesn't comply with the time limit specified in the question.

Ok. I opted for the routine that easily satisfies the speed requirement.

FRACTRAN: 156 98 fractions
Since this is a classic number theory problem, what better way to solve this than to use numbers!
37789/221 905293/11063 1961/533 2279/481 57293/16211 2279/611 53/559 1961/403 53/299 13/53 1/13 6557/262727 6059/284321 67/4307 67/4661 6059/3599 59/83 1/59 14279/871933 131/9701 102037079/8633 14017/673819 7729/10057 128886839/8989 13493/757301 7729/11303 89/131 1/89 31133/2603 542249/19043 2483/22879 561731/20413 2483/23701 581213/20687 2483/24523 587707/21509 2483/24797 137/191 1/137 6215941/579 6730777/965 7232447/1351 7947497/2123 193/227 31373/193 23533/37327 5401639/458 229/233 21449/229 55973/24823 55973/25787 6705901/52961 7145447/55973 251/269 24119/251 72217/27913 283/73903 281/283 293/281 293/28997 293/271 9320827/58307 9831643/75301 293/313 28213/293 103459/32651 347/104807 347/88631 337/347 349/337 349/33919 349/317 12566447/68753 13307053/107143 349/367 33197/349 135199/38419 389/137497 389/119113 389/100729 383/389 397/383 397/39911 397/373 1203/140141 2005/142523 2807/123467 4411/104411 802/94883 397/401 193/397 1227/47477 2045/47959 2863/50851 4499/53743 241/409 1/241 1/239

Takes in input of the form 2n × 193 and outputs 3a × 5b × 7c × 11d. Might run in 3 minutes if you have a really good interpreter. Maybe.

 ...okay, not really. This seemed to be such a fun problem to do in FRACTRAN that I had to try it. Obviously, this isn't a proper solution to the question as it doesn't make the time requirements (it brute forces) and it's barely even golfed, but I thought I'd post this here because it's not every day that a Codegolf question can be done in FRACTRAN ;)

Hint
The code is equivalent to the following pseudo-Python:
a, b, c, d = 0, 0, 0, 0

def square(n):
    # Returns n**2

def compare(a, b):
    # Returns (0, 0) if a==b, (1, 0) if a<b, (0, 1) if a>b

def foursquare(a, b, c, d):
    # Returns square(a) + square(b) + square(c) + square(d)

while compare(foursquare(a, b, c, d), n) != (0, 0):
    d += 1

    if compare(c, d) == (1, 0):
        c += 1
        d = 0

    if compare(b, c) == (1, 0):
        b += 1
        c = 0
        d = 0

    if compare(a, b) == (1, 0):
        a += 1
        b = 0
        c = 0
        d = 0

Python 3 (177)
N=int(input())
k=1
while N%4<1:N//=4;k*=2
n=int(N**.5)
R=range(int(2*n**.5)+1)
print([(a*k,b*k,c*k,d*k)for d in R for c in R for b in R for a in[n,n-1]if a*a+b*b+c*c+d*d==N][0])

After we reduce the input N to be not divisible by 4, it must be expressible as a sum of four squares where one of them is either the largest possible value a=int(N**0.5) or one less than that, leaving only a small remainder for sum of the three other squares to take care of. This greatly reduces the search space.
Here's a proof later this code always finds a solution. We wish to find an a so that n-a^2 is the sum of three squares. From Legendre's Three-Square Theorem, a number is the sum of three squares unless it is the form 4^j(8*k+7). In particular, such numbers are either 0 or 3 (modulo 4).
We show that no two consecutive a can make the leftover amount N-a^2 have such a shape for both consecutive values.. We can do so by simply making a table of a and N modulo 4, noting that N%4!=0 because we've extracted all powers of 4 out of N.
  a%4= 0123
      +----
     1|1010
N%4= 2|2121  <- (N-a*a)%4
     3|3232

Because no two consecutive a give (N-a*a)%4 in [0,3], one of them is safe to use. So, we greedily use the largest possible n with n^2<=N, and n-1. Since N<(n+1)^2, the remainder N-a^2 to be represented as a sum of three squares is at most (n+1)^2 -(n-1)^2, which equals 4*n. So, it suffices to check only values up to 2*sqrt(n), which is exactly the range R.
One could further optimize running time by stopping after a single solution, computing rather than iterating for the last value d, and searching only among values b<=c<=d. But, even without these optimizations, the worst instance I could find finished in 45 seconds on my machine.
The chain of "for x in R" is unfortunate. It can probably be shortened by string substitution or replacement by iterating over a single index that encodes (a,b,c,d). Importing itertools turned out not worth it.
Edit: Changed to int(2*n**.5)+1 from 2*int(n**.5)+2 to make argument cleaner, same character count.

This doesn't work for me... 5 => (2, 1, 0, 0)

Strange, it works for me: I get 5 => (2, 1, 0, 0) running on Ideone 3.2.3 or in Idle 3.2.2. What do you get?

@xnor BritishColour gets 5 => (2, 1, 0, 0). Did you even read the comment? (Now we have 3 comments in a row that have that code snippet. Can we keep the streak going?)

@Quincunx If we are to decipher 5 => (2, 1, 0, 0), it means 2^2 + 1^2 + 0^2 + 0^2 = 5.  So, yes, we can.

@Dr.Rebmu What are you saying about 5 => (2, 1, 0, 0)? I understood that's what that means. Did you mean that we can keep the streak going?

@Quincunx As with many things, when you ask about my quotation of 5 => (2, 1, 0, 0) a lot of the discernment of intention is left up to the reader.

Quincunx, I read @BritishColour's comment, and as far as I can see, 5 => (2, 1, 0, 0) is correct. The examples in the question consider 0^2=0 to be a valid square number. Therefore I interpreted (as I think xnor did) that British Colour got something else. British colour, as you hav not responded again, can we assume that you do in fact get 2,1,0,0?

Sorry, I ran it again and it worked, don't know what happened then, probably a Copy/Paste fail...

@xnor You can save 2 characters at the cost of "only" a factor of 1.5^3 slowdown by changing range(2*int(n**.5)+2) to range(3*int(n**.5)). Also, you could save a bunch of characters by switching to Python 2.7.

@isaacg The 3* improvement to the range finished at just over 3 minutes on my computer on the worst-case input of 999950883, so I'm not sure if I can include it.

CJam, 91 90 74 71 bytes
q~{W):W;:N4md!}gmqi257:B_**_{;)_[Bmd\Bmd]_N\{_*-}/mq_i@+\1%}g{2W#*}%`\;

Compact, but slower than my other approach.
Try it online! Paste the Code, type the desired integer in Input and click Run.
Background
This post started as a 99 byte GolfScript answer. While there was still room for improvement, GolfScript lacks built-in sqrt function. I kept the GolfScript version until revision 5, since it was very similar to the CJam version.
However, the optimizations since revision 6 require operators that are not available in GolfScript, so instead of posting separate explanations for both languages, I decided to drop the less competitive (and much slower) version.
The implemented algorithm computes the numbers by brute force:

For input m, find N and W such that m = 4**W * N.
Set i = 257**2 * floor(sqrt(N/4)).
Set i = i + 1.
Find integers j, k, l such that i = 257**2 * j + 257 * k + l, where k, l < 257.
Check if d = N - j**2 - k**2 - l**2 is a perfect square.
If it isn't, and go back to step 3.
Print 2**W * j, 2**W * k, 2**W * l, 2**W * sqrt(m).

Examples
$ TIME='\n%e s' time cjam lagrange.cjam <<< 999999999
[27385 103 15813 14]
0.46 s
$ TIME='\n%e s' time cjam lagrange.cjam <<< 805306368
[16384 16384 0 16384]
0.23 s

Timings correspond to an Intel Core i7-4700MQ.
How it works
q~              " Read and interpret the input. ";
{
  W):W;         " Increment “W” (initially -1). ";
  :N            " Save the integer on the stack in “N”. ';
  4md!          " Push N / 4 and !(N % 4). ";
}g              " If N / 4 is an integer, repeat the loop.
mqi             " Compute floor(sqrt(N/4)). ";
257:B_**        " Compute i = floor(sqrt(N)) * 257**2. ";
_               " Duplicate “i” (dummy value). ";
{               " ";
  ;)_           " Pop and discard. Increment “i”. ";
  [Bmd\Bmd]     " Push [ j k l ], where i = 257**2 * j + 257 * k + l and k, l < 257. ";
  _N\           " Push “N” and swap it with a copy of [ j k l ]. ";
  {_*-}/        " Compute m = N - j**2 - k**2 - l**2. ";
  mq            " Compute sqrt(m). ";
  _i            " Duplicate sqrt(m) and compute floor(sqrt(m)). ";
  @+\           " Form [ j k l floor(sqrt(m)) ] and swap it with sqrt(m). ";
  1%            " Check if sqrt(m) is an integer. ";
}g              " If it is, we have a solution; break the loop. ";
{2W#*}%         " Push 2**W * [ j k l sqrt(m) ]. ";
`\;             " Convert the array into a string and discard “i”. ";

C, 228
This is based on the algorithm on the Wikipedia page, which is an O(n) brute-force.
n,*l,x,y,i;main(){scanf("%d",&n);l=calloc(n+1,8);
for(x=0;2*x*x<=n;x++)for(y=x;(i=x*x+y*y)<=n;y++)l[2*i]=x,l[2*i+1]=y;
for(x=0,y=n;;x++,y--)if(!x|l[2*x+1]&&l[2*y+1]){
printf("%d %d %d %d\n",l[2*x],l[2*x+1],l[2*y],l[2*y+1]);break;}}

GolfScript, 133 130 129 bytes
~{.[4*{4/..4%1$!|!}do])\}:r~,(2\?:f;{{..*}:^~4-1??n*,}:v~)..
{;;(.^3$\-r;)8%!}do-1...{;;;)..252/@252%^@^@+4$\-v^@-}do 5$]{f*}%-4>`

Fast, but lengthy. The newline can be removed.
Try it online. Note that the online interpreter has a 5 second time limit, so it might not work for all numbers.
Background
The algorithm takes advantage of Legendre's three-square theorem, which states that every natural number n that is not of the form

can be expressed as the sum of three squares.
The algorithm does the following:

Express the number as 4**i * j.
Find the largest integer k such that k**2 <= j and j - k**2 satisfies the hypothesis of Legendre's three-square theorem.
Set i = 0.
Check if j - k**2 - (i / 252)**2 - (i % 252)**2 is a perfect square.
If it isn't, increment i and go back to step 4.

Examples
$ TIME='%e s' time golfscript legendre.gs <<< 0
[0 0 0 0]
0.02 s
$ TIME='%e s' time golfscript legendre.gs <<< 123456789
[32 0 38 11111]
0.02 s
$ TIME='%e s' time golfscript legendre.gs <<< 999999999
[45 1 217 31622]
0.03 s
$ TIME='%e s' time golfscript legendre.gs <<< 805306368
[16384 0 16384 16384]
0.02 s

Timings correspond to an Intel Core i7-4700MQ.
How it works
~              # Interpret the input string. Result: “n”
{              #
  .            # Duplicate the topmost stack item.
  [            #
    4*         # Multiply it by four.
    {          #
      4/       # Divide by four.
      ..       # Duplicate twice.
      4%1$     # Compute the modulus and duplicate the number.
      !|!      # Push 1 if both are truthy.
    }do        # Repeat if the number is divisible by four and non-zero.
  ]            # Collect the pushed values (one per iteration) into an array.
  )\           # Pop the last element from the array and swap it with the array.
}:r~           # Save this code block as “r” and execute it.
,(2\?          # Get the length of the array, decrement it and exponentiate.
:f;            # Save the result in “f”.
               # The topmost item on the stack is now “j”, which is not divisible by
               # four and satisfies that n = f**2 * j.
{              #
  {..*}:^~     # Save a code block to square a number in “^” and execute it.
  4-1??        # Raise the previous number to the power of 1/4.
               # The two previous lines compute (x**2)**(1/4), which is sqrt(abs(x)).
  n*,          # Repeat the string "\n" that many times and compute its length.
               # This casts to integer. (GolfScript doesn't officially support Rationals.)
}:v~           # Save the above code block in “v” and execute it.
)..            # Undo the first three instructions of the loop.
{              #
   ;;(         # Discard two items from the stack and decrement.
   .^3$\-      # Square and subtract from “n”.
   r;)8%!      # Check if the result satisfies the hypothesis of the three-square theorem.
}do            # If it doesn't, repeat the loop.
-1...          # Push 0 (“i”) and undo the first four instructions of the loop.
{              #
  ;;;)         # Discard two items from the stack and increment “i”.
  ..252/@252%  # Push the digits of “i” in base 252.
  ^@^@+4$\-    # Square both, add and subtract the result
  v^@-         # Take square root, square and compare.
}do            # If the difference is a perfect square, break the loop.
5$]            # Duplicate the difference an collect the entire stack into an array.
{f*}%          # Multiply very element of the array by “f”.
-4>            # Reduce the array to its four last elements (the four numbers).
`              # Convert the result into a string.

I didn't understand j-k-(i/252)-(i%252). From your comments (I cant actually read the code), it looks like you mean j-k-(i/252)^2-(i%252)^2. BTW, the equivalent of j-k-(i/r)^2-(i%r)^2 where r=sqrt(k) may save a few characters (and seems to work without problems even for k=0 in my C program.)

@steveverrill: Yes, I made a mistake. Thank you for noticing. It should be j-k^2-(i/252)^2-(i%252)^2. I'm still waiting for the OP to clarify if 0 is a valid input or not. Your program gives 1414 -nan 6 4.000000 for input 0.

I'm talking about my new program using Legendre's theorem, which I haven't posted yet. It looks like it never calls the code with the % or / when I have the equivalent of k=0, which is why it's not causing problems. You'll  see when I post it. Glad you got my old program running. Did you have the memory to build the full 2GB table in rev 1, and how long did it take?

Yeah, the C compiler can behave quite unexpectedly when optimizing. In GolfScript, 0/ => crash! :P I've run your rev 1 on my laptop (i7-4700MQ, 8 GiB RAM). On average, execution time is 18.5 seconds.

Wow is that 18.5 seconds including building the table? It takes over 2 minutes on my machine. I can see the problem is the Windows memory management. Rather than give the program the 2GB it needs straight away, it gives it in small chunks, so it must be doing a lot of unnecessary swapping until the full 2GB is allocated. Actually searching for the answer per user input is much faster, because by then the program doesn't have to go begging for memory.

@steveverrill: Yes, but that was with the default optimization level. With -Ofast, the entire process takes about 13 seconds. Fun fact: With -Ofast, you don't need -lm.

Rev B: C, 179
a,b,c,d,m=1,n,q,r;main(){for(scanf("%d",&n);n%4<1;n/=4)m*=2;
for(a=sqrt(n),a-=(3+n-a*a)%4/2;r=n-a*a-b*b-c*c,d=sqrt(r),d*d-r;c=q%256)b=++q>>8;
printf("%d %d %d %d",a*m,b*m,c*m,d*m);}

Thanks to @Dennis for the improvements. The rest of the answer below is not updated from rev A.
Rev A:  C,195
a,b,c,d,n,m,q;double r=.1;main(){scanf("%d",&n);for(m=1;!(n%4);n/=4)m*=2;a=sqrt(n);a-=(3+n-a*a)%4/2;
for(;(d=r)-r;q++){b=q>>8;c=q%256;r=sqrt(n-a*a-b*b-c*c);}printf("%d %d %d %d ",a*m,b*m,c*m,d*m);}

Much faster than my other answer and with much less memory!
This uses http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem. Any number not of the following form can be expressed as the sum of 3 squares (I call this the prohibited form):
4^a*(8b+7), or equivalently 4^a*(8b-1)
Note that all odd square numbers are of the form (8b+1) and all even square numbers are superficially of the form 4b. However this hides the fact that all even square numbers are of the form 4^a*(odd square)==4^a*(8b+1). As a result 2^x-(any square number < 2^(x-1)) for odd xwill always be of the prohibited form. Hence these numbers and their multiples are difficult cases, which is why so many of the programs here divide out powers of 4 as a first step.
As stated in @xnor's answer, N-a*a cannot be of the prohibited form for 2 consecutive values of a. Below I present a simplified form of his table. In addition to the fact that after division by 4 N%4 cannot equal 0, note that there are only 2 possible values for (a*a)%4.
(a*a)%4= 01
        +--
       1|10
  N%4= 2|21  <- (N-a*a)%4
       3|32

So, we want to avoid values of (N-a*a) that may be of the prohibited form, namely those where (N-a*a)%4 is 3 or 0. As can be seen this cannot occur for the same N with both odd and even (a*a).
So, my algorithm works like this:
1. Divide out powers of 4
2. Set a=int(sqrt(N)), the largest possible square
3. If (N-a*a)%4= 0 or 3, decrement a (only once)
4. Search for b and c such that N-a*a-b*b-c*c is a perfect square

I particularly like the way I do step 3. I add 3 to N, so that the decrement is required if (3+N-a*a)%4 = 3 or 2. (but not 1 or 0.) Divide this by 2 and the whole job can be done by a fairly simple expression.
Ungolfed code
Note the single for loop q and use of division/modulo to derive the values of b and c from it. I tried using a as a divisor instead of 256 to save bytes, but sometimes the value of a was not right and the program hung, possibly indefinitely. 256 was the best compromise as I can use >>8 instead of /256 for the division.
a,b,c,d,n,m,q;double r=.1;
main(){
  scanf("%d",&n);
  for(m=1;!(n%4);n/=4)m*=2;
  a=sqrt(n);
  a-=(3+n-a*a)%4/2;
  for(;(d=r)-r;q++){b=q>>8;c=q%256;r=sqrt(n-a*a-b*b-c*c);}
  printf("%d %d %d %d ",a*m,b*m,c*m,d*m);}

Output
An interesting quirk is that if you input a square number, N-(a*a)=0. But the program detects that 0%4=0 and decrements to the next square down. As a result square number inputs are always decomposed into a group of smaller squares unless they are of the form 4^x.
999999999
31621 1 161 294

805306368
16384 0 16384 16384

999950883
31621 1 120 221

1
0 0 0 1

2
1 0 0 1

5
2 0 0 1

9
2 0 1 2

25
4 0 0 3

36
4 0 2 4

49
6 0 2 3

81
8 0 1 4

121
10 1 2 4

Amazing! 0.003 s for every input! You can get those 5 chars back: 1. Declare m=1 before main. 2. Execute scanf in the for statement. 3. Use float instead of double. 4. n%4<1 is shorter than !(n%4). 5. There's an obsolete space in printf's format string.

A few more suggestions.

Thanks for the tips! n-=a*a doesn't work, because a can be modified afterwards (it gives some wrong answers and hangs on a small number of cases, like 100+7=107.) I included all the rest. It would be nice to something to shorten the printf but I think the only answer is to change the language. The key to speed is to settle on a good value for a quickly. Written in C and with a search space of less than 256^2, this is probably the fastest program here.

Right, sorry. Shortening the printf statement seems difficult without using a macro or an array, which would add bulk elsewhere. Changing languages seems the "easy" way. Your approach would weigh 82 bytes in CJam.

JavaScript - 175 191 176 173 168 chars
Brute force, but fast.
Edit
Fast but not enough for some nasty input. I had to add a first step of reduction by multiplies of 4.
Edit 2
Get rid of function, output inside the loop then force exit contition
Edit 3
0 not a valid input
Edit 4
Spurred by GalaxyCat105: a function again, but avoiding the length of return keyword
v=>(q=>{for(m=1;!(v%4);m+=m)v/=4;for(a=-~q(v);q&&a--;)for(w=v-a*a,b=-~q(w);q&&b--;)for(x=w-b*b,c=-~q(x);q&&c--;)(d=q(x-c*c))==~~d&&(q=0)})(Math.sqrt)||[m*a,m*b,m*c,m*d]

TEST

F=
v=>(q=>{for(m=1;!(v%4);m+=m)v/=4;for(a=-~q(v);q&&a--;)for(w=v-a*a,b=-~q(w);q&&b--;)for(x=w-b*b,c=-~q(x);q&&c--;)(d=q(x-c*c))==~~d&&(q=0)})(Math.sqrt)||[m*a,m*b,m*c,m*d]

test=[123456789, 805306368, 999999999, 265289728, 699731569];

console.log(test.map(x => x+' : '+F(x)))

Ungolfed:
v = prompt();

for (m = 1; ! (v % 4); m += m)
{
  v /= 4;
}
for (a = - ~Math.sqrt(v); a--;) /* ~ force to negative integer, changing sign lead to original value + 1 */
{
  for ( w = v - a*a, b = - ~Math.sqrt(w); b--;)
  {
    for ( x = w - b*b, c = - ~Math.sqrt(x); c--;)
    {
      (d = Math.sqrt(x-c*c)) == ~~d && prompt([m*a, m*b, m*c, m*d], a=b=c='') /* 0s a,b,c to exit loop */
    }
  }
}

You can remove the spaces in the array - v=(p=prompt)();for(m=1;!(v%4);m+=m)v/=4;for(a=-~(q=Math.sqrt)(v);a--;)for(w=v-a*a,b=-~q(w);b--;)for(x=w-b*b,c=-~q(x);c--;)(d=q(x-c*c))==~~d&&p([m*a,m*b,m*c,m*d],a=b=c='')

@GalaxyCat105 wow you spotted that mistake after 6 years! Thanx

Rev 1: C,190
a,z,m;short s[15<<26];p(){m=s[a=z-a];printf("%d %f ",m,sqrt(a-m*m));}
main(){m=31727;for(a=m*m;--a;s[z<m*m?z:m*m]=a%m)z=a/m*(a/m)+a%m*(a%m);scanf("%d",&z);for(;a*!s[a]||!s[z-a];a++);p();p();}

This is even more memory hungry than rev 0. Same principle: build a table with a positive value for all possible sums of 2 squares (and zero for those numbers which are not sums of two squares), then search it.
In this rev use an array of short instead of char to store the hits, so I can store the root of one of the pair of squares in the table instead of just a flag. This simplifies function p (for decoding the sum of 2 squares) considerably as there is no need for a loop.
Windows has a 2GB limit on arrays. I can get round that with short s[15<<26] which is an array of 1006632960 elements, enough to comply with the spec. Unfortunately, the total program runtime size is still over 2GB and (despite tweaking the OS settings) I have not been able to run over this size (though it is theoretically possible.) The best I can do is short s[14<<26] (939524096 elements.) m*m must be strictly less than this (30651^2=939483801.) Nevertheless, the program runs perfectly and should work on any OS that does not have this restriction.
Ungolfed code
a,z,m;
short s[15<<26];
p(){m=s[a=z-a];printf("%d %f ",m,sqrt(a-m*m));}
main(){
 m=31727;
 for(a=m*m;--a;s[z<m*m?z:m*m]=a%m)   //assignment to s[] moved inside for() is executed after the following statement. In this rev excessively large values are thrown away to s[m*m].
   z=a/m*(a/m)+a%m*(a%m);            //split a into high and low half, calculate h^2+l^2.
 scanf("%d",&z);
 for(;a*!s[a]||!s[z-a];a++);         //loop until s[a] and s[z-a] both contain an entry. s[0] requires special handling as s[0]==0, therefore a* is included to break out of the loop when a=0 and s[z] contains the sum of 2 squares.
 p();                                //print the squares for the first sum of 2 squares
 p();}                               //print the squares for the 2nd sum of 2 squares (every time p() is called it does a=z-a so the two sums are exchanged.)

Rev 0 C,219
a,z,i,m;double t;char s[1<<30];p(){for(i=t=.1;(m=t)-t;i++)t=sqrt(a-i*i);printf("%d %f ",i-1,t);}
main(){m=1<<15;for(a=m*m;--a;){z=a/m*(a/m)+a%m*(a%m);s[z<m*m?z:0]=1;}scanf("%d",&z);for(;1-s[a]*s[z-a];a++);p();a=z-a;p();}

This is a memory hungry beast. It takes a 1GB array, calculates all possible sums of 2 squares and stores a flag for each in the array. Then for the user input z, it searches the array for two sums of 2 squares a and z-a.
the function p then reconsitutes the original squares that were used to make the sums of 2 squares a and z-a and prints them, the first of each pair as an integer, the second as a double (if it has to be all integers two more characters are needed,t > m=t.)
The program takes a couple of minutes to build the table of sums of squares (I think this is due to memory management issues, I see the memory allocation going up slowly instead of jumping up as one might expect.) However once that is done it produces answers very quickly (if several numbers are to be calculated, the program from scanf onward can be put in a loop.
ungolfed code
a,z,i,m;
double t;
char s[1<<30];                              //handle numbers 0 up to 1073741823
p(){
 for(i=t=.1;(m=t)-t;i++)t=sqrt(a-i*i);      //where a contains the sum of 2 squares, search until the roots are found
 printf("%d %f ",i-1,t);}                   //and print them. m=t is used to evaluate the integer part of t.

main(){
 m=1<<15;                                   //max root we need is sqrt(1<<30);
 for(a=m*m;--a;)                            //loop m*m-1 down to 1, leave 0 in a
   {z=a/m*(a/m)+a%m*(a%m);s[z<m*m?z:0]=1;}  //split a into high and low half, calculate h^2+l^2. If under m*m, store flag, otherwise throw away flag to s[0]
 scanf("%d",&z);
 for(;1-s[a]*s[z-a];a++);                   //starting at a=0 (see above) loop until flags are found for sum of 2 squares of both (a) and (z-a)
 p();                                       //reconsitute and print the squares composing (a)
 a=z-a;                                     //assign (z-a) to a in order to...
 p();}                                      //reconsitute and print the squares composing (z-a)

Example output
The first is per the question. The second was picked as a difficult one to search for. In this case the program has to search as far as 8192^2+8192^2=134217728, but only takes a few seconds once the table is built.
123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Shouldn't you add a prototype for sqrt?

@edc65 I'm using GCC compiler (which is for Linux, but I have Cygwin Linux environment installed on my Windows machine.) This means I don't need to put #include <stdio.h> (for scanf/printf) or #include <math.h> (for sqrt.) The compiler links the necessary libraries automatically. I have to thank Dennis for that (he told me on this question http://codegolf.stackexchange.com/a/26330/15599) Best golfing tip I ever had.

I was already wondering why Hunt the Wumpus appeared in the linked questions. :) By the way, I don't know what GCC uses on Windows, but the GNU linker does not link the math library automatically, with or without the include. To compile on Linux, you need the flag -lm

@Dennis that's interesting, it does include stdio and several other libraries, but not math even with the include? By which I understand if you put the compiler flag, you don't need the include anyway? Well it's working for me, so I'm not complaining, thanks again for the tip. BTW I'm hoping to post a completely different answer taking advantage of Legendre's theorem (but it will still use a sqrt.)

-lm affects the linker, not the compiler. gcc opts to not require the prototypes for functions it "knows", so it works with or without the includes. However, the header filesprovide only function prototypes, not the functions themselves. On Linux (but not Windows, apparently), the math library libm is not part of the standard libraries, so you have to instruct ld to link to it.

Mathematica, 138 chars
So it turns out that this produces negative and imaginary results for certain inputs as pointed out by edc65 (e.g., 805306368), so this isn't a valid solution. I'll leave it up for now, and maybe, if I really hate my time, I'll go back and try to fix it.
S[n_]:=Module[{a,b,c,d},G=Floor@Sqrt@#&;a=G@n;b:=G[n-a^2];c:=G[n-a^2-b^2];d:=G[n-a^2-b^2-c^2];While[Total[{a,b,c,d}^2]!=n,a-=1];{a,b,c,d}]

Or, unsquished:
S[n_] := Module[{a, b, c, d}, G = Floor@Sqrt@# &;
 a = G@n;
 b := G[n - a^2];
 c := G[n - a^2 - b^2];
 d := G[n - a^2 - b^2 - c^2];
 While[Total[{a, b, c, d}^2] != n, a -= 1];
 {a, b, c, d}
]

I didn't look too hard at the algorithms, but I expect this is the same idea. I just came up with the obvious solution and tweaked it until it worked. I tested it for all numbers between 1 and one billion and... it works. The test only takes about 100 seconds on my machine.
The nice bit about this is that, since b, c, and d are defined with delayed assignments, :=, they don't have to be redefined when a is decremented. This saved a few extra lines I had before. I might golf it further and nest the redundant parts, but here's the first draft.
Oh, and you run it as S@123456789 and you can test it with {S@#, Total[(S@#)^2]} & @ 123456789 or # == Total[(S@#)^2]&[123456789]. The exhaustive test is
n=0;
AbsoluteTiming@ParallelDo[If[e != Total[(S@e)^2], n=e; Abort[]] &, {e, 1, 1000000000}]
n

I used a Print[] statement before but that slowed it down a lot, even though it never gets called. Go figure.

This is really clean! I'm surprised that it suffices to simply take every value but the the first as large as possible. For golfing, it's probably shorter to save n - a^2 - b^2 - c^2 as a variable and check that d^2 equals it.

Does it really work? What solution does it find for input 805306368?

S[805306368]={-28383, 536 I, 32 I, I}. Huh. That does produces 805306368 when you sum it, but obviously there is a problem with this algorithm. I guess I'll have to retract this for now; thanks for pointing that out...

The numbers that fail all seem to be divisible by large powers of 2. Specifically, they seem to be of the form a * 4^(2^k) for k>=2, having extracted out all powers of 4 so that a isn't a multiple of 4 (but could be even). Moreover, each a is either 3 mod 4, or twice such a number. The smallest one is 192.

Haskell 123+3=126
main=getLine>>=print.f.read
f n=head[map(floor.sqrt)[a,b,c,d]|a<-r,b<-r,c<-r,d<-r,a+b+c+d==n]where r=[x^2|x<-[0..n],n>=x^2]

Simple brute force over pre-calculated squares.
It needs the -O compilation option (I added 3 chars for this). It takes less than 1 minute for the worst case 999950883.
Only tested on GHC.

C: 198 characters
I can probably squeeze it down to just over 100 characters. What I like about this solution is the minimal amount of junk, just a plain for-loop, doing what a for-loop should do (which is to be crazy).
i,a,b,c,d;main(n){for(scanf("%d",&n);a*a+b*b-n?a|!b?a*a>n|a<b?(--a,b=1):b?++b:++a:(a=b=0,--n,++i):c*c+d*d-i?c|!d?c*c>i|c<d?(--c,d=1):d?++d:++c:(a=b=c=d=0,--n,++i):0;);printf("%d %d %d %d",a,b,c,d);}

And heavily prettified:
#include <stdio.h>

int n, i, a, b, c, d;

int main() {
    for (
        scanf("%d", &n);
        a*a + b*b - n
            ? a | !b
                ? a*a > n | a < b
                    ? (--a, b = 1)
                    : b
                        ? ++b
                        : ++a
                : (a = b = 0, --n, ++i)
            : c*c + d*d - i
                ? c | !d
                    ? c*c > i | c < d
                        ? (--c, d = 1)
                        : d
                            ? ++d
                            : ++c
                    : (a = b = c = d = 0, --n, ++i)
                : 0;
    );
    printf("%d %d %d %d\n", a, b, c, d);
    return 0;
}

Edit: It's not fast enough for all input, but I will be back with another solution. I'll let this ternary operation mess stay as of now.

Haskell, 117 bytes
w=floor.sqrt.fromIntegral
s 1 n=[[x]|x<-[w n],x^2==n]
s i n=[x:y|x<-[w n,w n-1..0],y<-s(i-1)(n-x^2),y/=[]]
f=head.s 4

Try it online!

Scala, 626 bytes
\$\text{626 bytes, it can be golfed much more.}\$
Golfed version. Try it online!
import scala.math.{sqrt,pow}

object Main {
  def isSq(n: Int) = pow(sqrt(n).toInt, 2) == n

  def decomp(n: Int): Seq[Int] = {
    for (i <- 1 until sqrt(n).toInt) {
      val r1 = n - pow(i, 2).toInt
      if (isSq(r1)) return Seq(i, sqrt(r1).toInt, 0, 0)
      for (j <- 1 until sqrt(r1).toInt) {
        val r2 = r1 - pow(j, 2).toInt
        if (isSq(r2)) return Seq(i, j, sqrt(r2).toInt, 0)
      }
    }
    Seq(sqrt(n).toInt, sqrt(n - pow(sqrt(n).toInt, 2)).toInt, 0, 0)
  }

  def main(args: Array[String]) = {
    val n = 123456789
    println(s"The four square numbers for $n are ${decomp(n).mkString(", ")}")
  }
}

Ungolfed version. Try it online!

Here is a basic implementation in Scala following the algorithm as described by Joseph Louis Lagrange. This algorithm is not the most efficient one but it should run within the required time constraints.
This algorithm uses Gauss's "Three-Square" theorem as an optimization. According to this theorem, a natural number can be expressed as the sum of three squares of integers if and only if it's not of the form 4^n(8m+7).
import scala.math.sqrt
import scala.math.pow

object Main {
  def isSquare(n: Int): Boolean = {
    val root = sqrt(n).toInt
    pow(root, 2) == n
  }

  def decompose(n: Int): Seq[Int] = {
    if (isSquare(n)) {
      Seq(n, 0, 0, 0)
    } else {
      for (i <- 1 until sqrt(n).toInt) {
        val remainder = n - pow(i, 2).toInt
        if (isSquare(remainder)) {
          return Seq(i, sqrt(remainder).toInt, 0, 0)
        }
        for (j <- 1 until sqrt(remainder).toInt) {
          val remainder2 = remainder - pow(j, 2).toInt
          if (isSquare(remainder2)) {
            return Seq(i, j, sqrt(remainder2).toInt, 0)
          }
          val k = sqrt(remainder2).toInt
          if (pow(k, 2) == remainder2) {
            return Seq(i, j, k, 0)
          }
        }
      }
      val root = sqrt(n).toInt
      Seq(root, sqrt(n - pow(root, 2)).toInt, 0, 0)
    }
  }

  def main(args: Array[String]): Unit = {
    val num = 123456789
    val result = decompose(num)
    println(s"The four square numbers for $num are ${result.mkString(", ")}")
  }
}

The above program tries to find 1 to 3 square numbers whose sum is equal to the input number. If the input number is a perfect square, it directly returns it along with three zeros. If not, it tries to find two square numbers whose sum is equal to the input number. If still not possible, it will look for three square numbers. For the remaining number, the fourth square number is calculated as the difference between the input number and the sum of squares of the first three numbers.
Please note that due to the limitation of the Scala language, the size of integers that can be processed by this program is limited. Numbers larger than Int.MaxValue (~2 billion) might cause overflow errors.

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

authors: АЛЕКСЕЙ ТРОФИМОВ


text:
Write a Boolfuck Compiler or Interpreter
From now on considering the word "fuck" to be non-profane would be a good idea
Make the smallest compiler or an interpreter for the Brainfuck variant, Boolfuck. It has to be executable in Windows, DOS, or a Unix operating system (so Linux, Mac, OpenBSD etc). Boolfuck is very similar to Brainfuck, and a complete definition can be found here.
Differences from Brainfuck

Works with bits not bytes (each cell is a single bit)
The output character is ; not .
The input is read in little endian mode and reads a single bit at a time.
Output is stored to a stream, then outputted in little endian mode

Rules

Use bits for the cells, not bytes
Support the complete list of operations from the official definition:
+  Flips the value of the bit under the pointer.

,  Reads a bit from the input stream, storing it under the pointer.
   The end-user types information using characters, though. Bytes
   are read in little-endian order; the first bit read from the
   character a, for instance, is 1, followed by 0, 0, 0, 0, 1, 1,
   and finally 0.

;  Outputs the bit under the pointer to the output stream. The bits
   get output in little-endian order, the same order in which they
   would be input. If the total number of bits output is not a
   multiple of eight at the end of the program, the last character
   of output gets padded with zeros on the more significant end. If
   the end-of-file character has been input, outputs a zero to the
   bit under the pointer.

<  This moves the pointer left by one bit.

>  This moves the pointer right by one bit.

[  If the value under the pointer is zero, jumps forward, just past
   the matching ] character.
]  Jumps back to the matching [ character.

Any character not in the official definition is ignored.
The program must start in the middle of the tape, not on an end.
The tape must be unbounded in size (infinite in both directions).

Here is a sample Boolfuck program; it outputs Hello, world!:
;;;+;+;;+;+;
+;+;+;+;;+;;+;
;;+;;+;+;;+;
;;+;;+;+;;+;
+;;;;+;+;;+;
;;+;;+;+;+;;
;;;;;+;+;;
+;;;+;+;;;+;
+;;;;+;+;;+;
;+;+;;+;;;+;
;;+;;+;+;;+;
;;+;+;;+;;+;
+;+;;;;+;+;;
;+;+;+;

Would a just-in-time compiler satisfy the requirement of 'compiler'?

In my sleep-deprived state, I was really excited when it occurred to me that it would be really simple to make a Brainfuck to Boolfuck compiler. Then I realized that that was the opposite of what was asked. :(

You could still do it... It wouldn't meet the requirements for the question though. It WOULD be very easy to do that, because the website I gave in the question gives all the conversion you need @millinon

@primo Yes it would

@OVER9000 how about an interpreter?

I believe you misunderstood what , does: "Reads a bit from the input stream, storing it under the pointer. The end-user types information using characters, though. Bytes are read in little-endian order; the first bit read from the character a, for instance, is 1, followed by 0, 0, 0, 0, 1, 1, and finally 0." That specifically says bit. Also, >,>,>,>,>,>,>,>,<<<<<<<< is listed as an equivalent of Brainf***'s , (which reads one 8 bit byte). Since there are 8 ,s, this means that each , reads one bit.

I edited the post to accept an interpreter @m.buettner

@Quincunx Yes indeed you are correct - edited.

Perl - 171 bytes
#!perl -0
%h=qw(< --$p > ++$p + $$p^=1 [ while($$p){ ] } ; $o|=$$p<<$g++;$g&=7or$o=!print+chr$o , $f&=7or$_=getc;$$p=1&ord>>$f++);$p=-1e9;eval<ARGV>=~s/./$h{$&};/gr;$g&&print+chr$o

Counting the shebang as 1 byte.
This translates the source into perl, and then evaluates the translation. The target program is taken as a command line argument, input sent to the program should be piped via stdin.
Sample usage:
$ perl boolf.pl hello.blf
Hello, world!

Details
Each of the commands is translated as follows:
< ⇒ --$p;
> ⇒ ++$p;
+ ⇒ $$p^=1;
[ ⇒ while($$p){
] ⇒ }
; ⇒ $o|=$$p<<$g++;$g&=7or$o=!print+chr$o;
, ⇒ $f&=7or$_=getc;$$p=1&ord>>$f++;
Negative integers are valid variable names, which is exploited here. $p, the cell pointer, is initialized as -1000000000, which is approximately at the center of the valid range, -2147483648 .. -1.
I've violated the specification in one regard: output is not printed all at once at the end, but rather as soon as a full byte has been collected. It's slightly longer this way, but I think it's more useful, as it allows for interactive scripts. For example, the Rot13 program below can be run in 'interactive mode', terminating the program with an eof char (^D on linux, ^Z on windows).

Edit - Demonstration of Endianness
There seems to be some doubt as to whether this interpreter reads bits in the proper (little-endian) order. I offer the following Boolfuck program, which reads a total of 8 bits, one at a time, and outputs each value read as its numeric ascii representation (0 outputs character 48, 1 outputs character 49). Each value is output entirely before the next bit is even read.
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;
,;[+];;;+;;+;;

Sample usage:
$ echo a | perl boolf.pl test.blf
10000110

From the problem description:
"Bytes are read in little-endian order; the first bit read from the character a, for instance, is 1, followed by 0, 0, 0, 0, 1, 1, and finally 0."
According to this description, the bits are read in the proper order.

Stress Testing
First, a script to translate Brainfuck programs into Boolfuck:
#!perl -n0

%h = (
  '+' => '>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<',
  '-' => '>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]<<<<<<<<<',
  '<' => '<<<<<<<<<',
  '>' => '>>>>>>>>>',
  ',' => '>,>,>,>,>,>,>,>,<<<<<<<<',
  '.' => '>;>;>;>;>;>;>;>;<<<<<<<<',
  '[' => '>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]',
  ']' => '>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]'
);

$o .= $h{$_} for /./g;
$o =~ s/<(?R)?>//g;
print $o

The Rot13 program from the Brainfuck Wiki Page translates thusly (with slight modification for null byte eof, rather than 255):
,>,>,>,>,>,>,>,>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]
<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<
[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]+<<<<<<<<+[
>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>
>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]
>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+
<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<
]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[
<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]>[>]+<[+<]>>>>>>>>>[+]>>
>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+<<<<<<<<+[>+]<[<]
>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>
>>>>>>>>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<
]>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>
>>>>>>>]<[+<]<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<
<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+
[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>+
<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>
>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>
>[+<<<<<<<<[>]+<[+<]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<
[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]
+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>
>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+<<<<<<<
<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>
>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>
>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>
>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>
>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<
[>]+<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[
+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]
>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>[>]+<
[+<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<[>]+<[+<]>>>>>>
>>>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>
>>>[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<<<
<<<<<<<<<<<<<<<<<<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<
[+<]>>>>>>>>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<
[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>
>>]<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<
<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<
<<<<<+[>+]<[<]>>>>>>>>>[+]<<<<<<<<<<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+
<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>
>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>
>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<
<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>>>>>>>>>>>>+
<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+
<]<<<<<<<<<<<<<<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>
[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+<<<<<<
<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<
<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<;>;>;>;>;>;>;>;>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<
<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+<
<<<<<<<+[>+]<[<]>>>>>>>>>[+]<<<<<<<<,>,>,>,>,>,>,>,>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]

On a 26 character string, this runs in approximately 0.17 seconds:
$ echo abcdefghijklmnopqrstuvwxyz | timeit perl boolf.pl rot13.blf
nopqrstuvwxyzabcdefghijklm

Elapsed Time:     0:00:00.171

The squares.b from Daniel Cristofani's site:
>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<
<<<[>]+<[+<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<
[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]
+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+
]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<
<<<<[>]+<[+<]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[
+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>+<<
<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<
<<<<<<<<<<<<<<[>]+<[+<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>
>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<[+<]>>>>>>
>>>[+]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<
+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>
>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<<<<<<<<<<[>]+<[+<]>>>>>
>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[
+<]>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]
>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>
>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[
>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+
<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[
+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<
<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<
<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<
<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>>>>>>>>>>+
<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<<<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<
<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<[>]+<[+<]>>
>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]
>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+
<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<
<<<<[>]+<[+<]>>>>>>>>>[+]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>
+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<<<<<<<<<;>;>;>;>;>;>;>;<<<<<<<<+<
<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>
>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]
>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+
<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+
<]>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+
<<<<<<<<[>]+<[+<]>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<
<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<
]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<
<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<
<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]
<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>
+]<[<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]
<[+<]>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]
<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[
+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>>[+]<<<<<<<<[>]+<[+<]>>>>>>>>
>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<
]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]+<<<<<<<<+[>+]<[<]>>>>>>>>
>[+]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>[>]+<[+<
]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]<<<<<<<<<<<<<<<<<<+<<<<<<<
<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]+<<<<<<<<+[>+]<[<]
>>>>>>>>>[+<<<<<<<<[>]+<[+<]>>>>>>>>>>[>]+<[+<]>>>>>>>>>[+]<<<<<<<<<+<<<<<<<<+[>+]<[<
]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>
>>>>>>]<[+<]<<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]<
<<<<<<<<+<<<<<<<<+[>+]<[<]>>>>>>>>>[+]+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]>>>>>>>>>+<<<<
<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]>;>;>;>;>;>;>;>;<<<<<<<;>;>;>;>;>;>;>;<<<<<<
<;>;>;>;>;>;>;>;>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]

Runs in about 2.4s:
$ timeit perl boolf.pl out.dat
0
1
4
9
16
25
36
49
64
81
100
121
144
169
196
225
256
289
324
361
400
441
484
529
576
625
676
729
784
841
900
961
1024
1089
1156
1225
1296
1369
1444
1521
1600
1681
1764
1849
1936
2025
2116
2209
2304
2401
2500
2601
2704
2809
2916
3025
3136
3249
3364
3481
3600
3721
3844
3969
4096
4225
4356
4489
4624
4761
4900
5041
5184
5329
5476
5625
5776
5929
6084
6241
6400
6561
6724
6889
7056
7225
7396
7569
7744
7921
8100
8281
8464
8649
8836
9025
9216
9409
9604
9801
10000

Elapsed Time:     0:00:02.421

Shame the tape isn't infinite (within memory constraints of course). Still: well done!

Also a shame the byte is read big-endian

@OVER9000 it isn't, though. $f&=7or$_=getc;$$p=1&ord>>$f++; The lowest bit is read first (if f is greater than 7, read a new byte and set f zero. Set the current cell to the ordinal value of the character just read, right shifted by f (post increment), mod 2.

@OVER9000 I've added a demonstration of endianness to the answer. Please review.

@primo yes you are correct, thanks for clearing that up*

Python (531 byte compiler)
This is a source-to-source compiler: it compiles Boolfuck to Python. Takes input from a file as its first argument, outputs Python code to file o in the current dict.
exec"eJxdUs1u8yAQvOcpEBez/lNS9RDZpS+C+CQ3xg5RwlKMmkTV9+5l7bZue2RmZ4YdGAJe2HSfmL14DDFvndy2o0RvnOjC+KZ2Gupgul5A20vFh7/zM2CjCRHxvMI+2UQ5mZh0KJXmuu2LpD9J9WKdwNCLO4Bqdk210wMGdmfWsSn21n3G6dbKS+eFdbE8HLskytWtyLZZfk53u0G1J9mNZCcNwO3AS54OIzPnyXCuN8QfiB8bynY5Z7wQPP6T7/4/zRdcysMyHuuu74UHOzBbe/QCCGYxXeaCbyYxS8C3AJOjJ/OoiWlXxldyR9DTD6hYoOcVuh7t2bDFoSFOrRxPy7tCPrJfMKtmQK/AdjN3ijLL6hOmhqivKYbyqzcEmOtHNTRDsdeMKhkoM3RuNGJbUpcI5T4l0ka47D7XLPbVzKr0QpB4H5LlV9LhGChBvKqGXrB8ACDvV/JGSK89/6A8w2sG9TWkDyKyzae2B/gAiGHBlw".decode('base64').decode('zip')

If you really really mind the cheating with the zip and all that, the actual code is 631 bytes long (longer than the original Brainfuck compiler, but shorter than the original False compiler):
from sys import*;n=0;g=open(argv[1]).read();d=["from sys import*;from itertools import*;p=0;t=set();o=[]"];d+=["j=[bin(ord(y))[:1:-1]for y in stdin.read()];i=map(int,chain(*[x+'0'*len(x)-8for x in j]))"if","in g else""]
for c in g:d+=[n*" "+("t^={p}"if"+"==c else"t.add(p)if i.pop()else t.remove(p)"if","==c else"o+=[p in t]"if";"==c else"p-=1"if"<"==c else"p+=1"if">"==c else"while p in t:"if"["==c else"")];n+=4 if"["==c else -4 if"]"==c else 0
d+=["o=''.join(map(str,map(int,o)));o=[o[f:f+8] for f in range(0,len(o),8)];o+=[o.pop()+'0'*(8-len(o[-1]))];print''.join(chr(int(q[::-1],2))for q in o)"];open(*'ow').write('\n'.join(d))

 This is possibly the fugliest python I've ever written, even though the open(*'ow') instead of open('o','w') made me feel really smart ;)

With a (slightly) better layout:
import sys

n = 0
g = open(sys.argv[1]).read()
d = ["from sys import*;from itertools import*;p=0;t=set();o=[]"]
d += [
    "j=[bin(ord(y))[:1:-1]for y in stdin.read()];i=map(int,chain(*["
    "x+'0'*len(x)-8for x in j]))" if "," in g else ""]
for c in g:
    d += [n * " " + (
        "t^={p}" if "+" == c
        else "t.add(p)if i.pop()else t.remove(p)" if "," == c
        else "o+=[p in t]" if ";" == c
        else "p-=1" if "<" == c
        else "p+=1" if ">" == c
        else "while p in t:" if "[" == c
        else "")]
    n += 1 if "[" == c else -1 if "]" else 0
d += [
    "o=''.join(map(str,map(int,o)));o=[o[f:f+8] for f in range(0,len(o),"
    "8)];o+=[o.pop()+'0'*(8-len(o[-1]))];print''.join(chr(int(q[::-1],"
    "2))for q in o)"]

open(*'ow').write('\n'.join(d))

And lastly, the example output for the "Hello World" test you gave:
from sys import*;from itertools import*;p=0;t=set();o=[]

o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]

t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]

o+=[p in t]
o+=[p in t]
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]

t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]

t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]
o+=[p in t]
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o+=[p in t]

o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
t^={p}
o+=[p in t]
o=''.join(map(str,map(int,o)));o=[o[f:f+8] for f in range(0,len(o),8)];o+=[o.pop()+'0'*(8-len(o[-1]))];print''.join(chr(int(q[::-1],2))for q in o)

A memfiller, +[>+] also works really well (and quick, consuming about 2 gigs in 10 seconds):
from sys import*;from itertools import*;p=0;t=set();o=[]

t^={p}
while p in t:
    p+=1
    t^={p}

o=''.join(map(str,map(int,o)));o=[o[f:f+8] for f in range(0,len(o),8)];o+=[o.pop()+'0'*(8-len(o[-1]))];print''.join(chr(int(q[::-1],2))for q in o)

Perl, 325 bytes
undef$/;open C,'<',pop;@I=split//,unpack'b*',<>;@C=split//,<C>;while($c<@C){$_=$C[$c];$D=\$D{$d};$$D=!$$D if/\+/;$$D=shift@I if/,/;/</&&$d--;/>/&&$d++;$O.=$$D?1:0if/;/;if(/\[/&&!$$D){$L=1;while($L){$L++if$C[++$c]eq'[';$L--if$C[$c]eq']'}}if(/]/){$L=1;while($L){$L++if$C[--$c]eq']';$L--if$C[$c]eq'['}$c--}$c++}print pack'b*',$O

The script boofuck.pl is an interpreter for Boolfuck. The argument is a file name that contains a Boolfuck program. It reads STDIN and outputs to STDOUT.
Remarks:

I have ignored the last sentence of the description for operator ;, because I haven't understood it:

If the end-of-file character has been input, outputs a zero to the
  bit under the pointer.

The script always reads STDIN, even if the boolfuck program does not contain a , operator.

Ungolfed:
# name conventions:
# "I" for input:
#         @I is array with "0" and "1" for the bits in STDIN
# "O" for output:
#         $O is output string with "O" and "1" for the output bits
# "C" for code:
#         @C is array of code byts
#         $c is code pointer
# "D" for data:
#         %D is "infinite" data band, hash is used to allow negative indexes
#         $d is data pointer
undef $/;                    # without input line separator <> reads the whole file.
open C, '<', pop;            # open program file, given as command line argument
@I = split //,             #/# split input bit string into array of (bit) characters
     unpack 'b*',            # convert input string into bit string with 0 and 1
     <>;                     # get STDIN
@C = split //,             #/# split program code string in byte array
     <C>;                    # get program code as string
while ($c < @C) {            # finished if the end of the program code is reached
    $_ = $C[$c];             # $_ contains operator
    $D = \$D{$d};            # $D is a reference to the data bit under the data pointer
# case '+':
    $$D = !$$D if /\+/;      # flip data bit value under the data pointer
# case ',':
    $$D = shift @I if /,/; #/# store input bit to data bit under the pointer
# case ';':
    $O .= $$D ? 1 : 0 if /;/;
                           #/# append current data bit to the output
# case '<':
    /</ && $d--;           #/# move data pointer to the left
# case '>':
    />/ && $d++;           #/# move data pointer to the right
# case '[':
    if (/\[/ && !$$D) {      # if current data bit is unset,
                             # jump to the matching ']'
        $L = 1;              # each '[' increments $L,
                             # each ']' decrements $L
        while ($L) {         # $L is zero in case of a matching ']'
            $L++ if $C[++$c] eq '[';
            $L-- if $C[$c]   eq ']'
        }
        # at the end of the outer loop, the code pointer is automatically increased
        # to get right after the matching ']'
    }
# case ']':
    if (/]/) {               # jump back to the matching '['
        $L = 1;
        while ($L) {
            $L++ if $C[--$c] eq ']';
            $L-- if $C[$c]   eq '['
        }
        $c--                 # at the end of the outer loop the code pointer
                             # is increased; therefore the code pointer is
                             # decreased to get right at the matching '['
    }
# else:
    # ignore other characters
# continue:
    $c++                     # increment code pointer
}
print pack'b*', $O           # convert output bits to bytes and print result

Examples:

"Hello World" example from the question
perl boolfuck.pl helloworld.boolfuck </dev/null

Result (with new line at the end):
Hello, world!

Example that reverses the input from the Boolfuck website:
>,>,>,>,>,>,>,>,<<<<<<<<
>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]
>>>>>>>>>
>,>,>,>,>,>,>,>,<<<<<<<<
>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]
<<<<<<<<<
>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>[+<<<<<<<<[>]+<[+<]
>;>;>;>;>;>;>;>;<<<<<<<<
<<<<<<<<<
>>>>>>>>>+<<<<<<<<+[>+]<[<]>>>>>>>>>]<[+<]

Running it:
echo -ne "\nHello World" perl boolfuck.pl reverse.boolfuck

Result (with new line at the end):
dlroW olleH

I also don't understand the ;.

source: Stack Exchange

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

sort: votes