#include <cassert>#include <complex>#include <cstddef>#include <iostream>#include <numbers>#include <valarray> using CD=std::complex<double>;using VA=std::valarray<CD>; // return all n complex roots out of a given complex number xVA root(CD x,unsigned n){constdouble mag=std::pow(std::abs(x),1.0/ n);constdouble step=2.0*std::numbers::pi/ n;double phase=std::arg(x)/ n;    VA v(n);for(std::size_t i{}; i!= n;++i, phase+= step)        v[i]=std::polar(mag, phase);return v;} // return n complex roots of each element in v; in the output valarray first// goes the sequence of all n roots of v[0], then all n roots of v[1], etc.VA root(VA v,unsigned n){    VA o(v.size()* n);    VA t(n);for(std::size_t i=0; i!= v.size();++i){        t= root(v[i], n);for(unsigned j=0; j!= n;++j)            o[n* i+ j]= t[j];}return o;} // floating-point numbers comparator that tolerates given rounding errorinlinebool is_equ(CD x, CD y,double tolerance=0.000'000'001){return std::abs(std::abs(x)- std::abs(y))< tolerance;} int main(){// input coefficients for polynomial x³ + p·x + qconst VA p{1,2,3,4,5,6,7,8};const VA q{1,2,3,4,5,6,7,8}; // the solverconst VA d=std::sqrt(std::pow(q/2,2)+std::pow(p/3,3));const VA u= root(-q/2+ d,3);const VA n= root(-q/2- d,3); // allocate memory for roots: 3 * number of input cubic polynomials    VA x[3];for(std::size_t t=0; t!=3;++t)        x[t].resize(p.size()); auto is_proper_root=[](CD a, CD b, CD p){return is_equ(a* b+ p/3.0,0.0);}; // sieve out 6 out of 9 generated roots, leaving only 3 proper roots (per polynomial)for(std::size_t i=0; i!= p.size();++i)for(std::size_t j=0, r=0; j!=3;++j)for(std::size_t k=0; k!=3;++k)if(is_proper_root(u[3* i+ j], n[3* i+ k], p[i]))                    x[r++][i]= u[3* i+ j]+ n[3* i+ k]; std::cout<<"Depressed cubic equation:   Root 1:\t\t Root 2:\t\t Root 3:\n";for(std::size_t i=0; i!= p.size();++i){std::cout<<"x³ + "<< p[i]<<"·x + "<< q[i]<<" = 0  "<<std::fixed<< x[0][i]<<"  "<< x[1][i]<<"  "<< x[2][i]<<std::defaultfloat<<'\n'; assert(is_equ(std::pow(x[0][i],3)+ x[0][i]* p[i]+ q[i],0.0));assert(is_equ(std::pow(x[1][i],3)+ x[1][i]* p[i]+ q[i],0.0));assert(is_equ(std::pow(x[2][i],3)+ x[2][i]* p[i]+ q[i],0.0));}}

Depressed cubic equation:   Root 1:              Root 2:                 Root 3:x³ + (1,0)·x + (1,0) = 0  (-0.682328,0.000000)  (0.341164,1.161541)  (0.341164,-1.161541)x³ + (2,0)·x + (2,0) = 0  (-0.770917,0.000000)  (0.385458,1.563885)  (0.385458,-1.563885)x³ + (3,0)·x + (3,0) = 0  (-0.817732,0.000000)  (0.408866,1.871233)  (0.408866,-1.871233)x³ + (4,0)·x + (4,0) = 0  (-0.847708,0.000000)  (0.423854,2.130483)  (0.423854,-2.130483)x³ + (5,0)·x + (5,0) = 0  (-0.868830,0.000000)  (0.434415,2.359269)  (0.434415,-2.359269)x³ + (6,0)·x + (6,0) = 0  (-0.884622,0.000000)  (0.442311,2.566499)  (0.442311,-2.566499)x³ + (7,0)·x + (7,0) = 0  (-0.896922,0.000000)  (0.448461,2.757418)  (0.448461,-2.757418)x³ + (8,0)·x + (8,0) = 0  (-0.906795,0.000000)  (0.453398,2.935423)  (0.453398,-2.935423)