A decision procedure for the Barendregt-Coppo-Dezani subtype relation on intersection types (“BCD subtyping”) is presented and formally verified in Coq. Types are extended with unary, covariant, distributing, preordered type constructors and binary products. A quadratic upper bound on the algorithm runtime is established. The formalization can be compiled to executable OCaml or Haskell code using the extraction mechanism of Coq.