title: Four Squares Together

text:
# Question
Title: 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

```

# Answer
# 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

1. 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...
2. 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

```
3. 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

```
4. 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+`

```
5. 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

1. 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`.
2. Set `j = p**a`.
3. Set `k = floor(sqrt(N - j**2) / A) * A`.
4. Set `l = floor(sqrt(N - j**2 - k**2) / A) * A`.
5. Set `m = floor(sqrt(N - j**2 - k**2 - l**2) / A) * A`.
6. 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

```
1. At only 4 iterations and 6.6 **micro**seconds per integer, primo's algorithm is incredibly fast.
2. 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.
3. 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.
4. 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.
5. 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.
> 1 votes

# Answer
# 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.
> 7 votes

# Answer
# 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

```
> 7 votes

# Answer
## 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 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.
> 7 votes

# Answer
# 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.
> 6 votes

# Answer
# 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:

1. For input `m`, find `N` and `W` such that `m = 4**W * N`.
2. Set `i = 257**2 * floor(sqrt(N/4))`.
3. Set `i = i + 1`.
4. Find integers `j, k, l` such that `i = 257**2 * j + 257 * k + l`, where `k, l < 257`.
5. Check if `d = N - j**2 - k**2 - l**2` is a perfect square.
6. If it isn't, and go back to step 3.
7. 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”. ";

```
> 5 votes

# Answer
# 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;}}

```
> 2 votes

# Answer
# 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:

1. Express the number as `4**i * j`.
2. Find the largest integer `k` such that `k**2 <= j` and `j - k**2` satisfies the hypothesis of Legendre's three-square theorem.
3. Set `i = 0`.
4. Check if `j - k**2 - (i / 252)**2 - (i % 252)**2` is a perfect square.
5. 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.

```
> 2 votes

# Answer
# 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 */
    }
  }
}

```
> 1 votes

# Answer
# 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

```
> 1 votes

# Answer
# ~~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.
> 1 votes

# Answer
# 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.
> 1 votes

# Answer
# ~~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.
> 1 votes

# Answer
# 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 `x`will 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

```
> 1 votes

# Answer
# 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!
> 0 votes

# Answer
# 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.
> 0 votes

---
Tags: code-golf, number-theory, restricted-time
---


# Question
Title: Find the nth number where the digit sum equals the number of factors
(This is OEIS A057531.)

### Your task

Given a positive integer, \\$n\\$, find the \\$n\\$th number where the digit sum equals the number of factors

### Explanation

For example, let's take 22:

Its factors are \\$\[1, 2, 11, 22\]\\$ (length: 4).

Its digit sum is 4.

This means that it is a number where the digit sum equals the number of factors.

### The series

The first few terms of this series are:

\\$\[1, 2, 11, 22, 36, 84, 101, 152, 156, 170\]\\$

### Test cases

Note: these are 1-indexed. You may use 0-indexing.

```
Input  Output
1      1
2      2
3      11
4      22
5      36
10     170
20     444
30     828
40     1111
50     1548
100    3588

```

### Clarifications

* **You may use either 0-indexing or 1-indexing**
* The sequence starts from 1, not from 0
* The factors of a number include 1 and the number itself
* **Default rules apply** \- you may output the first \\$n\\$ terms, or the infinite sequence, or something else
* **This is , so shortest answer in bytes wins!**

# Answer
# Vyxal, 7 bytes

```
KLn∑=)ȯ

```

Try it Online!

```
      ȯ # First n integers
-----)  # Where...
 L      # Length of
K       # Factor count
    =   # Equals
  n∑    # Digit sum

```
> 3 votes

# Answer
# Factor \+ `lists.lazy math.unicode project-euler.common`, 54 bytes

```
[ 1 lfrom [ dup >dec 48 v-n Σ swap tau* = ] lfilter ]

```

Try it online!

Returns an infinite lazy list of the sequence.

* `1 lfrom` The natural numbers
* `[ ... ] lfilter` Take the ones where...
* `>dec 48 v-n Σ` Digit sum
* `tau*` Divisor count
* `=` Are equal
> 6 votes

# Answer
# Raku, 39 bytes

```
grep {.comb.sum==grep $_%%*,1..$_},1..*

```

Try it online!

This is an expression for the lazy, infinite sequence of values.
> 6 votes

# Answer
# PARI/GP, 44 bytes

```
for(i=1,oo,sumdigits(i)-numdiv(i)||print(i))

```

Attempt This Online!

Prints the sequence forever.
> 5 votes

# Answer
# Jelly, 8 bytes

```
1DS=ÆdƲ#

```

Try it online!

Takes \\$n\\$ on STDIN, and outputs the first \\$n\\$ values.

## How it works

```
1DS=ÆdƲ# - Main link. Takes no arguments
      Ʋ  - Group the previous 4 links into a monad f(k):
 D       -   Digits of k
  S      -   Sum
    Æd   -   Divisor count of k
   =     -   Are they equal?
1      # - Read n from STDIN and find the first n k≥1 such that f(k) is true

```
> 4 votes

# Answer
# JavaScript (V8),  ~~60 59~~  58 bytes

Prints the sequence 'forever' (that is, until the call stack overflows).

```
for(n=0;g=d=>d&&!(n%d--)-~~(n+g)[d]+g(d);)g(++n)||print(n)

```

Try it online!

### Commented

We use a recursive function which simultaneously counts the number of divisors and subtracts the sum of the digits.

```
for(           // infinite loop:
  n = 0;       //   start with n = 0
  g = d =>     //   g is a recursive function taking a divisor d
  d &&         //   stop if d = 0
  !(n % d--) - //   add 1 if d is a divisor of n
  ~~(n + g)    //   coerce n to a string with some additional
               //   non-digit characters
  [d] +        //   subtract the d-th digit (0 if out of bounds)
  g(d);        //   add the result of a recursive call
)              //
g(++n)         // increment n; initial call to g with d = n
|| print(n)    // print n if the result is 0

```

---

# Python 3, 78 bytes

*-1 thanks to @Mukundan314*

A port of my initial JS code. Most probably sub-optimal.

```
f=lambda d,v:d and(n%d<1)-v%10+f(d-1,v//10)
n=1
while 1:f(n,n)or print(n);n+=1

```

Try it online!
> 4 votes

# Answer
# Ruby, ~~55~~ 51 bytes

-4 bytes thanks to G B

```
1.step{|a|a.digits.sum==(1..a).count{a%_1<1}&&p(a)}

```

Attempt This Online!
> 4 votes

# Answer
# Wolfram Language (Mathematica), 66 bytes

```
(k=t=0;While[k<#,If[Tr[1^Divisors@++t]==Tr@IntegerDigits@t,k++]])&

```

Try it online!
> 4 votes

# Answer
# Japt, 11 bytes

```
Èâ Ê¶Xìx}iU

```

Try it

```
Èâ Ê¶Xìx}iU     :Implicit input of integer U
È               :Function taking an integer X as argument
 â              :  Divisors
   Ê            :  Length
    ¶           :  Equal to
     Xì         :    Digit array
       x        :    Reduced by addition
        }       :End function
         iU     :Get the Uth integer >=0 that returns true

```
> 4 votes

# Answer
# Python 3, 69 bytes

```
n=1
while sum(n%~d//n+n//10**-~d*9for d in range(n))or[print(n)]:n+=1

```

Try it online!

Would print forever if I hadn't doctored the print function in the preamble.

Based on @Arnauld's answer.

This uses two additional little tricks:

1. sum d=1..n n%-d // n equals number of divisors - n
2. sum d=1.. n // 10d x 9 equals n - digit sum
> 4 votes

# Answer
# Python, ~~80~~ 77 bytes

```
i=1
while[sum(map(int,str(i)))-sum(i%-~j<1for j in range(i))or print(i)]:i+=1

```

Attempt This Online!

-3 bytes thanks to @Jonathan Allan
> 3 votes

# Answer
# Charcoal, 27 bytes

```
ＮθＷ‹№υ¹θ⊞υ⁼ΣＬυΣＥυ¬﹪Ｌυ⊕λＩ⊖Ｌυ

```

Try it online! Link is to verbose version of code. Outputs the `n`th term. Explanation:

```
Ｎθ

```

Input `n`.

```
Ｗ‹№υ¹θ

```

Repeat until `n` terms have been found.

```
⊞υ⁼ΣＬυΣＥυ¬﹪Ｌυ⊕λ

```

Determine whether the current length of the predefined empty list is a member of the sequence, and push the result to the predefined empty list. Note that this has been written in such a way that it thinks `0` is not a member of the sequence, although an alternative interpretation could have included it.

```
Ｉ⊖Ｌυ

```

Output the index of the last element in the list, which is therefore the `n`th term. (At a cost of 1 byte, this could be `Ｉ⌕Ａυ¹` to output the first `n` terms.)
> 3 votes

# Answer
# Retina, 68 bytes

```
K`
"$+"{`^
_
{`^
_$'¶$.'¶
\d
*
w`\b(_+)(?=\1*$)
_
_(¶_+)\1$

))0G`
_

```

Try it online! No test suite due to the way the program uses history. Explanation:

```
K`

```

Clear the buffer.

```
"$+"{`
)`

```

Repeat the original input number times.

```
^
_

```

Increment the value in the buffer.

```
{`
)`

```

Repeat until the buffer stops changing.

```
^
_$'¶$.'¶

```

Make three copies of the buffer; increment the first, and convert the second to decimal.

```
\d
*

```

Convert each digit separately back to unary, thus computing the sum of the digits.

```
w`\b(_+)(?=\1*$)
_

```

Count the number of factors in the current value.

```
_(¶_+)\1$

```

If the digit sum equals the factor count then decrement the first value again, causing the inner loop to exit.

```
0G`

```

Keep only the first value.

```
_

```

Convert the final value to decimal.
> 3 votes

# Answer
# Excel (ms365), 170 bytes

```
=INDEX(REDUCE(0,ROW(1:1048576),LAMBDA(a,b,IF(COUNT(REDUCE(0,SEQUENCE(b),LAMBDA(x,y,IF(MOD(b,y),x,VSTACK(x,y)))))-1=SUM(--MID(b,SEQUENCE(LEN(b)),1)),VSTACK(a,b),a))),A1+1)

```

Here `ROW(1:1048576)` will represent `n` and is the extend to which Excel can go (max rows). You'll notice that calculation wise the nesting will likely not compute but to be fair to the challenge of any possible positive integer I had to include this number. 'Much easier' to compute and able to calculate any of the given samples would be `ROW(1:3588)`:

What happens is:

* The outer `REDUCE()` will cycle through all possible integers in Excel-rows. It will `VSTACK()` any integer which follows the following rule;
* The inner `REDUCE()` will `VSTACK()` all integers in iteration if they equal zero when `MOD()` is applied **and** the `SUM()` of the numbers used in the integer equals the `COUNT()` of these stacked integers (I hope that is clear);
* `INDEX()` will then return a zero-based (hence the `+1)` integer from the `VSTACK()`. If you wish to return an array of the 1st n integers:

```
=DROP(TAKE(REDUCE(0,ROW(1:500),LAMBDA(a,b,IF(COUNT(REDUCE(0,SEQUENCE(b),LAMBDA(x,y,IF(MOD(b,y),x,VSTACK(x,y)))))-1=SUM(--MID(b,SEQUENCE(LEN(b)),1)),VSTACK(a,b),a))),A1+1),1)

```
> 3 votes

# Answer
# x86-64 machine code, 39 bytes

```
31 C0 FF C0 31 F6 6A 0A 59 50 99 F7 F9 29 D6 01 C2 75 F7 58 89 C1 50 F7 F9 58 83 EA 01 99 11 D6 E2 F4 19 D7 79 DC C3

```

Try it online!

Following the standard calling convention for Unix-like systems (from the System V AMD64 ABI), this takes a number `n` in EDI and returns the 0-indexed `n`th number where the digit sum equals the number of factors in EAX.

In assembly:

```
f:  xor eax, eax    # Set EAX to 0. EAX will be the number currently being tested.
nextnum:
    inc eax         # Increase EAX by 1.
    xor esi, esi    # Set ESI to 0.
    push 10; pop rcx# Set RCX (the 64-bit register containing ECX) to 10.
    push rax        # Save the value of RAX (containing EAX) onto the stack.
digitloop:
    cdq             # Set EDX to 0 by sign-extending EAX.
    idiv ecx        # Divide EDX:EAX by ECX (10).
                    # The quotient goes in EAX and the remainder goes in EDX.
    sub esi, edx    # Subtract EDX (the remainder; a digit of the number) from ESI.
    add edx, eax    # Add EAX to EDX.
    jnz digitloop   # Jump back if the result of the addition is nonzero.
            # This loop will exit when the quotient and remainder are both zero.
            # This performs one more iteration than otherwise necessary,
            # so that EDX will be 0 after the loop, which will be useful later.
    pop rax         # Restore the value of RAX from the stack.
    mov ecx, eax    # Set ECX to the same value, the number being tested.
                    #  (The higher bits of RCX are automatically zeroed.)
factorloop:
    push rax        # Save the value of RAX onto the stack.
    idiv ecx        # Divide EDX:EAX by ECX; quotient in EAX, remainder in EDX.
    pop rax         # Restore the value of RAX from the stack.
    sub edx, 1      # Subtract 1 from EDX. The carry flag CF becomes 1 iff EDX=0.
    cdq             # Set EDX to 0 by sign-extending EAX.
    adc esi, edx    # Add EDX+CF to ESI, increasing it by 1 with each factor found.
    loop factorloop # Decrease RCX by 1, and jump back if it's nonzero.
    sbb edi, edx    # Subtract EDX+CF from EDI.
                    # 1 is the last value checked, and it is always a factor.
                    # Thus, the last ADC added 1 to ESI, and CF=1 iff ESI went
                    #  from -1 to 0; iff the number satisfies the condition.
    jns nextnum     # Jump back if the result is nonnegative.
    ret             # Return.

```
> 3 votes

# Answer
# 05AB1E, 8 bytes

```
∞ʒÑgySOQ

```

Outputs the infinite sequence.

Try it online.

**Explanation:**

```
∞       # Push an infinite positive list: [1,2,3,...]
 ʒ      # Filter it by:
  Ñ     #  Get a list of divisors of the current integer
   g    #  Pop and push its length
  yS    #  Push the current integer again, and convert it to a list of digits
    O   #  Take the sum of those digits
     Q  #  Check if the amount of digits equals the digit-sum
        # (after which the filtered infinite list is output implicitly as result)

```
> 2 votes

# Answer
# MathGolf, 10 bytes

```
î·─£\Σ=╛p▲

```

Outputs the infinite sequence, with each term on a separated line.

Try it online.

**Explanation:**

```
         ▲  # Do-while true with pop:
î           #  Push the 1-based loop index
 ·          #  Quadruplicate it
  ─         #  Pop one, and get a list of its divisors
   £        #  Pop this list, and get its length
    \       #  Swap to put another 1-based index at the top of the stack
     Σ      #  Sum its digits
      =     #  Check if the amount of divisors and sum of digits are equal
       ╛    #  If this is truthy:
        p   #   Print a 1-based index with trailing newline

```
> 2 votes

# Answer
# , 53 bytes

```
while(F<-F+1)sum(F%/%10^(0:F)%%10,-!F%%1:F)||print(F)

```

Try it online!

Prints the sequence indefinitely (in theory). In practice, results over 308 are inaccurate as then `10^(0:F)` exceeds computing capabilities of R. We could counter this by using eg. `10^(0:nchar(F))`, but this costs bytes...

### Explanation:

```
while(F<-F+1)          # while increasing value of F (infinitely)
sum(F%/%10^(0:F)%%10,  # sum the digits of F (with a lot of unnecessary 0s)
-!F%%1:F)              # minus the count of divisors of F
||print(F)             # if the result is 0, print F

```
> 2 votes

# Answer
# Brachylog, 10 bytes

```
ℕ₁≜.ẹ+~l~f

```

Try it online!

This is a generator. The TIO link has a header to ask for the first 10 possible instantiations of the output.

### Explanation

```
ℕ₁≜.          The output is an instantiated natural integer
    .ẹ+        The sum of elements (i.e. digits) of the output…
       ~l      …is the length of a list…
         ~f    …which is the list of factors of the output

```
> 2 votes

# Answer
# Nibbles, 9 bytes (18 nibbles)

```
|,~~`^+`@~$,|,$~%@

```

```
|                   # filter
 ,~                 # 1..infinity
   ~                # for each number n that is falsy for
    `^              # bitwise xor of
      +             #   sum of
       `@ $         #   digits of n in base
         ~          #   10 (default)
                    # and
           ,        #   length of
             ,$     #   1..n
            |       #   filtered for each i that is
               ~    #   falsy for
                %@  #   n modulo i

```
> 2 votes

# Answer
# Excel, ~~112~~ 110 bytes

`LET(r,ROW(1:999),REDUCE("",r,LAMBDA(a,b,a&IF(SUM(MID(b,SEQUENCE(LEN(b)),1)+0)=SUM((MOD(b,r)<1)+0),b&",",""))))`
> 2 votes

# Answer
# Pyth, 15 bytes

```
.fqsjZTl{m*FdyP

```

Try it online!

Outputs the first \\$n\\$ terms as a list.

##### Explanation

```
.fqsjZTl{m*FdyP   | Full program
.fqsjZTl{m*FdyPZQ | with implicit variables
------------------+---------------------------------------------------------------
.f              Q | Find the first <input> integers Z (starting with 1) such that:
   s              |   the sum of
    jZT           |    the base 10 digits of Z
  q               |  equals
       l          |   the length of
                  |    the list of factors of Z, i.e.:
        {         |    the unique
         m*Fd     |     products of
             y    |      all subsets of
              PZ  |       the prime factors of Z

```
> 2 votes

# Answer
# PLIS `-i`, 78 bytes

```
1+@(H@({$0*$0}{([$0+$1]{$0%10}[$0/10](M*$0))@(spT@$0)-[$0+$1](H@($0%a))@$0}a))

```

The `-i` flag generates sequence entries infinitely, as such:

```
C:\Users\conorob\Programming\langs\PLIS (main -> origin)
λ plis -i -f example\digit-sum-equal-factor-count.min.plis
1 2 11 22 36 84 101 152 156 170 202 208 225 228 288 301 372 396 441 444 468 516 525 530
602 684 710 732 804 828 882 952 972 1003 1016 1034 1070 1072 1106 1111 1164 1236 ^C

```

PLIS is a programming language which primarily manipulates data using relevant OEIS sequences. While it has a lot of OEIS sequences implemented, PLIS does not aim to trivialize CGCC challenges by implementing random challenge sequences, hence the length.

Because OEIS sequences are the primary functions in this language, whenever a sequence of word characters are encountered, they are interpreted as a base 53 integer encoded using the alphabetic characters and `_`. (The OEIS sequence quoted in the OP, A057531, would be represented as `_ZU`, were it actually implemented in the language.)

Here is a commented, ungolfed explanation of what the code generally does.

```
# add 1 to each
1 +
# find indices where the sequence is nonzero
@(
  A000007@(
    # square each difference
    {$0*$0}
      # over each positive integer
      {
        # digit sum
        ([$0+$1] {$0%10} [$0/10] (A000012*$0))@(A055642@$0)
        -
        # factor count
        [$0+$1](A7@($0%A000027))@$0
      } A000027)
  )

```
> 2 votes

# Answer
# C (gcc), ~~84~~ 83 bytes

```
j,c;f(i){for(j=c=1;j<i;)c+=i%j++<1;for(;j;j/=10)c-=j%10;c||printf("%d ",i);f(i+1);}

```

Try it online!
> 1 votes

# Answer
# Fig, \\$10\log\_{256}(96)\approx\\$ 8.231 bytes

```
FmC'=SxLJk

```

Try it online!

Returns/prints the infinite list of such numbers.

```
FmC'=SxLJk
F          # Filter
 mC        # The infinite list of numbers
   '       # Where
     Sx    # The digit sum
    =      # Equals
       L   # The length
        Jk # Of its divisors + itself

```
> 1 votes

# Answer
# Pip Classic, ~~35~~ 18 bytes

```
W++xI$+x=0Nx%\,xPx

```

Try it online!

A lot of bytes saved thanks to @DLosc!
> 1 votes

# Answer
# Stax, 9 bytes

```
▐ë<of▌ëÇ@

```

Run and debug it

This is PackedStax, which unpacks to:

# Stax (unpacked), 11 bytes

```
VIf:d%_E|+=

```

Run and debug it

```
VI          # infinity
  f         # filtered by
   :d%      # the length of the divisors
      _E|+  # and the digit sum
          = # are equal

```
> 1 votes

# Answer
# Rockstar, 168 bytes

Rockstar can't do *anything* the challenge asks. So, of course, I'm going to do it anyway! On my phone!

```
listen to N
X's 0
while N
let X be+1
Y's 0
D's 0
while X-Y
let Y be+1
let M be X/Y
turn up M
let D be+M's X/Y

cut X in L
while L
let D be-roll L

let N be-D's 0

say X

```

Try it here (Code will need to be pasted in)

```
listen to N      :Read input string into variable N
X's 0            :Initialise X as 0
while N          :While N is not 0
let X be+1       :  Increment X
Y's 0            :  Initialise Y as 0
D's 0            :  Initialise D as 0
while X-Y        :  While Y is less than X
let Y be+1       :    Increment Y
let M be X/Y     :    Assign X/Y to variable M
turn up M        :    Round M up to nearest integer
let D be+        :    Increment D by
  M's X/Y        :      Is M equal to X/Y?
                 :  End while loop
cut X in L       :  Split X into a digit array and assign to variable L
while L          :  While L is not empty
let D be-        :    Decrement D by
  roll L         :      Pop the first element from L
                 :  End while loop
let N be-        :  Decrement N by
  D's 0          :    Is D equal to 0?
                 :End while loop
say X            :Output X

```
> 1 votes

# Answer
# Haskell, 73 bytes

```
s=[n|n<-[0..],length[d|d<-[1..n],0==mod n d]==(sum.map(read.pure).show)n]

```

Try it online!
> 1 votes

# Answer
# Scala, ~~83~~ 71 bytes

Saved 12 bytes thanks to the comment of @naffetS

---

Golfed version. Try it online!

```
for(a<-1 to l)if(s"$a".map(_-'0').sum==(1 to a).count(a%_<1))println(a)

```

Ungolfed version. Try it online!

```
object Main {
  def main(args: Array[String]): Unit = {
    val limit = 1000 // You can replace this with any limit you want
    for (a <- 1 to limit) {
      val sumOfDigits = a.toString.map(_.asDigit).sum
      val numFactors = (1 to a).count(a % _ == 0)
      if (sumOfDigits == numFactors) println(a)
    }
  }
}

```
> 1 votes

---
Tags: code-golf, number, sequence
---


title: Pendulum Encoding

text:
# Question
Title: Pendulum Encoding
Given a string as an input (which can be any acceptable/convenient format in your language), implement pendulum encoding. The test cases are split into individual items (which aren't quoted) for a visually appealing explanation.

## How do I do that?

The current **iteration index** starts at `0`.

* If the iteration index is even, append the current item onto the output string.
* If the iteration index is odd, prepend the current item onto the output string.

## An example

```
The input is [a b c d e f g].
Note that the letters a-g are individual one-character strings, to prevent confusion from the iteration index.
N: the iteration index

N:0 Out:      [a]
N:1 Out:    [b a]
N:2 Out:    [b a c]
N:3 Out:  [d b a c]
N:4 Out:  [d b a c e]
N:5 Out:[f d b a c e]
N:6 Out:[f d b a c e g]

```

The output should be `[f d b a c e g]`.

# Another example

```
The input is [u d l n u e m p].

N:0 Out:        [u]
N:1 Out:      [d u]
N:2 Out:      [d u l]
N:3 Out:    [n d u l]
N:4 Out:    [n d u l u]
N:5 Out:  [e n d u l u]
N:6 Out:  [e n d u l u m]
N:7 Out:[p e n d u l u m]

```

## Test cases

Here's a sample program doing this encoding.

Take note that the characters in the string aren't always unique.

```
Your output *has* to be flattened.

[a,b,c,d,e,f,g]   -> [f,d,b,a,c,e,g]
[]                -> []
[a]               -> [a]
[a,b,c,d]         -> [d,b,a,c]
[a,b]             -> [b,a]
[a,b,d]           -> [b,a,d]
[a,b,a,c,b,c]     -> [c,c,b,a,a,b]
[a,a,b,b,c,c]     -> [c,b,a,a,b,c]
[u,d,l,n,u,e,m,p] -> [p,e,n,d,u,l,u,m]

```

# Answer
# 05AB1E, 4 bytes

```
ι`Rì

```

I/O as a list of characters.

Try it online or verify all test cases.

**Explanation:**

```
ι     # Uninterleave the (implicit) input-list (into 2 parts by default for lists)
      #  i.e. ["u","d","l","n","u","e","m","p"] → [["u","l","u","m"],["d","n","e","p"]]
 `    # Push both parts separated to the stack
  R   # Reverse the second part
      #  → ["p","e","n","d"]
   ì  # And prepend it in front of the first
      #  → ["p","e","n","d","u","l","u","m"]
      # (after which the result is output implicitly)

```
> 7 votes

# Answer
# Python 3, 29 bytes

```
lambda l:l[1::2][::-1]+l[::2]

```

Try it online!

**Input**: A sequence
**Output**: The pendulum encoding of that sequence

**How**
Consider the sequence `[0,1,2,3,4,5]`, whose pendulum encoding is `[5,3,1,0,2,4]`. We can see that all even indices ended up in order on the right, and all odd indices are in reversed order on the left.

* `l[1::2][::-1]` takes all odd indices and reverses them, e.g `[5,3,1]`
* `l[::2]` takes all even indices, e.g `[0,2,4]`
> 21 votes

# Answer
# brainfuck, ~~24 21~~ 18 bytes

```
,[[<],[>],<]>>[.>]

```

Try it online!

Thanks to Jo King for -3 bytes

```
,[          while input
  [<],      add new character to start of memory
  [>],      add new character to end of memory
  <         go one back, so the loop will run another time, moving the pointer to the start of memory
]
>>[.>]      print memory

```
> 16 votes

# Answer
# APL (Dyalog Unicode), 11 bytes

I promised you'll see at least one interesting answer :)

```
{⍵[⍋-\⍳≢⍵]}

```

Try it online!

Uses my own tip about `-\⍳`, specifically the  variation, to generate the permutation needed for this challenge.

### How it works

`⍋-\⍳≢⍵` generates the target permutation for both even- and odd-length arrays:

```
⍋-\⍳≢⍵  ⍝ Length-7 vector  | Length-8 vector
    ≢   ⍝ Length
        ⍝ 7                | 8
   ⍳    ⍝ Range (1..n)
        ⍝ 1 2 3 4 5 6 7    | 1 2 3 4 5 6 7 8
 -\     ⍝ Cumulative alternating difference
        ⍝ 1 -1 2 -2 3 -3 4 | 1 -1 2 -2 3 -3 4 -4
⍋       ⍝ Grade up; permutation that will sort the input array
        ⍝ 6 4 2 1 3 5 7    | 8 6 4 2 1 3 5 7

```

Then `⍵[...]` arranges the original elements in that particular order.
> 9 votes

# Answer
# JavaScript, 36 bytes

Feels like forever since I posted a JS solution here!

Input as an array, output as a string. Handling the empty array cost 3 bytes.

```
a=>a.reduce((x,y,z)=>z%2?y+x:x+y,"")

```

Try it online!
> 7 votes

# Answer
# Haskell, ~~53~~ 37 bytes

```
([]#)
p#(a:b:s)=(b:p++[a])#s
p#l=p++l

```

Try it online!

*Golfed 16 bytes thanks to @xnor*
> 6 votes

# Answer
# JavaScript (ES6), 42 bytes

A recursive function taking and returning a string.

```
f=([c,...a],k,o='')=>c?f(a,!k,k?c+o:o+c):o

```

Try it online!
> 4 votes

# Answer
# Jelly,  ~~6~~  5 bytes

-1 thanks to Unrelated String! (Use an initial value for reduction.)

```
;ṭƭƒ“

```

**Try it online!**

A full program that accepts a string and prints the result.

### How?

```
;ṭƭƒ“ - Main Link: list of characters, S
    “ - get an empty list -> X
   ƒ  - reduce (X + S) using:
  ƭ   -   tie - this groups the last n (default 2) links together and calls them cyclically:
;     -     concatenate  (i.e. on the first, third, fifth,... we append)
 ṭ    -     tack         (     while on the second, fourth,... we prepend)
      - implicit, smashing print to STDOUT

```

It could be considered a bug that Jelly barfs if the non-initialised version of reduce, `/`, is provided with an empty list. If that were 'fixed' this becomes just the four-byte `ṭ;ƭ/`.
> 4 votes

# Answer
# Keg, `-lp`, `-ir`, `-hd`, 15 bytes

```
⑫&(&⑶+&)⒁2%[⒂|&

```

Try it online!

Purely a literal interpretation of the question plus a bit to ensure even length strings don't break.

## Explained

```
⑫&

```

First, we start by storing an empty string (pushed by ) in the register. This will be used as the final output, meaning that it needs to be initalised.

```
(&⑶+&)

```

Then, we enter the main for loop, which has the implicit condition `!` (take length of stack), as no explicit condition is provided. Now, there isn't anything on the stack at this point, so doing such a thing may seem pointless. But, by using the `-lp` flag (`--lengthpops`), we can have the `!` command take input if the stack is empty and push the length of the input. Also, the `-ir` (`--inputraw`) command ensures that the input word is taken as a series of letters, rather than a single string.

Inside the for loop, we push the contents of the register, reverse it ( reverses the top item of the stack) and then add whatever is next to the register. By doing so, we achieve the process behind the main algorithm, as consecutive letters are appended in the desired order.

```
⒁2%

```

At this point, the encoding has been fully completed. However, if the string is of an even length, the result will be reversed in the register. This requires us to push the length of the register () and check to see if it is even.

```
[⒂|&

```

If the register is of odd length, then we reverse the register and push it onto the stack (). Otherwise, we simply push the register. `-hd` will then ensure that only the top item on the stack is printed.
> 3 votes

# Answer
# Red, 52 bytes

```
func[a][append reverse extract next a 2 extract a 2]

```

Try it online!
> 3 votes

# Answer
# , 13 bytes

```
{~[:/:[:-/\#\

```

Try it online!

A  port of Bubbler's APL solution \- don't forget to upvote his answer!

```
           #\  length of successive prefixes
        -/\    cumulative alternating difference
      [:       function composition (caps the previous two verbs as a fork)
    /:         grade up
  [:           caps the fork
{~             use the list to index into the input (arguments reversed)

```
> 3 votes

# Answer
# Perl 5 `+p`, 30 bytes

```
$\=--$|%2?$\.$_:$_.$\for/./g}{

```

Try it online!
> 3 votes

# Answer
# Vyxal `r`, 3 bytes

```
yṘJ

```

Uninterleave, reverse the second argument, merge b + a. Requires the r flag because there's no way to merge to my knowledge the other way around in one command.

Try it Online!

And flagless from @AaronMiller if the input is taken as a string:

## Vyxal, 3 bytes

```
yṘp

```

Instead of J to merge, it prepends the element.

Try it Online!
> 3 votes

# Answer
# Japt `-P`, ~~7~~ ~~6~~ 5 bytes

I/O as a string.

```
ó ÔvÔ

```

Try it

```
ó ÔvÔ     :Implicit input of string    > "udlnuemp"
ó         :Uninterleave                > ["ulum","dnep"]
  Ô       :Reverse                     > ["dnep","ulum"]
   vÔ     :Reverse first element       > ["pend","ulum"]
          :Implicitly join and output  > "pendulum"

```
> 2 votes

# Answer
# Zsh, 26 bytes

```
for x y;a=($y $a $x)
<<<$a

```

Try it online!

Input is a list of characters.
> 2 votes

# Answer
# Retina, 20 bytes

```
*>^`.(.)?
$1
1,2,`.

```

Try it online! Link includes test suite. Explanation:

```
*>`

```

Execute this stage and immediately output the result without actually changing the working string.

```
^`

```

Before making the replacements, reverse their order.

```
.(.)?
$1

```

Keep only alternate characters.

```
1,2,`.

```

Delete alternate characters.

29 bytes as a reusable function:

```
,V,2,`
O^$`.((?=(..)*$))?
$#1

```

Try it online! Link includes test cases. Explanation:

```
,V,2,`

```

Reverse all the characters except alternate characters.

```
O^$`.((?=(..)*$))?
$#1

```

Sort all characters by the parity of their position from the end, and then reverse the result. This means that alternate characters end up reversed and sorted to the start, leaving the remaining characters at the end, although technically having been reversed twice.
> 2 votes

# Answer
# Charcoal, 13 bytes

```
Ｆ²«Ｐ✂θιＬθ²←↷⁴

```

Try it online! Link is to verbose version of code. Explanation:

```
Ｆ²«

```

Loop over the even and alternate characters.

```
Ｐ✂θιＬθ²

```

Print the current set of characters.

```
←↷⁴

```

Prepare to print the alternate characters backwards.
> 2 votes

# Answer
# Java (JDK), 60 bytes

```
a->{var r="";for(var c:a)r=r.length()%2<1?r+c:c+r;return r;}

```

Try it online!

## Alternative with Streams (62 bytes)

```
s->{int[]x={0};return s.reduce("",(a,b)->x[0]++%2<1?a+b:b+a);}

```

Try it online!
> 2 votes

# Answer
# Ruby, 42 bytes

```
f=->s,*w{s ?f[s[2..-1],s[1],*w,s[0]]:w*''}

```

Try it online!
> 2 votes

# Answer
# PHP, ~~49~~ 48 bytes

```
for(;$a=$argv[++$i];)$s=$i%2?$s.$a:$a.$s;echo$s;

```

Try it online!

Again, not a really great score for PHP..

EDIT: Thanks for @OlivierGrégoire for saving 1 byte
> 2 votes

# Answer
# Icon, 62 bytes

```
procedure f(s)
t:=""
t[k:=|(0|1)\*s:k]:=pop(s)&\z
return t
end

```

Try it online!

Takes the input as a list of chars.

# Icon, 66 bytes

```
procedure f(s)
t:=""
i:=1to*s&t[1-i%2:1-i%2]:=s[i]&\z
return t
end

```

Try it online!

Takes the input as a string.

`procedure`, `return` and `end` add a lot to the byte count :)

Icon's `slice` operator `:` can be used to insert substrings into strings, if the the two indices are equal. This `s[1:1]:="a"` prepends `s` with "a"; `s[0:0]:="b"` appends "b" to `s`. I start with an empty string `t`, scan the input string `s` and use the odd/even index `i` with the slice operator to prepend/append to `t`.
> 2 votes

# Answer
# C# (.NET Core), 93 bytes

```
public static string P(this string a){int z=0;return a.Aggregate("",(x,y)=>z++%2>0?y+x:x+y);}

```

Try it online!
> 2 votes

# Answer
# Pip, 10 bytes

```
ORVDQUW:aa

```

Try it online!

-4 from DLosc.

## Similar solution, 14 bytes

```
UW:a(RVa@1).@a

```

Try it online!
> 2 votes

# Answer
# Rockstar, ~~135~~ ~~125~~ ~~115~~ 96 bytes

```
listen to S
cut S
X's 0
O's ""
while S
let X be not X
let O be X and O+roll S or roll S+O

say O

```

Try it here (Code will need to be pasted in)
> 1 votes

# Answer
# Scala, 51 bytes

```
_./:(0,""){case(i->a,c)=>i+1->Seq(a+c,c+a)(i%2)}._2

```

Try it online
> 1 votes

# Answer
# Perl 5 `-p`, 26 bytes

```
s/(.)(.)/$\="$2$\$1";''/ge

```

Try it online!
> 1 votes

# Answer
# Factor, 47 bytes

```
[ dup <odds> ""reverse-as swap <evens> append ]

```

Try it online!

## Explanation

It's a quotation (anonymous function) that takes a string from the data stack as input and leaves a string on the data stack as output. Assuming `"udlnuemp"` is on top of the data stack when this quotation is called...

* `dup` Duplicate whatever's on top of the stack.

**Stack:** `"udlnuemp" "udlnuemp"`
* `<odds>` Get the odd-indexed elements as a virtual sequence.

**Stack:** `"udlnuemp" T{ odds f "udlnuemp" }`
* `""` Push an empty string to the stack.

**Stack:** `"udlnuemp" T{ odds f "udlnuemp" } ""`
* `reverse-as` Take a source sequence and an exemplar sequence, and create a new sequence where the source sequence is reversed, but with the type of the exemplar sequence.

**Stack:** `"udlnuemp" "pend"`
* `swap` Swap the top two objects on the data stack.

**Stack:** `"pend" "udlnuemp"`
* `<evens>` Get the even-indexed elements as a virtual sequence.

**Stack:** `"pend" T{ evens f "udlnuemp" }`
* `append` Append two sequences. The result of `append` has the same type as the first argument.

**Stack:** `"pendulum"`
> 1 votes

# Answer
# minigolf, 36 bytes

```
0,;i,unx2%!,,__a,n_nx2%,,__xo,;o_s,_

```

Attempt This Online!

## Explanation

```
0,;        Start with an empty list
i,         Foreach item in input:
  u            Store list to acc
  n            Push curr. item
  x2%!,        If iteration index is even:
    ,_             drop TOS
  _            end
  a            Push acc
  ,n_          Dump acc onto stack
  n            Push curr. item
  x2%,         If iteration index is odd:
    ,_             drop TOS
  _            end
  xo           Iteration index + 1
  ,;o          Wrap last n items of stack into array
_            End foreach
s,_          (For some reason we have an extra implicit input,
              so swap it upwards and drop it)

```
> 1 votes

# Answer
# Python 2, 25 bytes

```
lambda s:s[-1::-2]+s[::2]

```

Try it online!

Wrote this before I took a look at the submissions. The Python 3 answer is basically the same, but I improved the slicing to reverse the odd indices all in one.

But yeah, you can construct the pendulum encoding in two parts. First, the odd indices in reverse order. Second, the even indices in normal order.
> 0 votes

# Answer
# Lua, 56 bytes

```
s=''(...):gsub('(.)(.?)',load'a,b=...s=b..s..a')print(s)

```

Try it online!
> 0 votes

---
Tags: code-golf, string, sorting
---


title: Randomize points on a disc

text:
# Question
Title: Randomize points on a disc
I read about circles somewhere, and just now learned about discs (*it's actually a pretty common concept*) and thought about codegolf.

Your task is to randomize a point/several points on a disc with the radius 1.

**Rules:**

* All points must have an equal probability to be generated
* Floating point coordinates must be used; the minimum requirement is two decimals (e.g. the points `(0.12, -0.45)` or `(0.00, -1.00)` are valid)
* You get -20 bytes if your program actually displays the bounding circle and the point(s) generated in it. Coordinates still has to be valid but not displayed, and the image generated has to be at least 201 by 201 pixels in size
* You get -5 bytes if your program takes the number of points to be generated as input on stdin
* If you decide not to plot the bounding circle and the point(s), your program has to output the point(s) generated on the format `(x, y)` or `(x,y)` on stdout
* If you decide to take the number of generated points as input, but not to plot it - your program has to output all the randomized points on the format stated above with or without one space in between

Shortest submission in bytes wins!

# Answer
# Pyth, 26 - 5 = 21 bytes

```
VQp(sJ*@OZ2^.n1*yOZ.l_1)eJ

```

Takes the number of coordinates to generate on stdin, and outputs them on stdout like such:

```
(-0.5260190768964058, -0.43631187015380823)(-0.12127959509302746, -0.08556306418467638)(-0.26813756369750996, -0.4564539715526493)

```

Uses a strategy similar to @MartinBüttner, generating polar coordinates and radii, except it does it using complex exponentiation.
> 5 votes

# Answer
# CJam, ~~28~~ 27 bytes

```
PP+mr_mc\ms]1.mrmqf*"(,)".\

```

This solution is not rejection-based. I'm generating the points in polar coordinates, but with a non-uniform distribution of the radii to achieve a uniform density of the points.

Test it here.

## Explanation

```
PP+     e# Push 2π.
mr_     e# Get a random float between 0 and 2π, make a copy.
mc\     e# Take the cosine of one copy and swap with the other.
ms]     e# Take the sine of the other copy and wrap them in an array.
        e# This gives us a uniform point on the unit circle.
1.mr    e# Get a random float between 0 and 1.
mq      e# Take the square root. This is the random radius.
f*      e# Multiply x and y by this radius.
"(,)".\ e# Put the resulting numbers in the required format.

```

Why does it work? Consider a narrow annulus of radius `r` and (small) width `dr`. The area is approximately `2π*r*dr` (if the annulus is narrow, the inner and outer circumference are almost identical, and the curvature can be ignored, such that the area can be treated as that of a rectangle with side lengths of the circumference and the width of the annulus). So the area increases linearly with the radius. This means that we also want a linear distribution of the random radii, in order to achieve a constant density (at twice the radius, there is twice as much area to fill, so we want twice as many points there).

How do we generate a linear random distribution from 0 to 1? Let's look at the discrete case first. Say, we have a desired distribution of 4 values, like `{0.1, 0.4, 0.2, 0.3}` (i.e. we want `1` to be 4 times as common as `0`, and twice as common as `2`; we want `3` three times as common as `0`):

How can pick one of the four values with the desired distribution? We can stack them up, pick a uniformly random value between 0 and 1 on the y-axis and pick the segment at that point:

There's a different way to visualise this picking though. We could instead replace each value of the distribution with the accumulation of the values up to that point:

And now we treat the top line of this chart as a function `f(x) = y` and invert it to get a function `g(y) = f<sup>-1</sup>(y) = x`, which we can apply to a uniformly random value in `y ∈ [0,1]`:

Cool, so how can make use of this to generate a linear distribution of radii? This is the distribution we want:

The first step is to accumulate the values of the distribution. But the distribution is continuous, so instead of summing over all previous values, we take an integral from `0` to `r`. We can easily solve that analytically: `<sub>0</sub><sup>r</sup> r dr = 1/2 r<sup>2</sup>`. However, we want this to be normalised, i.e. to multiply it by a constant such that this gives `1` for the maximum value of `r`, so what we really want is `r<sup>2</sup>`:

And finally, we invert this to get a function we can apply to a uniform value in `[0,1]`, which we can again do analytically: it's just `r = √y`, where `y` is the random value:

This is a fairly useful technique that can often be used to generate simple distributions exactly (it works for any distribution, but for complicated ones the last two steps may have to be solved numerically). However, I wouldn't use it in this particular case in production code, because the square root, sine and cosine are prohibitively expensive: using a rejection-based algorithm is on average much faster, because it only needs addition and multiplication.
> 16 votes

# Answer
# Mathematica, ~~68~~ 44 - 20 = 24 bytes

*Many thanks for David Carraher for letting me know about `RandomPoint`, which saved 24 (!) bytes. Mathematica* does *have a built-in for everything.*

```
Graphics@{Circle[],Point@RandomPoint@Disk[]}

```

This plots the point and the bounding circle to qualify for the bonus:

The result is a vector image, so the size specification of 201x201 pixels doesn't really make sense, but by default it does render larger than that.
> 12 votes

# Answer
# CJam, ~~31~~ 26 bytes

```
{];'({2dmr(_}2*@mhi}g',\')

```

This works by repeatedly generating random points in a square of side length 2 and keeping the first that falls inside the unit disc.

*Thanks to @MartinBüttner for golfing off 3 bytes!*

Try it online in the CJam interpreter.

### How it works

```
{                  }g       Do:
 ];'(                         Clear the stack and push a left parenthesis.
     {      }2*               Do twice:
      2dmr                      Randomly select a Double between 0 and 2.
          (_                    Subtract 1 and push a copy.
               @              Rotate the copy of the first on top of the stack.
                mh            Compute the Euclidean norm of the vector consisting
                              of the two topmost Doubles on the stack.
                  i           Cast to integer.
                            If the result is non-zero, repeat the loop.
                     ',\    Insert a comma between the Doubles.
                        ')  Push a right parenthesis.

```
> 9 votes

# Answer
# iKe, ~~53~~ 51 bytes

Nothing particularly special, but I suppose we ought to have at least one graphical solution:

```
,(80+160*t@&{.5>%(x*x)+y*y}.+t:0N 2#-.5+?9999;cga;3)

```

Try it in your browser.

Edit: I can shave off two bytes by applying @MartinBüttner's approach for modifying the distribution of polar coordinates. I think it's also quite a bit more direct:

```
,(80*1+(%?c){x*(cos y;sin y)}'6.282*?c:9999;cga;3)

```
> 8 votes

# Answer
# Perl, 59 Bytes

```
while(($x=1-rand 2)**2+($y=1-rand 2)**2>1){};print"($x,$y)"

```

This is just a simple solution, generating points in a square and rejecting those too far away. My singular golfing trick is to include the assignments inside of the condition.

Edit: In the process of golfing, I found an interesting way to print random points on a *circle*.

```
use Math::Trig;$_=rand 2*pi;print"(",sin,",",cos,")"

```
> 7 votes

# Answer
# Octave, ~~24~~ 53 - 20 = 33 bytes

```
polar([0:2e-3:1,rand]*2*pi,[ones(1,501),rand^.5],'.')

```

Generates 501 equally-spaced theta values plus one random number and scales them all to \[0..2π\]. Then generates 501 1's for the radius of the circle, plus a random radius for the point and takes the square root to ensure uniform distribution over the disc. Then plots all the points as polar coordinates.

---

Here is a quick demonstration of the distribution (without the unit circle):

```
polar(2*pi*rand(99),rand(99).^.5,'.')

```
> 7 votes

# Answer
## Octave / Matlab, ~~74~~ 64 bytes

**Rejection method**, 64 bytes:

```
u=1;v=1;while u^2+v^2>1
u=rand;v=rand;end
sprintf('(%f,%f)',u,v)

```

**Direct method**, 74 bytes (thanks to Martin Büttner for helping me correct two errors):

```
t=rand*2*pi;r=1-abs(1-sum(rand(2,1)));sprintf('(%f,%f)',r*cos(t),r*sin(t))

```
> 5 votes

# Answer
# R, ~~99~~ ~~95~~ 81-20 = ~~79~~ ~~75~~ 61 Bytes

```
symbols(0,0,1,i=F,asp=1,ylim=c(-1,1));points(complex(,,,runif(9),runif(9,-1)*pi))

```

Use the complex number construction to build the x/y's from polar coordinates. Taking the input was a bit expensive and there is probably a better way to do this. The `ylim` ~~and `xlim`~~ is to ensure the entire circle is plotted and the `asp` ensures the points are shown under the circle symbol.

Thanks to @jbaums and @flodel for the savings

Try it here
> 5 votes

# Answer
# Processing /Java 141 bytes-20=121

the requirement for 201\*201 being the minimum size requires me to put in the `setup` method since Processing.org defaults to 200x200 :(

```
void setup(){noFill();size(201,201);}void draw(){float f=10,a=PI*2*random(),r=random();point(f+f*sin(a)*r,f+f*cos(a)*r);ellipse(f,f,f*2,f*2)}

```
> 4 votes

# Answer
# QBasic, 138 bytes - 20 - 5 = 113

```
INPUT n
r=200
SCREEN 12
RANDOMIZE TIMER
CIRCLE(r,r),r
PAINT(r,r)
FOR i=1TO n
DO
x=RND*r*2
y=RND*r*2
LOOP UNTIL POINT(x,y)
PSET(x,y),1
NEXT

```

Takes user input and draws the disc and points. Tested on QB64.

This is a pretty basic "throw at the dartboard and keep what sticks" strategy. The catch is that "what sticks" is not determined mathematically but graphically: a white disc is plotted on a black background, and then randomly generated points are rejected until they are not black. The points themselves are drawn in blue (though it's hard to tell when they're single pixels--click on the image to enlarge).
> 4 votes

# Answer
# awk - 95 - 5 = 90

```
{
    for(;$1--;printf"("(rand()<.5?x:-x)","(rand()<.5?y:-y)")")
        while(1<(x=rand())^2+(y=rand())^2);
}

```

Since i wasn't quite sure about the rand()\<.5 part, I did some distribution testing with this, using this script:

```
BEGIN{ srand() }
{
    split("0 0 0 0", s)
    split("0 0 0 0", a)

    for(i=$1; i--; )
    {
        while( 1 < r2 = ( x=rand() )^2 + ( y=rand() )^2 );

        x = rand()<.5 ? x : -x
        y = rand()<.5 ? y : -y

        ++s[ x>=0 ? y>=0 ? 1 : 4 : y>=0 ? 2 : 3 ]

        ++a[ r2>.75 ? 1 : r2>.5 ? 2 : r2>.25 ? 3 : 4]
    }

    print "sector distribution:"
        for(i in s) print "sector " i ": " s[i]/$1

    print "quarter area distribution:"
        for(i in a) print "ring " i ":   " a[i]/$1
}

```

which for an input of 1e7 gives me this result, after I sipped once or twice at my coffee:

```
1e7
sector distribution:
sector 1: 0.250167
sector 2: 0.249921
sector 3: 0.249964
sector 4: 0.249948
quarter area distribution:
ring 1:   0.24996
ring 2:   0.25002
ring 3:   0.250071
ring 4:   0.249949

```

which I think is quite alright.

A little explanation:
After scribbling about for a while it turned out, that if you want to divide the disc into four rings with equal area, the radii where you would have to cut are sqrt(1/4), sqrt(1/2) and sqrt(3/4). Since the actual radius of the point I test would be sqrt( x^2 + y^2), I can skip the square rooting all together. The 1/4, 2/4, 3/4 "coincidence" might be related to what M. Buettner pointed out earlier.
> 3 votes

# Answer
# HPPPL, 146 (171-20-5) bytes

```
EXPORT r(n)BEGIN LOCAL R,A,i,Q;RECT();Q:=118.;ARC_P(Q,Q,Q);FOR i FROM 1 TO n DO R:=√RANDOM(1.);A:=RANDOM(2*π);PIXON_P(G0,IP(Q+Q*R*COS(A)),IP(Q+Q*R*SIN(A)));END;FREEZE;END;

```

Example for 10000 points (including timing in seconds for the real device):

The function itself is called by `r(n)`. The rest in the image above is only for timing purposes.

Result (disc diameter is 236 pixels):

The version above does not store the point coordinates, so I wrote a version that takes two parameters `r(n,p)`. `n` is the number of points and `p=0` returns the points to the terminal, `p=1` plots the points and the disc), in case storing coordinates is mandatory. This version is 283 (308-20-5) bytes long:

```
EXPORT r(n,p)BEGIN LOCAL R,A,j,Q,x,y;Q:=118.0;CASE IF p==0 THEN print() END IF p==1 THEN RECT();ARC_P(Q,Q,Q) END END;FOR j FROM 1 TO n DO R:=√RANDOM(1.0);A:=RANDOM(2*π);x:=R*COS(A);y:=R*SIN(A);CASE IF p==0 THEN print("("+x+", "+y+")") END IF p==1 THEN PIXON_P(G0,IP(Q+Q*x),IP(Q+Q*y)) END END;END;FREEZE;END;

```

The ungolfed version:

```
EXPORT r(n,p)
BEGIN
LOCAL R,A,j,Q,x,y;
  Q:=118.0;
  CASE
    IF p==0 THEN print() END
    IF p==1 THEN RECT();ARC_P(Q,Q,Q) END
  END;
  FOR j FROM 1 TO n DO
    R:=√RANDOM(1.0);
    A:=RANDOM(2*π);
    x:=R*COS(A);
    y:=R*SIN(A);
    CASE
      IF p==0 THEN print("("+x+", "+y+")") END
      IF p==1 THEN PIXON_P(G0,IP(Q+Q*x),IP(Q+Q*y)) END
    END;
  END;
  FREEZE;
END;

```

Terminal output for `r(10,0)`:

`r(10,1)` shows the disc with the points, like shown above.
> 3 votes

# Answer
# JavaScript, 75 bytes

Rejection-based:

```
do x=(r=()=>4*Math.random()-2)(),y=r()
while(x*x+y*y>1)
alert(`(${[x,y]})`)

```

Direct method (80 bytes):

```
alert(`(${[(z=(m=Math).sqrt((r=m.random)()))*m.sin(p=m.PI*2*r()),z*m.cos(p)]})`)

```
> 2 votes

# Answer
# Python, 135 130 bytes

```
from random import*
def r():return uniform(-1,1)
p=[]
while not p:
    x,y=r(),r()
    if x**2+y**2<=1:p=x,y
print'(%.2f, %2f)'%p

```

Removed the `**0.5` thanks to @jimmy23013’s suggestion (because it is a unit circle, I am now checking whether the distance squared between (x, y) and (0, 0) is equal to 12. This is the same thing).

This also freed me to remove the parentheses.
> 2 votes

# Answer
# Python 3.8 (pre-release), ~~82~~ ~~79~~ 78 Bytes

```
from random import*
r=random
while(x:=r()*2-1)*x+(y:=r()*2-1)*y>1:1
print(x,y)

```

# -3 from TheThonnu by murdering whitespace and exponents.

# -1 from TheThonnu by showing that you can change function names by variable assignment.

### Verify that it actually works by running on TIO.

## The walrus likes to assign variables and return their value at the same time, and I don't understand why...
> 1 votes

# Answer
# Befunge-with-graphics, ~~224~~ 212 bytes (237 before deductions)

Takes input from stdin, draws image to the screen.

```
>99*3*:s000fc88*3*::f765+jv
v1j0-10::*:+29p02j50p01j50<
>0g29+:*-:*20g92+:*-:*+v
v:-1xg02g01_v#!`*:*:+92<
>!#v_       >94+2*1j
v
0vu<
6>z:1-#v_@
#^$$$$$<
>$$99*3*%r
0
|!:-1\<>0>~:52*-   #v_$r
\  > v*  ^+**25\-"0"<
>#v?v 4
  0123+
  >>>>^

```

## Explanation

Random coordinates are generated using the `?` operator, to build a large random number. Each of the 4 directions the `?` operator could go adds a new digit to a base-4 number which is then eventually modulo'd by the screen width to get a random coordinate on the screen. (bottom-most bit)

The distance from the centre of the circle is then calculated, and if it's outside the circle it generates a new pair of coordinates. Else, it draws a pixel and repeats, decrementing a counter. (top block of stuff)

Finally, the screen is updated and events are handled. (bit with all the `$$$$`s)

Circle with 10 points (to verify there are actually 10):

Circle with 10000 points:

I'm not certain this counts because of many factors

* Befunge-with-graphics is my "made up" language, but it predates my knowledge of this question.
* Since befunge doesn't support floating point, all coordinates are in pixels and may not have the required precision. The image, however, is greater than 201x201 pixels in size.
* Technically no points can appear directly on the edge of the circle.
> 0 votes

# Answer
# Desmos, 52 - 20 = 32 Bytes

```
r=1
a=1
b=1
(acosb,asinb)
a->random(),b->2πrandom()

```

I’ve never golfed in Desmos before. Click the bottom action to generate a random point. Try prettified version online. If anyone knows a shorter way to graph a polar coordinate point, please tell me, since I’m unaware of any builtin to graph a point in polar coordinates.
> 0 votes

---
Tags: code-golf, math, graphical-output, random
---