The Complexity of Divisibility

We address two sets of long-standing open questions in probability theory, from a computational complexity perspective: divisibility of stochastic maps, and divisibility and decomposability of probability distributions. We prove that finite divisibility of stochastic maps is an NP-complete problem, and extend this result to nonnegative matrices, and completely-positive trace-preserving maps, i.e. the quantum analogue of stochastic maps. We further prove a complexity hierarchy for the divisibility and decomposability of probability distributions, showing that finite distribution divisibility is in P, but decomposability is NP-hard. For the former, we give an explicit polynomial-time algorithm. All results on distributions extend to weak-membership formulations, proving that the complexity of these problems is robust to perturbations.


Introduction and overview
People have pondered divisibility questions throughout most of Western science and philosophy. Perhaps the earliest written mention of divisibility is in Aristotle's Physics in 350 BC, in the form of the Arrow paradox-one of Zeno of Elea's paradoxes (ca. 490-430 BC). Aristotle's lengthy discussion of divisibility (he devotes an entire chapter to the topic) was motivated by the same basic question as more modern divisibility problems in mathematics: can the behaviour of an object-physical or mathematical-be subdivided into smaller parts?
For example, given a description of the evolution of a system over some time interval t, what can we say about its evolution over the time interval t/2? If the system is stochastic, this question finds a precise formulation in the divisibility problem for stochastic matrices [19]: given a stochastic matrix P, can we find a stochastic matrix Q such that P = Q 2 ?
This question has many applications. For example, in information theory stochastic matrices model noisy communication channels, and divisibility becomes important in relay coding, when signals must be transmitted between two parties where direct end-toend communication is not available [23]. Another direct use is in the analysis of chronic disease progression [3], where the transition matrix is based on sparse observations of patients, but finer-grained time-resolution is needed. In finance, changes in companies' credit ratings can be modelled using discrete time Markov chains, where rating agencies provide a transition matrix based on annual estimates-however, for valuation or risk analysis, a transition matrix for a much shorter time periods needs to be inferred [17].
We can also ask about the evolution of the system for all times up to time t, i.e. whether the system can be described by some continuous evolution. For stochastic matrices, this has a precise formulation in the embedding problem: given a stochastic matrix P, can we find a generator Q of a continuous-time Markov process such that P = exp(Qt)? The embedding problem seems to date back further still, and was already discussed by Elfving in 1937 [10]. Again, this problem occurs frequently in the field of systems analysis, and in analysis of experimental time-series snapshots [7,22,27].
Many generalizations of these divisibility problems have been studied in the mathematics and physics literature. For example, the question of square-roots of (entry-wise) nonnegative matrices is an old open problem in matrix analysis [24]: given an entry-wise nonnegative matrix M, does it have an entry-wise nonnegative square-root? In quantum mechanics, the analogue of a stochastic matrix is a completely-positive trace preserving (cptp) map, and the corresponding divisibility problem asks: when can a cptp map T be decomposed as T = R • R, where R is itself cptp? The continuous version of this, whether a cptp can be embedded into a completely-positive semi-group, is sometimes called the Markovianity problem in physics [8]-the latter again has applications to subdivision coding of quantum channels in quantum information theory [26].
Instead of dynamics, we can also ask whether the description of the static state of a system can be subdivided into smaller, simpler parts. Once again, probability theory provides a rich source of such problems. The most basic of these is the classic topic of divisible distributions: given a random variable X, can it be decomposed into X = Y +Z where Y, Z are some other random variables? What if Y and Z are identically distributed? If we instead ask for a decomposition into infinitely many random variables, this becomes the question of whether a distribution is infinitely divisible.
In this work, we address two of the most long-standing open problems on divisibility: divisibility of stochastic matrices, and divisibility and decomposability of probability distributions. We also extend our results to divisibility of nonnegative matrices and completely positive maps. Surprisingly little is known about the divisibility of stochastic matrices. Dating back to 1962 [19], the most complete characterization remains for the case of a 2 × 2 stochastic matrix [14]. The infinite divisibility problem has recently been solved [8], but the finite case remains an open problem. Divisibility of random variables, on the other hand, is a widely-studied topic. Yet, despite first results dating back as far as 1934 [5], no general method of answering whether a random variable can be written as the sum of two or more random variables-whether distributed identically, or differently-is known.
We focus on the computational complexity of these divisibility problems. In each case, we show which of the divisibility problems have efficient solutions-for these, we give an explicit efficient algorithm. For all other cases, we prove reductions to the famous P = NP-conjecture, showing that those problems are NP-hard. This essentially implies that-unless P = NP-the geometry of the corresponding divisible and non-divisible is highly complex, and these sets have no simple characterization beyond explicit enumeration. In particular, this shows that any future concrete classification of these NP-hard problems will be at least as hard as answering P = NP.
The following theorems summarize our main results on maps. Precise formulations and proofs can be found in section 2.
Theorem 1. Given a stochastic matrix P, deciding whether there exists a stochastic matrix Q such that P = Q 2 is NP-complete.

Theorem 2. Given a cptp map B, deciding whether there exists a cptp map A such that B = A • A is NP-complete.
In fact, the last two theorems are strengthenings of the following result. The following theorems summarize our main results on distributions. Precise formulations and proofs can be found in section 3.
Theorem 4. Let X be a finite discrete random variable. Deciding whether X is n-divisible-i.e. whether there exists a random variable Y such that X = n i=1 Yis in P.
Theorem 5. Let X be a finite discrete random variable, and > 0. Deciding whether there exists a random variable Y -close to X such that Y is n-divisible, or that there exists such a Y that is nondivisible, is contained in P.
Theorem 6. Let X be a finite discrete random variable. Deciding whether X is decomposable-i.e. whether there exist random variables Y, Z such that X = Y + Z-is NP-complete.
Theorem 7. Let X be a finite discrete random variable, and > 0. Deciding whether there exists a random variable Y -close to X such that Y is decomposable, or that there exists such an -close Y that is indecomposable, is NP-complete.
It is interesting to contrast the results on maps and distributions. In the case of maps, the homogeneous 2-divisibility problems are already NP-hard, whereas finding an inhomogeneous decomposition is straightforward. For distributions, on the other hand, the homogeneous divisibility problems are efficiently solvable to all orders, but becomes NP-hard if we relax it to the inhomogeneous decomposability problem.
This difference is even more pronounced for infinite divisibility. The infinite divisibility problem for maps is NP-hard (shown in [8]), whereas the infinite divisibility and decomposability problems for distributions are computationally trivial, since indivisible and indecomposable distributions are both dense-see section 3.5.8 and 3.4.5.
The paper is divided into two parts. We first address stochastic matrix and cptp divisibility in section 2, obtaining results on entry-wise positive matrix roots along the way. Divisibility and decomposability of probability distributions is addressed in section 3. In both sections, we first give an overview of the history of the problem, stating previous results and giving precise definitions of the problems. We introduce the necessary notation at the beginning of each section, so that each section is largely self-contained.

Introduction
Mathematically, subdividing Markov chains is known as the finite divisibility problem. The simplest case is the question of finding a stochastic root of the transition matrix (or a cptp root of a cptp map in the quantum setting), which corresponds to asking for the evolution over half of the time interval. While the question of divisibility is rather simple to state mathematically, it is not clear a priori whether a stochastic matrix root for a given stochastic matrix exists at all. Historically, this has been a long-standing open question, dating back to at least 1962 [19]. Matrix roots were also suggested early on in other fields, such as economics and general trade theory, at least as far back as 1967 [31], to model businesses and the flow of goods. Despite this long history, very little is known about the existence of stochastic roots of stochastic matrices. The most complete result to date is a full characterization of 2 × 2 matrices, as given for example in [14]. The authors mention that ". . . it is quite possible that we have to deal with the stochastic root problem on a case-by-case basis." This already suggests that there might not be a simple mathematical characterization of divisible stochastic matrices-meaning one that is simpler than enumerating the exponentially many roots and checking each one for stochasticity.
There are similarly few results if we relax the conditions on the matrix normalization slightly, and ask for (entry-wise) nonnegative roots of (entry-wise) nonnegative matrices-for a precise formulation, see Definition 10 and 11. An extensive overview can be found in [24]. Following this long history of classical results, quantum channel divisibility recently gained attention in the quantum information literature. The foundations were laid in [33], where the authors first introduced the notion of channel divisibility. A divisible quantum channel is a cptp map that can be written as a nontrivial concatenation of two or more quantum channels.
A related question is to ask for the evolution under infinitesimal time steps, which is equivalent to the existence of a logarithm of a stochastic matrix (or cptp map) that generates a stochastic (resp. cptp) semi-group. Classically, the question is known as Elfving's problem or the embedding problem, and it seems to date back even further than the finite case, to 1937 [10]. In the language of Markov chains, this corresponds to determining whether a given stochastic matrix can be embedded into an underlying continuous time Markov chain. Analogously, infinite quantum channel divisibility-also known as the Markovianity condition for a cptp map-asks whether the dynamics of the quantum system can be described by a Lindblad master equation [21,12]. The infinite divisibility problems in both the classical and quantum case were recently shown to be NP-hard [8]. Formulated as weak membership problems, these results imply that it is NP-hard to extract dynamics from experimental data [7].
However, while related, it is not at all clear that there exists a reduction of the finite divisibility question to the case of infinite divisibility. In fact, mathematically, the infinite divisibility case is a special case of finite divisibility, as a stochastic matrix is infinitely divisible if and only if it admits an n th root for all n ∈ N [19].
The finite divisibility problem for stochastic matrices is still an open question, as are the nonnegative matrix and cptp map divisibility problems. We will show that the question of existence of stochastic roots of a stochastic matrix is NP-hard. We also extend this result to (doubly) stochastic matrices, nonnegative matrices, and cptp maps.
We start out by introducing the machinery we will use to prove Theorem 1 and 3 in section 2.2. A reduction from the quantum to the classical case can be found in section 2.4, from the nonnegative to the stochastic case in section 2.5 and the main result-in a mathematically rigorous formulation-is then presented as Theorem 20 in section 2.6.

Roots of matrices
In our study of matrix roots we restrict ourselves to the case of square roots. The more general case of p th roots of matrices remains to be discussed. We will refer to square roots simply as roots. To be explicit, we state the following definition.
Following the theory of matrix functions-see for example [15]-we remark that in the case of nonsingular M, √ M is nonempty and can be expressed in Jordan normal x of the Jordan block corresponding to the ith eigenvalue λ i , has the property that for all i ∈ N ≥0 , no more than one element of the sequence satisfies d i ∈ (2i, 2(i + 1)).
For a given matrix, the classification gives the set of all roots. If M is a real matrix, a similar theorem holds and there exist various numerical algorithms for calculating real square roots, see for example [15].

Roots of stochastic matrices
Remember the following two definitions.
Definition 11. A matrix Q ∈ K d×d is said to be stochastic if it is nonnegative and In contrast to finding a general root of a matrix, very little is known about the existence of nonnegative roots of nonnegative matrices-or stochastic roots of stochastic matrices-if d ≥ 3. For stochastic matrices and in the case d = 2, a complete characterization can be given explicitly, and for d ≥ 3, all real stochastic roots that are functions of the original matrix are known, as demonstrated in [14]. Further special classes of matrices for which a definite answer exists can be found in [16]. But even for d = 3, the general case is still an open question-see [20, ch. 2.3] for details.
Indeed, a stochastic matrix may have no stochastic root, a primary or nonprimary root-or both. To make things worse, if a matrix has a p th stochastic root, it might or might not have a q th stochastic root if p q-p is not a divisor of q-, q > p or q p, q < p.
A related open problem is the inverse eigenspectrum problem, as described in the extensive overview in [9]. While the sets Ω n ⊂ D-denoting all the possible valid eigenvalues of an n-dimensional stochastic matrix-can be given explicitly, and hence also Ω p n , almost nothing is known about the sets of valid eigenspectra. Any progress in this area might yield necessary conditions for the existence of stochastic roots.
In recent years, some approaches have been developed to approximate stochastic roots numerically, see the comments in [14, sec. 4]. Unfortunately, most algorithms are highly unstable and do not necessarily converge to a stochastic root. A direct method using nonlinear optimization techniques is difficult and depends heavily on the algorithm employed [20].
It remains an open question whether there exists an efficient algorithm that decides whether a stochastic matrix Q has a stochastic root.
In this paper, we will prove that this question is NP-hard to answer.

The Choi isomorphism
For the results on cptp maps, we will need the following basic definition and results. A map A which is completely positive and trace-preserving-i.e. tr(Aρ) = tr ρ ∀ρ ∈ H-is called a completely positive trace-preserving map, or short cptp map.
Remark 13. Let the notation be as in Definition 12 and pick a basis e 1 , . . . , e d of C d . Then A is completely positive if and only if the Choi matrix The condition of trace-preservation then translates to the following.

Remark 14.
A map A is trace-preserving if and only if tr 2 (C A ) = 1 d , where tr 2 denotes the partial trace over the second pair of indices.

Equivalence of computational questions
In the following we denote with S some arbitrary finite index set, not necessarily the same for all problems. We begin by defining the following decision problems.  Theorem 20. The reductions as shown in Fig. 1 hold.
Proof. The implication Stochastic Divisibility←−Stochastic Root needs one intermediate step. If P is not stochastic, the answer is negative. If it is stochastic, we can apply Stochastic Divisibility. The opposite direction holds for non-derogatory stochastic P: in this case we can enumerate all roots of P as a finite family which forms a valid instance for Stochastic Root.
The reduction Stochastic Root←−Nonnegative Root can be resolved by Lemma 25 and Lemma 26-we construct a family of matrices (Q s ) s∈S that contains a stochastic root iff (M s ) s∈S contains a nonnegative root. The result then follows from applying Stochastic Root. If our stochastic matrix P is irreducible, then any nonnegative root Q s : Q 2 s = P is stochastic, and in that case Stochastic Root←→Nonnegative Root-see [16, sec. 3] for details.
The link cptp Divisibility←−cptp Root again needs the following intermediate step. If A is not cptp, the answer is negative. If it is cptp, then we can apply cptp Divisibility. Similarly, if A is non-derogatory, the reduction works in the opposite direction as well.
The direction cptp Root←−Stochastic Root follows from Corollary 24. We start out with a family (Q s ) s∈S comprising all the roots of a stochastic matrix P. Then let (A s := emb Q s ) s∈S -this family then comprises all of the roots of B := A 2 k ≡ A 2 s ∀k, s. Furthermore, by Lemma 23, there exists a cptp A s if and only if there exists a stochastic Q s , and the reduction follows. Finally, we can extend our reduction to the programs Doubly Stochastic Root and Doubly Stochastic Divisibility as well as Nonnegative Divisibility, defined analogously, see our comment in Corollary 27 and the complete reduction tree in Fig. 1. 2 At this point, we observe the following fact. Proof. It is straightforward to come up with a witness and a verifier circuit that satisfies the definition of the decision class NP. For example in the cptp case, a witness is a matrix root that can be checked to be a cptp map using Remark 13 and squared in polynomial time, which is the verifier circuit. Both circuit and witness are clearly poly-sized and hence the claim follows. 2 By encoding an instance of 1-in-3sat into a family of nonnegative matrices (M s ) s∈S , we show the implication 1-in-3sat−→Nonnegative Root and 1-in-3sat−→(Doubly) Stochastic/cptp Divisibility, accordingly, from which NP-hardness of (Doubly) Stochastic/cptp Divisibility follows. The entire chain of reduction can be seen in Fig. 1.

Reduction of Stochastic Root to CPTP Root
This reduction is based on the following embedding.

Definition 22.
Let {e i } be an orthonormal basis of K d . The embedding emb is defined as We observe the following.
Lemma 23. We use the same notation as in Remark 13. Let A ∈ K d×d and B := emb A.
Proof. The first claim follows directly from the matrix representation of our operators. There, the Choi isomorphism is manifest as the reshuffling operation or partial transpose For more details, see e.g. [2].
The second statement follows from The final claim is trivial. 2 This remark immediately yields the following consequence.

Reduction of Nonnegative Root to Stochastic Root
The difference between Nonnegative Root and Stochastic Root is the extra normalization condition in the latter, see Definition 11. The following two lemmas show that this normalization does not pose an issue, so we can efficiently reduce the problem Nonnegative Root to Stochastic Root. Proof. We explicitly construct our family (Q s ) s∈S as follows. Pick an s ∈ S and denote M := M s . Let d be the dimension of M. We first pick a ∈ R + such that a max ij M ij = 1/2 1 and define where by sum of matrix M and scalar x we mean M + x1, 1 := (1) 1≤i,j≤d ∈ R d×d , and 1 The exact bound is a max ij M ij ≤ 43/81.
Observe that {A, B, C} form an orthogonal set-if one wishes, normalizing and pulling out the constant as eigenvalue to the corresponding eigenprojectors would work equally well. By construction, Q s is nonnegative if and only if M s is. Since the row sums of Q s are always 1, Q s is stochastic if and only if M s is nonnegative, and the claim follows. 2 Lemma 26. Let the notation be as in Lemma 25 and write √ N for the set of roots of N, Furthermore, the complement of (Q s ) s∈S in √ P does not contain any stochastic roots.
Proof. The first statement is obvious, since for all s ∈ S, and hence clearly (Q s ) s∈S ⊂ √ P. The last statement is not quite as straightforward-it is the main reason our carefully crafted matrix Q s has its slightly unusual shape. All possible roots of P are of the form It is easy to check that none of the other sign choices yields any stochastic matrix, so the claim follows. 2 2 Corollary 27. The results of Lemma 25 and 26 also hold for doubly stochastic matricesobserve that our construction of Q s is already doubly stochastic.

Reduction of 1-in-3 sat to Nonnegative Root
We now embed an instance of a boolean satisfiability problem, 1-in-3sat-see Definition 87 for details-into a family of matrices (M s ) s∈S in a way that there exists an s such that M s is nonnegative if and only if the instance of 1-in-3sat is satisfiable. The construction is inspired by [8].
We identify Denote with (m i1 , m i2 , m i3 ) ∈ {±1} 3 the three boolean variables occurring in the i th boolean clause, and let m i ∈ {±1} stand for the single i th boolean variable. Then 1-in-3sat translates to the inequalities Then there exists a family of matrices (M s ) s∈S such that ∃s : M s nonnegative iff the instance is satisfiable.
To prove this, we first need the following technical lemma.
Then there exists a family of matrices (C s ) s∈S such that ∃s: the first n c on-diagonal 4 ×4 blocks of C s are nonnegative iff the instance is satisfiable. In addition, we have s , and the complement contains no nonnegative root.
Proof. For every boolean variable m k , define a vector v k ∈ R d such that their first n c elements are defined via We will specify the dimension d later-obviously d ≥ n c , and the free entries are used to orthonormalize all vectors in the end. For now, we denote the orthonormalization region with o. We further define the vectors c 1 , c 2 ∈ R d to have all 1s in the first n c entries, i.e. c 1,2 = (1, . . . , 1, o 1,2 ). Let then The variables p k denote a specific rescaled choice of the boolean variables m i , which-in order to avoid degeneracy-have to be distinct, i.e. via The p ij are defined accordingly from the m ij and N ∈ N is large but fixed. Let further where we have used an obvious block notation to pad C s with zeroes, which will come into play later. Hatching signifies complex numbers, the colour scale is the same as in Fig. 3, i.e. light green denotes negative numbers, dark purple nonnegative entries.
The on-diagonal 2 × 2 blocks of C s then encode the 1-in-3sat inequalities from equation (2)-demanding nonnegativity-as the set of equations Note that we leave enough head space such that the rescaling in equation (4) does not affect any of the inequalities-see section 2.8 for details. Observe further that the eigenvalues corresponding to each eigenprojector in the last term of equation (3) necessarily have opposite sign, otherwise we create complex entries. We will later rescale C s by a positive factor, under which the inequalities are invariant, so the first claim follows.
We can always orthonormalize the vectors c 1,2 and v k using the freedom left in o, hence we can achieve that C 2 s = C 2 t ∀s, t. It is straightforward to check that no other sign choice for the eigenvalues of the first two terms yields nonnegative blocks-see Fig. 2 for details. From this, the last two claims follow. 2

Orthonormalization and handling the unwanted inequalities
As in [8], we have unwanted inequalities-the off-diagonal blocks in the first 4n c entries and the blocks involving the orthonormalization region o. We first deal with the off-diagonal blocks in favour of enlarging the orthonormalization region, creating more-potentially negative-entries in there, and then fix the latter.
Off-diagonal blocks. We begin with the following lemma.
Lemma 30. Let the family (C s ) s∈S be defined as in the proof of Lemma 29, and (n v , n c , m i , m ij ) the corresponding 1-in-3 sat instance. Then there exists a matrix E ∈ C d×d such that the top left 4n c × 4n c block of C s + E has at least one negative entry ∀s iff the instance is not satisfiable. Furthermore, imC s ⊥ imE ∀s, and C s + E has negative entries ∀s, ∀E ∈ √ Proof. Define Then E 1 has rank 1.
The variables t i are chosen close to 1 but distinct, e.g.
where M ∈ N large but fixed. Then E has rank n c + 1, and adding E to C s trivializes all unwanted inequalities in the upper left 4n c × 4n c block. By picking M large enough, the on-diagonal inequalities are left intact.
One can check that all other possible sign choices for the roots of E create negative entries in parts of the upper left block where C s is zero ∀s. Furthermore, C s and E have distinct nonzero eigenvalues by construction-the orthogonality condition is again straightforward, hence the last two claims follow. 2 Orthonormalization region.
Proof. Define One branch of an unsatisfiable instance of 1-in-3sat encoded into a matrix of total rank 19. The negative entries-two bright dots-in the upper left block in the combined matrix (d) indicate that this branch does not satisfy the given instance. By looking at all other blocks, one sees that none translates to a nonnegative matrix. Observe that in this naïve implementation the orthonormalization region is suboptimally large.
where 0 < a < 1 is used to orthonormalize Δ and E 2 , which is the case if By explicitly writing out the rank 2 matrix it is straightforward to check that D fulfils all the claims of the lemma-see Fig. 3 for an example. 2

Lifting singularities
The reader will have noted by now that even though we have orthonormalized all our eigenspaces, ensuring that the nonzero eigenvalues are all distinct, we have at the same time introduced a high-dimensional kernel in C s , E and D. The following lemma shows that this does not pose an issue. Proof. Take a matrix A ∈ (A s ) s∈S . We need to distort the zero eigenvalues {λ (0) i } slightly away from 0. Using notation from Definition 36, a conservative estimate for the required smallness without affecting positivity would be where we used the Jordan canonical form A = ZΛZ −1 for some invertible Z and Λ = diag(J 0 , J 1 ), such that J 0 collects all Jordan blocks corresponding to the eigenvalue 0, and J 1 collects the remaining ones. 2 This will lift all remaining degeneracies and singularities, without affecting our line of argument above. Observe that all inequalities in our construction were bounded away from 0 with enough head space independent of the problem size, so positivity in the lemma is sufficient.
We thus constructed an embedding of 1-in-3sat into non-derogatory and nondegenerate matrices, as desired. It is crucial to note that we do not lose anything by restricting the proof to the study of these matrices, as the following lemma shows.
Lemma 33. There exists a Karp reduction of the Divisibility problems when defined for all matrices to the case of non-degenerate and non-derogatory matrices.
Proof. As shown in Lemma 21, containment in NP for this problem is easy to see, also in the degenerate or derogatory case. Since 1-in-3sat is NP-complete, there has to exist a poly-time reduction of the Divisibility problems-when defined for all matricesto 1-in-3sat. Now embed this 1-in-3sat-instance with our construction. This yields a poly-time reduction to the non-degenerate non-derogatory case. 2

Complete embedding
We now finally come to the proof of Theorem 28.
Proof of Theorem 28. Construct the family (C s +E) s∈S using Lemma 29 and Lemma 30, ensuring that all orthonormalizing is done, which preliminarily fixes the dimension d. By Lemma 31, we now construct a mask D(δ) of dimension d + d , where d > 0 is picked such that we can also orthonormalize all previous vectors with respect to E 2 and δ.
By Lemmas 29, 30, 31 and 32, the perturbed family (M s ) s∈S := (C s + E + N D(δ)) swhere N and δ ∈ Q are chosen big enough so that all unwanted inequalities are trivially satisfied-fulfils the claims of the theorem and the proof follows. 2 We finalize the construction as follows. In Theorem 28, we have embedded a given 1-in-3sat instance into a family of matrices (M s ) s∈S , such that the instance is satisfiable if and only if at least one of those matrices is nonnegative.
By rescaling the entire matrix such that max ij (M s ) ij = 1/2, we could show that this instance of 1-in-3sat is satisfiable if and only if the normalized matrix family (Q s ) s∈S , which we construct explicitly, contains a stochastic matrix.
As shown in section 2.3, this can clearly be answered by Stochastic Divisibility, as the family (Q s ) s∈S comprises all the roots of a unique matrix P. If this matrix is not stochastic, our instance of 1-in-3sat is trivially not satisfiable. If the matrix is stochastic, we ask Stochastic Divisibility for an answer-a positive outcome signifies satisfiability, a negative one non-satisfiability.

Bit complexity of embedding
To show that our results holds for only polynomially growing bit complexity, observe the following proposition. Proof. We can ignore any construction that multiplies by a constant prefactor, for example Lemma 25 and Lemma 26. The renormalization for Lemma 25 to max ij M s,ij = 1/2 does not affect r either.
The rescaling in equation (4) and equation (5)  The only other place of concern is the orthonormalization region. Let us write a i for all vectors that need orthonormalization. In the n th step, we need to make up for O(n) entries with our orthonormalization, using the same amount of precision to solve the linear equations (a T i a n = 0) 1≤i<n . This has to be done with a variant of the standard Gauss algorithm, e.g. the Bareiss algorithm-see for example [1]-which has nonexponential bit complexity.
Together with the lifting of our singularities, which has polynomial precision, we obtain r(M s ) = O (poly(n v , n c )). Completing the embedding in section 2.9 changes the bit complexity by at most another polynomial factor, and hence the claim follows. 2

Introduction
Underlying stochastic and quantum channel divisibility, and-to some extent-a more fundamental topic, is the question of divisibility and decomposability of probability distributions and random variables. An illustrative example is the distribution of the sum of two rolls of a standard six-sided die, in contrast to the single roll of a twelve-sided die. Whereas in the first case the resulting random variable is obviously the sum of two uniformly distributed random variables on the numbers {1, . . . , 6}, there is no way to achieve the outcome of the twelve-sided die as any sum of nontrivial "smaller" dice-in fact, there is no way of dividing any uniformly distributed discrete random variable into the sum of non-constant random variables. In contrast, a uniform continuous distribution can always be decomposed 3 into two different distributions.
To be more precise, a random variable X is said to be divisible if it can be written as X = Y + Z, where Y and Z are non-constant independent random variables that are identically distributed (iid). Analogously, infinite divisibility refers to the case where X can be written as an infinite sum of such iid random variables.
If we relax the condition Y d = Z-i.e. we allow Y and Z to have different distributions-we obtain the much weaker notion of decomposability. This includes using other sources of randomness, not necessarily uniformly distributed.
Both divisibility and decomposability have been studied extensively in various branches of probability theory and statistics. Early examples include Cramer's theorem [6], proven in 1936, a result stating that a Gaussian random variable can only be decomposed into random variables which are also normally distributed. A related result on χ 2 distributions by Cochran [5], dating back to 1934, has important implications for the analysis of covariance.
An early overview over divisibility of distributions is given in [28]. Important applications of n-divisibility-the divisibility into n iid terms-is in modelling, for example of bug populations in entomology [18], or in financial aspects of various insurance models [30,29]. Both examples study the overall distribution and ask if it is compatible with an underlying subdivision into smaller random events. The authors also give various conditions on distributions to be infinitely divisible, and list numerous infinitely divisible distributions.
Important examples for infinite divisibility include the Gaussian, Laplace, Gamma and Cauchy distributions, and in general all normal distributions. It is clear that those distributions are also finitely divisible, and decomposable. Examples of indecomposable distributions are Bernoulli and discrete uniform distributions.
However, there does not yet exist a straightforward way of checking whether a given discrete distribution is divisible or decomposable. We will show in this work that the question of decomposability is NP-hard, whereas divisibility is in P. In the latter case, we outline a computationally efficient algorithm for solving the divisibility question. We extend our results to weak-membership formulations (where the solution is only required to within an error in total variation distance), and argue that the continuous case is computationally trivial as the indecomposable distributions form a dense subset.
We start out in section 3.2 by introducing general notation and a rigorous formulation of divisibility and decomposability as computational problems. The foundation of all our distribution results is by showing equivalence to polynomial factorization, proven in section 3.3. This will allow us to prove our main divisibility and decomposability results in section 3.4 and 3.5, respectively.

Discrete distributions
In our discussion of distribution divisibility and decomposability, we will use the standard notation and language as described in the following definition.
Definition 35. Let (Ω, F, p) be a discrete probability space, i.e. Ω is at most countably infinite and the probability mass function p : Ω −→ [0, 1]-or pmf, for short-fulfils x∈Ω p(x) = 1. We take the σ-algebra F to be maximal, i.e. F = 2 Ω , and without loss of generality assume that the state space Ω = N. Denote the distribution described by p with D. A random variable X : Ω −→ B is a measurable function from the sample space to some set B, where usually B = R.
For the sake of completeness, we repeat the following well-known definition of characteristic functions.
Definition 36. Let D be a discrete probability distribution with pmf p, and X ∼ D. It is well-known that two random variables with the same characteristic function have the same cumulative density function.
Definition 37. Let the notation be as in Definition 35. Then the distribution D is called finite if p(k) = 0 ∀k ≥ N for some N ∈ N.
Remark 38. Let D be a discrete probability distribution with pmf p. We will-without loss of generality-assume that p(0) = 0 and p(k) = 0 ∀k < 0 for the pmf p of a finite distribution. It is a straightforward shift of the origin that achieves this.

Continuous distributions
Definition 39. Let (X , A) be a measurable space, where A is the σ-algebra of X . The probability distribution of a random variable X on (X , A) is the Radon-Nikodym derivative f , which is a measurable function with P(X ∈ A) = A f dμ, where μ is a reference measure on (X , A).
Observe that this definition is more general than Definition 35, where the reference measure is simply the counting measure over the discrete sample space Ω. Since we are only interested in real-valued univariate continuous random variables, observe the following important remark.
Remark 40. We restrict ourselves to the case of X = R with A the Borel sets as measurable subsets and the Lebesgue measure μ. In particular, we only regard distributions with a probability density function f -or pdf, for short-i.e. we require the cumulative distribution function P(x) := P(X ≤ x) ≡ y≤x f (y)dy to be absolutely continuous.
Corollary 41. The cumulative distribution function P of a continuous random variable X is almost everywhere differentiable, and any piecewise continuous function f with R f (x)dx = 1 defines a valid continuous distribution.

Divisibility and decomposability of distributions
To make the terms mentioned in the introduction rigorous, note the two following definitions.
Definition 42. Let X be a random variable. It is said to be n-decomposable if X = Z 1 + . . . + Z n for some n ∈ N, where Z 1 , . . . , Z n are independent non-constant random variables. X is said to be indecomposable if it is not decomposable.

Definition 43. Let X be a random variable. It is said to be
If we are not interested in the exact number of terms, we also simply speak of decomposable and divisible. We will show in section 3.5.7 that-in contrast to divisibility-the question of decomposability into more than two terms is not well-motivated.
Observe the following extension of Remark 38.
Lemma 44. Let D be a discrete probability distribution with pmf p. If p obeys Remark 38, then we can assume that its factors do as well. In the continuous case, we can without loss of generality assume the same.
Proof. Obvious from positivity of convolutions in case of divisibility. For decomposability, we can achieve this by shifting the terms symmetrically. 2

Markov chains
To establish notation, we briefly state some well-known properties of Markov chains.
Remark 45. Take discrete iid random variables Y 1 , . . . , Y n ∼ D and write P(Y i = k) = p k := p(k) for all k ∈ N, independent of i = 1, . . . , n. Define further Then {X n , n ≥ 0} defines a discrete-time Markov chain, since This last property is also called stationary independent increments, i.e. we add an iid random variable at each step.
Remark 46. Let the notation be as in Remark 45. The transition probabilities of the Markov chain are then given by In matrix form, we write the transition matrix Working with transition matrices is straightforward-if the initial distribution is given by π := (1, 0, . . .), then obviously (πP) i = p i . Iterating P then yields the distributions of X 2 , X 3 , . . ., respectively-e.g. (πP 2 ) i = P(X 2 = i) ≡ P(Y 1 + Y 2 = i).
We know that X 2 is divisible-namely into X 2 = Y 1 + Y 2 , by construction-but what if we ask this question the other way round? We will show in the next section that there exists a relatively straightforward way to calculate if an (infinite) matrix in the shape of P has a stochastic root-i.e. if D is divisible. Observe that this is not in contradiction with Theorem 1, as the theorem does not apply to infinite operators.
In contrast, the more general question of whether we can write a finite discrete random variable as a sum of nontrivial, potentially distinct random variables will be shown to be NP-hard.

Equivalence to polynomial factorization
Starting from our digression in section 3.2.4 and using the same notation, we begin with the following definition.
Definition 47. Denote with S the shift matrix S ij := δ i+1,j . Then we can write Since S just acts as a symbol, we write and f ∼ g :⇔ f = cg, c > 0. We call f D the characteristic polynomial of D-not to be confused with the characteristic polynomial of a matrix. The equivalence space R defines the set of all characteristic polynomials, and can be written as We mod out the overall scaling in order to keep the normalization condition k p(k) = 1 implicit-if we write f D , we will always assume f D (1) = 1. An alternative way to define these characteristic polynomials is via characteristic functions, as given in Definition 36.
The reason for this definition is that it allows us to reduce operations on the transition matrix P or products of characteristic functions φ X to algebraic operations on f D . This enables us to translate the divisibility problem into a polynomial factorization problem and use algebraic methods to answer it. Observe that the R N are normed vector spaces, which we will make use of later. While this might seem obvious, it is worth clarifying, since this correspondence will allow us to directly translate results on polynomials to distributions.
The following lemma reduces the question of divisibility and decomposability-see Definition 42 and 43-to polynomial factorization.
Proof. Assume that D is n-divisible, i.e. that there exists a distribution D and random variables Z 1 , . . . , Z n ∼ D such that X = n i=1 Z i . Denote with Q the transition matrix of D , as defined in Remark 46, and write q for its probability mass function. Then as before. Write g D for the characteristic polynomial of D . By Definition 47, g n (S) ≡ f D (S), and hence g n D = f D . Observe that and hence i q(i) = 1 is normalized automatically. The other direction is similar, as well as the case of decomposability, and the claim follows. 2

Computational problems
We state an exact variant of the computational formulation of the question according to Definition 43-i.e. one with an allowed margin of error-as well as a weak membership formulation.
Instance. Finite discrete random variable X ∼ D.

Question. Does there exist a finite discrete distribution
Observe that this includes the case n = 2, which we defined in Definition 43.
Proof. By Lemma 51 it is enough to show that for a characteristic polynomial f ∈ R N , we can find a g ∈ R : g n = f in polynomial time. In order to achieve this, write (f ) 1/n as a Taylor expansion with rest, i.e.
If R ≡ 0, then g = p n-divides f , and then the distribution described by f is n-divisible. Since the series expansion is constructive and can be done efficiently-see [25]-the claim follows.
If the distribution coefficients are rational numbers, another method is to completely factorize the polynomial-e.g. using the LLL algorithm, which is known to be easy in this setting-sort and recombine the linear factors, which is also in O(poly(ordf )), see for example [13]. Then check if all the polynomial root coefficients are positive. 2 We collect some further facts before we move on.
Remark 55. Let p be the probability mass function for a finite discrete distribution D, and write supp p = {k : p(k) = 0}. If max supp p − min supp p =: w, then D is obviously not n-divisible for n > w/2, and furthermore not for any n that do not divide w, n < w/2. Indeed, D is not n-divisible if the latter condition holds for either max supp p or min supp p.
Remark 56. Let X ∼ D be an n-divisible random variable, i.e. ∃Z 1 , . . . , Z n ∼ D : Proof. This is clear, because R[x] is a unique factorization domain. 2

Divisibility with variation
As an intermediate step, we need to extend Theorem 54 to allow for a margin of error , as captured by the following definition.
Instance. Finite discrete random variable X ∼ D with pmf p X (k). Question. Do there exist finite random variables Z 1 , . . . , Z n ∼ D with pmfs p Z (k), such that p Z * . . . * p Z n times −p X ∞ < ?
Lemma 58. Distribution Divisibility n, is in P.
Proof. Let f (x) = N i=0 p i x i be the characteristic polynomial of a finite discrete distribution, and > 0. By padding the distribution with 0s, we can assume without loss of generality that N = deg f is a multiple of n. A polynomial root-if it exists-has the form g( Comparing coefficients in the divisibility condition f (x) = g(x) n , the latter translates to the set of inequalities a n 0 ∈ (p 0 − , p 0 + ) . . . h i (a 1 , . . . , a i−1 ))/na n−i 0 ). It is now easy to solve the system iteratively, keeping track of the allowed intervals I i for the a i .

Each term but the first one is of the form
If I i = ∅ for some i, we return No, otherwise Yes. We have thus developed an efficient algorithm to answer Distribution Weak Divisibility n, , and the claim of Lemma 58 follows. 2 Remark 59. Given a random variable X, the algorithm constructed in the proof of Lemma 58 allows us to calculate the closest n-divisible distribution to X in polynomial time.
Proof. Straightforward, e.g. by using binary search over . 2

Weak divisibility
For the weak membership problem, we reduce Weak Distribution Divisibility n, to Distribution Divisibility n, .
Proof. Let D be a finite discrete distribution. If Distribution Divisibility n, answers Yes, we know that there exists an n-divisible distribution -close to D. In case of No, D itself is not n-divisible, hence we know that there exists a non-n-divisible distribution close to D. 2

Continuous distributions
Let us briefly discuss the case of continuous distributions-continuous meaning a nondiscrete state space X , as specified in section 3.2.2. Although divisibility of continuous distributions is well-defined and widely studied, formatting the continuous case as a computational problem is delicate, as the continuous distribution must be specified by a finite amount of data for the question to be computationally meaningful. The most natural formulation is the continuous analogue of Definition 43 as a weak-membership problem. However, we can show that this problem is computationally trivial.
First observe the following intermediate result.
Lemma 61. Take  Proof. Due to symmetry, it is enough to show divisibility for f | A . Assuming f is divisible, we can write f = r * r, i.e. f (x) = R r(x − y)r(y)dy. It is straightforward to show that r(x) = 0 ∀x < 0. Definer We see that (r * r)(x) = 0 for x / ∈ A. For x ∈ A, the support of the integrand is contained in {y : y ∈ x −A/2 ∧y ∈ A/2} = x −A/2 ∩A/2 := S x , and hence we can write (r * r)(x) = The integrand r(x − y)r(y) = 0 ∀y < 0 ∨ y > x. The difference in the integration domains can be seen in Fig. 4. We get two cases.
Let x ∈ A. Assume ∃y ∈ (M/2, M ) such that r(x − y )r(y ) > 0. Let x := 2y . We then have r(y ) 2 = r(x − y )r(y ) > 0, and due to continuity f (x ) > 0, contradiction, because x ∈ (M, 2M ). Analogously Proof. It is enough to show the claim for functions f ∈ C + c,b with inf supp f ≥ 0. Let > 0, and M := sup supp f . Take j ∈ C c,b to be nondivisible with supp j ⊂ (2M, 3M ), and define By construction, f − g ∞ < , but g| (2M,3M ) ≡ j is not divisible, hence by Lemma 61 g is not divisible, and the claim follows. 2 Corollary 63. Let > 0. Let X be a continuous random variable with pdf p X (k). Then there exists a nondivisible random variable Y with pdf p Y (k), such that p X − p Y < .
and Proposition 62 finishes the claim. 2 Corollary 64. Any weak membership formulation of divisibility in the continuous setting is trivial to answer, as for all > 0, there always exists a nondivisible distribution close to the one at hand. Similar considerations apply to other formulations of the continuous divisibility problem.

Infinite divisibility
Let us finally and briefly discuss the case of infinite divisibility. While interesting from a mathematical point of view, the question of infinite divisibility is ill-posed computationally. Trivially, discrete distributions cannot be infinitely divisible, as follows directly from Theorem 54. A similar argument shows that neither the , nor the weak variant of the discrete problem is a useful question to ask, as can be seen from Lemma 58 and 60.
By the same arguments as in section 3.4.5, the weak membership version is thus easy to answer and therefore trivially in P.

Computational problems
We define the decomposability analogue of Definition 52 and 53 as follows.
Instance. Finite discrete random variable X ∼ D.
Question. Do there exist finite discrete distributions D , D : Definition 66 (Weak Distribution Decomposability ).
Instance. Finite discrete random variable X ∼ D with pmf p X (k). Question. If there exists a finite discrete random variable Y with pmf p Y (k), such that p X − p Y ∞ < and such that 1. Y is decomposable-return Yes 2. Y is indecomposable-return No.
In this section, we will show that Distribution Decomposability is NP-hard, for which we will need a series of intermediate results. Requiring the support of the first random variable Z 1 to have a certain size, i.e. | supp(p D )| = m, yields the following program. We then define Distribution Even Decomposability to be the case where the two factors have equal support.
The full reduction tree can be seen in Fig. 5. Analogous to Lemma 21, we state the following observation.
Lemma 68. All the above Decomposability problems in Definition 65 to 66 are contained in NP.
Proof. It is straightforward to construct a witness and a verifier that satisfies the definition of the decision class NP. For example in Definition 78, a witness is given by two tables of numbers which are easily checked to form finite discrete distributions. Convolving these lists and comparing the result to the given distribution can clearly be done in polynomial time. Both verification and witness are thus poly-sized, and the claim follows. 2

Even decomposability
We continue by proving that Distribution Even Decomposability is NP-hard. We will make use of the following variant of the well-known Subset Sum problem, which is NP-hard-see Lemma 92 for a proof. The interested reader will find a rigorous digression in section A.1.

Definition 69 ( Even Subset Sum).
Instance. Multiset S of reals with |S| even, l ∈ R.
Question. Does there exist a multiset T S with |T| = |S|/2 and such that | t∈T t − s∈S\T s| < l?
This immediately leads us to the following intermediate result.
Proof. Let (S, l) be an instance of Even Subset Sum. We will show that there exists a polynomial f ∈ R of degree 2|S| such that f is divisible into f = g · h with deg g = deg h iff (S, l) is a Yes instance. We will explicitly construct the polynomial f ∈ R. As a first step, we transform the Even Subset Sum instance (S, l), making it suited for embedding into f .
Let N := |S| and denote the elements in S with s 1 , . . . , s N . We perform a linear transformation on the elements s i via where a ∈ R >0 is a free scaling parameter chosen later such that |b i | < δ ∈ R + small.
The next step is to construct the polynomial f and prove that it is divisible into two polynomial factors f = g · h if and only if (B, al) is a Yes instance. We first define quadratic polynomials g(b i , x) := x 2 + b i x + 1 for i = 1, . . . , N , and set f T (x) := b∈T g(b, x) for T ⊂ B. Observe that for suitably small δ, the g(b i , x) are irreducible over R[x]. With this notation, we claim that f B (x) has the required properties.
In order to prove this claim, we first show that for sufficiently small scaling parameter a, a generic subset T ⊂ B with n := |T| and f T (x) =: sgn(c 1 ) = sgn(Σ) , c 2j > 0 for j = 1, . . . , |T| , sgn(c 2j+1 ) ≥ sgn(Σ) for j = 1, . . . , |T| − 1 , where Σ := t∈T t. Indeed, if then f B = g · h, where g, h ∈ R, then g = f B 1 and h = f B 2 for aforementioned subsets B 1 , B 2 B, and conversely if (B, al) is a Yes instance, then is a unique factorization domain, so all polynomials of the shape f T necessarily decompose into quadratic factors. By construction, c 0 = 1 and c 1 = nΣ, so the first two assertions follow immediately. To address equation (11) and (12), we further split up the even and odd coefficients into where c j,k is the coefficient of We thus have c j,k = O(δ k ) in the limit δ → 0we will implicitly assume the limit in this proof and drop it for brevity. Our goal is to show that the scaling in δ suppresses the combinatorial factors, i.e. that c j is dominated by its first terms c j,0 and c j,1 , respectively. In order to achieve this, we need some more machinery. First regard g(δ, x) = x 2 + δx + 1. It is immediate that for an expansion We will calculate the coefficients d j,k of g(δ, x) n explicitly and use them to bound the coefficients c j,k of f T (x). Using a standard Cauchy summation and the uniqueness of polynomial functions, we obtain With (n) l we denote the falling factorial, i.e. (n) l = n(n − 1)(n − 2) · · · (n − l + 1). By convention, (n) 0 = 1. Regarding even and odd powers of x separately, we can thus deduce that A straightforward estimate shows that for the even and odd case, we obtain the coefficient scaling which means that e.g. picking δ = O(1/n 2 ) is enough to exponentially suppress the higher order combinatorial factors.
We will now separately address the even and odd cases-equation (11) and (12).
Even case. As the constant coefficients c j,0 = O(1) in δ, it is the same as for g(δ, x) n and by equation (14), we immediately get Odd case. Note that if Σ < 0, we are done, so assume Σ > 0 in the following. A simple combinatorial argument gives so it remains to show that c j,1 > −c j,3 − . . . − c j,j . Analogously to the even case, by equation (14), we conclude which finalizes our proof. 2

m-Support decomposability
In the next two sections we will generalize the last result to Distribution Decomposability m . As a first observation, we note the following.
Lemma 71. Let f (n) be such that (f (n)β(f (n), n + 1 − f (n))) −1 = O(poly(n)). Then Proof. See proof of Theorem 54, and an easy scaling argument for n f (n) completes the proof. As in Remark 91, this symmetrically extends to Distribution Observe that f (n) = n/2 yields exponential growth, hence the remark is consistent with the findings in section 3.5.2.
We now regard the general case. As in the last section, we need variants of the Subset Sum problem, which are given in the following two definitions.
Instance. Multiset S of reals with |S| even, l ∈ R. Question. Does there exist a multiset T S with |T| = m and such that | t∈T t − s∈S\T s| < l?
Instance. Multiset S of positive integers or reals, x, y ∈ R : x ≤ y. Question. Does there exist a multiset T ⊂ S with |T| = m and such that x < t∈T t − s∈S\T s < y?
Both are shown to be NP-hard in Lemma 90 and 94, or by the following observation. In order to avoid having to take absolute values in the definition of Subset Sum m , we reduce it to multiple instances of Signed Subset Sum m , by using the following interval partition of the entire range (−l, l).
This finally leads us to the following result.
Proof. We will show the reduction Distribution Decomposability m ←− Subset Sum m . Let m be fixed. Let (S, l) be an Subset Sum instance. For brevity, we write Σ S := s∈S s. Without loss of generality, by Corollary 89, we again assume Σ S ≥ 0. Now define a := 2(|S|l + 2mΣ S − |S|Σ S )/(2m − |S|). Using Remark 74, pick a suitable subdivision of the interval (−l − 2a, l + 2a), such that where we chose c(m, i) = x i /(2m − |S|). The latter program we can answer using the same argument as for the proof of Lemma 70, and the claim follows. 2

One can verify that
As a side remark, this also confirms the following well-known fact.

General decomposability
We have already invented all the necessary machinery to answer the general case.
Proof. Follows immediately from Lemma 70, where we regard the special set of Subset Sum instances for which (S, l) is such that l = s∈S s. We show in Lemma 96 that Subset Sum(·, Σ · ) is still NP-hard, thus the claim follows. 2

Decomposability with variation
As a further intermediate result-and analogously to Definition 57-we need to allow for a margin of error .

Definition 78 (Distribution Decomposability ).
Instance. Finite discrete random variable X ∼ D with pmf p X (k). Question. Do there exist finite discrete random variables This definition leads us to the following result.
Proof. First observe that we can restate this problem in the following equivalent form. Given a finite discrete distribution D with characteristic polynomial f D , do there exist two finite discrete distributions D , D with characteristic polynomials Here, we are using the maximum norm from Definition 49, and assume without loss of generality that deg As f D is a polynomial, we can regard its Viète map v : C n −→ C n , where n = deg f D , which continuously maps the polynomial roots to its coefficients. It is a well-known factsee [32] for a standard reference-that v induces an isomorphism of algebraic varieties we again implicitly assume the limit → 0.
We continue by proving the reduction Distribution Divisibility ←− Subset Sum(·, Σ · + poly ), which is NP-hard as shown in Lemma 97. Let S = {s i } N i=1 be a Subset Sum multiset. We claim that it is satisfiable if and only if the generated characteristic function f S (x)-where we used the notation of the proof of Lemma 70-defines a finite discrete probability distribution and the corresponding random variable X is a Yes instance for Distribution Divisibility .
First assume f S is such a Yes instance. Then s∈S s ≥ 0, and there exist two characteristic polynomials g = i g i and h = i h i as above and such that f S − gh d < . We also know that if where T S and B (T) denotes an ball around the set T, and analogously for Regarding the linear coefficients, we thus have Now the case if f S is a No instance. Assume there exists a nontrivial multiset T S satisfying Then by construction t∈T t, s∈S\T s ≥ −O( ) and f T · f S\T = f S , contradiction, and the claim follows. 2

Weak decomposability
Analogously to section 3.4.4, we now regard the weak membership problem of decomposability.
Proof. In order to show the claim, we prove the reduction Weak Distribution Decomposability ←− Distribution Decomposability g( ) , where the function g = O( ). It is clear that the polynomial factor leaves the NP-hardness of the latter program intact. We use the same notation as in the proof of Lemma 79. Let f S be a Yes instance of Distribution Decomposability , and define S := {s + O( ) : s ∈ S}. From equation (15) it immediately follows that f S is a Yes instance of Distribution Decomposability g( ) , where we allow g = O( ). We have hence shown that there exists an O( ) ball around each Yes instance that solely contains Yes instances.
A similar argument holds for the No instances. It is clear that these cases can be answered using Weak Distribution Decomposability , and the claim follows. 2

Complete decomposability
Another interesting question to ask is for the complete decomposition of a finite distribution D into a sum of indecomposable distributions. We argue that this decomposition is not unique.
Proposition 81. There exists a family of finite distributions (D n ) n∈N with probability mass functions p n (k) : max supp p n (k) = 4n and such that, for each D n , there are at least n! distinct decompositions into indecomposable distributions.
Proof. We explicitly construct the family (D n ) n∈N . Let n ∈ N. We will define a set of irreducible quadratic polynomials {p k , n k for k = 1, . . . , n} such that n k are not positive, but p k n l are positive quartics ∀k, l-and thus define valid probability distributions. Since R[x] is a unique factorization domain the claim then follows.
Following the findings in the proof of Lemma 70, it is in fact enough to construct a set {a k , b k : 0 < |a k | < 2, −2 < b k < 0 for k = 1, . . . , n} ⊂ R 2n and such that a k + b l > 0 ∀k, l-then let p k := 1 + a k x + x 2 , n k := 1 + b k x + x 2 . It is straightforward to verify that e.g. Corollary 83. R is not a unique factorization domain.
Proposition 81 and Remark 82 show that an exponential number of complete decompositions-all of which have different distributions-do not give any further insight into the distribution of interest-indeed, as the number of positive indecomposable factors is not even unique. Asking for a non-maximal decomposition into indecomposable terms does therefore not answer more than whether the distribution is decomposable at all.
Indeed, the question whether one can decompose a distribution into indecomposable parts can be trivially answered with Yes, but if we include the condition that the factors have to be non-trivial, or for decomposability into a certain number of terms-say N ≥ 2 or the maximum number of terms-the problem is also obviously NP-hard by the previous results.
In short, by Theorem 77, we immediately obtain the following result.
Corollary 84. Let D be a finite discrete distribution. Deciding whether one can write D as any nontrivial sum of irreducible distributions is NP-hard.

Continuous distributions
Analogous to our discussion in section 3.4.5, the exact and variants of the decomposability question are computationally ill-posed. We again point out that answering the weak membership version is trivial, since the set of indecomposable distributions is dense, as the following proposition shows. Proof. We first extend Lemma 61, and again take f ∈ C + c,b : supp f ⊂ A ∪ B. While not automatically true that r(x), s(x) = 0 ∀x < 0, we can assume this by shifting r and s symmetrically. We also assume inf supp f = 0, and hence inf supp r = inf supp s = 0-see Lemma 44 for details. Since The integration domain difference is derived analogously, and can be seen in an example in Fig. 4. We again regard the two cases separately. Let x ∈ A. Assume ∃y ∈ (M − m, M ) such that r(x − y )s(y ) > 0. Then s(y ) > 0, contradiction. Now fix x ∈ (m, M ), and assume ∃y ∈ (0, x − m) : r(x − y )s(y ) > 0. Since x − y > x − x + m = m, r(x − y ) > 0 yields another contradiction.
The rest of the proof goes through analogously. 2 Corollary 86. Let > 0. Let X be a continuous random variable with pmf p X (k). Then there exists an indecomposable random variable Y with pmf p Y (k), such that p X − p Y < .

Conclusion
In section 2, we have shown that the question of existence of a stochastic root for a given stochastic matrix is in general at least as hard as answering 1-in-3sat, i.e. it is NP-hard. By Corollary 27, this NP-hardness result also extends to Nonnegative and Doubly Stochastic Divisibility, which proves Theorem 1. A similar reduction goes through for cptp Divisibility in Corollary 24, proving NP-hardness of the question of existence of a cptp root for a given cptp map.
In section 3, we have shown that-in contrast to cptp and stochastic matrix divisibility-distribution divisibility is in P, proving Theorem 4. On the other hand, if we relax divisibility to the more general decomposability problem, it becomes NP-hard as shown in Theorem 6. We have also extended these results to weak membership formulations in Theorem 5 and 7-i.e. where we only require a solution to within in the appropriate metric-showing that all the complexity results are robust to perturbation.
Finally, in section 3.4.5 and 3.5.8, we point out that for continuous distributionswhere the only computationally the only meaningful formulations are the weak membership problems or closely related variants-questions of divisibility and decomposability become computationally trivial, as the nondivisible and indecomposable distributions independently form dense sets.
As containment in NP for all of the NP-hard problems is easy to show (Lemma 21 and 68), these problems are also NP-complete. Thus our results imply that, apart for the distribution divisibility problem which is efficiently solvable, all other divisibility problems for maps and distributions are equivalent to the famous P = NP conjecture, in the following precise sense: A polynomial-time algorithm for answering any one of these questions-(Doubly) Stochastic, Nonnegative or cptp Divisibility, or either of the Decomposability variants-would prove P = NP. Conversely, solving P = NP would imply that there exists a polynomial-time algorithm to solve all of these Divisibility problems.
For Even Subset Sum, we generalize Corollary 89 to the following scaling property. Proof. Let S be the multiset of a Subset Sum instance (S, l), where we assume without loss of generality that all S s ≥ 0. Now first assume Σ S = 0. In that case the claim follows immediately, since the problems are identical.
If now Subset Sum(S , 0) = True, we know that there exists T S : | t∈T t − s∈S \T s| = 0. Now assume T contains both copies of −Σ S /2. Then clearly Finally observe the following extension of Lemma 96.