[Question]
  [
Consider the standard equilateral triangle, with nodes labeled using [barycentric coordinates](https://en.wikipedia.org/wiki/Barycentric_coordinate_system):
![](https://i.stack.imgur.com/URM7b.png)
We can turn this 3 node triangle into a 6 node triangle by adding a new line of 3 vertices (one more than was present on a side of the original 3 node triangle), remove any internal edges (but **not** internal nodes) and re-normalize the coordinates:
[![enter image description here](https://i.stack.imgur.com/cbF99.png)](https://i.stack.imgur.com/cbF99.png)
Repeating the process to go from a 6 node triangle to a 10 node triangle, add a line of 4 vertices (again, one more than was present on a side of the original 6 node triangle), remove any internal edges (but **not** internal nodes) and re-normalize the coordinates:
[![enter image description here](https://i.stack.imgur.com/Ff1Zt.png)](https://i.stack.imgur.com/Ff1Zt.png)
This process can be repeated indefinitely. The goal of this challenge is given an integer `N` representing how many times this process has been performed, output all the nodes for the associated triangle in barycentric coordinates.
# Input
Your program/function should take as input a single non-negative integer `N` representing how many times this process has been applied. Note that for `N=0`, you should output the original triangle with 3 nodes.
The input may come from any source (function parameter, stdio, etc.).
# Output
Your program/function should output all the nodes in normalized barycentric coordinates. The order of the nodes does not matter.
A number can be specified as a fraction (fraction reduction not required) or a floating point number. You may also output "scaled" vectors to specify a node. For example, all 3 of the following outputs are equivalent and allowed:
```
0.5,0.5,0
1/2,2/4,0
[1,1,0]/2
```
If using floating point output, your output should be accurate to within 1%. The output may be to any sink desired (stdio, return value, return parameter, etc.). Note that even though the barycentric coordinates are uniquely determined by only 2 numbers per node, you should output all 3 numbers per node.
# Examples
Example cases are formatted as:
```
N
x0,y0,z0
x1,y1,z1
x2,y2,z2
...
```
where the first line is the input `N`, and all following lines form a node `x,y,z` which should be in the output exactly once. All numbers are given as approximate floating point numbers.
```
0
1,0,0
0,1,0
0,0,1
1
1,0,0
0,1,0
0,0,1
0.5,0,0.5
0.5,0.5,0
0,0.5,0.5
2
1,0,0
0,1,0
0,0,1
0.667,0,0.333
0.667,0.333,0
0.333,0,0.667
0.333,0.333,0.333
0.333,0.667,0
0,0.333,0.667
0,0.667,0.333
3
1,0,0
0.75,0,0.25
0.75,0.25,0
0.5,0,0.5
0.5,0.25,0.25
0.5,0.5,0
0.25,0,0.75
0.25,0.25,0.5
0.25,0.5,0.25
0.25,0.75,0
0,0,1
0,0.25,0.75
0,0.5,0.5
0,0.75,0.25
0,1,0
```
# Scoring
This is code golf; shortest code in bytes wins. Standard loopholes apply. You may use any built-ins desired.
  
[Answer]
## CJam (22 bytes)
```
{):X),3m*{:+X=},Xdff/}
```
This is an anonymous block (function) which takes `N` on the stack and leaves an array of arrays of doubles on the stack. [Online demo](http://cjam.aditsu.net/#code=2%0A%0A%7B)%3AX)%2C3m*%7B%3A%2BX%3D%7D%2CXdff%2F%7D%0A%0A~%60)
### Dissection
```
{         e# Define a block
  ):X     e# Let X=N+1 be the number of segments per edge
  ),3m*   e# Generate all triplets of integers in [0, X] (inclusive)
  {:+X=}, e# Filter to those triplets which sum to X
  Xdff/   e# Normalise
}
```
[Answer]
# Haskell, 53 bytes
```
f n|m<-n+1=[map(/m)[x,y,m-x-y]|x<-[0..m],y<-[0..m-x]]
```
[Answer]
# Python 3, 87 bytes
This is actually supposed to be a comment to [the solution by TheBikingViking](https://codegolf.stackexchange.com/a/90801/58861) but I don't have enough reputation for comments.
One can save a few bytes by only iterating over the variables `i,j` and using the fact that with the third one they add up to `n+1`.
```
def f(n):d=n+1;r=range(n+2);print([[i/d,j/d,(d-i-j)/d]for i in r for j in r if d>=i+j])
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E/), 10 bytes
```
ÌL<¤/3ãDOÏ
```
**Explanation**
```
ÌL<          # range(1,n+2)-1
   ¤/        # divide all by last element (n+1)
     3ã      # cartesian product repeat (generate all possible triples)
       DO    # make a copy and sum the triples
         Ï   # keep triples with sum 1
```
[Try it online](http://05ab1e.tryitonline.net/#code=w4xMPMKkLzPDo0RPw48&input=Mw)
[Answer]
## Mathematica, ~~44~~ 43 bytes
```
Select[Range[0,x=#+1]~Tuples~3/x,Tr@#==1&]&
```
This is an unnamed function taking a single integer argument. Output is a list of lists of exact (reduced) fractions.
Generates all 3-tuples of multiples of `1/(N+1)` between 0 and 1, inclusive, and then selects those whose sum is 1 (as required by barycentric coordinates).
[Answer]
# [MATL](http://github.com/lmendo/MATL), 17 bytes
```
2+:qGQ/3Z^t!s1=Y)
```
[Try it online!](http://matl.tryitonline.net/#code=Mis6cUdRLzNaXnQhczE9WSk&input=Mw)
### Explanation
The approach is the same as in other answers:
1. Generate the array `[0, 1/(n+1), 2/(n+1), ..., 1]`, where `n` is the input;
2. Generate all 3-tuples with those values;
3. Keep only those whose sum is `1`.
More specifically:
```
2+     % Take input and add 2: produces n+2
:q     % Range [0 1 ... n+1]
GQ/    % Divide by n+1 element-wise: gives [0, 1/(n+1), 2/(n+1)..., 1]
3Z^    % Cartesian power with exponent 3. Gives (n+1)^3 × 3 array. Each row is a 3-tuple
t      % Duplicate
!s     % Sum of each row
1=     % Logical index of entries that equal 1
Y)     % Use that index to select rows of the 2D array of 3-tuples
```
[Answer]
## [Jellyfish](https://github.com/iatorm/jellyfish/blob/master/doc.md), ~~37~~ 33 bytes
*Thanks to Zgarb for saving 4 bytes.*
```
p
*%
# S
`
=E   S
`/
1+r#>>i
   3
```
[Try it online!](http://jellyfish.tryitonline.net/#code=cAoqJQojIFMKYAo9RSAgIFMKYC8KMStyIz4-aQogICAz&input=Mg)
Like my Mathematica and Peter's CJam answers, this generates a set of candidate tuples and then selects only those that sum to 1. I'm not entirely happy with the layout yet, and I wonder whether I can save some bytes with hooks or forks, but I'll have to look into that later.
[Answer]
# Perl 6: ~~50~~ 40 bytes
```
{grep *.sum==1,[X] (0,1/($_+1)...1)xx 3}
```
Returns a sequence of 3-element lists of (exact) rational numbers.
Explanation:
* `$_`  
Implicitly declared parameter of the lambda.
* `0, 1/($_ + 1) ... 1`  
Uses the sequence operator `...` to construct the arithmetic sequence
that corresponds to the possible coordinate values.
* `[X] EXPR xx 3`  
Takes the Cartesian product of three copies of EXPR, i.e. generates all possible 3-tuples.
* `grep *.sum == 1, EXPR`  
Filter tuples with a sum of 1.
[Answer]
# Ruby, 62
I'd be surprised if this can't be improved on:
```
->x{0.step(1,i=1.0/(x+1)){|a|0.step(1-a,i){|b|p [a,b,1-a-b]}}}
```
Taking the advice latent in the puzzle, this calculates the second node options based on the first, and the third node by subtracting the first two.
[Answer]
# [Brachylog](http://github.com/JCumin/Brachylog), 24 bytes
```
+:1f
yg:2j:eaL+?/g:Lz:*a
```
[Try it online!](http://brachylog.tryitonline.net/#code=KzoxZgp5ZzoyajplYUwrPy9nOkx6Oiph&input=MQ&args=Wg&debug=on)
[Answer]
# Python 3, 106 bytes
```
def f(n):r=range(n+2);print([x for x in[[i/-~n,j/-~n,k/-~n]for i in r for j in r for k in r]if sum(x)==1])
```
A function they takes input via argument and prints a list of lists of floats to STDOUT.
Python is not good at Cartesian products...
**How it works**
```
def f(n):                         Function with input iteration number n
r=range(n+2)                      Define r as the range [0, n+1]
for i in r for j in r for k in r  Length 3 Cartesian product of r
[i/-~n,j/-~n,k/-~n]               Divide each element of each list in the product
                                  by n+1
[x for x in ... if sum(x)==1]     Filter by summation to 1
print(...)                           Print to STDOUT
```
[Try it on Ideone](http://ideone.com/yuki0R)
[Answer]
# [Actually](https://github.com/Mego/Seriously), 15 bytes
This uses an algorithm similar to the one in [TheBikingViking's Python answer](https://codegolf.stackexchange.com/a/90801/47581). Golfing suggestions welcome. [Try it online!](http://actually.tryitonline.net/#code=dTt1cuKZgC8zQOKImWDOozE9YOKWkQ&input=Mw)
```
u;ur♀/3@∙`Σ1=`░
```
**Ungolfed:**
```
u;                Increment implicit input and duplicate.
  ur              Range [0..n+1]
    ♀/            Divide everything in range by (input+1).
      3@∙         Take the Cartesian product of 3 copies of [range divided by input+1]
         `Σ1=`    Create function that takes a list checks if sum(list) == 1.
              ░   Push values of the Cartesian product where f returns a truthy value.
```
[Answer]
# Ruby, ~~77~~ 74 bytes
Another answer using the algorithm in [TheBikingViking's Python answer](https://codegolf.stackexchange.com/a/90801/47581). Golfing suggestions welcome.
```
->n{a=[*0.step(1,i=1.0/(n+1))];a.product(a,a).reject{|j|j.reduce(&:+)!=1}}
```
Another 74-byte algorithm based on [Not that Charles's Ruby answer](https://codegolf.stackexchange.com/a/90828/47581).
```
->x{y=x+1;z=1.0/y;[*0..y].map{|a|[*0..y-a].map{|b|p [a*z,b*z,(y-a-b)*z]}}}
```
[Answer]
## JavaScript (Firefox 30-57), ~~88~~ 81 bytes
```
n=>[for(x of a=[...Array(++n+1).keys()])for(y of a)if(x+y<=n)[x/n,y/n,(n-x-y)/n]]
```
Returns an array of arrays of floating-point numbers. Edit: Saved 7 bytes by computing the third coordinate directly. I tried eliminating the `if` by calculating the range of `y` directly but it cost an extra byte:
```
n=>[for(x of a=[...Array(++n+1).keys()])for(y of a.slice(x))[x/n,(y-x)/n,(n-y)/n]]
```
]