title: Is it a factor of a polynomial?

text:
# Question
Title: Is it a factor of a polynomial?
A polynomial is divisible by a factor `(x-n)` if `f(n)=0` for a function `f`. Your job: to determine if a polynomial function `f(x)` is divisible by `(x-n)`.

# The input

The input is in the form of `(x-n), (Polynomial)`. **Remember, if n is negative, `(x-n)` will be in the input form of `(x+n)`**. For the polynomial, all exponents will be put in as `^`. Coefficients will be written next to the variable `x`. An example polynomial could be `2x^2 + x^1`. There will be no spaces between anything. The term `x` will be inputed as `x^1`. So what would "normally" look like `(x - 1)` will be `(x^1-1)`. The coefficients and powers will **always** be integers. **The coefficient one will be implicit if it is just `x`. I.e., `x` can be interpreted as `1x`**

# The output

A boolean value. Truthy, or Falsey.

*Thanks to @AlexA. For helping me clarify this!*

# Examples

```
Input:(x^1-1),(x^1-1)
Output: True

Input: (x^1+2),(2x^2+4x^1+2)
Output: False

Input: (x^1+7),(x^2-49)
Output: True

```

### Rules

* This is , so shortest code in bytes wins

Unfortunately, I don't know how to implement the snippet leaderboard. If anyone knows how, feel free to edit the post.

# Answer
# Pyth - 39 bytes

This is a monstrous combination of regexp and eval. I like the approach, but will try to improve the implementation.

It uses the Polynomial Remainder Theorem.

```
K_sPe:z"-|\+"3!v.ssXPtw,\^\x,"**""*K"\*

```

Doesn't work online because of eval use.
> 5 votes

# Answer
# Casio Basic, 19 bytes

```
judge(mod(b,a)=0

```

As it turns out, the fx-CP400 can do `mod` on algebraic expressions!

Polynomial and factor should be entered as expressions. 16 bytes for the code, 3 bytes to enter `a,b` into the parameter value box.
> 3 votes

# Answer
# MATLAB, ~~103 99 97 95~~ 93 bytes

I'm trying some different things, and got this to work to save a couple of bytes:

```
eval([regexprep(input(''),{'.+?1(.+)\),','(\d)x'},{'x=str2num(''$1'');disp(~','$1\*x'}) 41]);

```

If I can reduce that down even further, I will post an explanation.

---

Old code an explanation

```
t=sscanf(input(''),'(x^1%d),%s')';x=-t(1);disp(~eval(regexprep([t(2:end) ''],'(\d)x','$1\*x')))

```

This also works with **Octave**. You can try it online. I've save the program as a script named `isFactor.m`, so you can just enter `isFactor` at the prompt. \[Note: in Octave spits out a warning while running - MATLAB doesn't generate this\].

The input must be in the format `'(x^1+7),(x^2-49)'` as per the question. The quotes marks are added so MATLAB/Octave knows it is a string.

The output is either a `0` or a `1` depending on whether or not it is true or false.

---

So, the code works as follows. First we request an input, and then parse it. The parsing string extracts the signed number after the first `(x^1` in the string - this is our value of `n`. Then it continues to extract the string (`%s`) after the `),` in the input - this is our expression.

```
t=sscanf(input(''),'(x^1%d),%s')';

```

Next, we extract the value of `n` and set `x` equal to it - we are going to evaluate whether the expression equals zero when `n==x`, so this is why we store the value to x. Also we negate the extracted number, because of the minus sign when parsing.

```
x=-t(1);

```

We will then display the output which is a boolean

```
disp(

```

The output is basically the logical negation of our evaluated equation. If `f(x)` is zero, this will return 1, otherwise it will result in zero.

```
     ~eval(

```

We are evaluating the input expression, but to do this, we need to reformat it slightly so that MATLAB can understand. When we read the string, it is actually an array of `double` type, so we need to convert that to a character array. Before the conversion we also get rid of the first element as that is what we used for `n`. We then need to replace any occurrence of `x` which is preceded by a number (e.g `4x`) by the same thing but with a multiplication (`*`) sign in between so MATLAB can calculate it.

```
           regexprep(char([t(2:end) ''],'(\d)x','$1\*x')
     )
)

```
> 1 votes

# Answer
# VBScript, ~~118~~ 116 bytes

```
a=inputbox(""):for i=0 to 9:a=replace(a,i&"x",i&"*x"):next:b=split(a,","):x=-eval(b(0)):msgbox not cbool(eval(b(1)))

```

Since we know that the first part of the input is a linear polynomial, we only need to check if its root matches that of the second polynomial; and we need to prepare the term for `eval` by inserting `*` as needed.
> 1 votes

# Answer
# Axiom ~~77~~ 180 Bytes

```
f(a:UP(x,INT),b:UP(x,INT)):Boolean==(ground?(a)or ground?(b)=>false;p:=b;r:=a;if degree(a::POLY INT,x)>degree(b::POLY INT,x)then(p:=a;r:=b);(p rem r)$UP(x,FRAC INT)~=0=>false;true)

```

the previous solution

```
v(a,b)==(ground?(a) or ground?(b) or (b rem a)$UP(x,FRAC INT)~=0=>false;true)

```

was wrong because it assume degree(b)\>=degree(a) one bug i wrote... test and results

```
(3) -> f(x^1-1,x^1-1)
   (3)  true
                                                            Type: Boolean
(4) -> f(x^1+1,2*x^2+4*x^1+2)
   (4)  true
                                                            Type: Boolean
(5) -> f(x^1+2,2*x^2+4*x^1+2)
   (5)  false
                                                            Type: Boolean
(6) -> f(x^1+7,x^2-49)
   (6)  true
                                                            Type: Boolean
(7) -> f(1, 1)
   (7)  false
                                                            Type: Boolean
(8) -> f(1, x^2+1)
   (8)  false
                                                            Type: Boolean
(9) -> f(x^8-1, x^2-1)
   (9)  true

```
> 1 votes

# Answer
# Mathematica, 30 bytes

Try it online!

```
PolynomialRemainder[#2,#,x]<1&

```
> 0 votes

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


title: Is this polynomial a square?

text:
# Question
Title: Is this polynomial a square?
Given an integral polynomial \\$p\\$, determine if \\$p\\$ is a square of another integral polynomial.

An *integral polynomial* is a polynomial with only integers as coefficients.

For example, \\$x^2+2x+1\\$ should gives truthy, because \\$x^2+2x+1 = (x+1)^2\\$.

On the other hand, \\$2x^2+4x+2\\$ should gives falsy: \\$2x^2+4x+2 = (\sqrt{2}x+\sqrt{2})^2\\$. but \\$\sqrt{2}x+\sqrt{2}\\$ is not an integral polynomial.

## Input

A polynomial, in any reasonable format. For example, the polynomial \\$x^4-4x^3+5x^2-2x\\$ may be represented as:

* a list of coefficients, in descending order: `[1,-4,5,-2,0]`;
* a list of coefficients, in ascending order: `[0,-2,5,-4,1]`;
* a list of pairs of `(coefficient, degree)`, in any order: `[(1,4),(-4,3),(5,2),(-2,1),(0,0)]`;
* a map with degrees as keys and coefficient as values: `{4:1,3:-4,2:5,1:-2,0:0}`;
* a string representation of the polynomial, with a chosen variable, say `x`: `"x^4-4*x^3+5*x^2-2*x"`;
* a built-in polynomial object, e.g., `x^4-4*x^3+5*x^2-2*x` in PARI/GP.

## Output

A value representing whether the polynomial is a square. You can choose to

* output truthy/falsy using your language's convention (swapping is allowed), or
* use two distinct, fixed values to represent true (affirmative) or false (negative) respectively.

This is , so the shortest code in bytes wins.

## Testcases

Here I use coefficient lists in descending order:

### Truthy

```
[]
[25]
[1,2,1]
[1,2,1,0,0]
[1,4,0,-8,4]
[4,28,37,-42,9]
[4,0,-4,4,1,-2,1]
[1,-12,60,-160,240,-192,64]

```

### Falsy

```
[-1]
[24,1]
[1,111]
[2,4,2]
[1,2,1,0]
[1,3,3,1]
[1,-4,5,-2,0]
[4,0,-4,4,1,2,1]
[1,-9,30,-45,30,-9,1]

```

# Answer
> 6 votes
# JavaScript (Node.js), 102 bytes

```
a=>a.some(n=>(b[j=i]=b.map(m=>n-=m*~~b[j--])|b[0]?n/2/b[0]:n**.5)%(1/a[i+i++]?1:1/0)!=0,b=[],i=0)|1&~i

```

Try it online!

Input a list of coefficients, in descending order. The list should either contain only a zero or the first value in the list should be non-zero.

I *believe* it is correct. But I cannot *prove* it with my math knowledge. It at least passed all testcases.

# Answer
> 4 votes
# Python3, 289 bytes:

```
R=range
def m(a,b):
 r={}
 for x,y in a:
  for j,k in b:r[y+k]=r.get(y+k,0)+x*j
 return[(r[i],i)for i in r]
def c(d,k=[]):
 if[]==d:yield k;return
 for x in R(abs(d[0][0])+1):
  for e in R(d[0][1]+1):
   for i in[1,-1]:yield from c(d[1:],k+[(i*x,e)])
f=lambda x:any(x==m(i,i)for i in c(x))

```

Try it online!

Basic brute force.

# Answer
> 4 votes
# Vyxal, ~~33~~ 32 bytes

```
h√d£0¾ȦλȮḂÞ•-¥/J;RI÷?₂J†^¥J:⌊=JA

```

Try it Online!

–1 per emanresu A.

## How?

We find the square root of the polynomial using its power series, checking if the result has several integers followed by zeroes, plus additional checks.

```
           #                   stack (top ->)
h          # head              p0
 √         # square root       sqrt(p0)
  d        # double            2sqrt(p0)
   £       # set register    = 2k

0¾Ȧ        # top[0] = []       [[],p1,…,p(2n)]

λ          # λ(acc,cur)        acc cur
           #                 = [a1,…,an] p(n+1)
 Ȯ         # over              [a1,…,an] p(n+1) [a1,…,an]
  Ḃ        # dup and reverse   [a1,…,an] p(n+1) [a1,…,an] [an,…,a1]
   Þ•      # dot product       [a1,…,an] p(n+1) a1an+…+ana1
     -     # subtract          [a1,…,an] p(n+1)-a1an-…-ana1
      ¥    # push register     [a1,…,an] p(n+1)-a1an-…-ana1 2k
       /   # divide            [a1,…,an] (p(n+1)-a1an-…-ana1)/2k
           #                 = [a1,…,an] a(n+1)
        J; # join              [a1,…,an,a(n+1)]

R          # reduce            [a1,…,a(2n)]
 I         # halve             [[a1,…,an],[a(n+1),…,a(2n)]]
  ÷        # split to stack    [a1,…,an] [a(n+1),…,a(2n)]
   ?       # push input        [a1,…,an] [a(n+1),…,a(2n)] [p0,…,p(2n)]
    ₂      # length even?      [a1,…,an] [a(n+1),…,a(2n)] 0
     J     # join              [a1,…,an] [a(n+1),…,a(2n),0]
      †    # vectorized not    [a1,…,an] [1,…,1]
           # [if perfect square, then a(n+1) = … = a(2n) = 0,
           # and the length of input, 2n+1, would be odd]

^          # flip stack        [1,…,1] [a1,…,an]
 ¥         # push register     [1,…,1] [a1,…,an] 2k
  J        # join              [1,…,1] [a1,…,an,2k]
   :       # dup               [1,…,1] [a1,…,an,2k] [a1,…,an,2k]
    ⌊      # floor             [1,…,1] [a1,…,an,2k] [a1,…,an,2k]
     =     # equal?            [1,…,1] [1,…,1]
           # [if perfect square, then a0 = 2k, … , an would
           # all be integers]

J          # join              [1,…,1]
 A         # all truthy?       1

```

# Answer
> 3 votes
# Python NumPy, 103 bytes

```
lambda a:any(poly1d(rint((a[0]**.5*poly(sort(around(roots(a),7))[::2])).real))**2-a)
from numpy import*
```

Attempt This Online!

Inverted output: False for squares and True for non squares.

Direct approach using roots/linear factors. Floating point issues may lurk at some point (but not with the given test cases).

Basically, we take every other root, multiply corresponding linear factors and the square root of the leading coefficient, discard imaginary and fractional components, square and compare to the input.

For those not familiar with legacy numpy polynomials: `roots` takes a vector coefficients and returns tthe vector of roots. `poly` takes a vector of roots and returns the vector of coefficients of the monic polynomial with the given roots (with multiplicity). `poly1d` takes a vector of coefficients and returns a polynomial object. This changes semantics of basic algebraic operations.

## Python NumPy, 149 bytes

```
lambda l,a:all([any(poly1d(rint((fft.ifft(fft.fft(a)**.5*(1-2*int_([*bin(i)[3:]]))))[:l//2+1].real))**2-a)for i in r_[1<<l:2<<l]])
from numpy import*
```

Attempt This Online!

Inverted output. Takes the number of coefficients and the coefficient list.

This computes the square root directly in the Fourier domain. Unfortunately, for each Fourier coefficient we have to check both signs, making this approach more complicated than I had hoped.

# Answer
> 3 votes
# Wolfram Language (Mathematica), ~~58 68 64~~ 55 bytes

thanks @Steffan for the byte correction, and the suggestions
thanks @alephalpha

```
Tr[Last/@Mod[a=FactorList@#,2]]<2&&EvenQ[2√a[[1,1]]]&

```

This code factorizes the input polynomial and checks if the exponent of all factors is even and if the square root of the integer factor is an integer.

Try it online!

# Answer
> 2 votes
# Charcoal, 56 bytes

```
﹪Ｌθ²Ｆ⟦§θ⁰↨θΠ⊗⊕↔θ⟧«≔×⁰↨ι⁴ζＦＬζ«§≔ζκ¹§≔ζκ¬‹ιＸ↨²ζ²»¿⁻ιＸ↨²ζ²⎚

```

Try it online! Link is to verbose version of code. Does not support empty lists (use `[0]` instead). Outputs a Charcoal boolean, i.e. `-` if it thinks the polynomial is a square, nothing if not. Explanation: Since @tsh beat me to my original algorithm I've gone for this approach instead but I can only prove that it has no false negatives. If someone can come up with a counterexample I'll delete this answer pending a rewrite using @tsh's approach. Edit: A counterexample was found, but @Nitrodon's claim is that only polynomials of odd degree can be counterexamples, so I've spent 5 bytes checking that the polynomial is of even degree. Edit: Another counterexample was found, so I've spent 6 bytes checking that the leading coefficient is a perfect square.

```
﹪Ｌθ²

```

Output a `-` if the polynomial is of even degree.

```
Ｆ⟦§θ⁰↨θΠ⊗⊕↔θ⟧«

```

Calculate a very large integer given by taking the absolute values of all of the coefficients, incrementing and doubling them, and then taking the product. Evaluate the polynomial at this value. Loop over the leading coefficient of the polynomial and the result.

```
≔×⁰↨ι⁴ζＦＬζ«§≔ζκ¹§≔ζκ¬‹ιＸ↨²ζ²»

```

Perform a binary search for the nearest integer to the square root of this value.

```
¿⁻ιＸ↨²ζ²⎚

```

If its square does not equal the value, clear the output.

# Answer
> 2 votes
# Ruby, 122 bytes

```
->b{w=b.sum &:abs;(z=*-w..w).product(*[z]*k=b.size/2).any?{|c|x=-1;b==b.map{(0..x+=1).sum{|y|(a=c[x-y])&&c[y]?c[y]*a:0}}}}

```

Try it online!

Sheer brute force.

# Answer
> 2 votes
# Jelly, 26 bytes

```
Ærm2ÆṛJ}¹?×AṪ½Ɗ×þ`ŒdṙLƊṚ§⁼

```

A monadic Link that accepts a non-empty list of coefficients ordered by degree\* and yields `1` if the input represents the square of an integer polynomial, or `0` if not.

* i.e. the reverse order to the examples in the question, and \\$0\\$ being `[0]` rather than `[]`

**Try it online!**

### How?

Same approach as loopy walt's Python answer:

* get the polynomial with every other root of the polynomial defined by the input
* rescale using the square root of the coefficient of the highest degree in the input
* get the coefficients of the square of that polynomial
* check if that is equal to the input

Just a little fiddly when the input defines a constant since we find no roots whereupon finding coefficients would error (this is why `J}¹?` is there).

```
Ærm2ÆṛJ}¹?×AṪ½Ɗ×þ`ŒdṙLƊṚ§⁼ - Link: list of coefficients, C
Ær                         - get the roots of the polynomial P defined by C
  m2                       - modulo-2 slice
         ?                 - if...
        ¹                  - ...condition: no-op (falsey if P is just a constant)
    Æṛ                     - ...then: get polynomial coefficients
       }                   - ...else: use the right argument (implicitly C) with:
      J                    -            range of length (in this case that's [1])
                             (call this list NewCoefficients)
              Ɗ            - last three links as a monad - f(C):
           A               -   absolute values
            Ṫ              -   tail (absolute of coefficient of highest degree)
             ½             -   square root
          ×                - (NewCoefficients) multiplied by (that)
                 `         - use as both arguments of:
                þ          -   table of:
               ×           -     multiplication
                      Ɗ    - last three links as a monad - f(X=that)
                  Œd       -   anti-diagonals of X (starts with main, sweeps South-East)
                                 (has the effect of grouping product coefficients by degree)
                     L     -   length of X
                    ṙ      -   rotate (the anti-diagonals) left by (length of X)
                                 (ensures that the groups are sorted by degree descending)
                       Ṛ   - reverse
                        §  - sums (sum each group of like-degree product coefficients)
                         ⁼ - equals C?

```

# Answer
> 2 votes
# PARI/GP, 8 bytes

```
issquare
```

Attempt This Online!

A boring built-in. Return `1` for truthy and `0` for falsy.

---

# PARI/GP, 68 bytes

```
p->p&&for(i=!r=(abs(Vec(p)[1])^.5\/1)*x^d=#p\2,d,r+=(q=p-r^2)\r/2)+q
```

Attempt This Online!

Return `0` for truthy and `1` for falsy.

Basically a port of @tsh's JavaScript answer. Interestingly, the loop body `r+=(p-r^2)\r/2` looks very similar to the formula of calculating the square root of a real number using Newton's method.

# Answer
> 1 votes
# Mathematica, 144 bytes

Based on Factorization of polynomials.

The critical function used here is `FactorList[expr, GaussianIntegers -> False]`

---

Golfed version(It can be golfed much more). Try it online!

```
f[x_]:=Module[{a},If[IntegerQ@x,Return[!SquareFreeQ@x]];a=FactorList[x,GaussianIntegers->0];If[a[[1,1]]<2,a=a~Drop~1];AllTrue[a,EvenQ[#[[2]]]&]]

```

Ungolfed version. Try it online!

```
Clear["Global`*"];
f[expr_] := Module[{fl},
  If[IntegerQ@expr, Return[Not@SquareFreeQ[expr]]];
  fl = FactorList[expr, GaussianIntegers -> False];
  If[fl[[1, 1]] == 1, fl = Drop[fl, 1]];
  AllTrue[fl, EvenQ[#[[2]]] &]]
truthy = FromDigits[#, x] & /@ {{}, {25}, {1, 2, 1}, {1, 2, 1, 0,
     0}, {1, 4, 0, -8, 4}, {4, 28, 37, -42, 9}, {4, 0, -4, 4, 1, -2,
     1}, {1, -12, 60, -160, 240, -192, 64}};
falsy = FromDigits[#, x] & /@ {{-1}, {24, 1}, {1, 111}, {2, 4, 2}, {1,
      2, 1, 0}, {1, 3, 3, 1}, {1, -4, 5, -2, 0}, {4, 0, -4, 4, 1, 2,
     1}, {1, -9, 30, -45, 30, -9, 1}};
Print["Truthy tests:", f /@ truthy]
Print["Falsy tests:", f /@ falsy]

```

---
Tags: code-golf, math, decision-problem, polynomials
---


title: Cyclotomic polynomial

text:
# Question
Title: Cyclotomic polynomial
## Background (skip to definitions)

Euler proved a beautiful theorem about the complex numbers: \\$e^{ix} = \cos(x) + i \sin(x)\\$.

This makes de Moivre's theorem easy to prove:

$$ (e^{ix})^n = e^{i(nx)} \\ (\cos(x) + i\sin(x))^n = \cos(nx) + i\sin(nx) $$

We can plot complex numbers using the two-dimensional Euclidean plane, with the horizontal axis representing the real part and the vertical axis representing the imaginary part. This way, \\$(3,4)\\$ would correspond to the complex number \\$3+4i\\$.

If you are familiar with polar coordinates, \\$(3,4)\\$ would be \\$(5,\arctan(\frac 4 3))\\$ in polar coordinates. The first number, \\$r\\$, is the distance of the point from the origin; the second number, \\$\theta\\$, is the angle measured from the positive \\$x\\$-axis to the point, counter-clockwise. As a result, \\$3 = r \cos \theta\\$ and \\$4 = r \sin \theta\\$. Therefore, we can write \\$3+4i\\$ as \\$r \cos \theta + r i \sin \theta = r(\cos \theta + i \sin \theta) = re^{i\theta}\\$.

Let us solve the complex equation \\$z^n = 1\\$, where \\$n\\$ is a positive integer.

We let \\$z = re^{i\theta}\\$. Then, \\$z^n = r^ne^{in\theta}\\$. The distance of \\$z^n\\$ from the origin is \\$r^n\\$, and the angle is \\$n\theta\\$. However, we know that the distance of 1 from the origin is 1, and the angle is \\$\theta\\$. Therefore, \\$r^n=1\\$ and \\$n\theta=\theta\\$. However, if you rotate by \\$2π\\$ more, you still end up at the same point, because \\$2π\\$ is just a full circle. Therefore, \\$r=1\\$ and \\$n\theta = 2kπ\\$, giving us \\$z=e^{2ikπ / n}\\$.

We restate our discovery: **the solutions to \\$z^n=1\\$ are \\$z=e^{2ikπ / n}\\$**.

A polynomial can be expressed in terms of its roots. For example, the roots of \\$x^2-3x+2\\$ are 1 and 2, so \\$x^{2}-3x+2 = (x-1)(x-2)\\$. Similarly, from our discovery above:

$$z^n-1 = \prod^{n-1}\_{k=0} (z-e^{2ik\pi / n})$$

However, that product certainly contained roots of other n. For example, take \\$n=8\\$. The roots of \\$z^{4}=1\\$ would also be included inside the roots of \\$z^{8}=1\\$, since \\$z^{4}=1\\$ implies \\$z^{8} = (z^{4})^{2} = 1^{2} = 1\\$. Take \\$n=6\\$ as an example. If \\$z^{2}=1\\$, then we would also have \\$z^{6}=1\\$. Likewise, if \\$z^{3}=1\\$, then \\$z^{6}=1\\$.

If we want to extract the roots unique to \\$z^{n}=1\\$, we would need \\$k\\$ and \\$n\\$ to share no common divisor except \\$1\\$. Or else, if they share a common divisor \\$d\\$ where \\$d\>1\\$, then \\$z\\$ would be the \\$\frac k d\\$-th root of \\$z^{n / d}=1\\$. Using the technique above to write the polynomial in terms of its roots, we obtain the polynomial:

$$\prod\_{\substack{0 \le k \< n \\ \gcd(k,n) = 1}} (z - e^{2ik\pi / n})$$

Note that this polynomial is done by removing the roots of \\$z^{n / d}=1\\$ with d being a divisor of \\$n\\$. We claim that the polynomial above has integer coefficients. Consider the LCM of the polynomials in the form of \\$z^{n / d}-1\\$ where \\$d\>1\\$ and \\$d\\$ divides \\$n\\$. The roots of the LCM are exactly the roots we wish to remove. Since each component has integer coefficients, the LCM also has integer coefficients. Since the LCM divides \\$z^{n}-1\\$, the quotient must be a polynomial with integer coefficient, and the quotient is the polynomial above.

The roots of \\$z^{n}=1\\$ all have radius 1, so they form a circle. The polynomial represents the points of the circle unique to n, so in a sense the polynomials form a partition of the circle. Therefore, the polynomial above is the n-th cyclotomic polynomial. (cyclo- = circle; tom- = to cut)

## Definition 1

The n-th cyclotomic polynomial, denoted \\$\Phi\_n(x)\\$, is the unique polynomial with integer coefficients that divide \\$x^{n}-1\\$ but not \\$x^{k}-1\\$ for \\$k \< n\\$.

## Definition 2

The cyclotomic polynomials are a set of polynomials, one for each positive integer, such that:

$$x^n - 1 = \prod\_{k \mid n} \Phi\_k(x)$$

where \\$k \mid n\\$ means \\$k\\$ divides \\$n\\$.

## Definition 3

The \\$n\\$-th cyclotomic polynomial is the polynomial \\$x^{n}-1\\$ divided by the LCM of the polynomials in the form \\$x^{k}-1\\$ where \\$k\\$ divides \\$n\\$ and \\$k \< n\\$.

## Examples

1. \\$Φ\_{1}(x) = x - 1\\$
2. \\$Φ\_{2}(x) = x + 1\\$
3. \\$Φ\_{3}(x) = x^{2} + x + 1\\$
4. \\$Φ\_{30}(x) = x^{8} + x^{7} - x^{5} - x^{4} - x^{3} + x + 1\\$
5. \\$Φ\_{105}(x) = x^{48} + x^{47} + x^{46} - x^{43} - x^{42} - 2x^{41} - x^{40} - x^{39} + x^{36} + x^{35} + x^{34} + x^{33} + x^{32} + x^{31} - x^{28} - x^{26} - x^{24} - x^{22} - x^{20} + x^{17} + x^{16} + x^{15} + x^{14} + x^{13} + x^{12} - x^{9} - x^{8} - 2x^{7} - x^{6} - x^{5} + x^{2} + x + 1\\$

## Task

Given a positive integer \\$n\\$, return the \\$n\\$-th cyclotomic polynomial as defined above, in a reasonable format (i.e. e.g. list of coefficients is allowed).

## Rules

You may return floating point/complex numbers as long as they round to the correct value.

## Scoring

This is . Shortest answer in bytes wins.

## References

# Answer
> 12 votes
# Haskell, 120 bytes

```
import Data.Complex
p%a=zipWith(\x y->x-a*y)(p++[0])$0:p
f n=foldl(%)[1][cis(2*pi/fromInteger n)^^k|k<-[1..n],gcd n k<2]

```

Try it online!

Gives a list of complex floats that has entries like `1.0000000000000078 :+ 3.314015728506092e-14` due to float inacurracy. A direct method of multiplying out to recover the polynomial from its roots.

The `fromInteger` is a big concession to Haskell's type system. There's got to be a better way. Suggestions are welcome. Dealing with roots of unity symbolically might also work.

---

# Haskell, 127 bytes

```
(h:t)%q|all(==0)t=[]|1>0=h:zipWith(\x y->x-h*y)t q%q
f n=foldl(%)(1:(0<$[2..n])++[-1])[tail$f k++[0,0..]|k<-[1..n-1],mod n k<1]

```

Try it online!

No imports.

Uses the formula

Computes Φ\_n(x) by dividing the LHS by each of the other terms in the RHS.

The operator `%` does division on polynomials, relying on the remainder being zero. The divisor is assumed to be monic, and is given without the leading 1, and also with infinite trailing zeroes to avoid truncating when doing `zipWith`.

# Answer
> 7 votes
# Mathematica, ~~43~~ 41 bytes

```
Factor[x^#-1]/Times@@#0/@Most@Divisors@#&

```

Of course, we can always use the built-in, but if we don't, this divides *x**n*-1 by Φ*k*(*x*) (computed recursively) for every proper divisor *k* of *n*.

We use `Factor` to get a polynomial at the end. I think the reason this works is that `x^#-1` factors into all the cyclotomic polynomials of divisors of *n*, and then we divide out the ones we don't want.

*-2 bytes thanks to Jenny\_mathy, rewriting the `Factor` to only apply to the numerator.*

# Answer
> 5 votes
# MATL, ~~32~~ ~~31~~ ~~27~~ 25 bytes

```
Z\"@:@/]XHxvHX~18L*Ze1$Yn

```

Output may be non-integer due to floating-point inaccuracies (which is allowed). The footer does rounding for clarity.

**Try it online!**

# Answer
> 4 votes
# Haskell, ~~250 236 233 218~~ 216 bytes

This is a verbose version, (@xnor can do it in almost half the score) but it is guaranteed to work for any `n` as long as you have enough memory, but it does not use a builtin for generating the n-th cyclotomic polynomial. The input is an arbitrary size integer and the output is a polynomial type with (exact) rational coefficients.

The rough idea here is calculatin the polynomials recursively. For `n=1` or `n` prime it is trivial. For all other numbers this approach basically uses the formula from definition 2

solved for . Thanks @H.PWiz for quite a bunch of bytes!

```
import Math.Polynomial
import Data.Ratio
import NumberTheory
p=powPoly x
q=poly LE
c n|n<2=q[-1,1%1]|isPrime n=sumPolys$p<$>[0..n-1]|1>0=fst$quotRemPoly(addPoly(p n)$q[-1])$foldl1 multPoly[c d|d<-[1..n-1],n`mod`d<1]

```

For `n=105` this yields following polynomial (I tidied up all the `%1` denominators):

```
[1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1,1,1,1,0,0,-1,0,-1,0,-1,0,-1,0,-1,0,0,1,1,1,1,1,1,0,0,-1,-1,-2,-1,-1,0,0,1,1,1]

```

The polynomial for `n=15015` can be found here (the largest coefficient is 23).

Try it online!

# Answer
> 3 votes
# Jelly, 23 bytes

```
R÷
ÆḌÇ€FQœ-@Ç×ı2×ØPÆeÆṛ

```

Try it online!

Outputs as a list of coefficients.

Has floating point AND complex inaccuracies. Footer does rounding to make output prettier.

# Answer
> 3 votes
# , 36 bytes

```
1+//.@(*/)/@(,.~-)_1^2*%*i.#~1=i.+.]

```

Try it online!

Uses the formula

There are some floating-point inaccuracies, but it is allowed.

# Answer
> 3 votes
# Pari/GP, 8 bytes

A built-in.

```
polcyclo

```

Try it online!

---

# Pari/GP, 39 bytes, without built-in

```
f(n)=p=x^n-1;fordiv(n,d,d<n&&p/=f(d));p

```

Using the formula:

\\$ \Phi\_n(x)=\dfrac{x^{n}-1}{\prod\limits\_{\stackrel{d|n}{{}\_{d\<n}}}\Phi\_{d}(x)} \\$.

Try it online!

# Answer
> 3 votes
# APL(Dyalog Unicode), 32 bytes SBCS

```
{(0,⍵)-⍺×⍵,0}/⍤⌽1,¯1*2×÷×1⍸⍤=⊢∨⍳

```

Try it on APLgolf!

Almost direct port of miles' J answer. Returns a polynomial as a vector of coefficients wrapped once, with higher-order terms on the right side.

```
{(0,⍵)-⍺×⍵,0}/⍤⌽1,¯1*2×÷×1⍸⍤=⊢∨⍳    ⍵←n
                         1⍸⍤=⊢∨⍳    Numbers coprime to ⍵
                  ¯1*2×÷×    Divide by ⍵, double, raise -1 to the power of them
{(0,⍵)-⍺×⍵,0}/⍤⌽1,    Evaluate the polynomial with the roots of the above

```

# Answer
> 2 votes
# Mathematica, 81 bytes

```
Round@CoefficientList[Times@@(x-E^(2Pi*I#/k)&/@Select[Range[k=#],#~GCD~k<2&]),x]&

```

# Answer
> 2 votes
# , ~~176~~ ~~171~~ ~~139~~ 112 bytes

```
function(n){for(r in exp(2i*pi*(x=1:n)[(g=function(x,y)ifelse(o<-x%%y,g(y,o),y))(x,n)<2]/n))T=c(0,T)-r*c(T,0)
T}

```

Try it online!

Massively simpler version; uses a `for` loop rather than a `Reduce`.

# Answer
> 1 votes
## CJam (52 51 bytes)

```
{M{:K,:!W+K,0-{K\%!},{j($,\:Q,-[{(\1$Qf*.-}*;]}/}j}

```

Online demo. This is an anonymous block (function) which takes an integer on the stack and leaves a big-endian array of coefficients on the stack.

### Dissection

```
{                    e# Define a block
  M{                 e#   Memoised recursion with no base cases.
    :K,:!W+          e#     Store argument in K and build (x^K - 1)
    K,0-{K\%!},      e#     Find proper divisors of K
    {                e#     Foreach proper divisor D...
      j              e#       Recursive call to get Dth cyclotomic poly
      ($,\:Q,-       e#       The cleverest bit. We know that it is monic, and the
                     e#       poly division is simpler without that leading 1, so
                     e#       pop it off and use it for a stack-based lookup in
                     e#       calculating the number of terms in the quotient.
                     e#       Ungolfed this was (;:Q1$,\,-
                     e#       Store the headless divisor in Q.
      [              e#       Gather terms into an array...
        {            e#         Repeat the calculated number of times...
          (\         e#           Pop leading term, which goes into the quotient.
          1$Qf*.-    e#           Multiply Q by that term and subtract from tail.
        }*;          e#         Discard the array of Q,( zeroes.
      ]
    }/
  }j
}

```

# Answer
> 1 votes
# Maxima, 91 bytes

Using the formula:

\\$ \Phi\_n(x)=\dfrac{x^{n}-1}{\prod\limits\_{\stackrel{d|n}{{}\_{d\<n}}}\Phi\_{d}(x)} \\$.

---

Golfed version. Try it online!

```
f(n):=block(local(p,d),p:x^n-1,for d:1 thru n-1 do(if mod(n,d)=0 then p:p/f(d)),return(p))$

```

Ungolfed version. Try it online!

```
f(n) := block(
    local(p, d),
    p : x^n - 1,
    for d:1 thru n-1 do (
        if mod(n, d) = 0 then (
            p : p / f(d)
        )
    ),
    return(p)
)$

for n:1 thru 100 do (
    print(n, " -> ", ratsimp(f(n)))
)$

```

# Answer
> 0 votes
# JavaScript(ES6), ~~337~~ ~~333~~ ~~284~~ ... ~~252~~ ~~250~~ ~~245~~ 242 bytes

```
(v,z=[e=[1,u=0]],g=(x,y)=>y?g(y,x%y):x,h=Math,m=(l,x,p=h.cos(l),q=h.sin(l),i=0)=>x.map(()=>[(i&&(n=x[i-1])[0])-(w=x[i])[0]*p+w[1]*q,(i++&&n[1])-w[1]*p-w[0]*q]))=>{for(;++u<v;z=g(v,u)-1?z:[...m(h.PI*2*u/v,z),e]);return z.map(r=>h.round(r[0]))}

```

**Explanation (Selected):**

```
z=[e=[1,u=0]]

```

Initialize z=(1+0i) * x^0

```
g=(x,y)=>y?g(y,x%y):x

```

GCD calculation.

```
h=Math

```

Since I need to use Math functions quite a lot, I used another variable here.

```
m=(l,x,p=h.cos(l),q=h.sin(l),i=-1)=>blah blah blah

```

Polynomial multiplication.

```
for(;++u<v;z=g(v,u)-1?z:[...m(h.PI*2*u/v,z),e]);

```

The formula used is

```
return z.map(r=>h.round(r[0]))

```

Compress the output back to an integer array.

**Outputs:**

An array of integers, with the element at position i representing the coefficient of x^i.

One of the problem JS have is that since JS does not provide native library for calculations on polynomials and complex numbers, they needed to be implemented in an array-like way.

console.log(phi(105)) gives

```
Array(49)
 0:  1    1:  1    2:  1    3: -0    4: -0    5: -1    6: -1
 7: -2    8: -1    9: -1   10:  0   11: -0   12:  1   13:  1
14:  1   15:  1   16:  1   17:  1   18:  0   19: -0   20: -1
21:  0   22: -1   23: -0   24: -1   25:  0   26: -1   27: -0
28: -1   29:  0   30:  0   31:  1   32:  1   33:  1   34:  1
35:  1   36:  1   37: -0   38: -0   39: -1   40: -1   41: -2
42: -1   43: -1   44: -0   45: -0   46:  1   47:  1   48:  1
length: 49
__proto__: Array(0)

```

**337\>333(-4): Changed the code for checking undefined value**

**333\>284(-49): Changed the polynomial multiplication function because it can be simplified**

**284\>277(-7): Deleted some redundant code**

**277\>265(-12): Use 2 variables instead of a 2-element array to drop some bytes in array usage**

**265\>264(-1): Use Array.push() instead of Array.concat() to reduce 4 bytes, but added 3 for the for-loop braces and the z variable**

**264\>263(-1): Further golfed on the last amendment**

**263\>262(-1): Golfed on the for loop**

**262\>260(-2): Golfed out the if clause**

**260\>258(-2): Further combined the declarations**

**258\>252(-6): Golfed on reuse of array references**

**252\>250(-2): Replace some unary operators as binary operators**

**250\>245(-5): Move the increment in Array.map() to the last reference of counter to remove bytes**

**245\>242(-3): Use spread syntax instead of Array.push()**

---
Tags: code-golf, number-theory, polynomials
---


title: Self Referential Polynomials

text:
# Question
Title: Self Referential Polynomials
For every given degree \\$n\\$ it is possible to construct (at least one) an integral polynomial \\$p \in \mathbb Z\[X\]\\$ such that \\$p(k)\\$ (\\$p\\$ evaluated in \\$k\\$) is the coefficient of the term \\$x^k\\$ in the polynomial for all \\$0 \leqslant k \leqslant n\\$. To make them unique, we require the leading coefficient (the coefficient of \\$x^n\\$) to be positive and minimal.

These polynomials have some interesting properties, you can find some references in the thread that inspired me to do this challenge. You can also find those polynomials in https://oeis.org/A103423

One of the a priori unexpected properties is how the roots behave depending on \\$n\\$:

source (by /u/zorngov and /u/EpicSauceSc2)

### Task

Given a nonnegative integer \\$n\\$ output the self referential integral polynomial of degree \\$n\\$ with minimal positive leading coefficient.

### Details

The output can be in any human readable form, as string `x^2-x-1`, or also as a list of coefficients `[1,-1,-1]`. (The order of the coefficients can also be the other way around, it just needs to be consistent.)

### First few outputs

```
n=0: 1
n=1: x
n=2: x^2-x-1
n=3: 10*x^3-29*x^2-6*x+19
n=4: 57*x^4-325*x^3+287*x^2+423*x-19
n=5: 12813*x^5-120862*x^4+291323*x^3+44088*x^2-355855*x-227362

```

# Answer
> 4 votes
# Mathematica, 55 bytes

```
NullSpace@Table[x^c-Boole[r==c]/.x->r,{r,0,#},{c,0,#}]&

```

Output is the list coefficients, beginning from the constant term. Example:

```
In[1084] := Do[Print[%1077[n] // StandardForm], {n, 0, 7}]

{{1}}

{{0,1}}

{{-1,-1,1}}

{{19,-6,-29,10}}

{{-19,423,287,-325,57}}

{{-227362,-355855,44088,291323,-120862,12813}}

{{145991969,64989065,-123338281,-85635661,79841909,-18146731,1286795}}

{{-5958511844199,3384370785404,8437850634901,489428412300,-4499161007143,1776194531596,-258931801371,13131073916}}

```

---

This simply finds the vector such that `(A - I)v = 0`, similar to the MAPLE code in OEIS. The `NullSpace` method seems to always pick the minimal positive number for the last element, which matches the task description.

The `x^c-…/.x->r` indirection is to prevent having `0^0 == Indeterminate`.

# Answer
> 2 votes
# Sage, 74 bytes

```
lambda n:kernel(matrix(n+1,[j^-i-(-i==j)for i in[-n..0]for j in[0..n]])).0

```

The `-i` and `[-n..0]` could be `i` and `[0..n]`, if not for the positive leading coefficient requirement.

Try it on Sage Cell

# Answer
> 2 votes
# Pari/GP, 57 bytes

```
n->(a=matkerint(matrix(n++,,i,j,i--^j--)-1))*sign(a[n,1])

```

Try it online!

Returns the coefficients as a column vector.

# Answer
> 1 votes
# Julia, 172 bytes

Similar to the Mathematica code in OEIS.

---

Try it online!

```
using LinearAlgebra
function f(n)
n=n+1
V=nullspace([ifelse(k == 1, 1, (j - 1)^(k - 1)) for j in 1:n, k in 1:n] - I)[:,1]
if(V[end]<0) V=-V;end
(V/maximum(V))[end:-1:1]
end

```

---

Here, these coefficients are normalized by dividing by the maximum value of the coefficients.

**How to adjust the code?**

```
using LinearAlgebra

function VM(n)
    return [ifelse(k == 1, 1, (j - 1)^(k - 1)) for j in 1:n, k in 1:n]
end

function f(n)
    V = nullspace(VM(n + 1) - I)
    V = V[:, 1] # Get the first column
    if V[end] < 0
        V = -V
    end
    return (V / maximum(V))[end:-1:1] # Normalize and reverse
end

for n in 0:10
    print(f(n), "\n")
end

```

---
Tags: code-golf, math, number, linear-algebra, polynomials
---


title: Polynomialception

text:
# Question
Title: Polynomialception
Given two polynomials `f,g` of arbitrary degree over the integers, your program/function should evaluate the first polynomial in the second polynomial. `f(g(x))` (a.k.a. the *composition* `(fog)(x)` of the two polynomials)

### Details

Builtins are allowed. You can assume any reasonable formatting as input/output, but the input and output format should match. E.g. formatting as a string

```
x^2+3x+5

```

or as as list of coefficients:

```
[1,3,5] or alternatively [5,3,1]

```

Furthermore the input polynomials can be assumed to be fully expanded, and the outputs are also expected to be fully expanded.

### Examples

```
A(x) = x^2 + 3x + 5, B(y) = y+1
A(B(y)) = (y+1)^2 + 3(y+1) + 5 = y^2 + 5y + 9

A(x) = x^6 + x^2 + 1, B(y) = y^2 - y
A(B(y))= y^12 - 6y^11 + 15y^10 - 20y^9 + 15y^8 - 6y^7 + y^6 + y^4 - 2 y^3 + y^2 + 1

A(x) = 24x^3 - 144x^2 + 288x - 192, B(y) = y + 2
A(B(y)) = 24y^3

A(x) = 3x^4 - 36x^3 + 138x^2 - 180x + 27, B(y) = 2y + 3
A(B(y)) = 48y^4 - 96y^2

```

# Answer
> 11 votes
# Haskell, ~~86~~ 72 bytes

```
u!c=foldr1((.u).zipWith(+).(++[0,0..])).map c
o g=(0:)!((<$>g).(*))!pure

```

Defines a function `o` such that `o g f` computes the composition f ∘ g. Polynomials are represented by a nonempty list of coefficients starting at the constant term.

## Demo

```
*Main> o [1,1] [5,3,1]
[9,5,1]
*Main> o [0,-1,1] [1,0,1,0,0,0,1]
[1,0,1,-2,1,0,1,-6,15,-20,15,-6,1]
*Main> o [2,1] [-192,288,-144,24]
[0,0,0,24]
*Main> o [3,2] [27,-180,138,-36,3]
[0,0,-96,0,48]

```

## How it works

No polynomial-related builtins or libraries. Observe the similar recurrences

f(x) = a + f₁(x)x ⇒ f(x)g(x) = a g(x) + f₁(x)g(x)x,
f(x) = a + f₁(x)x ⇒ f(g(x)) = a + f₁(g(x))g(x),

for polynomial multiplication and composition, respectively. They both take the form

f(x) = a + f₁(x)x ⇒ W(f)(x) = C(a)(x) + U(W(f₁))(x).

The operator `!` solves a recurrence of this form for W given U and C, using `zipWith(+).(++[0,0..])` for polynomial addition (assuming the second argument is longer—for our purposes, it always will be). Then,

`(0:)` multiplies a polynomial argument by x (by prepending a zero coefficient);
`(<$>g).(*)` multiplies a scalar argument by the polynomial `g`;
`(0:)!((<$>g).(*))` multiplies a polynomial argument by the polynomial `g`;
`pure` lifts a scalar argument to a constant polynomial (singleton list);
`(0:)!((<$>g).(*))!pure` composes a polynomial argument with the polynomial `g`.

# Answer
> 10 votes
# Mathematica, 17 bytes

```
Expand[#/.x->#2]&

```

Example usage:

```
In[17]:= Expand[#/.x->#2]& [27 - 180x + 138x^2 - 36x^3 + 3x^4, 3 + 2x]

              2       4
Out[17]= -96 x  + 48 x

```

# Answer
> 7 votes
# TI-Basic 68k, 12 bytes

```
a|x=b→f(a,b)

```

The usage is straightforward, e.g. for the first example:

```
f(x^2+3x+5,y+1)

```

Which returns

```
y^2+5y+9

```

# Answer
> 6 votes
# Python 2, 138 ~~156 162~~ bytes

The inputs are expected to be integer lists with smallest powers first.

```
def c(a,b):
 g=lambda p,q:q>[]and q[0]+p*g(p,q[1:]);B=99**len(`a+b`);s=g(g(B,b),a);o=[]
 while s:o+=(s+B/2)%B-B/2,;s=(s-o[-1])/B
 return o

```

Ungolfed:

```
def c(a,b):
 B=sum(map(abs,a+b))**len(a+b)**2
 w=sum(B**i*x for i,x in enumerate(b))
 s=sum(w**i*x for i,x in enumerate(a))
 o=[]
 while s:o+=min(s%B,s%B-B,key=abs),; s=(s-o[-1])/B
 return o

```

In this computation, the polynomial coefficients are seen as digits (which may be negative) of a number in a very large base. After polynomials are in this format, multiplication or addition is a single integer operation. As long as the base is sufficiently large, there won't be any carries that spill over into neighboring digits.

-18 from improving bound on `B` as suggested by @xnor.

# Answer
> 4 votes
# Python + SymPy, ~~59~~ 35 bytes

```
from sympy import*
var('x')
compose

```

*Thanks to @asmeurer for golfing off 24 bytes!*

### Test run

```
>>> from sympy import*
>>> var('x')
x
>>> f = compose
>>> f(x**2 + 3*x + 5, x + 1)
x**2 + 5*x + 9

```

# Answer
> 3 votes
## Sage, 24 bytes

```
lambda A,B:A(B).expand()

```

As of Sage 6.9 (the version that runs on http://sagecell.sagemath.org), function calls without explicit argument assignment (`f(2) rather than f(x=2)`) causes an annoying and unhelpful message to be printed to STDERR. Because STDERR can be ignored by default in code golf, this is still valid.

This is very similar to Dennis's SymPy answer because Sage is a) built on Python, and b) uses Maxima, a computer algebra system very similar to SymPy in many ways. However, Sage is much more powerful than Python with SymPy, and thus is a different enough language that it merits its own answer.

Verify all test cases online

# Answer
> 2 votes
# PARI/GP, 19 bytes

```
(a,b)->subst(a,x,b)

```

which lets you do

```
%(x^2+1,x^2+x-1)

```

to get

> %2 = x^4 + 2\*x^3 - x^2 - 2\*x + 2

# Answer
> 1 votes
# MATLAB with Symbolic Toolbox, 28 bytes

```
@(f,g)collect(subs(f,'x',g))

```

This is an anonymous function. To call it assign it to a variable or use `ans`. Inputs are strings with the format (spaces are optional)

```
x^2 + 3*x + 5

```

Example run:

```
>> @(f,g)collect(subs(f,'x',g))
ans =
    @(f,g)collect(subs(f,'x',g))
>> ans('3*x^4 - 36*x^3 + 138*x^2 - 180*x + 27','2*x + 3')
ans =
48*x^4 - 96*x^2

```

# Answer
> 1 votes
# Pyth, ~~51~~ 34 bytes

```
AQsM.t*LVG.usM.t.e+*]0k*LbHN0tG]1Z

```

Test suite.

# Answer
> 1 votes
# Python 2, 239 232 223 bytes

```
r=range
e=reduce
a=lambda*l:map(lambda x,y:(x or 0)+(y or 0),*l)
m=lambda p,q:[sum((p+k*[0])[i]*(q+k*[0])[k-i]for i in r(k+1))for k in r(len(p+q)-1)]
o=lambda f,g:e(a,[e(m,[[c]]+[g]*k)for k,c in enumerate(f)])

```

A 'proper' implementation that does not abuse bases. Least significant coefficient first.

`a` is polynomial addition, `m` is polynomial multiplication, and `o` is composition.

# Answer
> 1 votes
## JavaScript (ES6), ~~150~~ 103 bytes

```
(f,g)=>f.map(n=>r=p.map((m,i)=>(g.map((n,j)=>p[j+=i]=m*n+(p[j]||0)),m*n+(r[i]||0)),p=[]),r=[],p=[1])&&r

```

Accepts and returns polynomials as an array a = \[a0, a1, a2, ...\] that represents a0 \+ a1\*x + a2\*x2 ...

Edit: Saved 47 bytes by switching from recursive to iterative polynomial multiplication, which then allowed me to merge two `map` calls.

Explanation: *r* is the result, which starts at zero, represented by an empty array, and *p* is *g**h*, which starts at one. *p* is multiplied by each *f**h* in turn, and the result accumulated in *r*. *p* is also multiplied by *g* at the same time.

```
(f,g)=>f.map(n=>            Loop through each term of f (n = f[h])
 r=p.map((m,i)=>(           Loop through each term of p (m = p[i])
  g.map((n,j)=>             Loop though each term of g (n = g[j])
   p[j+=i]=m*n+(p[j]||0)),  Accumulate p*g in p
  m*n+(r[i]||0)),           Meanwhile add p[i]*f[h] to r[i]
  p=[]),                    Reset p to 0 each loop to calculate p*g
 r=[],                      Initialise r to 0
 p=[1]                      Initialise p to 1
)&&r                        Return the result

```

# Answer
> 1 votes
# Ruby 2.4 \+ polynomial, 41+12 = 53 bytes

Uses flag `-rpolynomial`. Input is two `Polynomial` objects.

If someone outgolfs me in vanilla Ruby (no polynomial external library), I will be very impressed.

```
->a,b{i=-1;a.coefs.map{|c|c*b**i+=1}.sum}

```

# Answer
> 1 votes
# Maxima, 34 bytes

Try it online!

```
lambda([a,b],expand(subst(b,x,a)))

```

---
Tags: code-golf, math, sequence, polynomials, integer
---


title: The Chroma Key to Success

text:
# Question
Title: The Chroma Key to Success
The RGB color value `#00FF00` is a rather important one: it is used to make movies, TV shows, weather announcements, and more. It is the famous "TV green" or "green screen" color.

## The Challenge

Your task is to write a program that takes two input images, both in PNG format (or in your image library's image object type) and of the same dimensions. One image can be any old image. The other one is the image that will have a background of the color `#00FF00`. The output image will consist of the second image overlaid over the first, with no `#00FF00` color present (except in the first image). Input and output may be done with files, a GUI, etc. You are allowed to take an array of RGB values as input, as seen here. You may assume that an image has only pixels of full opacity.

## Basically...

Make a program that takes every `#00FF00` pixel in one image and replace it with the corresponding pixel in the background image.

## Test Cases

Generously provided by @dzaima: Background:
Foreground:
Output:

---

Of course, standard loopholes are *strictly forbidden*. This includes using an online resource to do it for you.
This is , so may the shortest code win and the best programmer prosper...

# Answer
> 14 votes
# x86-64 (and x86-32) machine code, ~~13~~ ~~15~~ 13 bytes

changelog:

1. Bugfix: the first version was only checking G=0xff, not requiring R and B to be 0. I changed to modifying the background in place so I could use `lodsd` on the foreground to have fg pixels in `eax` for the short-form `cmp eax, imm32` encoding (5 bytes), instead of `cmp dh,0xff` (3 bytes).
2. Save 2 bytes: noticed that modifying the bg in place allowed using a memory operand for `cmov`, saving a 2-byte `mov` load (and saving a register, in case that matters).

---

This is a function following the x86-64 System V calling convention, callable directly from C or C++ (on x86-64 non-Windows systems) with this signature:

```
void chromakey_blend_RGB32(uint32_t *background /*rdi*/,
                     const uint32_t *foreground /*rsi*/,
                  int dummy, size_t pixel_count /*rcx*/);

```

**The image format is RGB0 32bpp, with the green component at the 2nd lowest memory address within each pixel.** The ~~foreground~~ background image is modified in-place. `pixel_count` is rows\*columns. It doesn't care about rows/columns; it just chromekey blends however many dwords of memory you specify.

RGBA (with A required to be 0xFF) would require using a different constant, but no change in function size. Foreground DWORDs are compared for exact equality against an arbitrary 32-bit constant stored in 4 bytes, so any pixel-order or chroma-key colour can be easily supported.

The same machine code also works in 32-bit mode. To assemble as 32-bit, change `rdi` to `edi` in the source. All other registers that become 64-bit are implicit (lodsd/stosd, and loop), and the other explicit regs stay 32-bit. But note that you'll need a wrapper to call from 32-bit C, because none of the standard x86-32 calling conventions use the same regs as x86-64 SysV.

NASM listing (machine-code + source), commented for asm beginners with descriptions of what the more complex instructions do. (Duplicating the instruction reference manual is bad style in normal usage.)

```
 1                       ;; inputs:
 2                       ;; Background image pointed to by RDI, RGB0 format  (32bpp)
 3                       ;; Foreground image pointed to by RSI, RGBA or RGBx (32bpp)
 4          machine      ;; Pixel count in RCX
 5          code         global chromakey_blend_RGB32
 6          bytes        chromakey_blend_RGB32:
 7 address               .loop:                      ;do {
 8 00000000 AD               lodsd                   ; eax=[rsi], esi+=4. load fg++
 9 00000001 3D00FF0000       cmp    eax, 0x0000ff00  ; check for chromakey
10 00000006 0F4407           cmove  eax, [rdi]       ; eax = (fg==key) ? bg : fg
11 00000009 AB               stosd                   ; [rdi]=eax, edi+=4. store into bg++
12 0000000A E2F4             loop .loop              ;} while(--rcx)
13
14 0000000C C3               ret

##  next byte starts at 0x0D, function length is 0xD = 13 bytes

```

To get the original NASM source out of this listing, strip the leading 26 characters of each line with `<chromakey.lst cut -b 26- > chromakey.asm`. I generated this with
`nasm -felf64 chromakey-blend.asm -l /dev/stdout | cut -b -28,$((28+12))-` NASM listings leave more blank columns than I want between the machine-code and source. To build an object file you can link with C or C++, use `nasm -felf64 chromakey.asm`. (Or `yasm -felf64 chromakey.asm`). To strip the addresses and line numbers as well, `cut -b16-30,41-` works well.

**untested**, but I'm pretty confident that the basic idea of load / load / cmov / store is sound, because it's so simple.

I could save 3 bytes if I could require the caller to pass the chroma-key constant (0x00ff00) as an extra arg, instead of hard-coding the constant into the function. I don't think the usual rules allow writing a more generic function that has the caller set up constants for it. But if it did, the 3rd arg (currently `dummy`) is passed in `edx` in the x86-64 SysV ABI. Just change `cmp eax, 0x0000ff00` (5B) to `cmp eax, edx` (2B).

---

With SSE4 or AVX, you might do this faster (but larger code size) with `pcmpeqd` and `blendvps` to do a 32-bit element size variable-blend controlled by the compare mask. (With `pand`, you could ignore the high byte). For packed RGB24, you might use `pcmpeqb` and then 2x `pshufb`+`pand` to get TRUE in bytes where all 3 components of that pixel match, then `pblendvb`.

(I know this is code-golf, but I did consider trying MMX before going with scalar integer.)

# Answer
> 13 votes
# Mathematica ~~57~~ 35 bytes

update: by default, a green background is removed using `RemoveBackground`. The first submission included the unnecessary second parameter, \`{"Background",Green}".

---

```
#~ImageCompose~RemoveBackground@#2&

```

Removes the background of image 2 and composes the result with image 1.

---

**Example**

The following, in prefix rather than infix form, shows more clearly how the code works.

# Answer
> 11 votes
# Python 3 \+ numpy, 59 bytes

```
lambda f,b:copyto(f,b,'no',f==[0,255,0])
from numpy import*

```

Try it online!

Input is given in the format of a `numpy` array, with integer triplets representing pixels (where `#00FF00` in hex color code is equivalent to `[0, 255, 0]`). The input array is modified in place, which is allowed per meta.

## Example Images

### Input (from the question)

Background:

Foreground:

Foreground image after running the function:

## Reference Implementation (uses `opencv` to read image files)

```
g = lambda f,b:copyto(f,b,'no',f==[0,255,0])
from numpy import*

import cv2

f = cv2.imread("fg.png")
b = cv2.imread("bg.png")

g(f, b)

cv2.imshow("Output", f)
cv2.imwrite("out.png", f)

```

Displays the image to the screen and writes it to an output file.

# Answer
> 9 votes
# Processing, ~~116~~ 99 bytes

```
PImage f(PImage b,PImage f){int i=0;for(int c:f.pixels){if(c!=#00FF00)b.pixels[i]=c;i++;}return b;}

```

Unfortunately, processing doesn't support java 8 stuff, like lambdas.

Example implementation: (saves image as `out.png` and also draws it on screen)

```
PImage bg;
void settings() {
  bg = loadImage("bg.png");
  size(bg.width,bg.height);
}
void setup() {
  image(f(bg, loadImage("fg.png")), 0, 0);
  save("out.png");
}
PImage f(PImage b,PImage f){int i=0;for(int c:f.pixels){if(c!=#00FF00)b.pixels[i]=c;i++;}return b;}

```

# Answer
> 6 votes
## Bash + ImageMagick, 45 bytes

```
convert $1 $2 -transparent lime -composite x:

```

Takes two images as arguments and displays the output on the screen. Change `x:` to `$3` to write to a third file argument instead. The method is simple: read the "background" image; read the "foreground" imagek; reinterpret the color "lime" (#00ff00) as transparency in the second image; then composite the second image onto the first and output.

## ImageMagick: 28 bytes?

I could have submitted this as an ImageMagick answer but it's not clear how to deal with the arguments. If you want to posit that ImageMagick is a stack-based language (which is kinda sorta not really true but almost... it's weird) then `-transparent lime -composite` is a function that expects two images on the stack and leaves one merged image on the stack... maybe that's good enough to count?

# Answer
> 3 votes
# MATL, ~~40~~ ~~37~~ 31 bytes

```
,jYio255/]tFTF1&!-&3a*5M~b*+3YG

```

Example run with the offline interpreter. The images are input by their URL's (local filenames could also be provided).

### Explanation

```
,        % Do this twice
  j      %   Input string with URL or filename
  Yi     %   Read image as an M×N×3 uint8 array
  o      %  Convert to double
  255/   %   Divide by 255
]        % End
t        % Duplicate the second image
FTF      % Push 1×3 vector [0 1 0]
1&!      % Permute dimensions to give a 1×1×3 vector
-        % Subtract from the second image (M×N×3 array), with broadcast
&3a      % "Any" along 3rd dim. This gives a M×N mask that contains
         % 0 for pure green and 1 for other colours
*        % Mulltiply. This sets green pixels to zero
5M       % Push mask M×N again
~        % Negate
b        % Bubble up the first image
*        % Multiply. This sets non-green pixels to zero
+        % Add the two images
3YG      % Show image in a window

```

# Answer
> 3 votes
# Pyth, 27 bytes

```
M?q(Z255Z)GHG.wmgVhded,V'E'

```

It takes quoted input. The input are the two paths of the image files. Output a file `o.png` Unfortunately that cannot be tested on the online interpreter for safety reason (`'` is disabled on it). You will need to get Pyth on your computer to test it.

**Explanation**

```
M?q(Z255Z)GHG                  # Define a function g which takes two tuples G and H and returns G if G != (0, 255, 0), H otherwise
                       V'E'    # Read the images. They are returned as lists of lists of colour tuples
                      ,        # Zip both images
               m  hded         # For each couple of lists in the zipped list...
                gV             # Zip the lists using the function g
             .w                # Write the resulting image to o.png

```

# Answer
> 2 votes
# Matlab 2016b and Octave, ~~62~~ 59 bytes

Input: A = MxNx3 unit8 foreground matrix, B = MxNx3 unit8 background matrix.

```
k=sum(A(:,:,2)-A(:,:,[1 3]),3)==510.*ones(1,1,3);A(k)=B(k);

```

Output: A = MxNx3 unit8 matrix

Sample use:

```
A = imread('foreground.png');
B = imread('backgroundimg.png');

k=sum(A(:,:,2)-A(:,:,[1 3]),3)==510.*ones(1,1,3);A(k)=B(k);

imshow(A)

```

# Answer
> 1 votes
# C++, 339 bytes

This uses CImg, and it can take files in other formats, too. The result is displayed in a window.

```
#include<CImg.h>
using namespace cimg_library;
int main(int g,char** v){CImg<unsigned char> f(v[1]),b(v[2]);for(int c=0;c<f.width();c++){for(int r=0;r<f.height();r++){if((f(c,r)==0)&&(f(c,r,0,1)==255)&&(f(c,r,0,2)==0)){f(c,r)=b(c,r);f(c,r,0,1)=b(c,r,0,1);f(c,r,0,2) = b(c,r,0,2);}}}CImgDisplay dis(f);while(!dis.is_closed()){dis.wait();}}

```

Compile with `g++ chromakey.cpp -g -L/usr/lib/i386-linux-gnu -lX11 -o chromakey -pthread`.

# Answer
> 1 votes
# R, 135 bytes

```
function(x,y,r=png::readPNG){a=r(x);m=apply(a,1:2,function(x)all(x==0:1));for(i in 1:4)a[,,i][m]=r(y)[,,i][m];png::writePNG(a,"a.png")}

```

Anonymous function, takes 2 png file paths as arguments and output a png picture called `a.png`.

Slightly ungolfed, with explanations:

```
function(x,y){
    library(png)
    # readPNG output a 3D array corresponding to RGBA values on a [0,1] scale:
    a = readPNG(x)
    # Logical mask, telling which pixel is equal to c(0, 1, 0, 1),
    # i.e. #00FF00 with an alpha of 1:
    m = apply(a, 1:2, function(x) all(x==0:1))
    # For each RGB layer, replace that part with the equivalent part of 2nd png:
    for(i in 1:4) a[,,i][m] = readPNG(y)[,,i][m]
    writePNG(a,"a.png")
}

```

# Answer
> 1 votes
# SmileBASIC, 90 bytes whats the key

```
DEF C I,J
DIM T[LEN(I)]ARYOP.,T,I,16711936ARYOP 2,T,T,T
ARYOP 6,T,T,0,1ARYOP 5,I,I,J,T
END

```

`I` is the foreground and output, `J` is the background. Both are integer arrays of pixels, in 32 bit ARGB format.

## Ungolfed

```
DEF C IMAGE,BACKGROUND 'function
 DIM TEMP[LEN(IMAGE)]  'create array "temp"
 ARYOP #AOPADD,TEMP,IMAGE,-RGB(0,255,0)    'temp = image - RGB(0,255,0)
 ARYOP #AOPCLP,TEMP,TEMP,-1,1              'temp = clamp(temp, -1, 1)
 ARYOP #AOPMUL,TEMP,TEMP,TEMP              'temp = temp * temp
 ARYOP #AOPLIP,IMAGE,IMAGE,BACKGROUND,TEMP 'image = linear_interpolate(image, background, temp)
END

```

## Explanation:

ARYOP is a function that applies a simple operation to every element in an array.
It is called like `ARYOP mode, output_array, input_array_1, input_array_2, ...`

First, to determine which pixels in the image are green, `-16711936` (the RGBA representation of the green color) is subtracted from each pixel in the foreground image. This gives an array where `0` represents green pixels, and any other number represents non-green pixels.

To convert all nonzero values to `1`, they are squared (to remove negative numbers), then clamped to between `0` and `1`.

This results in an array with only `0`s and `1`s.
`0`s represent green pixels in the foreground image, and should be replaced with pixels from the background.
`1`s represent non-green pixels, and those will need to be replaced with pixels from the foreground.

This can easily be done using linear interpolation.

# Answer
> 0 votes
# PHP, 187 bytes

```
for($y=imagesy($a=($p=imagecreatefrompng)($argv[1]))-1,$b=$p($argv[2]);$x<imagesx($a)?:$y--+$x=0;$x++)($t=imagecolorat)($b,$x,$y)-65280?:imagesetpixel($b,$x,$y,$t($a,$x,$y));imagepng($b);

```

assumes 24bit PNG files; takes file names from command lines arguments, writes to stdout.
Run with `-r`.

**breakdown**

```
for($y=imagesy(                                 # 2. set $y to image height-1
        $a=($p=imagecreatefrompng)($argv[1])    # 1. import first image to $a
    )-1,
    $b=$p($argv[2]);                            # 3. import second image to $b
    $x<imagesx($a)?:                            # Loop 1: $x from 0 to width-1
        $y--+$x=0;                              # Loop 2: $y from height-1 to 0
        $x++)
            ($t=imagecolorat)($b,$x,$y)-65280?:     # if color in $b is #00ff00
                imagesetpixel($b,$x,$y,$t($a,$x,$y));   # then copy pixel from $a to $b
imagepng($b);                                   # 5. output

```

# Answer
> 0 votes
# JavaScript (ES6), 290 bytes

```
a=>b=>(c=document.createElement`canvas`,w=c.width=a.width,h=c.height=a.height,x=c.getContext`2d`,x.drawImage(a,0,0),d=x.getImageData(0,0,w,h),o=d.data,o.map((_,i)=>i%4?0:o[i+3]=o[i++]|o[i++]<255|o[i]?255:0),x.drawImage(b,0,0),createImageBitmap(d).then(m=>x.drawImage(m,0,0)||c.toDataURL()))

```

Takes input as two `Image` objects (in currying syntax), which can be created with an HTML `<image>` element. Returns a Promise that resolves to the Base64 data URL of the resulting image, which can be applied to the `src` of an `<image>`.

The idea here was to set the alpha value for each `#00FF00` pixel to `0` and then paint the foreground, with its background keyed out, on top of the background.

## Test Snippet

Including the foreground and background by their data URLs was too large to post here, so it was moved to CodePen:

Try it online!

# Answer
> 0 votes
# OSL, 83 bytes

```
shader a(color a=0,color b=0,output color c=0){if(a==color(0,1,0)){c=b;}else{c=a;}}

```

Takes two inputs. The first is the foreground, and the second, the background.

# Answer
> 0 votes
# C (gcc), ~~248~~ 245 bytes

*-3 bytes thanks to @ceilingcat*

```
#define F(N)for(i=0;i<N;++i)
#define V x,y,z,w
*S="P%d %d %d %d ";V;v[6];main(i,f)int**f;{F(2)fscanf(f[i]=fopen(f[i+1],"r"),S,&x,&y,&z,&w);for(printf(S,V);~*v;x=0){F(3)x||putchar(i[!(*v|v[1]-255|v[2])*3+v]);F(3)v[i+3]=getc(*f),v[i]=getc(f[1]);}}

```

Try it online!

---
Tags: code-golf, graphical-output, image-processing
---


# Question
Title: Tips for golfing in TypeScript's type system
What general tips do you have for golfing in TypeScript's type system? I'm looking for ideas which can be applied to code-golf problems and which are also at least somewhat specific to TypeScript's type system (e.g. "remove comments" is not an answer).

Please post one tip per answer.

If you're not familiar with programming in TypeScript's type system, the TypeScript Handbook and TypeScript Discord Server are both great resources, as well as the TypeScript Type System chatroom.

If you post any answers in TypeScript's type system, make sure to check out this deadlineless 50-100 rep bounty.

# Answer
# Embrace the type errors

Type parameter constraints can often be replaced by a well-placed `//@ts-ignore` comment:

```
// 34 bytes
type Head<T extends any[]> = T[0];

```

```
// 33 bytes
//@ts-ignore
type Head<T> = T[0];

```
> 9 votes

# Answer
# Use `{}` instead of `any` or `unknown` where applicable

The type `{}` counterintuitively means "anything that is not `null` or `undefined`, so it can be useful for shorter constraints, if you don't care about `null`/`undefined`:

```
// 34 bytes
type Head<T extends any[]> = T[0];

```

```
// 32 bytes
type Head<T extends{}[]> = T[0];

```
> 5 votes

# Answer
# Use indexed access types as switch statements

Instead of

```
type Foo<T extends "a" | "b" | "c"> =
    T extends "a"
    ? A
  : T extends "b"
    ? B
  : T extends "c"
    ? C
  : never

```

use

```
type Foo<T extends "a" | "b" | "c"> = {
  a: A,
  b: B,
  c: C,
}[T]

```

Note that in some cases this won't work for recursive types, as it might give the error `'foo' is referenced directly or indirectly in its own type annotation`.
> 3 votes

# Answer
# Use unions as sets

You can use `T | U | V` as a set, or `[T] | [U] | [V]` if you want to prevent merging (e.g. `number | 1` becomes `number`).

Set operations:

* Union: `A | B`
* Intersection: `A & B` for sets of primitives, or `Extract<A, B>` more more complex types
* Difference: `Exclude<A, B>`
* Map: `A extends A ? Foo<A> : 0` (if `A` is a type parameter)
* Map: `A extends infer B ? B extends A ? Foo<B> : 0 : 0` (if `A` is not a type parameter)
* Map: `A extends infer B ? [B]|[B] extends [A] ? Foo<B, A>` (if `A` is a type parameter and you need to preserve the reference to the full set)
* Empty set: `never`
* Size:

```
//@ts-ignore
type CountUnion<T,N=[]>=(T extends T?(x:()=>T)=>0:0)extends(x:infer U)=>0?U extends()=>infer V?CountUnion<Exclude<T,V>,[...N,0]>:N:0;

```
> 3 votes

# Answer
# Use `any` or shorter types for template literals

```
type EndsInPeriod<I> = I extends `${any}.` ? true : false;

```

Saving around 3 bytes (instead of the alternative `${string}.`
> 3 votes

# Answer
# Use type parameter defaults

Instead of

```
type RecursiveThing<T> = /* something recursive */
type Foo = RecursiveThing<Bar>;

```

Use

```
type Foo<T=Bar> = /* something recursive */

```

Instead of

```
type SomeConstant = /* something */;
type Foo<T> = /* something using that type */

```

Use

```
type Foo<T,SomeConstant = /*something*/>= /* something using that type */

```
> 2 votes

# Answer
# Use tuples to do non-negative integer arithmetic

```
type IntToTuple<T, A extends any[] = []> = T extends A["length"] ? A : IntToTuple<T, [...A, 0]>;
type TupleToInt<T extends any[]> = T["length"];
type AddTuples<A extends any[], B extends any[]> = [...A, ...B];
type SubTuples<A extends any[], B extends any[]> = A extends [...B, ...infer C] ? C : never /* overflow */;
type MulTuples<A extends any[], B extends any[], C extends any[] = []> = A extends [0, ...infer A] ? MulTuples<A, B, [...C, ...B]> : C
type Add<A, B> = TupleToInt<AddTuples<IntToTuple<A>, IntToTuple<B>>>;
type Sub<A, B> = TupleToInt<SubTuples<IntToTuple<A>, IntToTuple<B>>>;
type Mul<A, B> = TupleToInt<MulTuples<IntToTuple<A>, IntToTuple<B>>>;

type Answer = Mul<Add<5, 2>, Sub<9, 3>>;
//   ^? - 42

```
> 1 votes

# Answer
# Use stringification of number literal types to parse them

``${42}`` becomes `"42"`, which can be worked with more easily:

```
type Mod2<N extends number> = `${N}` extends `${string}${"1" | "3" | "5" | "7" | "9"}` ? 1 : 0;

```
> 1 votes

# Answer
# Use homomorphic mapped types to map tuples

```
type Foo<T> = [[[T]]];
type MapFoo<T> = {
  [K in keyof T]: Foo<T[K]>
}

type X = MapFoo<[1, 2, 3]>;
//   ^? - [Foo<1>, Foo<2>, Foo<3>]

```
> 1 votes

# Answer
# Use this large collection of golfed utility types
> 1 votes

# Answer
# A couple minor byte-saves

## Use `A extends{length:N}` instead of `A["length"]extends N`

This saves 1 byte off of decimal to unary conversion, among a few other things.

## When taking the first element and the rest of a tuple, use `{}` and indexing instead of `infer`

This saves 2 bytes if you only use the inferred value once. For example, this

```
A extends[infer B,...infer C]?[B,C]:0

```

becomes

```
A extends[{},...infer C]?[A[0],C]:0

```

---

I will add more if I think of / encounter any.
> 1 votes

# Answer
# Combine multiple `extends` clauses into one

If you have multiple `extends` clauses where you're just extracting part of a tuple/string/whatever and know that the match will succeed and don't care about the failure case, you can combine many clauses into one:

```
// Before:
type F<A,B>=A extends[infer D,...infer K]?B extends[5,6,infer F,...infer L]?[D,K,F,L]:0:0
// After:
type F<A,B>=[A,B]extends[[infer D,...infer K],[5,6,infer F,...infer L]]?[D,K,F,L]:0

```

In a case like this where you're just selecting part of each input, you can instead take the inputs as a tuple:

```
// Even better:
type F<A>=A extends[[infer D,...infer K],[5,6,infer F,...infer L]]?[D,K,F,L]:0

```

This example is purposely useless, but this can come in handy in real challenges. One use case is when zipping two tuples:

```
type F<A>=A extends[[infer X,...infer G],[infer Y,...infer H]]?[[X,Y],...F<A>]:A

```

I've used this technique in two submissions of mine: Add two really big numbers and Order by Earliest Lower Digit.
> 1 votes

---
Tags: code-golf, tips
---