ABSTRACT We consider a class of collocation matrices associated with the Newton basis of the space of polynomials of degree at most , evaluated at a set of nodes. In the most general setting, we allow of these nodes to either coincide with or differ from those defining the Newton basis. We establish necessary and sufficient conditions for to be strictly totally positive (STP), based solely on an appropriate ordering of the nodes. We derive explicit formulas for the bidiagonal decomposition of and present an efficient algorithm for computing it to high relative accuracy (HRA). This factorization enables accurate and fast computation of solutions to linear systems with coefficient matrix , as well as the evaluation of its eigenvalues, singular values, and inverse, all to HRA. Furthermore, we apply our framework to the least squares approximation problem in the Newton basis. We show that the STP property of guarantees both uniqueness of the solution and, under mild conditions, its computation to an accuracy comparable to HRA. Extensive numerical experiments confirm the theoretical properties and demonstrate the superior accuracy and performance of the proposed algorithms when compared to standard numerical methods, especially in the presence of severe ill‐conditioning.