C, finding root, Newton Method

newton_method.c

// implementation of Newton Method (Simple Iteration)
// 2017, Copyleft lion137
// Znajduje pierwiastek wielomianu Metodą Newtona: 
// https://en.wikipedia.org/wiki/Newton%27s_method czyli:
// x0 = wartość poczatkowa (zgadywana - nie powinna za daleko odbiegać od pierwiastka - tu jest problem tej metody,
// ale wszystkie mają podobne problemy, patrz Wikipedia: https://en.wikipedia.org/wiki/Root-finding_algorithm , chociaz 
// 100 dalej nie jest za daleko od pierwiastka -0.169... patrz przykład w main)
// i iteracja:
// x1 = x0 - f(x0) / f'(x0)
// x0 = x1 -> podstawienie

// kompilowane:
// konmpilator: GNU g++ 5.4.0
// kompilacja: g++ -Wall -std=c++11 <nazwa_pliku>

#include <iostream>
#include <cmath> // Potrzebna funkcja podnoszenia do potegi (pow)
#include <limits> // NaN
#include <random>

using namespace std;

//deklaracje funkcji:
double root(double ar[], int size, int degree, double init_guess, double accuracy);
double apply(double tab[], double x, int size, int degree);



int main(){
	
	int degree = 3;
	int size = 4;
	double ar[size] = {5, -7, 3, 0.5}; // opisuje wielomian: 4x^3 + x^2 + 3x + 0.5 
	double accuracy = 0.01;
	double guess = -1;
	cout << root(ar, size, degree, guess, accuracy)<<endl; // -> -0.16975:
	//https://www.wolframalpha.com/input/?i=4+*+x%5E3+%2B+x%5E2+%2B+3*x+%2B+0.5+%3D+0
	return 0;
	}




//definicje funkcji:
double root(double ar[], int size, int degree, double init_guess, double accuracy){
	double der[size];  // tablica do wpisania wielomianu pochodnych
	int degree_der = degree - 1; // stopień pochodnej wejścia - zawsze o 1 mniejszy
	int ind_der = size - 1;  // ustawia na stopień najwyższej potęgi x
	
	// policzenie pochodnych:
	for (int i = 0; i < size; ++i){ 
		der[i] = ar[i] * ind_der;  // jak we wzorze na pochodną
		ind_der -= 1;  // dekrementujemy - czyli schodzimy do niższej potęgi w wielomianie
		}
		
	// Iteracyjne liczenie pierwiastka:
	double x0 = init_guess; //punkt startowy
	double x1; // kolejna iteracja
	
	for (int i = 0;i < 300 ; ++i){ // pętla ograniczona na wypadek rozbiezności (wtedy zwróci NaN)
		//krok iteracyjny:
		x1 = x0 - (apply(ar, x0, size, degree) / apply(der, x0, size, degree_der));
		//sprawdzenie dokładności:
		if ( (x0 - x1 <= accuracy) && (x0 - x1 >= -accuracy)) return x1; 
		x0 = x1; //podstawienie iteracyjne
		}
	
	return numeric_limits<double>::quiet_NaN(); // w razie rozbieżności lub braku zadanej dokładności
	// zwraca Nan (not a number)
	}

double apply(double tab[], double x, int size, int degree){
	// Wylicza wartość wielomianu w punkcie x
	double s = 0;  // wartosć początkowa do ewaluacji, oczywiście zero
	
	for (int i = 0; i < size; ++i){
		if (tab[i] == 0) s += 0;
		s += tab[i] * (pow(x, degree)); // podstawia x do kolejnego członu wielomianu i sumuje
		degree -= 1;  // jak w funkcji root
		}
		
	return s;
	}