Newton-Raphson/newton.cpp at master · Salil03/Newton-Raphson · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <bits/stdc++.h>
using namespace std;
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

void derivative(double polynomial[], double slope[], int deg)
{
    for (int i = 0; i < deg; i++)
    {
        slope[i] = polynomial[i + 1] * (i + 1);
    }
}

complex<double> poly_value(double polynomial[], int deg, complex<double> input)
{
    complex<double> value(0, 0);
    for (int i = 0; i <= deg; i++)
    {
        value += pow(input, i) * polynomial[i];
    }
    return value;
}

double max_radius(double polynomial[], int deg)
{
    double radius = 0;
    for (int i = 0; i < deg; i++)
    {
        radius += abs(polynomial[i] / polynomial[deg]);
    }
    return radius;
}

int main()
{
    srand((unsigned int)time(NULL));
    cout << "Enter degree of polynomial: \n";
    int deg;
    cin >> deg;
    double polynomial[deg + 1], slope[deg] = {0};
    cout << "Enter coefficients from largest power to smallest: \n";
    for (int i = deg; i >= 0; i--)
    {
        cin >> polynomial[i];
    }
    derivative(polynomial, slope, deg);
    double d = (double)deg;
    int circles = max(ceil(0.26632 * log2(d)), 2.0);
    int points = max(ceil(8.32547 * d * log2(d)), 2.0);
    double start_radius = max_radius(polynomial, deg);

    ofstream rad;
    rad.open("radius.txt");
    rad << start_radius;
    rad.close();

    complex<double> inputs[circles * points];
    for (int v = 1; v <= circles; v++)
    {
        for (int j = 0; j < points; j++)
        {
            inputs[(v - 1) * points + j] = polar(start_radius * (1 + sqrt(2)) * pow((d - 1) / d, (2 * v - 1) / (4 * circles)), 2 * M_PI * j / points);
        }
    }

    ofstream start_points;
    start_points.open("start_points.txt");

    vector<complex<double>> roots;
    for (int start = 0; start < circles * points; start++)
    {
        double tolerance = 0.0000001;
        int cycle_lim = ceil(d * log2((1 + sqrt(2)) / tolerance));
        int cycle_curr = 0;
        complex<double> input = inputs[start];
        complex<double> output = poly_value(polynomial, deg, input);
        start_points << real(input) << " " << imag(input) << "\n";
        while (cycle_curr < cycle_lim && abs(output) > tolerance)
        {
            input = input - output / poly_value(slope, deg - 1, input);
            output = poly_value(polynomial, deg, input);
            cycle_curr++;
        }
        bool flag = true;
        for (int index = 0; index < (int)roots.size(); index++)
        {
            if (abs(roots[index] - input) <= 0.001)
            {
                if (abs(poly_value(polynomial, deg, input)) < abs(poly_value(polynomial, deg, roots[index])))
                {
                    roots[index] = input;
                }

                flag = false;
                break;
            }
        }
        if (flag && abs(poly_value(polynomial, deg, input)) < tolerance)
        {
            roots.push_back(input);
        }
    }
    start_points.close();
    ofstream sols;
    sols.open("roots.txt");
    cout << "Roots of this polynomial are: \n";
    for (auto root : roots)
    {
        cout << fixed << setprecision(10) << root << "\n";
        sols << fixed << setprecision(10) << real(root) << " " << imag(root) << "\n";
    }
    sols.close();
}