SLAU

#include <iostream>
#include <omp.h>
#include <cmath>
#include <complex>
#define pi 3.14159265
#define eps 0.000001
using namespace std;

complex<double> gamma(double n, const double& K)
{
    complex<double> res;
    if (K * K >= pi * pi * n * n)
    {
        res = sqrt(K * K - pi * pi * n * n);
    }
    else
    {
        res = 1i * sqrt(pi * pi * n * n - K * K);
    }
    return res;
}
complex<double> Jnm(double n, double m, const double& a0, const double& a1)
{
    complex<double> res;
    if (n == m)
    {
        res = a1 - a0 - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
    }
    else
    {
        res = (sin((n - m) * pi * a1) - sin((n - m) * pi * a0)) / ((n - m) * pi) - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
    }
    return res;
}
complex<double> Inm(double n, double m, const double& a0, const double& a1)
{
    complex<double> res;
    if (n == m)
    {
        res = 1.0 - Jnm(n, m, a0, a1);
    }
    else
    {
        res = -Jnm(n, m, a0, a1);
    }
    return res;
}
complex<double> sum(complex<double>* x0, int k, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
{
    int i, j;
    complex<double>s;
    complex<double>h;
    for (i = 0; i < N; i++)
    {
        h = (1.0 - exp(2.0 * 1i * c * gamma((double)i + 1, K)));
        for (j = 0; j < N; j++)
        {
            s = s + (x0[i] * h * gamma((double)j + 1, K) * Jnm((double)i + 1, (double)j + 1, a0, a1) * Inm((double)j + 1, (double)k + 1, a0, a1)) / (1.0 - exp(2.0 * 1i * gamma((double)j + 1, K) * c));
        }
    }
    return s;
}
complex<double>* parallel_area(complex<double>* x0, const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
{

    int i;
    complex<double>* x = new complex<double>[N];
    #pragma omp parallel num_threads(n)
    {
        #pragma omp for
        for (i = 0; i < N; i++)
        {
            x[i] = (sum(x0, i, N, c, l, a0, a1, K) + (gamma(l, K) * Inm(l, (double)i + 1, a0, a1))) / gamma((double)i + 1, K);
        }
    }
    return x;
    //Последовательно
    /*
    //complex<double>s1, s2, s3;
    for (i = 0; i < N; i++)
    {
        //s1 = sum(x0, i, N, c, l, a0, a1, K);
        //s2 = gamma(l, K) * Inm(l, (double)i + 1, a0, a1);
        //s3 = gamma((double)i + 1, K);
        x[i] = (sum(x0, i, N, c, l, a0, a1, K) + gamma(l, K) * Inm(l, (double)i + 1, a0, a1)) / gamma((double)i + 1, K);
        //x[i] = (s1 + s2) / s3;
    }
    return x;
    */
}
bool mensheeps(complex<double>* x, complex<double>* x0, int N)
{
    double max = abs(x[0] - x0[0]);
    double temp;
    for (int i = 1; i < N; i++)
    {
        temp = abs(x[i] - x0[i]);
        if (temp > max)
            max = temp;
    }
    if (max <= eps)
        return true;
    else return false;
}
void show(complex<double>* x, int N)
{
    for (int i = 0; i < N; i++)
        cout <<abs(x[i])<< "\n";
    cout << "\n";
}
void Yacobi(const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
{
    int i;
    //Задаем вектор нач. приближения
    complex<double>* x0 = new complex<double>[N]; //по умолчанию иниц. нулями
    for (i = 0; i < N; i++)
        x0[i] = 1.0;
    complex<double>* x = parallel_area(x0, n, N, c, l, a0, a1, K);
    while (!mensheeps(x, x0, N))
    {
        for (i = 0; i < N; i++)
            x0[i] = x[i];
        x = parallel_area(x0, n, N, c, l, a0, a1, K);
    }
    show(x, N);
    delete[]x0;
    delete[]x;
}
int main()
{
    setlocale(LC_ALL, "RU");
    double c, l, a0, a1, K;
    cout << "Введите параметры с, l, a0, a1, к cоответственно: \n";
    cin >> c >> l >> a0 >> a1 >> K;
    int n, N;
    cout << "Введите размерность системы N: \n";
    cin >> N;
    cout << "Задайте количество потоков: \n";
    cin >> n;
    Yacobi(n, N, c, l, a0, a1, K);
    /*
    ofstream out("den.txt");
    for (int i = 1; i <= 10; i++)
    {
        //cout << "n = " << i << " gamma = " << gamma(i, K) << "\n";
        for (int j = 1; j <= 10; j++)
        {
            out << "n = " << i << " m = " << j << "\n"
                << "Jnm = " << Jnm(i, j, a0, a1) << " Inm = " << Inm(i, j, a0, a1) << " \n";
        }
    }
    out.close();
    */
    return 0;
}