A multiple of n in every base!

Question

In November 2019, Alon Ran published a particularly lovely sequence in the OEIS, A329126:

\$a(n)\$ is the lexicographically earliest string of digits which yields a multiple of \$n\$ when read in any numeric base. 1, 110, 101010, 111100, 100010001000100010, 1111110, 10000010000010000010000010000010000010, 11111111000, 1010101010101010100, ...

For example, \$a(3) = 101010\$ because this is the smallest number that is a multiple of \$3\$ in every integer base $$ \begin{alignat}{2} 101010_1 &= 1^5 + 1^3 + 1^1 = 3 &&= 3\cdot1 \\ 101010_2 &= 2^5 + 2^3 + 2^1 = 42 &&= 3\cdot14 \\ 101010_3 &= 3^5 + 3^3 + 3^1 = 273 &&= 3\cdot91 \\ 101010_4 &= 4^5 + 4^3 + 4^1 = 1092 &&= 3\cdot364 \\ 101010_5 &= 5^5 + 5^3 + 5^1 = 3255 &&= 3\cdot1085 \\ &\hspace{0.4em}\vdots \end{alignat} $$ Another way of conceptualizing this is that the polynomial \$f(n) = n^5 + n^3 + n\$ is divisible by \$3\$ for every integer \$n\$, which you can check manually with modular arithmetic.

---

We can generalize this idea further by looking at polynomials with restricted coefficients. For example, if we allow the coefficients of the polynomial to be \$\{-1,0,1\}\$, then \$g(n) = n^4-n^2\$ is divisible by \$4\$ for all integers \$n\$.

---

The challenge!

In any reasonable format of your choosing, you will be given a list of allowed integer coefficients, coefficients (which always contains both \$0\$ and \$1\$), and a positive integer, k.

Your job is to find the degree, d, of the smallest monic polynomial as measured by degree, $$ p(n) = n^d + c_{d-1}n^{d-1} + \dots + c_2n^2 + c_1n + c_0, $$ such that each \$c_i\$ is in coefficients and \$p(n)\$ is a multiple of \$k\$ for all integers \$n \in \mathbb{Z}\$.

In order to avoid the simplest form of brute-forcing, your program must be able to compute each individual test case from the test data on TIO (or, if TIO does not support your language, on my machine within 60 seconds). This is code-golf, so shortest code wins.

Test data

coefficients | k | d  | example
-------------+---+----+--------------------------------------------
[0,1]        | 1 | 0  | 1
[0,1]        | 2 | 2  | x^2  + x
[0,1]        | 3 | 5  | x^5  + x^3  + x
[0,1]        | 4 | 5  | x^5  + x^4  + x^3  + x^2
[0,1]        | 5 | 17 | x^17 + x^13 + x^9  + x^5  + x
[0,1]        | 6 | 6  | x^6  + x^5  + x^4  + x^3  + x^2  + x
[0,1]        | 7 | 37 | x^37 + x^31 + x^25 + x^19 + x^13 + x^7 + x
[-1,0,1]     | 4 | 4  | x^4  - x^2
[-1,0,1]     | 5 | 5  | x^5  - x
[-1,0,1]     | 6 | 3  | x^3  - x
[-1,0,1]     | 7 | 7  | x^7  - x
[-1,0,1]     | 8 | 5  | x^5  - x^3
[-1,0,1]     | 9 | 8  | x^8  - x^6  - x^4 + x^2  
[ 0,1,2]     | 4 | 4  | x^4  + x^2  + 2x
[ 0,1,3]     | 5 | 9  | x^9  + x^5  + 3x
[ 0,1,2,3]   | 6 | 3  | x^3  + 3x^2 + 2x

Answer 1

Python 3.8+, 385 372 bytes

from itertools import*
E=enumerate
def f(c,k):
 p=[range(k)]
 while(n:=[a*b%k for a,b in E(p[-1])])not in p:p+=n,;S=[{0}];M=[D:=0]*len(p)
 for d in range(k*k):
  if[]==p[D:]:D=p.index(n)
  if any(all(sum(s)%k<1for s in zip(*T))for T in product(*[[[P*(C+(i==D))for P in p[i]]for C in S[r]]for i,r in E(M)])):return(k!=1)+d
  S+=[{a+b for b in S[-1]for a in c}];M[D]+=1;D+=1

Try it online!

This code is surprisingly fast and runs in only about 1.5s on TIO.

-12 bytes thanks to @ovs

Python 3.8+ is used for the walrus operator (:=), but it costs three bytes because Python 3 doesn't allow mixed tabs and spaces.

How

A naive way to go about the challenge is to test all possible coefficient combinations at a given degree, then proceed to the next degree. This is viable for all of the test cases except for [0,1] | 7, which requires degree d=37. This would require brute-forcing more than 2^37 combinations, even if the order is chosen intelligently.

My next approach was to work out what possible sums could be obtained for each n then work backwards on sets of coefficients which work for those sums, but that ran into much of a similar problem.

Since only that one test case fails, it would be possible to optimize for coeffs=[0,1] specifically, but that seems too cheaty / against the spirit of the challenge.

The method used in this answer relies on the property that the powers of n (mod k) are periodic, at least after a certain point. As an example for k=4:

Table of n^i mod k for k=4
i=0 corresponds to the constant term

n\i 0 1 2 3 4 5 6 7 8 9
0   1 0 0 0 0 0 0 0 0 0
1   1 1 1 1 1 1 1 1 1 1
2   1 2 0 0 0 0 0 0 0 0
3   1 3 1 3 1 3 1 3 1 3

Thus the powers mod 4 overall have a period of length 2 starting from i=1. (Note that i=0 is never a solution unless k=1, and the constant term should always be divisible by k, so it should normally be zero).

We use this periodicity to our advantage. For example, in the [0,1] | 7 | 37 case, we don't need to consider 37 exponents. The powers have a period of length 6 starting from i=0, so there are only 6 distinct tuples of residues of of powers.

Table of n^i mod k for k=7
i=0 corresponds to the constant term

n\i 1 2 3 4 5 6  7 8 9
0   0 0 0 0 0 0  0 0 0
1   1 1 1 1 1 1  1 1 1
2   2 4 1 2 4 1  2 4 1
3   3 2 6 4 5 1  3 2 6
4   4 2 1 4 2 1  4 2 1
5   5 4 6 2 3 1  5 4 6
6   6 1 6 1 6 1  6 1 6

We only need to consider the coefficients with which to weight the six columns i∊[1...6]. Now we have to be careful. The coefficients are not just chosen from [0,1]. As we proceed to higher degrees, we can reuse the coefficients, so for example the n term is the same as the n^7 term is the same as the n^13 term; these all have powers in the 1 column. Thus we can use 7 copies of the 1 column to obtain the desired polynomial.

Ungolfed Code

import itertools


def min_degree(coeffs, k):
    if k == 1:
        return 0
    # powers[n-1] will be the tuple of [1...k]^m mod k
    # e.g. for k=4, powers=[(1,2,3),(1,0,1),(1,0,3)]; next power will be (1,0,1) which repeats
    powers = [tuple(range(1, k))]
    # keep computing more powers until two are the same,
    # (we know the sequence of powers must be periodic eventually)
    while (
        next_power := tuple([a * b % k for a, b in enumerate(powers[-1], 1)])
    ) not in powers:
        powers.append(next_power)
    # period_start_index is 0-indexed
    # for k=4, period_start_index=1 and period_length=2
    period_start_index = powers.index(next_power)
    period_length = len(powers) - period_start_index
    # summed_coeffs[i] is all possibilities of sums of i coeffs with replacement.
    # for example for i=2 and coeffs=[0,1,3], summed_coeffs[i]={0,1,2,3,4,6}={0+0,0+1,1+1,0+3,1+3,3+3}
    summed_coeffs = [{0}]
    for i in range(k * k):
        # mod by k here is unnecessary but can theoretically improve speed
        summed_coeffs.append({(a + b) % k for b in summed_coeffs[-1] for a in coeffs})
    max_repeats = [0 for _ in powers]
    # A simple upper bound for d is d ≤ k*φ(k)
    # because Σ_{m=1}^{k}(n^(mφ(k)) always works
    # Since φ(k) ≤ k-1, we know d ≤ k*φ(k) ≤ k*(k-1) < k*k
    for d in range(1, k * k):
        # degree_index is the index of [1...n]^(d-1) mod k in the list `powers`
        # we use d-1 because the `d`-th coefficient must be 1 ∵ monic polynomial
        degree_index = (
            d - 1
            if d - 1 < period_start_index
            else period_start_index + (d - 1 - period_start_index) % period_length
        )
        possible_column_totals = [
            [
                tuple(
                    [p * (coeff + (1 if i == degree_index else 0)) for p in powers[i]]
                )
                for coeff in summed_coeffs[r]
            ]
            for i, r in enumerate(max_repeats)
        ]
        max_repeats[degree_index] += 1
        for column_totals in itertools.product(*possible_column_totals):
            if all(sum(s) % k == 0 for s in zip(*column_totals)):
                return d

Answer 2

Haskell, 221 bytes

import Data.Set
k%g|w<-fromList.fmap(\q->(`mod`k).sum.zipWith(*)q.(<$>[1..]).(^)<$>[1..k]).mapM id=[n|n<-[0..],t<-[div n 2],n>(k-1)^2||size(intersection(w$(0-)<$>g<$[1..t])$w$([0]<$[1..t])++(g<$[1..n-t-1])++[[1]])>0]!!0

Use compiler flags -O2 (I didn't test if -O is enough. It might be)

To get the result, evaluate k%digits, e.g, 6%[0,1]. Here's a link to TIO with a main function that runs all the tests, except 7%[0,1]:

Try it online!

The code works by generating all the negations of the possible first halves of the polynomial, and all the possible second halves of the polynomial. Then it evaluates all of them mod k at all possible inputs, puts the results in two sets, and checks if they have any intersection.

It also has two spot optimizations: the constant coefficient is always zero, and a polynomial of degree k^2-2k+2 is always possible.

The code doesn't run to completion on the TIO site. however, it reached 30s on my computer. So I'm waiting for the OP to decide if this program is ok. Also, there is a slightly shorter but less performant version that also might make the cut:

import Data.Set
k%g|w<-fromList.fmap(\q->(`mod`k).sum.zipWith(*)q.(<$>[1..]).(^)<$>[1..k]).mapM id=[n|n<-[0..],t<-[div n 2],size(intersection(w$(0-)<$>g<$[1..t])$w$([0]<$[1..t])++(g<$[1..n-t-1])++[[1]])>0]!!0

that removes the k^2-2k+2 optimization.

Answer 3

JavaScript (ES6),  297 290 286  284 bytes

Expects (coefficients)(k).

c=>k=>(g=(i,q,d,A,m,h)=>i<n?c.map(c=>c*q&&i>=r|!m&&g(i,q-1,(D=m-~i)>d?D:d,A.map((v,j)=>v+B[i][j]*c),m+n-r,D>d?c:h),g(i+1,q,d,A,0,h)):M=h^!A.some(x=>x%k)|M<d?M:d)(0,k,0,(M=A=>(r=M[A+A])?A.fill(!r--):M(A.map((v,i)=>v*i%k),M[A+A]=n=B.push(A)))([...Array(k)].map((_,i)=>i),B=[]),0)&&M^k<2

Try it online!

Commented

I think this is very close to @fireflame241's solution, although I came up with this independently. The first part was working fine from the beginning. But for some reason, I got the second part poorly written and very slow on my early attempts.

Step 1

We store in \$B[\:]\$ all the distinct tuples \$(0,1,\dots,k-1)\$ raised to increasing powers, modulo \$k\$. At the end of the process, \$n\$ is the number of entries in \$B[\:]\$ and \$r\$ is the 0-based index at which the sequence restarts.

( M = A =>            // M is a recursive function taking a vector A[]
  (r = M[A + A]) ?    // if A[] was already encountered at index r:
    A.fill(!r--)      //   stop the recursion, clear all values in A[]
                      //   and decrement r
  :                   // else:
    M(                //   do a recursive call:
      A.map((v, i) => //     update each value in A[] ...
        v * i % k     //       ... to the next power modulo k
      ),              //     end of map()
      M[A + A] = n =  //     store A[] in M, set n = length of B[] ...
      B.push(A)       //     ... once A[] has been appended to B[]
    )                 //   end of recursive call
)(                    // initial call to M with:
  [...Array(k)]       //   an array filled with all values from 0 to k - 1
  .map((_, i) => i),  //
  B = []              //   initialize B[] to an empty array
)                     //

Step 2

We now use a recursive function \$g\$ to find a combination of tuples from \$B[\:]\$ multiplied by some coefficients from \$c[\:]\$ that sum up to \$0\$'s modulo \$k\$. We eventually return the lowest degree \$M\$ of the corresponding polynomials, ignoring those that are not monic.

By using any tuple that repeats, it's always possible to build a solution with exactly \$k\$ terms that lead to the same result modulo \$k\$. That's why we use \$k\$ as an upper bound for the number of terms in the polynomial.

The initial call to \$g\$ is g(0, k, 0, A, 0) and the returned value is M ^ k < 2 to support the edge case \$k=1\$.

g = (                        // g is a recursive function taking:
  i,                         //   i = pointer in B[]
  q,                         //   q = max. number of terms in the polynomial
  d,                         //   d = degree
  A,                         //   A[] = current vector, initialized to the
                             //         array of k 0's, returned by step 1
  m,                         //   m = number of times B[i] has been used
  h                          //   h = coefficient of the leading term
) =>                         //
  i < n ?                    // if i is less than n:
    c.map(c =>               //   for each coefficient c in c[]:
      c * q &&               //     if both c and q are not 0
      i >= r | !m &&         //     and B[i] can be repeated or m = 0,
      g(                     //     do a recursive call:
        i,                   //       leave i unchanged
        q - 1,               //       decrement q
        (D = m - ~i)         //       D is the degree of the new term
        > d ? D : d,         //       update d to D if D > d
        A.map((v, j) =>      //       update A[]:
          v + B[i][j] * c),  //         add B[i][j] * c to each A[i]
        m + n - r,           //       add the period to m
        D > d ? c : h        //       update the coeff. of the leading term
      ),                     //     end of recursive call
      g(                     //     do another recursive call with:
        i + 1, q, d, A, 0, h //       i + 1, m = 0
      )                      //       all other parameters unchanged
    )                        //   end of map()
  :                          // else (i = n):
    M = h ^                  //   if the leading coeff. is not equal to 1
        !A.some(x => x % k)  //   or some value in A[] is not 0 modulo k
        | M < d ? M          //   or M is less than d, leave M unchanged
                : d          //   else update it do d