# Question
Title: Fastest count of certain hypercubes with labeled vertices
**CHALLENGE**

This is a fastest-code challenge.

Count how many n-dimensional hypercubes with n=1,2,3,4 exist, with vertices labeled with either 1 or 0, such that there does not exist any rectangle formed by 4 vertices such that one diagonal is 1-1 and the other diagonal is 0-0. To be more clear, the forbidden rectangles have the form:

```
10 or 01
01    10

```

Note that we are talking about all rectangles, not just the faces of the hypercube. For example in the cube:

```
   0*---1*
  /|   /|
 0----0 |
 | 0--|-0
 |/   |/
 1*---0*

```

there is an invalid rectangle formed by the four vertices marked with a star (\*), although all faces are valid.

The list of all known expected results is at OEIS A000609. This is actually a conjecture, but a probable one. For more information and some "hints", see the BACKGROUND section.

The count for n=1,2,3,4 must be repeated in a loop 40 times (just to have more resolution on the measured time), and the results for n=1,2,3,4 (4,14,104,1882 respectively) printed on the first round of the loop only.

You should provide a link for tio.run.

Here is my best in C (about 30 s).

**BACKGROUND (optional)**

I have asked a question related to this challenge at mathoverflow. See the links within that question and especially the one related to the generation of the rectangles.

Another possible idea is using the inclusion-exclusion principle: see this question.

**SCORE**

I will make a first rough placement list based on the time shown in the tio.run link provided, in the "Debug" section, "User time". Some time later a more precise one with times tested on my machine, since performance on tio.run may vary over time.

**RECORDS**

There will be another record list if someone is able to compute the case n = 5 or more in a reasonable time (let's say less than 10 minutes for n = 5 and 3 days for n = 6). In this case it may be that I won't double check, so I will trust your claim. I personally computed the case n = 5 in about 5 minutes with a slightly modified version of the above linked C software (which took 8 minutes).

**STANDINGS**

| FASTEST CODE | USER | LANGUAGE | BEST TIME ON MY MACHINE |
| --- | --- | --- | --- |
| 1 | 138 Aspen | Rust | 0.1364747 s |
| 2 | gsitcia | Python | 0.1462836 s |
| 3 | user1502040 | Python | 1.0249261 s |
| 4 | Level River St | Ruby | 1.3259241 s |
| 5 | Me | C | 25.888948 s |

| "n RECORD" | USER | LANGUAGE | RECORD and TIME MY PC |
| --- | --- | --- | --- |
| 1 | gsitcia | Python | n = 7 in 1h 19m 26s |
| 2 | gsitcia | Python | n = 6 in 1.946 s |
| 3 | 138 Aspen | Rust | n = 6 in 2h 45m 24s |
| 4 | user1502040 | Python | n = 6 in 22h 6m |
| 5 | gsitcia | Python | n = 5 in 0.156 s |
| 6 | 138 Aspen | Rust | n = 5 in 1.947 s |
| 7 | user1502040 | Python | n = 5 in 25.372 s |
| 8 | Level River St | Ruby | n = 5 in 27.900 s |
| 9 | Me | C | n = 5 in ~480 s |

The time for "n RECORD" is for a single computation of the specified value for n.

# Answer
> 3 votes
# Python + python-sat (about 1 second on my machine)

```
import itertools
from pysat.solvers import Glucose4

def count_cubes(n):
    clauses = []
    for (x, y, z) in itertools.combinations(range(1 << n), 3):
        e0_p = y & (~x)
        e0_n = x & (~y)
        e0 = e0_p | e0_n
        e1_p = z & (~x)
        e1_n = x & (~z)
        e1 = e1_p | e1_n
        if e0 & e1:
            continue
        w = x ^ (e0 | e1)
        x, y, z, w = x + 1, y + 1, z + 1, w + 1
        clauses.append([x, w, -y, -z])
        clauses.append([-x, -w, y, z])
    if not clauses:
        assert n == 1
        clauses.append([-1, 1, -2, 2])
    with Glucose4(bootstrap_with=clauses, incr=True) as m:
        count = 0
        while m.solve():
            count += 1
            model = m.get_model()
            m.add_clause([-v for v in model])
    return count

for t in range(40):
    for n in range(1, 5):
        r = count_cubes(n)
        if t == 0:
            print(n, r)

```

# Answer
> 3 votes
# Ruby, n=5 in 25 seconds on TIO

This code has been updated to enable `n=6` to be reached without memory issues.

Specifically `[0,(1<<2**n/2)-1].product(a,a).each` has been modified to two nested loops. Facility for counting the iterations of the outer loop `#p y` has been added but commented out. The number of iterations is equal to the number of solutions for the previous `n` as explained below.

Timing facilities have also been added. On my machine this program uses about 33% of my 6 core processor and for `n=6` it handles about 1000 iterations of the outer loop per hour, which means it would take 4 days to complete 94572 iterations for `n=6`. I'm glad to see others have used the ideas in this answer. There are some minor optimisations that could be made, but a different approach such as taking advantage of symmetry (see gsitcia's answer) is needed for major progress.

```
N=5
starttime=Time.now
f=->n=1,a=[0]{return if n>N
#generate rectangles
u=[]
(6**(n-1)).times{|i|
  j=i*6
  p=q=r=s=0
  t=0
  n.times{
    k=j%6
    p=p<<1|(k&1); q=q<<1|(41>>k&1); r=r<<1|(25>>k&1); s=s<<1|(26>>k&1)
    t=k if k>3
    j/=6
  }
  u<<[p,q,r,s] if t>4
}
#check colourings
b=[]
[0,(1<<2**n/2)-1].product(a).each_with_index{|z,y|c,d=z;#p y
  a.each{|e|
    e=(e^c)|d<<2**n/2
    h=false
    u.each{|w|p,q,r,s=w
      g=e>>p&1
      h = g!=e>>q&1 && g==e>>r&1 && g!=e>>s&1
      break if h
    }
  b<<e unless h
  }
}
p [n,b.size*2]
f[n+1,b]
}
f[]
p [starttime, Time.now, Time.now-starttime]

```

Try it online!

# Ruby, n=4 repeated 40 times in 1.3 seconds on TIO

```
N=4
printoutput=true
f=->n=1,a=[0]{return if n>N
#generate rectangles
u=[]
(6**(n-1)).times{|i|
  j=i*6
  p=q=r=s=0
  t=0
  n.times{
    k=j%6
    p=p<<1|(k&1); q=q<<1|(41>>k&1); r=r<<1|(25>>k&1); s=s<<1|(26>>k&1)
    t=k if k>3
    j/=6
  }
  u<<[p,q,r,s] if t>4
}
#check colourings
b=[]
[0,(1<<2**n/2)-1].product(a,a).each{|z|c,d,e=z
  e=(e^c)|d<<2**n/2
  h=false
  u.each{|w|p,q,r,s=w
    g=e>>p&1
    h = g!=e>>q&1 && g==e>>r&1 && g!=e>>s&1
    break if h
  }
  b<<e unless h
}
p b.size*2 if printoutput
f[n+1,b]
}
40.times{f[];printoutput=false}

```

Try it online!

# Explanation

This is a recursive function which takes the solution for the previous `n`, makes a cartesian product of two copies, calculates all rectangles that bridge between them, then checks the colourings of these bridging rectangles. Because the solution for the previous `n` is known to be free of illegal colourings, it is only necessary to check rectangles that span the new dimension. This approach means there are only `1882 x 1882 = 3541924` colourings to check for `n=5` instead of `2**32 = 4 billion.` Also we only calculate colourings where the most significant bit is `0` which halves the number of calculations; the remaining colourings are the same but with all the 0's and 1's flipped. Precalculating the rectangles before checking colourings also helps.

**Generating rectangles**

Paraphrasing this post referenced in the question, each rectangle can be represented as 4 points `[p, q=p+x, r=p+x+y, s=p+y]` where `x` and `y` are n-dimensional vectors with values of `-1`, `0` or `+1` in each dimension. For the shape to be a rectangle `x` and `y` must be orthogonal, which means that for every dimension `i`, the product `x[i]*y[i]` must be zero. `p[i]` can be `0` or `+1` but if either `x[i]` or `y[i]` is nonzero `p[i]` is constrained to only one possibility. This gives the following possibilities:

```
 k    =     0   1   2   3   4   5
 p[i] =     0   1   0   1   0   1
 x[i] =    +1  -1   0   0   0   0
 y[i] =     0   0   0   0  +1  -1

```

The value of `j` is treated as a base 6 number with `n` digits, the last of which is zero. Because the least significant digit of `j` generates the most significant bit of `p`,`q`,`r` and `s` (representing the new dimension `n`) this ensures that we only generate rectangles where `x[n]` is `+1` ensuring that the generated rectangles span the new dimension.

To avoid double counting, we keep track of the last value of `k` in `t` where `k>3` (hence `y[i]<>0`) and only accept rectangles where `t==5` meaning the last nonzero value of `y[i]` was negative (not positive.) We also reject cases where `y[i]` is zero for all dimensions, since if `y` is zero, the shape is a line, not a rectangle.

**Commented code**

This is the same as the code for n=5 in 25 second above. The comments do not have significant impact on runtime.

```
# Set a constant N for recursion depth, and record  the start time
N=5
starttime=Time.now

# Define a recursive function f where n is the number of dimensons
# a is the set of solutions for the previous recursion. Each solution is represented by an integer
# Initially this is [0,1] for n=0 but we consider only colourings where the last bit is zero.

f=->n=1,a=[0]{return if n>N

# generate rectangles
u=[]                         # Make an array to store rectangles as arrays of 4 integers corresponding to indices of vertices
(6**(n-1)).times{|i|         # Iterate through all base 6 numbers with n-1 digits
  j=i*6                      # Multiply 6 by so that the last digit of j is zero. This ensures we only generate rectangles where x=+1 in the new dimension
  p=q=r=s=0                  # Set up variabes for the 4 vertices
  t=0                        # Set up variable to record a code for value of y
  n.times{
    k=j%6                    # Extract the last digit of j and decode the values for the last bit of p,q,r,s shifting any existing bits left
    p=p<<1|(k&1)
    q=q<<1|(41>>k&1)
    r=r<<1|(25>>k&1)
    s=s<<1|(26>>k&1)
    t=k if k>3               # If y<>0, record the value of k in t: 4 for positive, 5 for negative
    j/=6                     # Shift j one base 6 digit right
  }
  u<<[p,q,r,s] if t>4        # If y is nonzero in at least one dimension, and the last nonzero dimension was negative, add the rectangle to the list
}

#check colourings
b=[]                         # Make an array to store solutions for dimension n
[0,(1<<2**n/2)-1].product(a).# Consider all possible solutions for n-1 from a with 0 in the most significant bit. Consider also an XOR mask for flipping bits
  each_with_index{|z,y|      # Iterate through all possibilities. y gives the index for use in monitoring progress
  c,d=z                      # c is the XOR mask for bit flipping. d is the n-1 cube for the high order bits
  #p y                       # (Uncomment to display iteration progress)
  a.each{|e|                 # Iterate through all possible solutions for n-1 with 0 in the most significant bit. e is the n-1 cube for the low order bits.
    e=(e^c)|d<<2**n/2        # Proposed solution is e | d<<2**n/2 . c is either 0 or (1<<2**n/2)-1 which is used to flip the low order bits.
    h=false                  # h=false to indicate no illegal rectangles found
    u.each{|w|p,q,r,s=w      # Iterate through rectangles where w is an array of integers [p,q,r,s] corresponding to indices of vertices
      g=e>>p&1               # Find the value g representing the colour (0 or 1) of vertex p
      h = g!=e>>q&1 &&       # If g is the same as the diagonally opposite vertex r and different from the adjacent vertices p and s...
       g==e>>r&1 && g!=e>>s&1
      break if h             # ...h becomes true showing the rectangle is illegal. Break.
    }
  b<<e unless h              # If no illegal rectangles were found add e to the list of solutions
  }
}
p [n,b.size*2]               # Print the value of n and the number of solutions found
f[n+1,b]                     # (double the size of b since b only contains solutions where the last bit is zero)
}
f[]                          # Call the function, then print the start time, end time and elapsed time.
p [starttime, Time.now, Time.now-starttime]

```

Try it online!

# Answer
> 3 votes
# Rust, n=4 repeated 40 times in 0.6 seconds on TIO

A port of @Level River St's Ruby answer in Rust.

---

Try it online!

```
const N: usize = 4;

fn main() {
    let mut print_output = true;
    let mut a: Vec<u64> = vec![0];
    rec_function(1, &mut a, &mut print_output);

    for _ in 0..40 {
        print_output = false;
        rec_function(1, &mut a, &mut print_output);
    }
}

fn rec_function(n: usize, a: &mut Vec<u64>, print_output: &mut bool) {
    if n > N {
        return;
    }

    let mut u: Vec<(u64, u64, u64, u64)> = vec![];

    for i in 0..6_u64.pow((n - 1) as u32) {
        let mut j = i * 6;
        let (mut p, mut q, mut r, mut s) = (0, 0, 0, 0);
        let mut t = 0;
        for _ in 0..n {
            let k = j % 6;
            p = p << 1 | (k & 1);
            q = q << 1 | (41 >> k & 1);
            r = r << 1 | (25 >> k & 1);
            s = s << 1 | (26 >> k & 1);
            if k > 3 {
                t = k;
            }
            j /= 6;
        }
        if t > 4 {
            u.push((p, q, r, s));
        }
    }

    let mut b: Vec<u64> = vec![];
    let pairs = [0, (1 << (1 << (n - 1))) - 1];

    for &c in pairs.iter() {
        for &d in a.iter() {
            for &e in a.iter() {
                let e = (e ^ c) | d << (1 << (n - 1));
                let mut h = false;
                for &(p, q, r, s) in &u {
                    let g = e >> p & 1;
                    h = g != (e >> q & 1) && g == (e >> r & 1) && g != (e >> s & 1);
                    if h {
                        break;
                    }
                }
                if !h {
                    b.push(e);
                }
            }
        }
    }

    if *print_output {
        println!("{}", b.len() * 2);
    }
    rec_function(n + 1, &mut b, print_output);
}

```

# Answer
> 3 votes
# Python 3 (+Numpy) (n=4 repeated 40 times in 0.15 seconds, n=6 in 2.5 seconds on TIO)

```
import itertools
import functools
from collections import Counter
import time
import numpy as np

@functools.lru_cache(None)
def gen_rectangles(n):
    x = []
    for t in itertools.product([
        (0,0,0),
        (1,1,1),
        (0,0,1),
        (1,1,0),
        (0,1,0),
        (1,0,1),
    ], repeat=n):
        a = b = c = 0
        for a1, b1, c1 in t:
            a, b, c = a*2+a1, b*2+b1, c*2+c1
        if a < b < c:
            x.append((a, b, c, a^b^c))
    return x

def mask(n, k):
    return ((1<<(1<<n))-1)//((1<<(1<<k))+1)

def swap_parts(n, k1, k2):
    m1 = mask(n, k1)
    m2 = mask(n, k2)
    m12 = m1 ^ m2
    ma =  m1 & ~m2
    mb = ~m1 &  m2
    d = 2**k2 - 2**k1
    return m12, ma, mb, d

def swap(x, m12, ma, mb, d):
    xs = x & m12
    a = x & ma
    b = x & mb
    return x ^ xs ^ (a>>d) ^ (b<<d)

def flip_parts(n, k):
    m = mask(n, k)
    k1 = 1 << k
    return m, k1

def flip(x, m, k1):
    a = x & m
    b = (x >> k1) & m
    return (a << k1) | b

@functools.lru_cache(None)
def parts(n):
    return [flip_parts(n, k) for k in range(n)], [swap_parts(n, k1, k2) for k2 in range(n) for k1 in range(k2)]

def reduce(x, n):
    fp, sp = parts(n)
    for p in fp:
        x1 = flip(x, *p)
        if x1 < x:
            return x1, 0, p
    for p in sp:
        x1 = swap(x, *p)
        if x1 < x:
            return x1, 1, p

def make_order(d):
    d = {a:b[0] for a,b in d.items()}
    C = Counter(d.values())
    empty = d.keys() - d.values()
    l = []
    while empty:
        i = empty.pop()
        j = d.get(i, None)
        if j is None: continue
        l.append(i)
        k = C[j] - 1
        if k == 0:
            empty.add(j)
        else:
            C[j] = k
    return l

@functools.lru_cache(None)
def gen_cubes(n):
    if n == 0:
        return [0, 1]
    if n == 1:
        return [0, 1, 2, 3]
    prev = gen_cubes(n-1)
    half = 2 ** (n-2)
    k = 2 ** (n-1)
    full = 2 ** k - 1
    n1 = n - 1
    bad = {}
    good = []
    for x in prev:
        bc = bin(x).count('1')
        if bc > half:
            bad[x] = full ^ x, 2, full
            continue
        t = reduce(x, n1)
        if t:
            bad[x] = t
            continue
        good.append(x)
    order = make_order(bad)
    compatible = {}
    masks = []
    for a,b,c,d in gen_rectangles(n):
        if b >= k or c < k: continue
        a, b, c, d = 1 << a, 1 << b, 1 << (c-k), 1 << (d-k)
        masks.append((a|b, c|d, a, d, b, c))
    num_good = sum(1 for b in good for t in masks if b & t[1] in (t[3], t[5]))
    prev = np.array(prev, dtype=np.uint64)
    for b in good:
        l = prev
        for ma, mb, pa1, pb1, pa2, pb2 in masks:
            q = b & mb
            if q == pb1:
                pa = pa1
            elif q == pb2:
                pa = pa2
            else:
                continue
            l = l[l&ma != pa]#[a for a in l if a&ma != pa]
        compatible[b] = list(map(int, l))
    for x in reversed(order):
        pre, disc, p = bad[x]
        l = compatible[pre]
        if disc == 0:
            l = [flip(i, *p) for i in l]
        elif disc == 1:
            l = [swap(i, *p) for i in l]
        else:
            l = [p ^ i for i in l]
        compatible[x] = l
    out = []
    for x, l in compatible.items():
        x <<= k
        out.extend(x|i for i in l)
    return out

def count_cubes(n):
    if n == 0:
        return 2
    if n == 1:
        return 4
    prev = gen_cubes(n-1)
    half = 2 ** (n-2)
    k = 2 ** (n-1)
    full = 2 ** k - 1
    n1 = n - 1
    bad = {}
    good = []
    for x in prev:
        bc = bin(x).count('1')
        if bc > half:
            bad[x] = full ^ x, 2, full
            continue
        t = reduce(x, n1)
        if t:
            bad[x] = t
            continue
        good.append(x)
    order = make_order(bad)
    compatible = {}
    masks = []
    for a,b,c,d in gen_rectangles(n):
        if b >= k or c < k: continue
        a, b, c, d = 1 << a, 1 << b, 1 << (c-k), 1 << (d-k)
        masks.append((a|b, c|d, a, d, b, c))
    num_good = sum(1 for b in good for t in masks if b & t[1] in (t[3], t[5]))
    prev = np.array(prev, dtype=np.uint64)
    for b in good:
        l = prev
        for ma, mb, pa1, pb1, pa2, pb2 in masks:
            q = b & mb
            if q == pb1:
                pa = pa1
            elif q == pb2:
                pa = pa2
            else:
                continue
            l = l[l&ma != pa]#[a for a in l if a&ma != pa]
        compatible[b] = len(l)
    for x in reversed(order):
        pre, disc, p = bad[x]
        compatible[x] = compatible[pre]
    return sum(compatible.values())

if __name__ == '__main__':
    t = time.time()
    for i in range(7):
        print(i, count_cubes(i))
    print(f'{time.time() - t:.1f} seconds 0-6')
    t = time.time()
    for _ in range(40):
        gen_cubes.cache_clear()
        parts.cache_clear()
        gen_rectangles.cache_clear()
        for i in range(5):
            gen_cubes(i)
    print(f'{time.time() - t:.3f} seconds 0-4 x40')

```

Try it online!

Based on @Level River St's Ruby answer. The main difference is that this answer exploits additional symmetry.

As a simple example, consider a single square of the hypercube. There are 3 cases:

* All four vertices are colored the same (2 cases like this (flip))
* A single vertex is colored differently (8 cases like this (flip/rotate))
* 2 adjacent vertices are colored 1 (4 cases like this (rotate only))

This answer considers the symmetries of the n-1 hypercube, drastically reducing the amount of work necessary (in the case of n=6, it reduces 94572 cases to only 63 cases)

---
Tags: math, fastest-code
---


title: Derivative of a product

text:
# Question
Title: Derivative of a product
In calculus, the derivative of a mathematical function defines the rate at which it changes. The derivative of a function `f(x)` can be marked as `f'(x)`, and these can sometimes be abbreviated to `f` and `f'`.

The derivative of the product of two functions `fg` is `f'g + g'f`. The derivative of three `fgh` is `f'gh + fg'h + fgh'`. In general, the derivative of a product of any number of functions is the sum of the product of all but one, multiplied by the derivative of the remaining one, for each individual function.

Your challenge is to take a string of distinct alphabetical characters and transform it into its derivative. No spacing or simplification is required, and terms separated by `+` may be in any order. The string will contain at least two characters.

This is , so the shortest code wins!

## Testcases

```
abc -> a'bc+ab'c+abc'
de -> d'e+de'
longabcdef -> l'ongabcdef+lo'ngabcdef+lon'gabcdef+long'abcdef+longa'bcdef+longab'cdef+longabc'def+longabcd'ef+longabcde'f+longabcdef'
short -> s'hort+sh'ort+sho'rt+shor't+short'

```

# Answer
# Pip, 15 bytes

```
aRL#aJ'+<>#aJ''

```

Attempt This Online!

### Explanation

```
aRL#aJ'+<>#aJ''
a                Command-line arg       "abc"
   #a            Length of arg          3
 RL              Repeat-list            ["abc"; "abc"; "abc"]
     J'+         Join on "+"            "abc+abc+abc"
        <>#a     Groups of size len(a)  ["abc"; "+ab"; "c+a"; "bc"]
            J''  Join on "'"            "abc'+ab'c+a'bc"

```
> 25 votes

# Answer
# Python, 54 bytes

```
lambda s:((l:=len(s))*(l*"'"+-~l*(s+"+")))[-~l*~l::-l]
```

Attempt This Online!

# Python, 56 bytes

```
lambda s:((l:=-~len(s))*(l*(s+"+")+-~l*"'"))[l:-3*l:-~l]
```

Attempt This Online!

## Old Python, 58 bytes

```
lambda s:((l:=-~len(s))*(l*(s+"+")+l*"'"+"+"))[l:-3*l:-~l]
```

Attempt This Online!

## How?

```
 a b c d +
(a)b c d +
 a(b)c d +
 a b(c)d +
 a b c(d)+
 ' ' ' '(')'
 a b c d(+)
 a b c d +
(a)b c d +
 a(b)c d +
 a b(c)d +
 ' ' '(')' '
 a b c(d)+
 a b c d(+)
 a b c d +
(a)b c d +
 a(b)c d +
 ' '(')' ' '
 a b(c)d +

   etc.

```
> 12 votes

# Answer
# C (clang), ~~69~~ ~~64~~ 59 bytes

-10 thanks to @Noodle9 and @ceilingcat

```
l;f(char*a){for(l=0;a[l];)printf("+%.*s'%s"+!l,++l,a,a+l);}

```

Try it online!
> 10 votes

# Answer
# Brainfuck, ~~87~~ 84 bytes

```
>,>---[>+<+++++++]>++[>,]<[<]<[<]>[[.>]>.>[.[-<<+>>]>[.>]<[<]<++++.----[->+<]<[<]]>]

```

*Saved 3 bytes thanks to Nitrodon.*

Try it online!

**Explanation**:

```
memory progress for abc:
null a null 'bc
null ab null 'c
null abc' null

string is split to two parts by null character

>,>                load first character
---[>+<+++++++]>++ load ' to memory; taken form esolangs wiki
[>,]               load rest of input
<[<]<[<]>          go to beginning of memory
[
 [.>]              print first part
 >.>               print '
 [                 check if there is text after '
  .                print first char of second part
  [-<<+>>]         move char to the first part
  >[.>]            print rest of second part
  <[<]<            go to position of '
  ++++.----        use ' to print plus
  [->+<]           move '
  <[<]             go before beginning
 ]
 >                 move back to beginning
                   if there wasn't text after ' move after end instead
]

```
> 8 votes

# Answer
# Python, 47 bytes

```
lambda s:"+".join(s.replace(c,c+"'")for c in s)

```

Attempt This Online!

How it works:

```
a  b  c  d --> a' b  c  d
a  b  c  d --> a  b' c  d
a  b  c  d --> a  b  c' d
a  b  c  d --> a  b  c  d'

then join with '+'

```
> 8 votes

# Answer
# Factor, 53 bytes

```
[ dup length [1,b] [ cut "'"glue ] with map "+"join ]

```

Try it online!

If the order of each product doesn't matter:

# Factor, 43 bytes

```
[ all-rotations [ 39 suffix ] map "+"join ]

```

Try it online!

## Explanation

```
                  ! "short"
all-rotations     ! { "short" "horts" "ortsh" "rtsho" "tshor" }
[ 39 suffix ] map ! { "short'" "horts'" "ortsh'" "rtsho'" "tshor'" }
"+"join           ! "short'+horts'+ortsh'+rtsho'+tshor'"

```
> 6 votes

# Answer
# JavaScript (Node.js), 46 44 bytes

```
l=>[...l].map(s=>l.replace(s,s+"'")).join`+`

```

Try it online!

-2 bytes thanks to emanresu A
> 6 votes

# Answer
# Jelly, 14 bytes

```
J‘œṖ€j€”'j”+”'

```

Try It Online!

```
J‘œṖ€j€”'j”+”'  Full Program
J               [1, 2, ..., len]
 ‘              [2, 3, ..., len + 1]
  œṖ€           For each of these, partition (*)
     j€”'       Join each partition on ' (**)
         j”+    Join on +
            ”'  An unparseable nilad forces the current string to be printed

(And then the ' is printed at the end)

(*)  partition takes a list on the left but since we give it a single number,
     it auto-wraps it so we just get all needed splits
(**) the last partition only has one sublist so join doesn't work properly here

```
> 5 votes

# Answer
# Wolfram Language (Mathematica), ~~45~~ 36 bytes

```
aRiffle[a/.#->#<>"'"&/@a,"+"]<>""

```

Try it online!

Input `[{characters...}]`. Assumes all characters in input are unique. Returns a string.

---

### Wolfram Language (Mathematica), 35 bytes

```
ToString@Expand[Dt[1##]/._@a_:>a']&

```

Try it online!

Input `[symbols...]`. Assumes symbols are unique and not defined. Returns a spaced string.

`Dt` treats strings as constants.
> 5 votes

# Answer
# Retina, ~~15~~ 13 bytes

```
L$|+`.
$>`'$'

```

Try it online!

-2 thanks to Neil.

Match every character. For each one, take the text before the separator after the matched string, followed by an apostrophe, followed by the text after the matched string. Print as a list with separator `+`.
> 5 votes

# Answer
# K (ngn/k), 21 bytes

Most straightforward K solution.

```
{"+"/?[x;;"'"]'1+!#x}

```

Try it online!

### Explanation

```
{"+"/?[x;;"'"]'1+!#x}
     ?[x;;"'"]          / insert "'" in string at
              '         / each position
               1+!#x    / 1..len of string
 "+"/                   / join the resulting strings with "+"

```
> 5 votes

# Answer
# JavaScript (Node.js), ~~83~~ \\$\cdots\\$ ~~72~~ 69 bytes

```
s=>g=(r="",i=0)=>s[i]?g(r+("+"+s).slice(!i,++i+1)+"'"+s.slice(i),i):r

```

Try it online!

*Saved ~~7~~ 8 bytes thanks to emanresu A!!!*
*Saved 3 bytes thanks to Samathingamajig!!!*
> 5 votes

# Answer
# Charcoal, 10 bytes

```
⪫⪪⪫Ｅθθ+Ｌθ'

```

Try it online! Link is to verbose version of code. Explanation: Port of @DLosc's Pip answer.

```
    θ       Input string
   Ｅ        For each character
     θ      Input string
  ⪫         Join on
      +     Literal string `+`
 ⪪          Split into substrings of length
        θ   Input string
       Ｌ    Length
⪫           Join on
         '  Literal string `'`
            Implicitly print

```
> 4 votes

# Answer
# Husk, ~~14~~ 12 bytes

```
J''CL¹J'+SRL

```

Try it online!

Same approach as DLosc's Pip answer.
> 4 votes

# Answer
# BQN, 28 bytesSBCS

```
{1↓⥊'''∾˘⟨∘,≠⟩⥊«⥊'+'∾˘⊢⌜˜}

```

Run online!

returns the required string with a trailing space. removing that is +3 bytes

-3 from ovs, who also made this tacit.
> 4 votes

# Answer
# BQN, 21 bytesSBCS

```
1↓⟜∾'+'∾¨1↓↑∾⟜'''⊸∾¨↓

```

Run online!

gives prefixes,  suffixes. `∾⟜'''⊸∾¨` joins matching prefixes and suffixes with a single quote.

`1↓` removes the first value (empty prefix ∾ `'` ∾ full suffix).

`'+'∾¨` prepends a + to each string.

`1↓⟜∾` joins the string and removes the leading +.
> 4 votes

# Answer
# Vyxal, 8 bytes

```
ż\'vṀ\+j

```

Try it Online!

No clue why I didn't vectorise insertion back in May 2022.

## Explained

```
żƛ?n\'Ṁ;\+j
żƛ           # For each item P in the range 1...len(input):
  ?n         #   at position P in the input,
    \'Ṁ      #   insert a '
       ;     # end map
        \+j  # join the result of that on "+" and output

```
> 3 votes

# Answer
# Bash, 83 bytes

```
v=$1\';until [ $v = \'$1 ];do a+=+$v;v=$(sed "s/\(.\)'/'\1/"<<<$v);done;echo ${a#+}

```

Try it online!
> 3 votes

# Answer
# C (gcc), 151 bytes

```
char*f(char*c,int i){char*r=malloc(i*i+i);char*o=r,k=0;for(;k<i*i;){if(k++%(i+1)==1)*o++='\'';*o++=*c++;if(!*c)c-=i,*o++='+';}*o--=0;*o='\'';return r;}

```

Try it online!
> 3 votes

# Answer
# Julia, 40 bytes

```
a->join([replace(a,i=>"$i'") for i=a],+)
```

Attempt This Online!

only works with all letters distinct

#### Why...how...what...

* `[... for i=a]` iterates over the string `a` and store the results in a list. `i` will be a `Char`
* `replace(a, i=>"$i'")` replaces all characters equal to `i` in `a` with `i` followed by `'`
* `join(..., +)` joins the list with `+`. `+` is a function, so `string(+) == "+"` is used.
> 3 votes

# Answer
# , ~~81~~ 75 bytes

*Edit: -6 bytes thanks to Giuseppe*

```
function(x,`[`=substring)paste0(x[1,s<-1:nchar(x)],"'",x[s+1],collapse="+")

```

Try it online!

Takes advantage of the vectorization of `substring` over a vector of start (& optional end) positions, and the vectorization & argument-recycling of `paste` over vectors of strings to paste-together.
> 3 votes

# Answer
# PARI/GP, 58 bytes

```
s->j=strjoin;j([j(strsplit(s,c),Str(c"'"))|c<-Vec(s)],"+")
```

Attempt This Online!

A port of @att's Mathematica answer.

Takes input as a string. Requires that all characters are unique.

---

## PARI/GP, 64 bytes

```
v->i=1;j=strjoin;j([j([j(v[1..i]),j(v[i++..#v])],"'")|c<-v],"+")
```

Attempt This Online!

Take input as a list of characters.
> 2 votes

# Answer
# Retina 0.8.2, 14 bytes

```
.
$`$&'$'+
.$

```

Try it online! Link includes test cases. Explanation:

```
.
$`$&'$'+

```

Replace each variable with its derivative multiplied by the other terms plus a trailing `+` to sum the products together.

```
.$

```

Delete the final trailing `+`.
> 2 votes

# Answer
# Jelly, 10 bytes

```
ṛ€j”+sLj”'

```

Try it online!

Port of DLosc's Pip solution.

```
 €            For each element of the input,
ṛ             the input.
  j”+         Join the copies on +,
     sL       split into slices the length of the original input,
       j”'    and join those on '.

```

# Jelly, 11 bytes

```
p”'$ṛ¦ⱮJj”+

```

Try it online!

The implementation of `¦` is a bit questionable, so we can just abuse Cartesian product to make up for it.

```
p”'            Pair each element of the input with ',
   $ṛ¦         and substitute elements of that into the input at
      ⱮJ       each index, individually.
        j”+    Join the results on +.

```

Just for fun, without preserving the order within each term, `ṭ-Ƥ”'ż@j”+` is 10.
> 2 votes

# Answer
# 05AB1E, 10 bytes

```
ε''«Nǝ}'+ý

```

Try it online!

```
ε       map each character (of the implicit input string)
 ''     the ' character
 «      append to the current character
 N      iteration index
 ǝ      replace the input string's Nth character with the character with ' appended
}
'+      + character
ý       join by it

```
> 2 votes

# Answer
# J-uby, 52 bytes

```
~(:-&(:& &:sub|(~:%&(~:+&?')|:& &:*)|:|&A)).D|~:*&?+

```

Attempt This Online!
> 2 votes

# Answer
# Japt, 10 bytes

```
¬£i''YÄÃq+

```

Try it

```
¬£i''YÄÃq+     :Implicit input of string
¬              :Split
 £             :Map each element at 0-based index Y
  i''          :  Insert a "'" in U at index
     YÄ        :  Y+1
       Ã       :End map
        q+     :Join with "+"

```
> 2 votes

# Answer
# Raku, 33 bytes

This 33-byte version works on my local recent-ish version of Raku:

```
{join '+',.comb.map:{S/$^a/$a'/}}

```

But not on TIO, where a few more bytes are needed:

```
{join '+',.comb.map:{$^a;S/$a/$a'/}}

```

Try it online!

* `.comb` breaks the input string up into characters.
* `.map: { S/$^a/$a'/ }` returns a copy of the input string for each character, in which that character has had an apostrophe appended to it.
* `join '+'` joins those strings together with a plus sign.
> 2 votes

# Answer
# Ruby, 36 bytes

```
->s{s.chars.map{s.sub _1,'\&\''}*?+}

```

Attempt This Online!
> 2 votes

# Answer
# Rust, 123 bytes

```
|s|{s.chars().map(|_|s).enumerate().map(|(i,w)|{let mut t=w.repeat(1);t.insert(i+1,'\'');t}).collect::<Vec<_>>().join("+")}

```

Attempt This Online!
> 2 votes

# Answer
# Nibbles, 8.5 bytes

```
*"+".,,@::`<$_"'"

```

Attempt This Online!

```
*"+".,,@::`<$_"'"
*                 Join
 "+"               with plus signs
    .              for each n in
     ,              range
      ,              length
       @              input
        :           join
         :           join
          `<$         first n characters of
             _         input
              "'"      apostrophe
                      the remaining characters

```
> 2 votes

# Answer
# Burlesque, 20 bytes

```
JiTjiS{''IC}Z[[-'+IC

```

Try it online!

```
J      # Dup
iT     # Tails
j      # Swap
iS     # Heads
{
 ''IC  # Intercalate "'"
}Z[    # Zip together then map
[-     # Drop head
'+IC   # Intercalate "+"

```
> 1 votes

# Answer
# SNOBOL4 (CSNOBOL4), 99 bytes

```
	I =INPUT
N	X =X + 1
	I TAB(X) . L REM . R	:F(O)
	O =O "+" L "'" R	:(N)
O	O TAB(1) REM . OUTPUT
END

```

Try it online!
> 1 votes

# Answer
# Go, 150 bytes

```
import."strings"
type S=string
func f(s S)S{b:=[]S{}
for i:=range s{k:=""
for j,d:=range s{k+=S(d)
if j==i{k+="'"}}
b=append(b,k)}
return Join(b,"+")}

```

Attempt This Online!

### Using `strings.Builder`, 153 bytes

```
import."strings"
func f(s string)string{b:=Builder{}
w:=b.WriteRune
for i:=range s{w('+')
for j,d:=range s{w(d)
if j==i{w('\'')}}}
return b.String()[1:]}

```

Attempt This Online!

Surprised how the string builder solution is only 3 bytes logner.
> 1 votes

# Answer
# Pyth, 10 bytes

```
j\+XRQ\'Sl

```

Test suite

##### Explanation:

```
j\+XRQ\'Sl  | Full code
j\+XRQ\'SlQ | with implicit variables
------------+--------------------------------------
    R   SlQ | For each 1-index d of the input,
   X Q\'    |  Insert "'" into the input at index d
j\+         | Join with "+"

```
> 1 votes

# Answer
# Python 3, 78 bytes

```
d=lambda s:'+'.join([s[:i+1]+'\''+s[i+1:] for i in range(len(s)) if i<len(s)])

```

Try it online!
> 1 votes

# Answer
# Arturo, ~~49~~ 47 bytes

```
$[s][join.with:"+"map..1size s=>[insert s&"'"]]

```

Try it
> 1 votes

# Answer
# Swift, 104 bytes

Swift is not a good language for string manipulation.

```
{var a=""
for i in $0.indices{var c=$0
c.insert("'",at:c.index(after:i))
a+=c+"+"}
a.popLast()
return a}

```

Type is `(String) -> String` or `(Substring) -> String`.

Try it online!

## Swift, 120 bytes

Not as short, but I think the way that `a` is initialized is funny.

```
{var a=[$0][..<0]
for i in $0.indices{var c=$0
c.insert("'",at:c.index(after:i))
a+=[c]}
return a.joined(separator:"+")}

```

If `Array` is a pair `(pointer, length)`\*, indexing operations are like this:

* `arr[i]` = `*(pointer + i)`
* `arr[i...]` = `(pointer + i, length - i)`
* `arr[..<i]` = `(pointer, i)`
* `arr[...i]` = `(pointer, i + 1)`
* `arr[i..<j]` = `(pointer + i, j - i)`
* `arr[i...j]` = `(pointer + i, j - i + 1)`
* `arr[...]` = `(pointer, length)`\**

So, `arr[..<0]` creates an empty collection no matter what the array is. Here, I'm applying it to the array literal `[$0]`, which you would never do in practice but is technically legal.

\*It's actually a triplet `(pointer, length, capacity)`, but this gets the point across.

\**Yes, this exists for some reason.
> 1 votes

# Answer
# Java 8, 93 bytes

```
s->s.chars().mapToObj(e->s.replace((char)e+"",(char)e+"'")).reduce("",(a,b)->a==""?b:a+"+"+b)

```

Try it online!

---

Another close solution not using Streams (**95 bytes**) :

```
s->{String r="";for(String c:s.split("")){r+=(s+"+").replaceAll(c+"\\+{0,1}",c+"'");}return r;}

```
> 1 votes

# Answer
# Thunno, \\$ 16 \log\_{256}(96) \approx \\$ 13.17 bytes

```
LDIA*'+sjsAp''sj

```

Attempt This Online!

Port of DLosc's Pip answer.

#### Explanation

```
LDIA*'+sjsAp''sj  # Implicit input
LD                # Length of input and duplicate
  IA*             # Wrap input in a list and repeat length times
     '+sj         # Join this list by "+"
         sAp      # Split into chunks of length of input
            ''sj  # Join this list by "'"
                  # Implicit output

```
> 1 votes

# Answer
# Scala, 72 bytes

Golfed version. Try it online!

```
def f(l:String)=l.indices.map(i=>l.patch(i,s"${l(i)}'",1)).mkString("+")

```

Ungolfed version

```
object Main extends App {
  def f(l: String): String = {
    l.indices.map(i => l.patch(i, s"${l(i)}'", 1)).mkString("+")
  }

  val test = List("abc", "de", "longabcdef", "short")
  val expected = List(
    "a'bc+ab'c+abc'",
    "d'e+de'",
    "l'ongabcdef+lo'ngabcdef+lon'gabcdef+long'abcdef+longa'bcdef+longab'cdef+longabc'def+longabcd'ef+longabcde'f+longabcdef'",
    "s'hort+sh'ort+sho'rt+shor't+short'"
  )
  val checker = List("✖", "✓")

  for ((t, e) <- test zip expected) {
    val s = f(t)
    val checkMark = checker(if (s == e) 1 else 0)
    println(s"$t -> $s $checkMark")
  }
}

```
> 1 votes

# Answer
# Uiua SBCS, 21 bytes

```
/$"_+_"≡⍜↻(⊂@')+1⇡⊃⧻¤

```

Try it!
> 1 votes

# Answer
# Javascript (Node.js), 114 bytes

```
function a(b){l=b.length;s="";for(i=1;i<=l;i++){if(i>1)s+='+';s+=b.substring(0,i)+"'"+b.substring(i,l);}return s;}

```

Try it online!

I forgot to add "+"

+10 characters to remove a trailing '+'
> 0 votes

---
Tags: code-golf, string, calculus
---


title: Recursive Z-matrix

text:
# Question
Title: Recursive Z-matrix
Originally from a CMC I proposed for the last BMG event

## Challenge

Given a non-negative integer \\$n\\$, create a 2D array of size \\$2^n × 2^n\\$ which is generated in the following manner:

1. Divide the matrix into four quadrants of size \\$2^{n-1} × 2^{n-1}\\$.
2. Visiting order of the four quadrants is defined to be the Z-shape (top-left, top-right, bottom-left, then bottom-right).
3. Recursively apply the ordering (steps 1-2) to each quadrant, until the ordering is defined for each cell in the matrix.
4. Visit each cell in the defined order, sequentially writing down 0, 1, 2, 3, ... to each cell.

You can output 1-based instead of 0-based (add 1 to all cells in the examples below).

Standard rules apply. The shortest code in bytes wins.

## Examples

```
n = 0:
[[0]]

n = 1:
[[0, 1],
 [2, 3]]

n = 2:
[[0,  1,  4,  5],
 [2,  3,  6,  7],
 [8,  9,  12, 13],
 [10, 11, 14, 15]]

n = 3:
[[0,  1,  4,  5,  16, 17, 20, 21],
 [2,  3,  6,  7,  18, 19, 22, 23],
 [8,  9,  12, 13, 24, 25, 28, 29],
 [10, 11, 14, 15, 26, 27, 30, 31],
 [32, 33, 36, 37, 48, 49, 52, 53],
 [34, 35, 38, 39, 50, 51, 54, 55],
 [40, 41, 44, 45, 56, 57, 60, 61],
 [42, 43, 46, 47, 58, 59, 62, 63]]

```

Brownie points for beating or tying with my ~~9~~ 6 bytes in Jelly or 19 bytes in J.

# Answer
# Jelly, 5 bytes

```
4ṗ>2Ġ

```

Try it online! Neil saved a byte. Thanks!

```
4ṗ       n-th Cartesian power of [1,2,3,4]
  >2     Replace 1,2,3,4 with 0,0,1,1
    Ġ    Group indices by values

```

For n=2, `4ṗ>2` generates a list of 16 pairs:

```
 1: [0, 0]
 2: [0, 0]
 3: [0, 1]
 4: [0, 1]
 5: [0, 0]
 6: [0, 0]
 7: [0, 1]
 8: [0, 1]
 9: [1, 0]
10: [1, 0]
11: [1, 1]
12: [1, 1]
13: [1, 0]
14: [1, 0]
15: [1, 1]
16: [1, 1]

```

Then `Ġ` turns this into a list of lists:

* all indices where the value is `0 0`: namely `[ 1, 2, 5, 6]`,
* all indices where the value is `0 1`: namely `[ 3, 4, 7, 8]`,
* all indices where the value is `1 0`: namely `[ 9, 10, 13, 14]`,
* all indices where the value is `1 1`: namely `[11, 12, 15, 16]`.

The result is a Z-matrix.

## Why does this work?

Suppose we bump the input up to \\$n=3\\$. How does our output compare to that for \\$n=2\\$?

Let's call `4ṗ>2` a function \\$f\\$. The list of 16 pairs above is the output of \\$f(2)\\$. By definition of the Cartesian power, \\$f(3) = \[0,0,1,1\] \times f(2)\\$ =

```
  [[0] + x for x in f(2)]     (indices  1..16)
+ [[0] + x for x in f(2)]     (indices 17..32)
+ [[1] + x for x in f(2)]     (indices 33..48)
+ [[1] + x for x in f(2)]     (indices 49..64)

```

There are now 8 values for `Ġ` to group the indices by, ranging from `[0,0,0]` to `[1,1,1]`.

All indices in `1..32` end up in the first four groups, and because we repeated `[[0] + x for x in f(2)]` twice, `a[i+16]` is always in the same group as `a[i]`. Then, all indices in `33..64` end up in the last four groups, and their relative order is the same as the first half of the list (the only difference is they have `1` instead of `0` in front). Within each quadrant, the order is the same as it was in `f(2)`. So we get:

```
[0, 0, 0]:  1  2  5  6 17 18 21 22
[0, 0, 1]:  3  4  7  8 19 20 23 24
[0, 1, 0]:  9 10 13 14 25 26 29 30
[0, 1, 1]: 11 12 15 16 27 28 31 32

[1, 0, 0]: 33 34 37 38 49 50 53 54
[1, 0, 1]: 35 36 39 40 51 52 55 56
[1, 1, 0]: 41 42 45 46 57 58 61 62
[1, 1, 1]: 43 44 47 48 59 60 63 64

```
> 19 votes

# Answer
# Charcoal, ~~27~~ ~~23~~ 18 bytes

```
≔ＥＸ²Ｎ↨⁴↨ι²υＩＥ⊗υ⁺ιυ

```

Try it online! Link is to verbose version of code. Explanation: Now a port of @tsh's Python answer.

```
≔ＥＸ²Ｎ↨⁴↨ι²υ

```

Map the numbers up to 2ⁿ by converting them to base 2 and then interpreting the result as base 4.

```
ＩＥ⊗υ⁺ιυ

```

Vectorised add the doubled array to itself as a matrix.

Previous ~~27~~ 23 byte answer uses a different approach that is interesting in its own right because some golfing languages have built-ins for things like group index by value or sort array by key.

```
≔ＥＸ⁴Ｎ↨²÷↨ι⁴¦²υＩＥ⊕⌈υ⌕Ａυι

```

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

```
≔ＥＸ⁴Ｎ↨²÷↨ι⁴¦²υ

```

For each element in the final diagram, convert it to base 4, halve the digits, then convert from base 2. This gives the row number for that element.

```
ＩＥ⊕⌈υ⌕Ａυι

```

For each row list the elements assigned to that row.
> 6 votes

# Answer
# , ~~65~~ ~~63~~ ~~59~~ ~~57~~ 55 bytes

```
for(i in seq(l=scan()))T=T^4%x%2^rbind(0:1,2:3);log2(T)

```

Try it online!

Relies on the observation that each next iteration `N` relates to the previous iteration `P` by Kronecker product-like operation with the matrix `M` for `n=1`, but using summation instead of multiplication:

```
N = kronecker(4*P, M, FUN = "+")

```

Here, we transform the values with powers and logarithms in order to make use of the actual Kronecker product function, which is conveniently abbreviated as `%x%` operator.
> 6 votes

# Answer
# , 21 19 bytes

Found Bubbler's solution. :-)

```
2(*+/])4#.2#:@i.@^]

```

Try it online!

The first row is A000695, the first column is the same but times two. In the comments Marc LeBrun notes that this is rebasing `n` from base 2 into base 4. So that's helpful with J:

```
2(*+/])4#.2#:@i.@^]
          2   i.@^]  0 1 … 2^n
           #:        digits in base 2
                      0 0 0
                      0 0 1
                          …
                      1 1 0
                      1 1 1
       4#.           interpret as digits in base 4
                      0 1 4 5 16 17 20 21
2(*+/])              addition table of list*2 with list

```
> 6 votes

# Answer
# Python 2, 81 bytes

```
R=[int(bin(x)[2:],4)for x in range(2**input())]
for y in R:print[x+y*2for x in R]

```

Try it online!
> 5 votes

# Answer
# Wolfram Language (Mathematica), ~~49~~ 32 bytes

```
PositionIndex@Tuples[0|0|1|1,#]&

```

Try it online!

Returns an `Association` where the values correspond with rows of the matrix (associations are ordered). For a traditional list-of-lists output, prepend `List@@` (+6 bytes).

Port of Lynn's Jelly solution.

---

Older solution:

```
Outer[#+##&,t=#~FromDigits~4&/@{0,1}~Tuples~#,t]&

```

Try it online!
> 4 votes

# Answer
# , 33 bytes

```
,~@^~&2$4(/:[:<.@-:4#.inv])@i.@^]

```

Try it online!

The true insight of this idea is due to Neil's excellent answer.

I'll mention a couple of other insights I had which didn't help with bytes but might be useful to others:

1. The deltas along each row follow a repeating pattern whose unique elements are http://oeis.org/A007583: 1, 3, 11, 43, 171...
2. Similarly the deltas going down each column follow http://oeis.org/A047849: 1, 2, 6, 22, 86, 342...

## how

* `4...@i.@^]` Produces `0, 1, ... 4^n`, where n is the input.
* `/:[:<.@-:4#.inv]` Sort those according to row number, which we get by converting to base 4, halving, and flooring.
* `,~@^~&2$` Now shape that into a square matrix of side length `2^n`.
> 4 votes

# Answer
# Jelly, 9 bytes

```
4*Ḷb4:2ḄĠ

```

Try it online!

Pretty much a port of Neil's Charcoal answer, so go upvote them too.

```
4*Ḷb4:2ḄĠ    Main Link; take x
4*           4 ** x
  Ḷ          [0, 1, 2, ..., 4 ** x - 1]
   b4        Convert to base 4
     :2      Floor-divide by 2
       Ḅ     Convert from binary
        Ġ    Group the indices by their corresponding values

```

This answer can be ported directly into yuno.

# yuno, 9 bytes

```
4*ʀb4:2ɮG

```

Try it online!

---

This approach is probably what Bubbler's 6-byter is going for, but a) I need to find how to get `0,0,1,1` in 3 bytes, and b) I need to figure out how to make the format valid without the `s2*$` being that long:

# Jelly, 11 bytes

```
Ø.x2¤ṗỤs2*$

```

Try it online!

# Jelly, 16 bytes

```
,Uz0$FU
2*ḶBçþ`Ḅ

```

Try it online!
> 4 votes

# Answer
# , ~~57~~ ~~52~~ 49 bytes

*Edit: -8 bytes by changing to iterative approach with `scan()` for input (instead of recursive function)*

```
for(i in seq(l=2*scan())-1)F=t(cbind(F,F+2^i));+F

```

Try it online!

**Explanation**

~~Recursive function that~~ Each loop `i`creates two side-by-side copies of the columns of a matrix, adding ~~the maximum value of the first +1~~ `2^i` to the second copy, and then transposes it all.

This creates the top left & top right copies of the last matrix (`F`), transposed: so, every second loop creates a full 'Z' of four copies of the last matrix, untransposed.

So we need to loop twice the value of `n` times.
> 4 votes

# Answer
# ARM A32 machine code, 48 bytes

```
e3a0c000 e1a0215c e1b03152 112fff1e
e02c3112 ec432b10 f2800e00 f2810590
f2211110 f480180d e28cc001 eafffff4

```

Following the AAPCS, this takes in r0 the address of an *array* of *array*s of 32-bit integers to be filled with the result, and takes n in r1.

Assembly:

```
.section .text
.code 32
.global zm
zm:
    mov r12, #0               @ Initialise counter to 0
loop:
    asr r2, r12, r1           @ Shift counter right by n, to obtain the top half
    asrs r3, r2, r1           @ Shift that right by n again, to check for ending (2^2n)
    bxne lr                   @ If so, return
    eor r3, r12, r2, lsl r1   @ Exclusive-or cancels out the top half, to obtain the bottom half
    vmov d0, r2, r3           @ Place the two halves into a 64-bit SIMD&FP register
    vmull.p8 q0, d0, d0       @ Multiply by itself as polynomials over GF(2);
               @ because this is a ring of characteristic 2, (a+b)^2 = a^2 + b^2 holds,
               @ thus squaring doubles each exponent in the polynomial, spreading out the bits.
               @ Results go into Q0, which is also D0 and D1.
    vshl.i64 d0, #1           @ Shift the squared top half left by 1 ...
    vorr d1, d0               @ ... and combine by OR with the bottom half
    vst1.32 {d1[0]}, [r0]!    @ Place the result into memory, advancing the pointer
    add r12, #1               @ Increment the counter
    b loop                    @ Repeat

```

(This should work on v7 or later from what I can tell, but in the testing setup I'm using (similar to here), it gives 'instruction not supported' errors on versions lower than v8.4, which doesn't seem right.)
> 4 votes

# Answer
# Jelly,  ~~14 12~~  9 bytes

### 9 byter

```
Ø.ṗḅ4+þḤ$

```

**Try it online!**

##### How?

Makes a table under \\$x+2y\\$ where both \\$x\\$ and \\$y\\$ are a prefix of the Moser–De Bruijn sequence.

```
Ø.ṗḅ4+þḤ$ - Link: n
Ø.        - bits -> [0,1]
  ṗ       - Cartesian power (n) -> binary representations of [0..2^n-1]
   ḅ4     - convert from base four -> first 2^n terms of the Moser–De Bruijn sequence
        $ - last two links as a monad, f(S):
       Ḥ  -   double -> D=[2x for x in S]
      þ   -   (s in S) table with (d in D):
     +    -     (s) add (d)

```

(Still no clue for the 6 byte code though!)

---

### 12 Byter

-1 thanks to caird coinheringaahing (we can use 1-based values)

```
Ġ;"FL+ƊZƊ⁸¡⁺

```

A monadic Link that accepts a non-negative integer and yields the matrix (using the 1-based values option).

**Try it online!** (Footer formats as a Python list for clarity.)

##### How?

Starts with `[[1]]` and repeatedly applies a function that creates the required new entries to the right (and then below, via the transposed matrix) by adding the number of current elements to each element.

```
Ġ;"FL+ƊZƊ⁸¡⁺ - Link: non-negative integer, n
Ġ            - group indices (of implicit [n]) by value -> [[1]] (our initial M)
          ¡  - repeat...
         ⁸   - ...number of times: chain's left argument -> n
        Ɗ    - ...what: last three links as a monad, f(M):
      Ɗ      -   last three links as a monad, g(M):
   F         -     flatten
    L        -     length
     +       -     add to elements of M
  "          -   zip (rows of M) with (rows of that) and apply:
 ;           -     concatenation
       Z     -   transpose
           ⁺ - repeat last link (i.e. repeat f() n more times)

```
> 3 votes

# Answer
# Python 3 with numpy, ~~91~~ ~~87~~ ~~74~~ ~~70~~ ~~69~~ 67 bytes

```
lambda n:c_[s:=(c_[:1<<n]&(b:=1<<r_[:n]))@b]*2+s
from numpy import*
```

Attempt This Online!

`c_[:1<<n]` produces an array of the numbers 0 (inclusive) to 2n (exclusive), vertically. `r_[:n]` produces 0 (inc.) to n (exc.), horizontally, and the left shift makes the corresponding powers of 2, saved in `b`. The bitwise AND, with broadcasting, produces a table of values, splitting the numbers into their constituent bits. The matrix product (`@`) is equivalent to taking the dot product of each broken-down number with `b`; this squares each nonzero bit, changing powers of 2 into powers of 4, and adds them back up, producing a horizontal sequence 0, 1, 4, 5, 16, 17, 20, 21, ..., which is saved in `s`. `c_` makes a vertical version of that sequence, which is doubled and added to the original sequence, which (by broadcasting) produces a table of sums, which is the required result.

---

Previously:

```
from numpy import*
lambda n:(s:=vander([4],8)@unpackbits(mat(r_[:1<<n],uint8),0))+s.T*2

```

`r_[:1<<n]` produces 0 (inc.) to 2n (exc.). It is passed to `mat` ("Interpret the input as a matrix"), which in particular makes it 2-dimensional, while also changing its type to `uint8`, which is necessary for the next function `unpackbits`, which breaks the numbers down into their binary representations, across the added dimension. `vander([4],8)` produces an array of powers of 4, which in the `@` (matrix multiplication) gets dot-producted with each binary vector, producing a similar `s` (transposed differently), and the rest is similar.
> 3 votes

# Answer
# Haskell, 61 bytes

```
f n|r<-iterate(>>= \x->[4*x,4*x+1])[0]!!n=[[x+2*y|x<-r]|y<-r]

```

Try it online!

Same idea as tsh's Python, but the list `r` (e.g. `[0,1,4,5,16,17,20,21]`) is created in a Haskell-y way.
> 2 votes

# Answer
# Julia 1.0, 42 bytes

```
!x=(x-=1)<0||(~n=n*2^2x.+!x;[~0 ~1;~2 ~3])

```

Try it online! (one-indexed)
> 2 votes

# Answer
# JavaScript (ES6),  ~~89 86~~  79 bytes

*Saved 7 bytes thanks to @tsh*

```
n=>[...Array(1<<n)].map((_,y,a)=>a.map(g=(k=1,x)=>k>>n?0:x&k|(y&k|g(k*2,x))*2))

```

Try it online!
> 1 votes

# Answer
# Ruby, ~~90...~~ 64 bytes

```
->n{(k=0...2**n).map{|y|k.map{|x|k.sum{|z|[x,y][z%2][z/2]<<z}}}}

```

Try it online!

### Quick explanation:

Every number can be expressed by interleaving the binary digits of the number of row and column.
> 1 votes

# Answer
# Octave/MATLAB, 67 bytes

```
function x=f(n);x=1;for i=1:n;k=4^i/4;x=[x,k+x;2*k+x,3*k+x];end;end

```

Try it online!
> 1 votes

# Answer
# Haskell, 88 bytes

```
f 0=[[1]]
f n=let(#)=zipWith(++);m=n-1;a=f m;[b,c,d]=map(map(+4^m))<$>[a,b,c]in a#b++c#d

```

Try it online!
> 1 votes

# Answer
# Pip, 19 bytes

```
$+{TBg}FB4*^21MC Ea

```

Try it here! Or, here's a 20-byte equivalent in Pip Classic: Try it online!

### Explanation

Inspired by Neil's comment about interleaving bit patterns, we can simply map the x and y coordinates to the right value (example x = 5, y = 9):

```
5    101    0 1 0 1
9   1001   1 0 0 1    10010011   147

```

We're not going to use Pip's `WV` operator to do the interleaving, since it works left-to-right and doesn't fill in with zeros; instead, observe that an equivalent formula is to convert to base 4, multiply the y value by 2, and sum:

```
5    101      1 0 1      10001 (2) =  101 (4) =  17 (10)
9   1001   1 0 0 1    10000010 (2) = 2002 (4) = 130 (10)   147

```

So:

```
$+{TBg}FB4*^21MC Ea
                 Ea   2^input
              MC      Map this function to a grid of coordinate pairs of that size:
  {TBg}                Convert each coordinate to binary
       FB4             Treat as base 4 and convert back to decimal
          *            Multiply by
           ^21         list containing 2 and 1
                       (multiplying the row by 2 and the column by 1)
$+                     Sum

```
> 1 votes

# Answer
# Python 3 \+ numpy, 66 bytes

```
from numpy import*
f=lambda m:c_[m and(f(m-.5),4**m//2+f(m-.5))].T

```

Try it online!

Mostly straight-forward recursion. Splits the 2x2 update into 2 1x2 updates with intermittent transpose.

# Python 3 \+ numpy, 68 bytes

```
from numpy import*
f=lambda m:m and log2(kron(4**f(m-.5),c_[1:3])).T

```

Try it online!

*Just noticed that @Kirill L.'s R answer uses the exact same idea.*

By componentwise exponentiating adding in the original problem becomes multiplying so friend Kronecker can do the heavy lifting for us.

Downside: will get numerically inaccurate pretty fast.
> 1 votes

# Answer
# BQN, 18 bytes

```
{(⍉∾⍉+·≠⥊)⍟2⍟≍≍0}

```

Try it here.

How it works:

```
{(⍉∾⍉+·≠⥊)⍟2⍟≍≍0}
              ≍≍0  # Starting from a rank 2 array containing 0
            ⍟     # apply the preceding function  times where  is the input
  ⍉∾               # The transposed array joined to
    ⍉+·≠⥊          # the transposed array plus the number of elements
 (       )⍟2       # applied twice

```
> 1 votes

# Answer
# APL(Dyalog Unicode), ~~10~~ 9 bytes SBCS

```
⊢⌸2<,⍳⎕⍴4

```

Try it on APLgolf!

```
⊢⌸2<,⍳⎕⍴4
      ⎕⍴4  # 4 repeated n times
     ⍳     # create a n-dimensional array with shape [4, 4, ..., 4] where value equals index
    ,      # flatten
  2<       # greater than 2 (vectorized)
⊢⌸         # group by value and return indices for each group

```

Created with the help of Luminespire.
> 1 votes

# Answer
# flax, 5 bytes

```
G>1ṗ4

```

Port of @lynn's Jelly answer.

Since the lastest version isn't on ATO, here is an image of the program running (in the REPL).

## Explanation

It is same as the jelly answer.

```
G>1ṗ4
G      " group by values
 >1    "   greater than 1
   ṗ4  "   nth cartesian power of [0,1,2,3]

```
> 1 votes

# Answer
# Python 3 \+ numpy, 113 bytes

A very literal (and long) implementation of the spec.

```
lambda n:sum(tile(R(R(c_[:2,2:4].T,1<<i,0),1<<i,1)<<2*i,(1<<n+~i,)*2)for i in r_[:n])
from numpy import*
R=repeat

```

Try it online!
> 0 votes

# Answer
# Japt, 14 bytes

Port of tsh's Python solution.

```
2pU Ç¤n4
£Ë+XÑ

```

Try it
> 0 votes

# Answer
# JavaScript (ES6), 72 bytes

```
f=(n,a=[0])=>n?f(n-1,a.flatMap(x=>[x*=4,x+1])):a.map(y=>a.map(x=>x+2*y))

```

Try it online!

Use the idea from Lynn's Haskell answer which saves many bytes on generating the list for iterate. Using base conversion as what my Python answer do will make it much longer (79 bytes).
> 0 votes

# Answer
# Python 3.8 + numpy, ~~85~~ 84 bytes

```
import numpy;z=lambda n:n and numpy.block([[m:=z(n-1),m+(k:=4**~-n)],[m+2*k,m+3*k]])

```

Not as clever as some of the other python ones, but shorter than the other python 3 ones. Also 0 returns a scalar instead of a list, not sure if that's allowed.
> 0 votes

# Answer
# Pari/GP, 62 bytes

```
n->matrix(2^n,,i,j,(g=k->fromdigits(binary(k-1),4))(i)*2+g(j))

```

Try it online!
> 0 votes

---
Tags: code-golf, array, matrix
---


# Question
Title: How compactly can your language perform accurate numerical integration?
**WARNING: This challenge may need 128 bit floats.1**

The task is to perform numerical integration. Consider the following three functions.

\\$ f(x) = cx^{c - 1}e^{-x^c} \\$

\\$ g\_1(x) = 0.5e^{-x} \\$

\\$ g\_2(x) = 5 e^{-10 x} \\$

We will have that \\$c \geq 0.2\\$. Your code should be correct for any value between 0.2 and 1.0.

The task is to perform the integral:

\\$ \int\_{0}^{\infty} f(x) \ln\left(\frac{f(x)}{g\_1(x)+g\_2(x)}\right) \mathrm{d}x. \\$

Your output should be correct to within plus or minus \\$10^{-10}\\$ percent of the correct answer.

**Input**

A single real number \\$c\\$.

**Output shown to 15 significant places**

```
c = 0.2. Output: 119.2233798550179
c = 0.3. Outout: 8.077346771987397
c = 0.4. Output: 2.07833944013377
c = 0.5. Output: 0.8034827042913205
c = 0.6. Output: 0.3961739639940003
c = 0.7. Output: 0.2585689391629249
c = 0.8. Output: 0.2287758419066705
c = 0.9. Output: 0.2480070283065089
c = 1.0. Output: 0.2908566108890576

```

**Timing**

Your code must complete on TIO before it times out.

As is usual, this is a competition *per* language. In other words, Java coders should not be put off by the single character solution in a language they can't believe exists. You may use any language/library you choose.

**Questions and notes**

* Does this integral actually exist/converge?

Yes. It is the Kullback Leibler divergence between two probability distributions.

* Show me any code that can compute these integrals.

See this Wolfram Alpha example.

* Why such high precision that we need 128 bit floats?

It is quite easy to solve the integral to low precision using simple ideas. However, to get high precision seems to need careful thought both about the integral itself but also the quadrature method that will be used. I felt this made the question more challenging and interesting.

1. For lots of languages this is easy. E.g Fortran has real128, gcc/clang have \_Float128 for C/C++, Python has Decimal (and many other options), Dyalog APL has 128 bit (decimal) floats, C# has the Decimal type, Haskell has rounded, PHP has decimal, MATLAB has the Multiprecision Computing Toolbox, Julia has BigFloat, Java has BigDecimal, JavaScript has big.js, Ruby has BigDecimal, R has Rmpfr, Perl has Math::BigFloat, SOGL has 200 digit decimal precision, etc. If your favourite golfing language doesn't, maybe consider adding it?

# Update

Alexey Burdin has shown that

\\$ \int\_{0}^{\infty} f(x) \ln\left(\frac{f(x)}{g\_1(x)+g\_2(x)}\right) \mathrm{d}x = \int\limits\_0^1 \left(\ln c+\left(\ln \left(-\ln t\right)\right)\cdot\frac{c-1}{c}+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +(-\ln{t})^{1/c}\right)dt\\$

This potentially greatly simplifies the problem.

# Answer
> 7 votes
## python3, 3982 2183 1866 1821 1677, ~14s

```
from functools import*
#let's start with the comments. I'll exclude them from total count
#from functools I need only reduce
#I've been coding in python2 for 5 years so I'm used to built-in reduce much
#and it's very handy here
g=range
#just keep it in mind while reading) ranGe -> g
#and we're getting rid of classes in favor of functions, that gives
#us ~100 bytes more. It's how it was:
#class P:
#it's a Polynomial, self.c are coefficients and self.d is a divisor. So all polynomial computations are exact
#we need polynomials to implement https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature
#   c=[0]
#   d=1
#don't even know if these initial values are necessary
#   def __init__(self,c,d=1):
#       self.c=list(c)
#       self.d=d
#n'th derivative
#   def dn(self,n=1):
#       c=self.c[:]
#       for _ in g(n):
#           c=[c[i]*i for i in g(1,len(c))]
#       d=reduce(gcd,c,self.d)
        #gcd is imported from math along with exp,log and gamma
#       return P([i//d for i in c],self.d//d)
#   def m(self,x):
#the value of poly at point x. __call__ may look better
#       return sum((x**i)*self.c[i] for i in range(len(self.c)-1,-1,-1))/self.d
        #the reverse order was for floats and really needed there for precision,
        #but python3 have exact Fractions, so maybe it's to be golfed too
#Now let's convert this class to functions
#derivative=lambda c:reduce(lambda c,_:[c[i]*i for i in g(1,len(c))],g(n),c)
#the_gcd=lambda c,d:reduce(gcd,c,d)
#the_result=lambda c,d,m:([i//m for i in c],d//m)
#where c is now the derivative
#and combined:
v=lambda c,d,n=1:(lambda c:
(lambda c,d,m:([i//m for i in c],d//m))(c,d,reduce(gcd,c,d))
)\
(reduce(lambda c,_:[c[i]*i for i in g(1,len(c))],g(n),c))
#the computing of a plynomial's value:
m=lambda c,d,x:sum((x**i)*c[i] for i in g(len(c)-1,-1,-1))/d
d={}
#https://en.wikipedia.org/wiki/Combination written recursively with Pascal's triangle, with the dict d for speed
def c(n,k):
    if (n,k) not in d:
        d[(n,k)]=1 if k==0 or k==n or n==0 else c(n-1,k)+c(n-1,k-1)
    return d[(n,k)]
#p=30
#it's the floating point precision for Decimal
#it was needed for computing the polynomial roots, you can see the raw version
#maybe I'll get rid of it if it's used in the only place
from math import*
from fractions import Fraction as F
n=61
l=v(sum(([c(n,i)*(-1)**(n-i),0] for i in g(n+1)),[]),
        2**n*reduce(lambda x,y:x*y,g(1,n+1),1),n)
#Legendre polynomial of order n, pretty straightforward, -3bytes for space
#you can't see these bytes in the tio version, there is 1 line
d=v(*l)
#and the derivative for weights
#see we can use c and d now 'cause we don't need them more
c=2**100
#maybe the hardest-to-read part, the roots of l
#taking too long to compute (2-3 mins) so hard-stored
#b(b'GKSmRyJoJz...') is the 390-bytes byte array,
#representing the <0 roots (the function is odd, so -x and x are both roots and we have to store only half)
#30 13-byte(100-bit, for 10^{-30} precision) integers are stored in lsbf order of bytes
#and c is common multiplier for all computed roots' denominators
from base64 import b64decode as b
#extract numbers from base-64,
r=(lambda b:[sum(2**(n*8)*i for n,i in enumerate(b[j:j+13])) for j in g(0,391,13)])(b(b'GKSmRyJoJzZNd978D7hBJ/SAXF6k3Q+D7w9eWlMN1hwU7v+nhNcPP3k+PPY3ALm75/G0DwxkhnBW7vYeTLvhhw8XTAWijBbPeSIxclAPxUQBLJIlU0r8S8gODz6W/ntGb0p0WeUPww6bS5oF+R0NzuqOe20OpkOWAx5DBcFecEQODqjUMBgWzmCaCiGqpQ2xekr44w56Xl998jMNN0/S8UF6TE04eGm5DL0i6Wm5wA/nHOhgNgxOR3kgsJNHFpZQMKsLdK2AMa1xRv65pzQYCxSYft+TX736FBjQfQqY8CuipXP9Bxm/adwJEKYLfykKeTs/aG00CVuQ6lnO2j04CUVLhggLyUIxC/eIqRKid9IHdSe6aSQRS8hkmWoZB2rhfvNKb9bHP8KfWwb379JeWbyiuY/elZkFsIKVLmMYZwdEhs7TBAvB47snQ/UMwtDNCgRb9soEqlNzn6z8GT8DTfSlRAVDCG07FjtxAheAQ01VIxgWXpy6oQEELmiRzDIqr+ckI9EA'))
#and make fractions
#0 is in r for magically mistaken 391 for 390 )
#maybe +[0]+r[::-1] is 1 byte shorter (or longer?)
#but at least more readable )
#roots computation code is in the raw below
r=[F(i,c) for i in [-j for j in r]+r[-2::-1]]
#the weights
w=[2/((1-x**2)*(m(*d,x))**2) for x in r]
from decimal import*
D=Decimal
setcontext(Context(prec=30, rounding=ROUND_HALF_DOWN))
#that was float=
#it's only converting Fraction to Decimal. no built-in
f=lambda r:(lambda n,d:D(n)/D(d))(*str(r).split('/'))
#number of sub-intervals of (0,1)
d=2**8
s=input().strip()
c=D(s)
#the inverse of c, not the math.exp(1) ))
e=D(1)/c
#pretty straightforward, I is the integrator,
I=lambda h:sum(sum(f(v)*h(f(((x+1)/2+i)/d))/(2*d) for v,x in zip(w,r))
        for i in range(d))
#the rest follows (at math:)
print(D(gamma(1/float(s)+1))+(c.ln())-I(lambda t:(D('0.5')+5*((-9*((-((t).ln()))**(e))).exp())).ln())-D('0.5772156649015328606065120900824')*(1-e)-1)

```

Try it online! -- takes input as a valid Decimal, like `0.2`
Math:
0. We do google numerical integration.
1. Let's see provocative \\$f(x)=c\cdot x^{c-1}\cdot e^{-x^c}= -\left(e^{-x^c}\right)'\\$
2. So, we have \\$-\int\limits\_0^\infty \ln\left(\frac{c\cdot x^{c-1}\cdot e^{-x^c}}{0.5e^{-x}+5e^{-10x}}\right)\left(e^{-x^c}\right)'dx\\$,
we change the variable to \\$t=e^{-x^c}\\$, \\$x^c=-\ln(t)\\$, \\$x=\left(-\ln t\right)^\frac1c\\$,
\\$t(0)=1\\$, \\$t(+\infty)=0\\$ and the integral is
\\$\int\limits\_0^1 \ln\left(\frac{c\cdot x(t)^{c-1}\cdot t}{\left(0.5+5e^{-9x(t)}\right)e^{-x(t)}}\right)dt=\\$
\\$\int\limits\_0^1 \left(\ln c+\left(\ln x(t)\right)\cdot(c-1)+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +x(t)\right)dt=\\$
\\$\int\limits\_0^1 \left(\ln c+\left(\ln \left(-\ln t\right)\right)\cdot\frac{c-1}{c}+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +x(t)\right)dt\\$
3. Next we wolframalpha these things:
\\$\int\limits\_0^1\ln(-\ln t)dt\\$, \\$\int\limits\_0^1 (-\ln t)^a dt\\$, \\$\int\limits\_0^1 \ln t\\,dt\\$ and do the rest numerically (the latter thing maybe can be done with trapeziums or rectangles or Simpson -- the idea for short submissions).
The raw thing:

```
from copy import copy
from functools import reduce
class Poly:
    c=[0]
    d=1
    def __init__(self,c,d=1):
        self.c=list(c)
        self.d=d
    def dn(self,n=1):
        c=copy(self.c)
        for _ in range(n):
            c=[c[i]*i for i in range(1,len(c))]
        d=reduce(gcd,c,self.d)
        return Poly([i//d for i in c],self.d//d)
    def m(self,x):
        return sum((x**i)*self.c[i] for i in range(len(self.c)-1,-1,-1))/self.d
def gcd(x,y):
    return abs(gcd(y,x%y)) if y else x
cdict={}
def c(n,k):
    if (n,k) not in cdict:
        cdict[(n,k)]=1 if k==0 or k==n or n==0 else c(n-1,k)+c(n-1,k-1)
    return cdict[(n,k)]
def root(x,f,f_):
    x_old=0
    while x!=x_old:# and abs(f(x))>1e-15:
        print(x,f(x),f_(x))
        x_old=x
        x=x-f(x)/f_(x)
    print('===')
    return x
precision=30
def b(b,e,f,eps=1e-17):
    mold,m=0,1
    while e-b>eps and mold!=m:
        mold=m
        m=(e+b)/2
        if (l.m(m)>0)!=(l.m(b)>0):
            e=m
        else:
            b=m
        #print(*map(float,[b,e,f(b),f(e)]))
    print('.',sep='',end='')
    return b
##class Frac:
##    n=0
##    d=1
##    def __init__(self,n,d):
##        self.n=n
##        self.d=d
##    def
lcm=lambda x,y:x*y//gcd(x,y)
import math
from fractions import Fraction
degree=61
n=degree
l=Poly(sum(([c(n,i)*(-1)**(n-i),0] for i in range(n+1)),[]),
       2**n*reduce(lambda x,y:x*y,range(1,n+1),1)).dn(n)
ld=l.dn()
#print(l.coeffs,l.d)
#r=[root(math.cos(math.pi*(4*i-1)/(4*n+2)),l.m,ld.m) for i in range(1,n+1)]
##d=2**8#4096
##while True:
##    i_=[(b,e) for b,e in [(Fraction(i,d)-1,Fraction(i+1,d)-1)
##                          for i in range(2*d)] if (l.m(b)>0)!=(l.m(e)>0)]
##    print(d,len(i_))
##    if len(i_)==n:
##        break
##    d*=2
##r=[b(*i,l.m,Fraction(1,10**precision)) for i in i_]
##print('\n',sep='',end='')
cm=1267650600228229401496703205376
#cm=reduce(lambda x,i:lcm(int((str(i).split('/')+[1])[1]),x),r,1);cm
#[str(i*cm) for i in r]
r=[Fraction(int(i),cm) for i in ['-1266681605106811426160584860697', '-1262547799267707740863710314936', '-1255122086390021131342967954015', '-1244422184933863703004183427391', '-1230475806835976018691325584397', '-1213319288075032955237836344344', '-1192997371846808208190052844742', '-1169563069074221666262667138623', '-1143077514053466175428649175963', '-1113609802945881973204995949479', '-1081236812728978305124643427497', '-1046043000289730543348383120049', '-1008120181921745998647915597624', '-967567293702390287305984778942', '-924490133333634877286205966159', '-879001084101412705885157895540', '-831218821664589525108846925845', '-781268004433877062952267018393', '-729278948345971539276192261648', '-675387286880006687633990914140', '-619733617202494872284783888651', '-562463133363442075271048079222', '-503725247500289525435137450347', '-443673200037730123766618583031', '-382463659900222496924931293872', '-320256315780122313563780137228', '-257213459527711400732339664476', '-193499562749975057795207525453', '-129280847722706546595631169560', '-64724853735360852478581550597', '0', '64724853735360852478581550596', '129280847722706546595631169559', '193499562749975057795207525452', '257213459527711400732339664475', '320256315780122313563780137227', '382463659900222496924931293871', '443673200037730123766618583030', '503725247500289525435137450346', '562463133363442075271048079221', '619733617202494872284783888650', '675387286880006687633990914139', '729278948345971539276192261647', '781268004433877062952267018392', '831218821664589525108846925844', '879001084101412705885157895539', '924490133333634877286205966158', '967567293702390287305984778941', '1008120181921745998647915597623', '1046043000289730543348383120048', '1081236812728978305124643427496', '1113609802945881973204995949478', '1143077514053466175428649175962', '1169563069074221666262667138622', '1192997371846808208190052844741', '1213319288075032955237836344343', '1230475806835976018691325584396', '1244422184933863703004183427390', '1255122086390021131342967954014', '1262547799267707740863710314935', '1266681605106811426160584860696']]
#r=[Fraction(-72002513042367095, 72057594037927936), Fraction(-35883766692782835, 36028797018963968), Fraction(-35672715239059617, 36028797018963968), Fraction(-141474422995755015, 144115188075855872), Fraction(-69944451686068079, 72057594037927936), Fraction(-137938433007892963, 144115188075855872), Fraction(-135628098615462501, 144115188075855872), Fraction(-66481963419518119, 72057594037927936), Fraction(-32488216960856169, 36028797018963968), Fraction(-126602776952709289, 144115188075855872), Fraction(-61461197033678349, 72057594037927936), Fraction(-59460660411885473, 72057594037927936), Fraction(-7163124720376473, 9007199254740992), Fraction(-109999665903886841, 144115188075855872), Fraction(-105102359763536115, 144115188075855872), Fraction(-99930853605361871, 144115188075855872), Fraction(-23624659822433771, 36028797018963968), Fraction(-88819888837164975, 144115188075855872), Fraction(-41454708727199633, 72057594037927936), Fraction(-19195661220692463, 36028797018963968), Fraction(-70455555169530143, 144115188075855872), Fraction(-63944654967081315, 144115188075855872), Fraction(-28633465234423143, 72057594037927936), Fraction(-6304975386764939, 18014398509481984), Fraction(-10870271009372985, 36028797018963968), Fraction(-36408927801417385, 144115188075855872), Fraction(-29241784833142379, 144115188075855872), Fraction(-21998353389560073, 144115188075855872), Fraction(-14697530755564293, 144115188075855872), Fraction(-7358363943167303, 144115188075855872), Fraction(0, 1), Fraction(3679181971583651, 72057594037927936), Fraction(3674382688891073, 36028797018963968), Fraction(2749794173695009, 18014398509481984), Fraction(14620892416571189, 72057594037927936), Fraction(4551115975177173, 18014398509481984), Fraction(43481084037491939, 144115188075855872), Fraction(50439803094119511, 144115188075855872), Fraction(57266930468846285, 144115188075855872), Fraction(31972327483540657, 72057594037927936), Fraction(35227777584765071, 72057594037927936), Fraction(76782644882769851, 144115188075855872), Fraction(82909417454399265, 144115188075855872), Fraction(44409944418582487, 72057594037927936), Fraction(94498639289735083, 144115188075855872), Fraction(49965426802680935, 72057594037927936), Fraction(52551179881768057, 72057594037927936), Fraction(13749958237985855, 18014398509481984), Fraction(114609995526023567, 144115188075855872), Fraction(118921320823770945, 144115188075855872), Fraction(122922394067356697, 144115188075855872), Fraction(15825347119088661, 18014398509481984), Fraction(129952867843424675, 144115188075855872), Fraction(132963926839036237, 144115188075855872), Fraction(33907024653865625, 36028797018963968), Fraction(68969216503946481, 72057594037927936), Fraction(139888903372136157, 144115188075855872), Fraction(70737211497877507, 72057594037927936), Fraction(142690860956238467, 144115188075855872), Fraction(143535066771131339, 144115188075855872), Fraction(144005026084734189, 144115188075855872)]
#r=[Fraction(-144104924687379029, 144115188075855872), Fraction(-144061113633417193, 144115188075855872), Fraction(-71991152681636629, 72057594037927936), Fraction(-143868501524294767, 144115188075855872), Fraction(-35929931908622757, 36028797018963968), Fraction(-143536019344824829, 144115188075855872), Fraction(-71658710569760757, 72057594037927936), Fraction(-143063986084641675, 144115188075855872), Fraction(-142775775751233995, 144115188075855872), Fraction(-35613215044775229, 36028797018963968), Fraction(-142095317851501427, 144115188075855872), Fraction(-141703235672459401, 144115188075855872), Fraction(-141276708943911307, 144115188075855872), Fraction(-70407920670826367, 72057594037927936), Fraction(-140320744889653401, 144115188075855872), Fraction(-34947884983137447, 36028797018963968), Fraction(-139228355106226913, 144115188075855872), Fraction(-138631327306447751, 144115188075855872), Fraction(-69000300827755959, 72057594037927936), Fraction(-68668165733468573, 72057594037927936), Fraction(-136638678208163739, 144115188075855872), Fraction(-135907811461286037, 144115188075855872), Fraction(-67571954440908669, 72057594037927936), Fraction(-134347156155496359, 144115188075855872), Fraction(-16689718369143059, 18014398509481984), Fraction(-132655882883583863, 144115188075855872), Fraction(-131761773444627501, 144115188075855872), Fraction(-65417817986076227, 72057594037927936), Fraction(-64938847793634331, 72057594037927936), Fraction(-32222046285398959, 36028797018963968), Fraction(-127867345160661637, 144115188075855872), Fraction(-31703855946358703, 36028797018963968), Fraction(-62866338356003711, 72057594037927936), Fraction(-15577420891180097, 18014398509481984), Fraction(-123475765655790143, 144115188075855872), Fraction(-61151075136161875, 72057594037927936), Fraction(-121098806255952059, 144115188075855872), Fraction(-119866026109883705, 144115188075855872), Fraction(-118604109492524957, 144115188075855872), Fraction(-117313363144639965, 144115188075855872), Fraction(-115994100814789491, 144115188075855872), Fraction(-28661660795766563, 36028797018963968), Fraction(-28317829445786347, 36028797018963968), Fraction(-111868458922669023, 144115188075855872), Fraction(-110438407601984251, 144115188075855872), Fraction(-27245377857813575, 36028797018963968), Fraction(-26874531136490757, 36028797018963968), Fraction(-105988607520833283, 144115188075855872), Fraction(-26113331820545011, 36028797018963968), Fraction(-102892657018719659, 144115188075855872), Fraction(-50653488045428427, 72057594037927936), Fraction(-49848334969235775, 72057594037927936), Fraction(-98062129987227927, 144115188075855872), Fraction(-96403753553428489, 144115188075855872), Fraction(-94721943747436259, 144115188075855872), Fraction(-23254277343922147, 36028797018963968), Fraction(-45644832420663197, 72057594037927936), Fraction(-89540030043462969, 144115188075855872), Fraction(-10971078784389609, 18014398509481984), Fraction(-85975896119833679, 144115188075855872), Fraction(-42081131673510875, 72057594037927936), Fraction(-20582043201506859, 36028797018963968), Fraction(-40237035159487741, 72057594037927936), Fraction(-39300203286200335, 72057594037927936), Fraction(-19176909251924261, 36028797018963968), Fraction(-74796221710411347, 144115188075855872), Fraction(-72866625298407583, 144115188075855872), Fraction(-35459658404464939, 72057594037927936), Fraction(-68954769584591107, 144115188075855872), Fraction(-33486730579157493, 72057594037927936), Fraction(-32487936568629803, 72057594037927936), Fraction(-31481245542875313, 72057594037927936), Fraction(-60933804407252565, 144115188075855872), Fraction(-58890306225406907, 144115188075855872), Fraction(-56832493264165915, 144115188075855872), Fraction(-54760865727051299, 144115188075855872), Fraction(-26337963587783547, 72057594037927936), Fraction(-25289092203398143, 72057594037927936), Fraction(-48468147330210951, 144115188075855872), Fraction(-46346328843725859, 144115188075855872), Fraction(-11053311177256411, 36028797018963968), Fraction(-21034706713097939, 72057594037927936), Fraction(-19957678053844251, 72057594037927936), Fraction(-4718949543956531, 18014398509481984), Fraction(-35578660114658865, 144115188075855872), Fraction(-33397075583856055, 144115188075855872), Fraction(-15603686524289111, 72057594037927936), Fraction(-29010084771446223, 144115188075855872), Fraction(-1675359053686717, 9007199254740992), Fraction(-24594889131807827, 144115188075855872), Fraction(-22378054994346827, 144115188075855872), Fraction(-20155781304247947, 144115188075855872), Fraction(-4482152060343901, 36028797018963968), Fraction(-15697077176510773, 144115188075855872), Fraction(-6730865269878557, 72057594037927936), Fraction(-5611555844344803, 72057594037927936), Fraction(-4490882388138869, 72057594037927936), Fraction(-6738234618615355, 144115188075855872), Fraction(-2246533281244659, 72057594037927936), Fraction(-2246806352819167, 144115188075855872), Fraction(0, 1), Fraction(1123403176409583, 72057594037927936), Fraction(4493066562489317, 144115188075855872), Fraction(3369117309307677, 72057594037927936), Fraction(8981764776277737, 144115188075855872), Fraction(11223111688689605, 144115188075855872), Fraction(13461730539757113, 144115188075855872), Fraction(3924269294127693, 36028797018963968), Fraction(17928608241375603, 144115188075855872), Fraction(10077890652123973, 72057594037927936), Fraction(11189027497173413, 72057594037927936), Fraction(12297444565903913, 72057594037927936), Fraction(26805744858987471, 144115188075855872), Fraction(14505042385723111, 72057594037927936), Fraction(31207373048578221, 144115188075855872), Fraction(16698537791928027, 72057594037927936), Fraction(2223666257166179, 9007199254740992), Fraction(37751596351652247, 144115188075855872), Fraction(39915356107688501, 144115188075855872), Fraction(42069413426195877, 144115188075855872), Fraction(44213244709025643, 144115188075855872), Fraction(23173164421862929, 72057594037927936), Fraction(24234073665105475, 72057594037927936), Fraction(50578184406796285, 144115188075855872), Fraction(52675927175567093, 144115188075855872), Fraction(27380432863525649, 72057594037927936), Fraction(28416246632082957, 72057594037927936), Fraction(29445153112703453, 72057594037927936), Fraction(15233451101813141, 36028797018963968), Fraction(62962491085750625, 144115188075855872), Fraction(64975873137259605, 144115188075855872), Fraction(66973461158314985, 144115188075855872), Fraction(34477384792295553, 72057594037927936), Fraction(70919316808929877, 144115188075855872), Fraction(36433312649203791, 72057594037927936), Fraction(37398110855205673, 72057594037927936), Fraction(76707637007697043, 144115188075855872), Fraction(78600406572400669, 144115188075855872), Fraction(80474070318975481, 144115188075855872), Fraction(82328172806027435, 144115188075855872), Fraction(84162263347021749, 144115188075855872), Fraction(42987948059916839, 72057594037927936), Fraction(87768630275116871, 144115188075855872), Fraction(11192503755432871, 18014398509481984), Fraction(91289664841326393, 144115188075855872), Fraction(93017109375688587, 144115188075855872), Fraction(47360971873718129, 72057594037927936), Fraction(12050469194178561, 18014398509481984), Fraction(49031064993613963, 72057594037927936), Fraction(99696669938471549, 144115188075855872), Fraction(101306976090856853, 144115188075855872), Fraction(51446328509359829, 72057594037927936), Fraction(104453327282180043, 144115188075855872), Fraction(52994303760416641, 72057594037927936), Fraction(107498124545963027, 144115188075855872), Fraction(108981511431254299, 144115188075855872), Fraction(55219203800992125, 72057594037927936), Fraction(55934229461334511, 72057594037927936), Fraction(113271317783145387, 144115188075855872), Fraction(114646643183066251, 144115188075855872), Fraction(57997050407394745, 72057594037927936), Fraction(29328340786159991, 36028797018963968), Fraction(29651027373131239, 36028797018963968), Fraction(14983253263735463, 18014398509481984), Fraction(60549403127976029, 72057594037927936), Fraction(122302150272323749, 144115188075855872), Fraction(61737882827895071, 72057594037927936), Fraction(124619367129440775, 144115188075855872), Fraction(125732676712007421, 144115188075855872), Fraction(126815423785434811, 144115188075855872), Fraction(31966836290165409, 36028797018963968), Fraction(128888185141595835, 144115188075855872), Fraction(129877695587268661, 144115188075855872), Fraction(130835635972152453, 144115188075855872), Fraction(32940443361156875, 36028797018963968), Fraction(66327941441791931, 72057594037927936), Fraction(133517746953144471, 144115188075855872), Fraction(67173578077748179, 72057594037927936), Fraction(135143908881817337, 144115188075855872), Fraction(33976952865321509, 36028797018963968), Fraction(68319339104081869, 72057594037927936), Fraction(137336331466937145, 144115188075855872), Fraction(138000601655511917, 144115188075855872), Fraction(69315663653223875, 72057594037927936), Fraction(4350886097069591, 4503599627370496), Fraction(139791539932549787, 144115188075855872), Fraction(17540093111206675, 18014398509481984), Fraction(140815841341652733, 144115188075855872), Fraction(70638354471955653, 72057594037927936), Fraction(17712904459057425, 18014398509481984), Fraction(71047658925750713, 72057594037927936), Fraction(142452860179100915, 144115188075855872), Fraction(71387887875616997, 72057594037927936), Fraction(71531993042320837, 72057594037927936), Fraction(143317421139521513, 144115188075855872), Fraction(35884004836206207, 36028797018963968), Fraction(143719727634491027, 144115188075855872), Fraction(71934250762147383, 72057594037927936), Fraction(143982305363273257, 144115188075855872), Fraction(18007639204177149, 18014398509481984), Fraction(36026231171844757, 36028797018963968)]
w=[2/((1-x**2)*(ld.m(x))**2) for x in r]
from decimal import Decimal,Context,setcontext,ROUND_HALF_DOWN
setcontext(Context(prec=30, rounding=ROUND_HALF_DOWN))
float=lambda r:(lambda x:Decimal(x[0])/Decimal(x[1]))(str(r).split('/'))
d=2**8
s='{120.000000000000000000000000000,9.26052826812554737559992621571,3.32335097044784255118406403126,2.00000000000000000000000000000,1.50457548825155601882809780906,1.26582350605728327736200759570,1.13300309631934634747833911121,1.05218372089129332722399903882,1.00000000000000000000000000000}'
l=list(map(Decimal,s[1:-1].split(',')))
gamma=Decimal('0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495')
for (c,c_1,real),l_ in zip([
    (Decimal('0.1')*i,Decimal('10')/i,10/i) for i in range(2,11)],l):
    #let t=e^{-x^c},\ln t=-x^c,x=(-\ln t)^{\frac1c}
    #t(0)=1,t(\infty)=0
    #dt=-e^{-x^c}\cdot x^{c-1}\cdot c
    #f(x)=e^{-x^c}\cdot x^{c-1}\cdot c=t\cdot (-\ln t)^{\frac{c-1}c}\cdot c
    xt=lambda t:(-((t).ln()))**(c_1)
##    f=lambda x:c*(x**(c-1))#*(math.exp(-(x**c)))
    fln=lambda t:((-((t).ln())).ln())#+(t.ln())
##    g=lambda x:Decimal('0.5')+5*((-9*x).exp())#*math.exp(-x)
    g=lambda t:Decimal('0.5')+5*((-9*xt(t)).exp())#*math.exp(-x)
##    h=lambda x:f(x)*((-(x**c)).exp())*\
##       ((f(x)).ln()-(x**c)-((g(x)).ln()-x))
##    #0,inf->0,1: x=t/(t-1), xt-x=t, (x-1)t=x, t=x/(x-1)
##    h1=lambda t:h(t/(1-t))/((1-t)**2) # 0,1
##    h2=lambda v:h1((v+1)/2)/2
    h1=lambda t:g(t).ln()#f(t).ln()-((g(t).ln())-xt(t))
    int_=lambda h1:sum(
          sum(float(w_)*h1(float(((x+1)/2+i)/d))/(2*d) for w_,x in zip(w,r))
          for i in range(d))
    print(c,Decimal(math.gamma(real+1))+(c.ln())-int_(h1)-gamma*(1-c_1)-1)
h1_=r'''0.2 0.476044892154085268532810576285
0.3 0.326045243249124607083205564208
0.4 0.194544295792217104304557844994
0.5 0.0805857800502670858642542661562
0.6 -0.0176136562407463311404194282348
0.7 -0.102042234943715186336125809639
0.8 -0.174612380676135222634071174702
0.9 -0.237048749195042896134197564780
1.0 -0.290856610889057599214369459143'''
x=[i.split(': ') for i in r'''c = 0.2. Output: 119.2233798550179
c = 0.3. Outout: 8.077346771987397
c = 0.4. Output: 2.07833944013377
c = 0.5. Output: 0.8034827042913205
c = 0.6. Output: 0.3961739639940003
c = 0.7. Output: 0.2585689391629249
c = 0.8. Output: 0.2287758419066705
c = 0.9. Output: 0.2480070283065089
c = 1.0. Output: 0.2908566108890576'''.split('\n')]
y=[i.split(' ') for i in r'''0.2 119.223379855017945799292478451
0.3 8.07734677198740392689540304712
0.4 2.07833944013376923348715605315
0.5 0.80348270429132046532502570246
0.6 0.39617396399400051391328713871
0.7 0.25856893916292307961043114604
0.8 0.22877584190665494099938058972
0.9 0.24800702830645810748432393277
1.0 0.29085661088905759921436945914'''.split('\n')]
#[(lambda l:(l and min(l)) or 0)([n for n,(d1,d2) in enumerate(zip(i[1],j[1])) if d1!=d2]) for i,j in zip(x,y)]
#[0, 14, 15, 17, 17, 16, 15, 14, 17]

```

# Answer
> 5 votes
# Pari/GP, ~~88~~ 87 bytes, \< 1 sec

Great work from @AlexeyBurdin!

I also saw the "provocative" way in which \\$f\\$ was a derivative of an easier function, but only had the idea to use this for partial integration, without success.

For golfing reasons I recombined all the logarithms, except the one giving \\$-1\\$ because it makes no difference for the code length and subtracting one seemed slightly more efficient than multiplying the integrand by \\$t\\$. I also took the logarithms of \\$1/t\\$ instead of \\$t\\$ to be able to get rid of the otherwise needed parentheses for the minus sign. The `2` in the last argument to `intnum` reduces the integration step size. That is needed to achieve the required precision.

```
l=log;i(c)=intnum(t=0,1,l(c*l(1/t)^(1-1/c)/(.5+5*exp(-9*l(1/t)^c^-1)))+l(1/t)^c^-1,2)-1

```

Try it online!

# Answer
> 5 votes
# Python 3, ~~150~~ ~~133~~ 122 bytes

*(thanks to Delta and Alexey Burdin for further golfing this)*

SymPy's mpmath module is a great way to handle this problem. All math functions (`exp`, `log`, `euler`, `gamma` and `quad` for integration) are run in high precision.

```
lambda c:gamma(1/c+1)+l(c)-euler*(1-1/c)-1-quad(lambda t:l(.5+5*exp(-9*((-l(t))**(1/c)))),[0,1])
from mpmath import*
l=log

```

Try it online!

# Answer
> 4 votes
# Wolfram Language (Mathematica), 105 bytes ~7.5 seconds

Mathematica is good at automatically determining the best method needed to perform the integration. Setting the working precision high gets us to results that match the stated answer. The time stated is to perform all the example integrals.

```
NIntegrate[((f=#*x^(#-1))*Log[f/(E^x^#*(5/E^(10*x)+0.5/E^x))])/E^x^#,{x,0,-Log@0},WorkingPrecision->128]&

```

Try it online!

# Answer
> 3 votes
# Desmos, ~~62~~ 59 bytes

```
a=ln(1/t)^{1/c}
f(c)=∫_0^1(ln(a^c/a/e(.5+5e^{-9a})c)+a)dt

```

This is a port of Christian Siever's solution with some modifications for numeric precision.

Try it on Desmos!

*-3 bytes thanks to Aiden Chow*.

The final column in the table shows the relative error, which reaches 5.71×10^-13 as the worst case at c=0.2. This barely slips under the 10^-12 relative error requirement.

As for the speed requirement, this run in under 1500ms on my computer.

## Precision notes

Desmos uses 64-bit floats for everything, so it's tough to get the desired precision, considering the machine precision is on the order of 10^-16. However, Desmos's quadrature routine is pretty good for general integrals because it works by recursively splitting the integration interval when it needs more precision.

`ln(1/t)` can be changed to `(-lnt)` for -1 byte, but this costs about two orders of magnitude of relative precision.

# Answer
> 1 votes
# APL(NARS), 929 chars, time 1' 7''

```
⎕fpc←128⋄grv←gerr←0⋄ge←q←0v⋄funOne←{f←q×(⍵*q-1v)×*-⍵*q⋄f×⍟f÷(0.5v×*-⍵)+5×*-10v×⍵}

r←(f S)w
gerr←1⋄r←¯1J1⋄→0×⍳3≠≢w⋄→0×⍳w[3]>⎕fpc÷3.322⋄gerr←grv←0⋄ge←÷10v*w[3]⋄r←f It w[1]w[2]
→0×⍳∼gerr≠0⋄r←¯1J1

r←(f I)w;a;b;n;nn;h;s1;s2;h1;h2;i;t;p;m
r←⍬⋄→0×⍳gerr≠0⋄h2←h1←s2←s1←0v⋄(a b n)←w⋄→2×⍳n≤0⋄h←(2×nn←2v×n)÷⍨b-a⋄i←0v⋄→3
r←⍬⋄gerr←1⋄→0
t←f a+h×2×i+←1⋄s1+←t⋄→4×⍳∼0=2∣i⋄h1+←t⋄→5
h2+←t
→3×⍳i<nn-1⋄i←0v⋄→2×⍳∞=∣s1
s2+←f a+h×¯1+2×i+←1⋄→6×⍳i<nn⋄→2×⍳∞=∣s2⋄p←f¨a b⋄→2×⍳∞=+/∣p⋄t←+/p
p←(3÷⍨h×t+(4×s2)+2×s1)(3÷⍨h×2×t+(4×h2)+2×h1)
m←(h×s2+s1+t÷2)(h× 2×s1+t÷2)⋄→2×⍳(∞=∣m[1])∨∞=∣p[1]
s2←∣-/p⋄s1←∣-/m⋄→10×⍳∼s2>s1⋄r←m,s1⋄→0
r←p,s2

r←(f It)w;t;y;m
r←0v⋄→2×⍳0≠gerr⋄→3×⍳∼grv>250⋄gerr←1
r←0v⋄→0
grv+←1⋄t←f I w,1⋄→2×⍳0≠gerr⋄→6×⍳t[3]≥1⌊ge×1E5⋄→5×⍳ge>t[3]⋄t←f I w,19⋄→2×⍳0≠gerr
→6×⍳∼ge>t[3]
y←f I w,23⋄→2×⍳0≠gerr⋄→6×⍳ge≤∣y[1]-t[1]⋄grv-←1⋄r←↑y⋄→0
m←2v÷⍨-/⌽w⋄r←(f It w[1],w[1]+m)+(f It(w[1]+m),w[2])⋄grv-←1

:for q :in 0.1v+10÷⍨⍳9⋄q,(11⍕ funOne S 1E¯64v (7E8v) 11)⋄:endfor

```

82+ 8+82+18+2 +39+74+13+40+5+25+63+44+50+37+6+10+ 15+35+7+79+12+54+58+6 +65=929

"grv" is one global variable used for make not too deep the recursive function

"gerr" is one global variable used for find some possible error

"ge" is one global variable used for store epsilon of precision

"funOne" would be the function use in exercise where the "c" is here "q".

"S" is the function to call for integration of one continue function in one interval, it has arguments on the left the function, in the right 3 numbers, 2 are a b for the integration interval \[a, b\], and one for the precision of calculation. It is possible that S function fail (return ¯1J1 in that case) or not show the error in some integral even if the number result is wrong (example some integral of periodic functions). It is possible S give the wrong value to the function that stop the computation with "DOMAIN ERROR".

"I" is a function for calculate the integral in 2 differents way the trapezes and Simpson formula, return the two sums have less difference.

"It" is a recursive function from the model of qsort, that use function "I" and change partition for recalculate the result each little interval it find ok for to be a little more sure of result even if that means 2x slow (so it is possible one silent fail).

```
  :for q :in 0.1v+10÷⍨⍳9⋄q,(11⍕ funOne S 1E¯64v (7E8v) 11)⋄:endfor
0.2  119.22337985502
0.3  8.07734677200
0.4  2.07833944015
0.5  0.80348270430
0.6  0.39617396400
0.7  0.25856893917
0.8  0.22877584191
0.9  0.24800702831
1  0.29085661089

```

1E¯64v would be one value not 0 sufficient little, the value for +oo that not make wrong results in function funOne, it seems to be 7E8v.It gets 1' 7'' for all value of c (here it is q variable) from 0.1 to 1. The function integrate of sys (∫) is more fast, but the right result it seems to be only the first

```
     q←0.2v⋄30⍕  1E¯70v funOne ∫7E8v
 119.223379855016584445335168065350
     q←0.3v⋄30⍕  1E¯70v funOne ∫7E8v
 0.253368882611387716254867856476
     q←0.4v⋄30⍕  1E¯70v funOne ∫7E8v
 0.016195181203412657409325335505
     q←1v⋄30⍕  1E¯70v funOne ∫7E8v
 0.001677768739402358961260938889

```

For some function f, +/f¨a b is not the same +/f a b. One has to eliminate all the places where is used "=" or use ⎕ct←0 in each function use "=".

In the function "I" these are the 4 sums, the first (p\[1\]) the more precise Simpson sum, than the second (p\[2\]) less precise Simpson the 3rd the more precise use trapezes (m\[1\]), the less precise use trapezes (m\[2\]).

```
p←(3÷⍨h×t+(4×s2)+2×s1)(3÷⍨h×2×t+(4×h2)+2×h1)
m←(h×s2+s1+t÷2)(h× 2×s1+t÷2)

```

# Answer
> 1 votes
# C (gcc), 1933 bytes, 3''

```
    #include <quadmath.h>
    #include <stdio.h>
    #define C signed char
    #define U unsigned
    #define I int
    #define LD __float128
    #define G goto
    #define P printf
    #define R return
    #define FM FLT128_MAX
    I PP(C*a,C*fmt,LD v){I r;C arr[1010];arr[0]=0;r=quadmath_snprintf(arr,1000,fmt,v);if(r<0)R r;R P("%s%s",a,arr);}
    I grv, gerr; LD ge,t[3],y[3],gq=0.2q;I Er(LD x){if(x<0)R-x>=FM;R x>=FM;}
    LD*In(LD*rf,LD(*f)(LD),LD a,LD b,U n){LD  h,h1,h2,s1,s2,t,x,p[2],m[2],r;U nn,i;if(gerr)R rf;h2=h1=s2=s1=0;if(n==0||n>1e6){L0:gerr=1;R rf;}nn=2*n;h=(b-a)/(2*nn);i=0;L1:++i;s1+=t=f(a+h*2*i);if(Er(t))G L0;if(i%2==0)h1+=t;else h2+=t;if(i<nn-1)G L1;i=0;if(Er(s1))G L0;L2:++i;s2+=t=f(a+h*(2*i-1));if(Er(t))G L0;if(i<nn)G L2;if(Er(s2))G L0;t=f(a);if(Er(t))G L0;x=f(b);if(Er(x))G L0;x=f(b);t+=x;p[0]=h*(t+4*s2+2*s1)/3;p[1]=2*h*(t+4*h2+2*h1)/3;m[0]=h*(s2+s1+t/2);m[1]=2*h*(s1+t/2);if(Er(p[0])||Er(m[0]))G L0;s1=fabsq(p[0]-p[1]);s2=fabsq(m[0]-m[1]);if(s1<s2){rf[0]=p[0];rf[1]=p[1];rf[2]=s1;}else{rf[0]=m[0];rf[1]=m[1];rf[2]=s2;}R rf;}
    LD It(LD(*f)(LD),LD a,LD b){LD r,m;if(gerr){L3: r=0;R r;}if(grv>250){gerr=1;G L3;}++grv;In(t,f,a,b,1);if(gerr)G L3;if(t[2]>=1||t[2]>=ge*1e5)G L6;if(ge>t[2])G L4;In(t,f,a,b,19);if(gerr)G L3;if(ge>t[2]){L4:In(y,f,a,b,23);if(gerr)G L3;if(ge<=fabsq(y[0]-t[0]))G L6;--grv;R y[0];}L6:m=(b-a)/2;r=It(f,a,a+m)+It(f,a+m,b);--grv;R r;}
    LD Si(LD(*f)(LD),LD a,LD b,U p){LD r; gerr=1;if(p>38)R FM;gerr=grv=0;ge=powq(10,-(LD)p);r=It(f,a,b);R gerr?FM:r;}
    LD fN(LD x){LD f,h,z;f=logq(gq)+(gq-1.0q)*logq(x)-powq(x,gq)+x-logq(0.5q+expq(-9.0q*x));return gq*powq(x,gq-1.0q)*expq(-powq(x,gq))*f;}
    LD fO(LD w){LD f=gq*powq(w,gq-1.0q)*expq(-powq(w,gq)),z;z=0.5q*expq(-w)+5.0q*expq(-10.0q*w);R f*logq(f/z);}
    LD funOne(LD x){R x<1000?fO(x):fN(x);}
    I main(){LD r[3],q,i,f,a,b;for(gq=0.2q;gq<=1;gq+=0.1q){q=Si(funOne,1e-64,7e8,11);PP("","%.2Qf",gq);PP(", ","%.30Qf",q);P("\n");}R 0;}

```

Try it online!

ungolfed:

```
// > gcc this.c -lquadmath
#include <quadmath.h>
#include <stdio.h>

#define C  signed char
#define U  unsigned
#define I  int
#define LD __float128
#define G   goto
#define P   printf
#define R   return
#define FM  FLT128_MAX

I PP(C*a,C*fmt,LD v)
{I r;C arr[1010];
 arr[0]=0;r=quadmath_snprintf(arr,1000,fmt,v);if(r<0)R r;
 R P("%s%s",a,arr);
}

I grv, gerr; LD ge,t[3],y[3],gq=0.2q;

I Er(LD x){if(x<0)R-x>=FM;R x>=FM;}

LD*In(LD*rf,LD(*f)(LD),LD a,LD b,U n)
{LD  h,h1,h2,s1,s2,t,x,p[2],m[2],r;U nn,i;
     if(gerr)R rf;
     h2=h1=s2=s1=0;if(n==0||n>1e6){L0:gerr=1;R rf;}
     nn=2*n;h=(b-a)/(2*nn);i=0;
 L1: ++i;s1+=t=f(a+h*2*i);if(Er(t))G L0;if(i%2==0)h1+=t;else h2+=t;if(i<nn-1)G L1;i=0;if(Er(s1))G L0;
 L2: ++i;s2+=t=f(a+h*(2*i-1));if(Er(t))G L0;if(i<nn)G L2;if(Er(s2))G L0;t=f(a);if(Er(t))G L0;x=f(b);if(Er(x))G L0;x=f(b);t+=x;
     p[0]=h*(t+4*s2+2*s1)/3;p[1]=2*h*(t+4*h2+2*h1)/3;
     m[0]=h*(s2+s1+t/2);m[1]=2*h*(s1+t/2);
     if(Er(p[0])||Er(m[0]))G L0;
     s1=fabsq(p[0]-p[1]);s2=fabsq(m[0]-m[1]);
     if(s1<s2){rf[0]=p[0];rf[1]=p[1];rf[2]=s1;}
     else     {rf[0]=m[0];rf[1]=m[1];rf[2]=s2;}
     R rf;
}

LD It(LD(*f)(LD),LD a,LD b)
{LD r,m;
    if(gerr){L3:  r=0;  R r;}
    if(grv>250){gerr=1;G L3;}
    ++grv;In(t,f,a,b,1);if(gerr)G L3;if(t[2]>=1||t[2]>=ge*1e5)G L6;if(ge>t[2])G L4;
    In(t,f,a,b,19);if(gerr)G L3;
    if(ge>t[2])
      {L4:In(y,f,a,b,23);if(gerr)G L3;if(ge<=fabsq(y[0]-t[0]))G L6;--grv;R y[0];}
L6: m=(b-a)/2;r=It(f,a,a+m)+It(f,a+m,b);--grv;
    R r;
}

LD Si(LD(*f)(LD),LD a,LD b,U p){LD r; gerr=1;if(p>38)R FM;gerr=grv=0;ge=powq(10,-(LD)p);r=It(f,a,b);R gerr?FM:r;}

LD fN(LD x)
{LD f,h,z;
 f=logq(gq)+(gq-1.0q)*logq(x)-powq(x,gq)+x-logq(0.5q+expq(-9.0q*x));
 return gq*powq(x,gq-1.0q)*expq(-powq(x,gq))*f;
}
LD fO(LD w)
{LD f=gq*powq(w,gq-1.0q)*expq(-powq(w,gq)),z;
 z=0.5q*expq(-w)+5.0q*expq(-10.0q*w);
 R f*logq(f/z);
}
LD funOne(LD x){R x<1000?fO(x):fN(x);}

I main()
{LD r[3],q,i,f,a,b;
 for(gq=0.2q;gq<=1;gq+=0.1q)
       {q=Si(funOne,1e-64,7e8,11);
        PP("","%.2Qf",gq);PP(", ","%.30Qf",q);P("\n");
       }
 R 0;
}

```

This would be one port of APL NARS answer to C....

"grv" is one global variable used for make not too deep the recursive function

"gerr" is one global variable used for find some possible error

"ge" is one global variable used for store epsilon of precision

"funOne" would be the function use in exercise where the "c" is here "q".

"Si" is the function to call for integration of one continue function in one interval, it has arguments on the left the function, in the right 3 numbers, 2 are a b for the integration interval \[a, b\], and one for the precision of calculation. It is possible that Si function fail (return FLT128\_MAX in that case) or not show the error in some integral even if the number result is wrong (example some integral of periodic functions). It is possible S give the wrong value to the function that stop the computation.

"In" is a function for calculate the integral in 2 differents way the trapezes and Simpson formula, return the two sums have less difference.

"It" is a recursive function from the model of qsort, that use function "In" and change partition for recalculate the result each little interval it find ok for to be a little more sure of result even if that means 2x slow (so it is possible one silent fail).

```
0.20, 119.223379855017788373126945485560
0.30, 8.077346772003659028828166522766
0.40, 2.078339440146599411449845152080
0.50, 0.803482704300521900602181696767
0.60, 0.396173964001614436175399637559
0.70, 0.258568939171353519263491445983
0.80, 0.228775841912709730320218436446
0.90, 0.248007028314673262106569598459
1.00, 0.290856610889722950747717956860

```

1E¯64 would be one value not 0 sufficient little, the value for +oo that not make wrong results in function funOne, it seems to be 7E8.It gets 3'' for all value of c (here it is q variable) from 0.1 to 1.

There was some problem but one see the different behaviors of functions one know what to do: Cut and paste 2 function in 1...

```
fO                fN
args: 1..10*114000
inf, 7.146180987711778877830383715838e-05
inf, 1.783936268160517786280124237027e-05
inf, 7.145191389859739510802225203435e-06
inf, 3.548880901445954254961267159123e-06
inf, 1.999831511094571465956022639528e-06
inf, 1.225727400360527391484147042523e-06
inf, 7.980629270833437955446020180540e-07
inf, 5.439501158831436480969842592496e-07
inf, 3.843424608223386417573364622201e-07

fO                                             fN
args: 1..10*10
3.675246391455555827226018131679e-02, 3.675246391455555827226018131679e-02
4.382838598960554328826955735253e-02, 4.382838598960554328826955735253e-02
4.457247966758306754921782013739e-02, 4.457247966758306754921782013739e-02
4.397010507150066869496325771966e-02, 4.397010507150066869496325771966e-02
4.298551944407423573690051874755e-02, 4.298551944407423573690051874755e-02
4.190240448327716144966191450276e-02, 4.190240448327716144966191450276e-02
4.082071320292452338988660547692e-02, 4.082071320292452338988660547692e-02
3.977794012002491268811702806721e-02, 3.977794012002491268811702806721e-02
3.878734785898533961995102431946e-02, 3.878734785898533961995102431946e-02

fO                                           fN
args: 1/1..10
-8.243106934757116432656018641562e-01, -1.344494569228624722492452318334e-01
-1.290429743446295599947997195038e+00, 2.015992534869343736972957346392e-01
-1.460826278910182993294146550044e+00, 7.738815018537929006207034760483e-01
-1.455905510003998243340646950725e+00, 1.487514582923200556163817499455e+00
-1.330519622929357039421233297650e+00, 2.297630406711373536179161098540e+00
-1.115252743630621385523306314288e+00, 3.179047478956116831216159688529e+00
-8.293358937016342749872510061477e-01, 4.115949553797663332308038109261e+00
-4.858491347654902257043425759574e-01, 5.097618989802348559994390744177e+00
-9.419553179254290987076726693566e-02, 6.116387165599517173497262275172e+00

```

---
Tags: code-golf, math
---