/*==============================================================================
  Лабораторна робота №4. Завдання 1. Варіант 19.
  Тема: рекурентні співвідношення, степеневі ряди.

  Умова (таблиця 4.2, варіант 19): обчислити значення функції y, розвинувши
  функцію sin(x) у ряд Тейлора. Визначити похибку обчислення.

              | sin^2(x) - sin(x),     0 < x <= 1
          y = |
              | sin^3(x) + sin(2*x),  -2 <= x <= 0

  Вимоги методичних вказівок:
    - параметр функції змінюється від заданого з клавіатури мінімального
      значення до максимального із заданим кроком;
    - точність розвинення вводиться з клавіатури в діапазоні 1e-2 ... 1e-6;
    - для розвинення створити власну функцію, що обчислює суму ряду за
      рекурентним співвідношенням;
    - функції обчислення факторіалу та степеня НЕ використовувати;
    - похибка — різниця абсолютних значень наближеного та стандартного
      обчислення функції;
    - стандартне значення обчислювати бібліотечною функцією.

  Виконав: Одарчук Олексій, КНУ імені Тараса Шевченка, ФІТ, група ІПЗ-11.

  Компілятор: gcc -std=c17
==============================================================================*/

#include <stdio.h>
#include <stdbool.h>
#include <math.h>

/* Межі області визначення функції за умовою варіанта. */
const double X_MIN = -2.0;
const double X_MAX = 1.0;

/* Страховка від зациклення, якщо ряд не досягає заданої точності. */
const int MAX_TERMS = 1000;

/*------------------------------------------------------------------------------
  utf8Width — ширина рядка в символах, а не в байтах.

  У кодуванні UTF-8 кирилична літера займає два байти, тому специфікатор
  формату виду %-16s, який рахує байти, вирівнює заголовки таблиць
  неправильно. Функція рахує лише початкові байти символів: у продовжувальних
  байтах старші біти дорівнюють 10 у двійковій системі.

  Параметри: s [вхідний] — рядок у кодуванні UTF-8.
  Повертає : кількість символів рядка.
------------------------------------------------------------------------------*/
int utf8Width(const char *s)
{
    int width = 0;

    for (const unsigned char *p = (const unsigned char *)s; *p != '\0'; ++p)
        /* Продовжувальний байт UTF-8 має вигляд 10xxxxxx: маска 0xC0 лишає
           два старші біти, і якщо вони не дорівнюють 10, це початок символу. */
        if ((*p & 0xC0) != 0x80)
            ++width;

    return width;
}

/*------------------------------------------------------------------------------
  printPadded — вивести рядок, доповнивши його пропусками до заданої ширини.

  Параметри:
      s     [вхідний] — рядок, що виводиться;
      width [вхідний] — потрібна ширина поля в символах;
      left  [вхідний] — true: вирівнювання ліворуч; false: праворуч.

  Локальні змінні:
      pad — кількість пропусків, які треба додати.
------------------------------------------------------------------------------*/
void printPadded(const char *s, int width, bool left)
{
    const int pad = width - utf8Width(s);

    if (!left)
        for (int i = 0; i < pad; ++i)
            putchar(' ');

    printf("%s", s);

    if (left)
        for (int i = 0; i < pad; ++i)
            putchar(' ');
}

/*------------------------------------------------------------------------------
  sinTaylor — обчислити sin(t) розвиненням у ряд Тейлора (ряд Маклорена).

  Ряд:  sin(t) = t - t^3/3! + t^5/5! - t^7/7! + ...

  Рекурентне співвідношення. Позначимо a(n) — доданок з номером n (n = 0, 1, 2…):
        a(n) = (-1)^n * t^(2n+1) / (2n+1)!
  Тоді
        a(n+1)             -t^2
        ------  =  ------------------------
         a(n)        (2n+2) * (2n+3)
  Початковий доданок a(0) = t.

  Отже ані факторіал, ані піднесення до степеня окремо не обчислюються:
  кожен наступний доданок отримується з попереднього сталою кількістю
  множень і одним діленням, незалежно від номера доданка.

  Параметри:
      t     [вхідний]  — аргумент синуса;
      eps   [вхідний]  — точність: підсумовування триває, доки модуль
                         доданка не стане меншим за eps;
      terms [вихідний] — адреса лічильника доданків; може бути NULL,
                         якщо кількість доданків не потрібна.

  Повертає: наближене значення sin(t).

  Локальні змінні:
      sum  — накопичувана сума ряду;
      term — поточний доданок ряду;
      n    — номер поточного доданка;
      t2   — квадрат аргументу (обчислюється один раз).
------------------------------------------------------------------------------*/
double sinTaylor(double t, double eps, int *terms)
{
    const double t2 = t * t;

    double sum = t; /* a(0) = t */
    double term = t;
    int n = 0;
    int count = 1;

    while (fabs(term) >= eps && count < MAX_TERMS) {
        /* Рекурентний перехід a(n) -> a(n+1). */
        term *= -t2 / ((2.0 * n + 2.0) * (2.0 * n + 3.0));
        sum += term;
        ++n;
        ++count;
    }

    if (terms != NULL)
        *terms = count;

    return sum;
}

/*------------------------------------------------------------------------------
  yApprox — обчислити значення кусково-заданої функції y(x), використовуючи
            наближене (рядом Тейлора) значення синуса.

  Параметри:
      x     [вхідний]  — аргумент функції;
      eps   [вхідний]  — точність розвинення в ряд;
      terms [вихідний] — адреса лічильника доданків ряду (сумарно за виклик).

  Повертає: значення y(x).

  Передумова: x належить області визначення [-2; 1]; перевірку виконує
              функція inDomain() перед викликом.

  Локальні змінні:
      s   — наближене значення sin(x);
      s2x — наближене значення sin(2x);
      n1, n2 — кількість доданків у відповідних розвиненнях.
------------------------------------------------------------------------------*/
double yApprox(double x, double eps, int *terms)
{
    int n1 = 0, n2 = 0;

    if (x > 0.0) {
        /* Гілка 0 < x <= 1:  y = sin^2(x) - sin(x) */
        const double s = sinTaylor(x, eps, &n1);
        *terms = n1;
        return s * s - s;
    }

    /* Гілка -2 <= x <= 0:  y = sin^3(x) + sin(2x) */
    const double s = sinTaylor(x, eps, &n1);
    const double s2x = sinTaylor(2.0 * x, eps, &n2);
    *terms = n1 + n2;
    return s * s * s + s2x;
}

/*------------------------------------------------------------------------------
  yExact — обчислити те саме значення y(x) через бібліотечну функцію sin().
           Використовується як еталон для визначення похибки.

  Параметри: x [вхідний] — аргумент функції.
  Повертає : значення y(x), обчислене стандартними засобами.
------------------------------------------------------------------------------*/
double yExact(double x)
{
    if (x > 0.0) {
        const double s = sin(x);
        return s * s - s;
    }
    const double s = sin(x);
    return s * s * s + sin(2.0 * x);
}

/*------------------------------------------------------------------------------
  inDomain — чи належить аргумент області визначення функції [-2; 1].

  Параметри: x [вхідний] — аргумент функції.
  Повертає : true — функцію визначено в точці x.

  Порівняння виконується з невеликим допуском, щоб точки на межі діапазону
  не відкидались через похибку округлення при обчисленні аргументу.
------------------------------------------------------------------------------*/
bool inDomain(double x)
{
    const double eps = 1e-12;
    return x >= X_MIN - eps && x <= X_MAX + eps;
}

/*------------------------------------------------------------------------------
  readDouble — прочитати дійсне число з контролем коректності введення.

  Параметри:
      prompt [вхідний]  — текст запрошення;
      value  [вихідний] — адреса змінної для введеного числа.
  Повертає : true — число прочитано; false — вхідні дані вичерпано.
------------------------------------------------------------------------------*/
bool readDouble(const char *prompt, double *value)
{
    for (;;) {
        printf("%s", prompt);
        const int scanned = scanf("%lf", value);

        if (scanned == EOF)
            return false;
        if (scanned == 1)
            return true;

        int c;
        /* очищення буфера клавіатури */
        while ((c = getchar()) != '\n' && c != EOF) {
        }
        printf("Помилка: очікується дійсне число. Повторіть.\n");
    }
}

/*------------------------------------------------------------------------------
  printTableHeader — вивести заголовок таблиці значень функції: аргумент,
                     наближене та стандартне значення, похибка, доданки ряду.
------------------------------------------------------------------------------*/
void printTableHeader(void)
{
    printf("  ");
    printPadded("x", 10, false);
    printf(" | ");
    printPadded("y (ряд Тейлора)", 18, false);
    printf(" | ");
    printPadded("y (бібліотечна)", 18, false);
    printf(" | ");
    printPadded("похибка", 14, false);
    printf(" | ");
    printPadded("доданків", 8, false);
    printf("\n");
    printf("  -----------+--------------------+--------------------"
           "+----------------+---------\n");
}

/*------------------------------------------------------------------------------
  Головна функція. Читає межі та крок зміни аргументу, точність розвинення,
  табулює функцію і виводить таблицю значень з похибками.

  Локальні змінні:
      xStart, xEnd, step — межі та крок зміни аргументу;
      eps                — точність розвинення в ряд;
      x                  — поточне значення аргументу;
      approx, exact      — наближене та еталонне значення функції;
      error              — похибка обчислення;
      terms              — кількість доданків ряду;
      steps              — номер поточного кроку табулювання;
      maxError           — найбільша похибка за всю таблицю.
------------------------------------------------------------------------------*/
int main()
{
    printf("Лабораторна робота №4, завдання 1 (варіант 19)\n");
    printf("Виконав: студент групи ІПЗ-11 Одарчук Олексій\n");
    printf("Обчислення y(x) з розвиненням sin(x) у ряд Тейлора\n\n");
    printf("    y = sin^2(x) - sin(x),    якщо  0 < x <= 1\n");
    printf("    y = sin^3(x) + sin(2x),   якщо -2 <= x <= 0\n\n");
    printf("Область визначення: [-2; 1].\n\n");

    double xStart = 0.0, xEnd = 0.0, step = 0.0, eps = 0.0;

    if (!readDouble("Уведіть початкове значення аргументу: ", &xStart))
        return 1;
    if (!readDouble("Уведіть кінцеве значення аргументу:   ", &xEnd))
        return 1;

    do {
        if (!readDouble("Уведіть крок зміни аргументу (> 0): ", &step))
            return 1;
        if (step <= 0.0)
            printf("Помилка: крок має бути додатним.\n");
    } while (step <= 0.0);

    do {
        if (!readDouble("Уведіть точність (від 1e-6 до 1e-2): ", &eps))
            return 1;
        if (eps < 1e-6 || eps > 1e-2)
            printf("Помилка: точність має бути в діапазоні 1e-6 ... 1e-2.\n");
    } while (eps < 1e-6 || eps > 1e-2);

    if (xStart > xEnd) {
        printf("\nПочаткове значення більше за кінцеве — таблиця порожня.\n");
        return 2;
    }

    printf("\nТочність розвинення: %.1e\n\n", eps);
    printTableHeader();

    double maxError = 0.0;
    int steps = 0;

    /* Цикл табулювання. Лічильник кроків використано замість накопичення
       x += step, щоб похибка додавання не накопичувалась від кроку до кроку. */
    for (;;) {
        const double x = xStart + steps * step;
        if (x > xEnd + 1e-12)
            break;

        if (!inDomain(x)) {
            printf("  %10.4f | %s\n", x,
                   "функцію не визначено (аргумент поза межами [-2; 1])");
        } else {
            int terms = 0;
            const double approx = yApprox(x, eps, &terms);
            const double exact = yExact(x);

            /* Похибка за умовою варіанта — різниця абсолютних значень
               наближеного та стандартного обчислення функції. */
            const double error = fabs(fabs(approx) - fabs(exact));

            if (error > maxError)
                maxError = error;

            printf("  %10.4f | %18.12f | %18.12f | %14.3e | %8d\n", x, approx, exact,
                   error, terms);
        }

        ++steps;
    }

    printf("\nОброблено точок: %d\n", steps);
    printf("Найбільша похибка за таблицею: %.3e\n", maxError);
    printf("Задана точність розвинення:    %.3e\n", eps);

    return 0;
}