APR2

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

/*
Metody Numeryczne zadanie APR2
Przygotował: Dawid Grzeszuk
Pomoc: Radosław Szkorla
Kierunek: Informatyka 3 semestr
Środowisko: Visual Studio 2017 RC
Język programowania: C#
*/

namespace APR2
{
    internal static class Details
    {
        static public double eps = 1e-5;
        static public int precision = 10;
    }
    internal static class DisplayMethods
    {
        public static void Show(double[,] tab)
        {
            for (int x = 0; x < tab.GetLength(0); x++)
            {
                for (int y = 0; y < tab.GetLength(1); y++)
                {
                    if (y < tab.GetLength(1) - 1)
                        Console.Write($"{Math.Round(tab[x, y], Details.precision)} ");
                    else
                        Console.Write($"| {Math.Round(tab[x, y], Details.precision)} ");
                }
                Console.WriteLine();
            }
        }
    }

    internal static class MathFunctions
    {
        #region basic functions
        static double f(double x) => x;
        static double g(double x, int j) => f(x) * Math.Pow(x, j);
        public static double h(int k, double a, double b) => (b - a) / Math.Pow(2, k);
        public static double x(long i, int k, double a, double b) => a + i * h(k, a, b);
        public static double y(long i, int k, double a, double b, int j) => g(x(i, k, a, b), j);
        public static double delta(double a, double b, int i) => (Math.Pow(b, i) - Math.Pow(a, i));
        #endregion

        #region additional simpson functions
        public static double P(double[,] M, double x)
        {
            double wynik = 0;
            for (int i = 0; i < M.GetLength(0); i++)
            {
                wynik += M[i, M.GetLength(1) - 1] * Math.Pow(x, i);
            }
            return wynik;
        }
        public static double[,] GlobalM;
        public static double f2(double x) => Math.Pow(f(x) - P(GlobalM, x), 2);
        public static double h2(int k, double a, double b) => (b - a) / Math.Pow(2, k);
        public static double x2(long i, int k, double a, double b) => a + i * h(k, a, b);
        public static double y2(long i, int k, double a, double b, int j) => f(x(i, k, a, b));
        #endregion

        #region gauss method
        public static bool GaussMethod(int n, ref double[,] tab)
        {
            int r;
            double max;
            for (int k = 0; k < n; k++)
            {
                max = tab[k, k];
                r = k;
                for (int x = k; x < n; x++)
                {
                    if (Math.Abs(tab[x, k]) > Math.Abs(max))
                    {
                        max = tab[x, k];
                        r = x;
                    }
                }
                if (Math.Abs(max) < Details.eps)
                {
                    return false;
                }
                for (int x = k; x < n + 1; x++)
                {
                    double tmp = tab[k, x];
                    tab[k, x] = tab[r, x];
                    tab[r, x] = tmp;
                }
                for (int x = k; x < n + 1; x++)
                {
                    tab[k, x] = tab[k, x] / max;
                }
                for (int x = 0; x < n; x++)
                {
                    if (x == k)
                        continue;
                    double t = tab[x, k];
                    for (int y = k; y < n + 1; y++)
                    {

                        tab[x, y] = tab[x, y] - t * tab[k, y];
                    }
                }
            }
            return true;
        }
        #endregion

        #region m detection
        public static double[,] mDetection(double a, double b, int k)
        {
            double[,] M = new double[k + 1, k + 2];

            for (int i = 0; i < k + 1; i++)
            {
                for (int j = 0; j < k + 1; j++)
                {
                    M[i, j] = delta(a, b, i + j + 1) / (i + j + 1.0);
                }
            }
            for (int i = 0; i < M.GetLength(0); i++)
            {
                double S1 = 0;
                double S2 = 0;
                int counter = 1;
                S1 = SimpsonSums.SimpsonFunction1(counter, a, b, i);
                do
                {
                    counter++;
                    S2 = S1;
                    S1 = SimpsonSums.SimpsonFunction1(counter, a, b, i);

                } while (Math.Abs(S1 - S2) > Details.eps);
                M[i, M.GetLength(1) - 1] = S1;
            }
            return M;
        }
        #endregion

        #region s detection
        public static double S(double a, double b, int i)
        {
            double S1 = 0;
            double S2 = 0;
            int counter = 1;
            S1 = SimpsonSums.SimpsonFunction2(counter, a, b, i);
            do
            {
                counter++;
                S2 = S1;
                S1 = SimpsonSums.SimpsonFunction2(counter, a, b, i);

            } while (Math.Abs(S1 - S2) > Details.eps);
            return S1;
        }

        #endregion

        #region p detection
        public static void P(double[,] M)
        {
            Console.WriteLine();
            Console.Write("P(x)=");
            for (int i = 0; i < M.GetLength(0); i++)
            {
                if (i == M.GetLength(0) - 1)
                {
                    Console.Write($"{M[i, M.GetLength(1) - 1]}x^{i} ");
                    break;
                }
                Console.Write($"{M[i, M.GetLength(1) - 1]}x^{i} +");
            }
            Console.WriteLine();
        }
        #endregion
    }

    internal static class SimpsonSums
    {
        public static double SimpsonSum1(int k, double a, double b, int j)
        {
            double temp = 0;
            for (int i = 1; i < Math.Pow(2, k); i += 2)
            {
                temp += MathFunctions.y(i, k, a, b, j);
            }
            return temp;
        }
        public static double SimpsonSum2(int k, double a, double b, int j)
        {
            double temp = 0;
            for (int i = 2; i < Math.Pow(2, k) - 1; i += 2)
            {
                temp += MathFunctions.y(i, k, a, b, j);
            }
            return temp;
        }
        public static double SimpsonSum3(int k, double a, double b, int j)
        {
            double temp = 0;
            for (int i = 1; i < Math.Pow(2, k); i += 2)
            {
                temp += MathFunctions.y2(i, k, a, b, j);
            }
            return temp;
        }
        public static double SimpsonSum4(int k, double a, double b, int j)
        {
            double temp = 0;
            for (int i = 2; i < Math.Pow(2, k) - 1; i += 2)
            {
                temp += MathFunctions.y2(i, k, a, b, j);
            }
            return temp;
        }
        public static double SimpsonFunction1(int k, double a, double b, int j) => (MathFunctions.h(k, a, b) / 3) * (MathFunctions.y(0, k, a, b, j) + MathFunctions.y(Convert.ToInt64(Math.Pow(2, k)), k, a, b, j) + 4 * SimpsonSum1(k, a, b, j) + 2 * SimpsonSum2(k, a, b, j));
        public static double SimpsonFunction2(int k, double a, double b, int j) => (MathFunctions.h2(k, a, b) / 3) * (MathFunctions.y2(0, k, a, b, j) + MathFunctions.y2(Convert.ToInt64(Math.Pow(2, k)), k, a, b, j) + 4 * SimpsonSum3(k, a, b, j) + 2 * SimpsonSum4(k, a, b, j));

    }

    public class Execute
    {
        #region variables
        readonly int k = 0;
        double[,] tab;
        #endregion

        public Execute()
        {
            Console.Write("Podaj k: ");
            k = Convert.ToInt32(Console.ReadLine());

            Console.WriteLine("\nPodaj przedział");
            Console.Write("od: ");
            double a = Convert.ToDouble(Console.ReadLine());
            Console.Write("do: ");
            double b = Convert.ToDouble(Console.ReadLine());
            Console.Clear();

            tab = MathFunctions.mDetection(a, b, k);

            Console.WriteLine("Macierz przed rozwiązaniem");
            DisplayMethods.Show(tab);
            Console.WriteLine();

            try
            {
                if (MathFunctions.GaussMethod(k + 1, ref tab))
                {
                    Console.WriteLine("Macierz po rozwiązaniu");
                    DisplayMethods.Show(tab);
                    MathFunctions.GlobalM = tab;
                    Console.WriteLine();
                    MathFunctions.P(tab);
                    Console.WriteLine($"S(k)= {MathFunctions.S(a, b, k)}");
                }
                else throw new Exception();
            }
            catch
            {
                Console.WriteLine("Macierzu układu osobliwa");
            }
            Console.WriteLine();
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            var program = new Execute();
        }
    }
}