/******************************************************************************
 * This program computes the roots of cubic equations
 * of the form Ax^3 + Bx^2 + Cx + D = 0
 * 
 * Copyright © 2020 Richard Lesh.  All rights reserved.
 *****************************************************************************/

#undef NDEBUG

#include <clocale>
#include <cmath>
#include <codecvt>
#include <iostream>
#include <string>

std::locale utf8loc(std::locale(), new std::codecvt_utf8<wchar_t>);
using namespace std;

static double const A = 1.0;
static double const B = -6.0;
static double const C = 11.0;
static double const D = -6.0;

int main(int argc, char **argv) {
    setlocale(LC_ALL, "en_US.UTF-8");
    wcout.imbue(utf8loc);
    wcin.imbue(utf8loc);

/******************************************************************************
 * General cubic can be changed to a depressed cubic 
 * t^3 + pt + q = 0
 * by substituting
 * x = t - B / 3A
 *****************************************************************************/

    double p = (3. * A * C - B * B) / (3. * A * A);
    double q = (2. * B * B * B - 9. * A * B * C + 27. * A * A * D) / (27. * A * A * A);

/******************************************************************************
 * Roots of the depressed cubic are...
 * tk = a cos(b - 2πk/3) for k = 0, 1, 2
 * with a = 2 sqrt(-p/3)
 * b = (1/3) arccos(3q/2p * sqrt(-3/p))
 *****************************************************************************/
    double a = 2. * sqrt(-p / 3.);
    double b = (1. / 3.) * acos(3. * q / 2. / p * sqrt(-3. / p));

    double t1 = a * cos(b);
    double t2 = a * cos(b - 2. * M_PI / 3.);
    double t3 = a * cos(b - 4. * M_PI / 3.);
    double root1 = t1 - B / (3. * A);
    double root2 = t2 - B / (3. * A);
    double root3 = t3 - B / (3. * A);
    wcout << L"root #1: " << root1 << endl;
    wcout << L"root #2: " << root2 << endl;
    wcout << L"root #3: " << root3 << endl;
    return 0;
}