using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace QuadFormula
{
    class Program
    {

        static void Main(string[] args)
        {
            var qf = new QF(1.5, 3.3, 3, true);
            var result = qf.Calculate(P, F);
            Console.WriteLine(result);

        }
        public static double P(double x)
        {
            double a = 1.5;
            double b = 3.3;
            double alfa = (double)1 / 3;
            double beta = 0;
            var z = Math.Pow((x - a), -alfa) * Math.Pow((b - x), -beta);
            if (double.IsInfinity(z))
                return 0;
            return z;
        }
        public static double F(double x)
        {
            return 2 * Math.Cos(2.5 * x) * Math.Exp(x / 3) + 4 * Math.Sin(3.5 * x) * Math.Exp(-3 * x) + x;
        }

    }
    public class QF
    {
        public delegate double Function(double x);
        public bool Debug;
        public double[] Nodes;
        public double[] Moments;
        public double[,] SystemA;
        public double A;
        public double B;
        public double[] SystemResultA;
        public int N;
        public QF(double a, double b, int n, bool debug)
        {
            this.A = a;
            this.B = b;
            this.N = n;
            this.Debug = debug;
        }
        public double Calculate(Function P, Function F)
        {
            if (this.Debug)
                Console.Write("Генерируем узлы");
            this.GenerateNodes();

            if (this.Debug)
                Console.WriteLine("Считаем моменты...");
            this.CalculateMoments(P);
            if (this.Debug)
                Console.WriteLine("Составляем систему...");
            this.GenerateSystem();

            if (this.Debug)
                Console.WriteLine("Решаем систему методом Гаусса...");
            this.SolveSystem();
            if (this.Debug)
                Console.WriteLine("Считаем интеграл по квадратурной формуле...");
            return CalculateIntegral(F);
        }
        private void GenerateNodes()
        {
            Nodes = new double[N];
            var dx = (B - A) / (N - 1);
            var node = A;
            for (int i = 0; i < N; i++)
            {
                Nodes[i] = node;
                node += dx;
            }
            if (this.Debug)
                Console.WriteLine(" с шагом {0}, количество {1}", dx, N);
        }
        private void CalculateMoments(Function P)
        {
            Moments = new double[N];
            for (int i = 0; i < N; i++)
            {
                Moments[i] = RombergIntegration(x => P(x) * Math.Pow(x, i), A, B, 3000, 1E-15);
                if (this.Debug)
                    Console.WriteLine("\t {0}-момент = {1}", i, Moments[i]);
            }
        }
        private void GenerateSystem()
        {
            SystemA = new double[N, N];
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < N; j++)
                {
                    SystemA[i,j] =Math.Pow(Nodes[j], i);
                    Console.Write("{0,10}", SystemA[i, j]);
                }
                if (this.Debug)
                    Console.WriteLine();
            }
        }
        private void SolveSystem()
        {
            SystemResultA = GaussSolve(SystemA, Moments);
        }
        private double CalculateIntegral(Function f)
        {
            double result = 0;
            for (int i = 0; i < N; i++)
            {
                result += SystemResultA[i] * f(Nodes[i]);
            }
            return result;
        }
        private static double RombergIntegration(Function f, double a, double b, int maxIterations, double tolerance)
        {
            int n = 2;
            double h = b - a;
            double sum = 0.0;
            int j = 0;
            double[,] R = new double[maxIterations, maxIterations];

            R[1, 1] = h * (f(a) + f(b)) / 2.0;
            h = h / 2;
            R[2, 1] = R[1, 1] / 2 + h * f(a + h);
            R[2, 2] = (4 * R[2, 1] - R[1, 1]) / 3;

            for (j = 3; j <= maxIterations; j++)
            {
                n = 2 * n;
                h = h / 2;
                sum = 0.0;
                for (int k = 1; k <= n; k += 2)
                {
                    sum += f(a + k * h);
                }
                R[j, 1] = R[j - 1, 1] / 2 + h * sum;
                double factor = 4.0;
                for (int k = 2; k <= j; k++)
                {
                    R[j, k] = (factor * R[j, k - 1] - R[j - 1, k - 1]) / (factor - 1);
                    factor = factor * 4.0;
                }
                if (Math.Abs(R[j, j] - R[j, j - 1]) < tolerance * Math.Abs(R[j, j]))
                {
                    sum = R[j, j];
                    return sum;
                }
            }
            sum = R[n, n];
            return sum;
        }
        private static double[,] ResizeArray(double[,] original, int rows, int cols)
        {
            var newArray = new double[rows, cols];
            int minRows = Math.Min(rows, original.GetLength(0));
            int minCols = Math.Min(cols, original.GetLength(1));
            for (int i = 0; i < minRows; i++)
                for (int j = 0; j < minCols; j++)
                    newArray[i, j] = original[i, j];
            return newArray;
        }
        private static double[] GaussSolve(double[,] A, double[] B)
        {
            // размерность квадратной матрицы A
            int n = A.GetUpperBound(0) + 1;
            var M = ResizeArray(A, n, n + 1);
            for (int i = 0; i < n; i++)
            {
                M[i, n] = B[i];
            }
            for (int col = 0; col + 1 < n; col++) if (M[col, col] == 0)
                // проверка на нулевые коэффициенты
                {
                    // поиск ненулевого коэффициента
                    int swapRow = col + 1;
                    for (; swapRow < n; swapRow++) if (M[swapRow, col] != 0) break;

                    if (M[swapRow, col] != 0) // найден ненулевой коэффициент
                    {
                        // смена значений местами
                        double[] tmp = new double[n + 1];
                        for (int i = 0; i < n + 1; i++)
                        {
                            tmp[i] = M[swapRow, i];
                            M[swapRow, i] = M[col, i];
                            M[col, i] = tmp[i];
                        }
                    }
                    else return null;
                }
            // исключение
            for (int sourceRow = 0; sourceRow + 1 < n; sourceRow++)
            {
                for (int destRow = sourceRow + 1; destRow < n; destRow++)
                {
                    double df = M[sourceRow, sourceRow];
                    double sf = M[destRow, sourceRow];
                    for (int i = 0; i < n + 1; i++)
                        M[destRow, i] = M[destRow, i] * df - M[sourceRow, i] * sf;
                }
            }
            // обратный ход гаусса
            for (int row = n - 1; row >= 0; row--)
            {
                double f = M[row, row];
                if (f == 0) return null;

                for (int i = 0; i < n + 1; i++) M[row, i] /= f;
                for (int destRow = 0; destRow < row; destRow++)
                { M[destRow, n] -= M[destRow, row] * M[row, n]; M[destRow, row] = 0; }
            }
            // формирование ответа
            var x = new double[n];
            for (int i = 0; i < n; i++)
            {
                x[i] = M[i, n];
            }
            return x;
        }

    }
}