/*
 *   Program to calculate roots
 *
 *   Copyright (c) 1988-1995 Shamus Software Ltd.
 *   Modified as a test driver by D. Corbit.
 */

#include <iostream.h>
#include <flash.h>

miracl          precision = 50;

Flash           SqrtAB03(Flash x, Flash c);
Flash           CbrtAB03(Flash x, Flash c);
int             main()
{                               /* Find roots */
    Flash           x,
                    y;
    Flash           t;
    Flash           crappy_estimate;
    const Flash     OnePointFive = "1.5";
    const Flash     One = "1.";
    const Flash     Three = "3.";
    const Flash     Third = One / Three;
    int             n;
    RPOINT = ON;
    cout << "enter number\n";
    cin >> x;
    cout << "to the power of 1/2\n";
    y = sqrt(x);
    cout << y;
    x = y * y;
    cout << "\nto the power of 2 = \n";
    cout << x << endl;
    cout << "Guess is too high\n";
    crappy_estimate = y * OnePointFive;
    t = SqrtAB03(crappy_estimate, x);
    cout << "Estimate for t (Pass 1)" << t << endl;
    t = SqrtAB03(t, x);
    cout << "Estimate for t (Pass 2)" << t << endl;
    t = SqrtAB03(t, x);
    cout << "Estimate for t (Pass 3)" << t << endl;

    cout << "Guess is too low\n";
    crappy_estimate = y / OnePointFive;
    t = SqrtAB03(crappy_estimate, x);
    cout << "Estimate for t (Pass 1)" << t << endl;
    t = SqrtAB03(t, x);
    cout << "Estimate for t (Pass 2)" << t << endl;
    t = SqrtAB03(t, x);
    cout << "Estimate for t (Pass 3)" << t << endl;

    cout << "enter number\n";
    cin >> x;
    cout << "to the power of 1/3\n";
    y = pow(x, Third);
    cout << y;
    x = y * y * y;
    cout << "\nto the power of 3 = \n";
    cout << x << endl;
    cout << "Guess is too high\n";
    crappy_estimate = y * OnePointFive;
    t = CbrtAB03(crappy_estimate, x);
    cout << "Estimate for t (Pass 1)" << t << endl;
    t = CbrtAB03(t, x);
    cout << "Estimate for t (Pass 2)" << t << endl;
    t = CbrtAB03(t, x);
    cout << "Estimate for t (Pass 3)" << t << endl;

    cout << "Guess is too low\n";
    crappy_estimate = y / OnePointFive;
    t = CbrtAB03(crappy_estimate, x);
    cout << "Estimate for t (Pass 1)" << t << endl;
    t = CbrtAB03(t, x);
    cout << "Estimate for t (Pass 2)" << t << endl;
    t = CbrtAB03(t, x);
    cout << "Estimate for t (Pass 3)" << t << endl;

    return 0;
}

/********************************************************************/
/* This is Schroder's Method(B, Lamda=0, Omega = 3) for x^2 - c = 0 */
/********************************************************************/
Flash           SqrtAB03(Flash x, Flash c)
{
    static const Flash Four = "4.";
    static const Flash Half = ".5";
    static const Flash Six = "6.";
    Flash x2 = x * x;
    Flash ans0 = Four * x * c * (x2 + c) / (x2 * (x2 + Six * c) + c * c);
    return (ans0 + c / ans0) * Half;
}
Flash           CbrtAB03(Flash x, Flash c)
{
    static const Flash Four = "4.";
    static const Flash Half = ".5";
    static const Flash Sixteen = "16.";
    static const Flash Nineteen = "19.";
    static const Flash Ten = "10.";
    static const Flash One = "1.";

    Flash x2 = x * x;
    Flash x6 = x2 * x2 * x2;
    Flash c2 = c * c;
    Flash x3c = x2 * x * c;
    Flash x3c16 = x3c * Sixteen;
    Flash tenx6 = Ten * x6;
    Flash temp = x3c16 + tenx6 + c2;
    Flash ans0 = x * (Four * x6 + Nineteen * x3c + Four * c2) / (temp);
    Flash ans1 = One / x2 * c * (temp) / (x3c16 + Ten * c2 + x6);
    return (ans0 + ans1) * Half;
}