#include <iostream>
#include <math.h>

using namespace std;
typedef double(*userFuncPtr)(double);

/* FUNKCJE UZYTKOWNIKA */
double user_function1(double x){
       return pow(x-2.0f,2) - 25.0f;
}

double user_function2(double x){
       return (x-5)*(x-3)*x*(x+10);
}

double user_function3(double x){
       return pow(x, 4) - 625.0f;
}


double derivateIn(double x0, double h, int nNable, userFuncPtr func, int n);
void metoda_siecznych(double a, double b, double eps, userFuncPtr func, double *out);
void metoda_stycznych(double a, double b, double eps, userFuncPtr func, double *out);
void metoda_bernoulliego(int n, const double * a_i_vec, double eps, userFuncPtr func);

int main(){
    const double A = -10.0f;
    const double B = -1.0f;
    const double EPS = 0.0001;
    double out;

    cout<<"Metoda siecznych:"<<endl;
    metoda_siecznych(A, B, EPS, &user_function3, &out);

    cout<<"\n\nMetoda stycznych:"<<endl;
    metoda_stycznych(A, B, EPS, &user_function3, &out);


    double a_i_vec[] = {1.0f, -18.0f, -65.f, 150.0f, 0.0f};
    cout<<"\n\nMetoda Bernoulliego:"<<endl;
    metoda_bernoulliego(4, a_i_vec, EPS, &user_function2);

    system("pause");
    return 0;
}

void metoda_bernoulliego(int n, const double *a_i_vec, double eps, userFuncPtr func){

     double y_i;
     double *y_b = new double[n*1000];
     memset(y_b, 1, sizeof(double)*n);

     double fac = 0;
     /*
     do{
        for(int i = 0; i<(n-1); i++)
                //fac += a_i_vec[i+1]

        y_i = (-1.0f/a_i_vec[0]);
     }while( (1.0f) > eps);
     */
     //horner


     delete [] y_b;
}

void metoda_stycznych(double a, double b, double eps, userFuncPtr func, double *out){
     if(func(a)*func(b) > 0.0f){
                        cout<<"W tym przedziale nie ma pierwiastkow!\n";
                        return;
                        }

     double deriv2A, deriv2B; // 2 pochodne w a i b
     double x = 0;
     int counter = 0;

     deriv2A = derivateIn(a, 0.001, 3, func, 2);
     deriv2B = derivateIn(b, 0.001, 3, func, 2);

     cout<<"Drugie pochodne:"<<deriv2A<<"  "<<deriv2B<<endl;

     double x_i;

     if( (func(a)*deriv2A) > 0.0f ){
         x_i = a;
     }else if( (func(b)*deriv2B) > 0.0f ){
         x_i = b;
     }

     do{
        x_i = x_i - func(x_i)/(derivateIn(x_i, 0.001, 3, func, 1));
        cout<<"iteracja["<<counter<<"]: x_i = "<<x_i<<"    wartosc func:"<<func(x_i)<<endl;
        counter++;
     }while( fabs(func(x_i)) > eps);

}



 /* ------------ NIELONIOWE ROWNANIA ------------*/



void metoda_siecznych(double a, double b, double eps, userFuncPtr func, double *out){
     if(func(a)*func(b) > 0.0f){
                        cout<<"W tym przedziale nie ma pierwiastkow!\n";
                        return;
                        }

     double deriv2A, deriv2B; // 2 pochodne w a i b
     double x = 0;
     int counter = 0;

     deriv2A = derivateIn(a, 0.001, 3, func, 2);
     deriv2B = derivateIn(b, 0.001, 3, func, 2);

     cout<<"Drugie pochodne:"<<deriv2A<<"  "<<deriv2B<<endl;

     double x_const, x_i;

     if( (func(a)*deriv2A) > 0.0f ){
         x_const = a;
         x_i = b;
     }else if( (func(b)*deriv2B) > 0.0f ){
         x_const = b;
         x_i = a;
     }

     do{
        x_i = x_i - ((func(x_i)*(x_const-x_i))/(func(x_const)-func(x_i)));
        cout<<"iteracja["<<counter<<"]: x_i = "<<x_i<<"    wartosc func:"<<func(x_i)<<endl;
        counter++;
     }while( fabs(func(x_i)) > eps);

}






 /* ------------ ROZNICZKOWANIE ------------*/
/* ------------ CORE FUNCTIONS ------------*/

void initVectors(double *vec1, double *vec2, int size){
        memset(vec2, 0x00, size*sizeof(double));
        memset(vec1, 0x00, size*sizeof(double));

        for (size_t i = 0; i < size; i++)
                        vec1[i] = 1.0f / (i + 1);

        for (size_t i = 0; i < size; i++)
                for (size_t j = 0; j < size; j++){
                        if ((i+j) < size)
                                vec2[i+j] += vec1[i] * vec1[j];
                }

}

void substractionVec(double x0, double h, userFuncPtr func, double *out_vec, int size){
        memset(out_vec, 0x00, size*sizeof(double));

        double *vecOne, *vecTwo, *temp;

        vecOne = new double[size+1];
        vecTwo = new double[size+1];

        //tworzmy wektor wartosciami funkcji
        for (size_t i = 0; i <= size; i++)
                vecOne[i] = func(x0 - (i*h));


        for(size_t i = 0; i < (size); i++)
        {
                for (size_t j = 0; j < (size-i); j++)
                {
                        vecTwo[j] = vecOne[j] - vecOne[j + 1];
                }

                out_vec[i] = vecTwo[size-(i+1)];

                temp = vecTwo;
                vecTwo = vecOne;
                vecOne = temp;
        }

        delete[] vecOne;
        delete[] vecTwo;
}


/*Funkcja obliczajšca pochodnš funkcji
        parametry:
        x0 - punkt w który obliczana jest pochodna
        h - krok różniczkowania
        nNable - dokładnoœć różniczkowania (iloœć operatorów nabla dla n=1)
        func - badana funkcja
        n - rzšd pochodnej
*/

double derivateIn(double x0, double h, int nNable, userFuncPtr func, int n){
        //pochodne wyższego rzędu niż 2 nie sš zaimplementowane
        if (n > 2)
                return 0.0f;

        double *vec1 = new double[nNable];
        double *vec2 = new double[nNable];
        double *subVec = new double[nNable];
        double outPut = 0;

        initVectors(vec1, vec2, nNable);
        substractionVec(x0, h, func, subVec, nNable);


        if (n == 1){
                for (size_t i = 0; i < nNable; i++)
                {
                        outPut += vec1[i] * subVec[i];
                }

                outPut /= h;
        }

        if (n == 2){
                for (size_t i = 0; i < (nNable-1); i++)
                {
                        outPut += vec2[i] * subVec[i+1];
                }
                outPut /= pow(h,2);
        }


        delete [] vec1;
        delete [] vec2;
    delete [] subVec;

        return outPut;
}