неділю, 13 червня 2021 р.

Розв'язок задачі SubmultiplesOfN з Топкодер.

У цьому пості розгляну як розв'язати задачу SubmultiplesOfN. Це задача середньої складності з першого дивізіону, що була частиною SRM 805. У ній потрібно було знайти скільки разів зустрічаються всі можливі кратні числу $N$ в числі $B$. Під зустрічаються мається на увазі чи є число представлене в текстовому вигляді в десятковій системі підпослідовністю такого ж представлення іншого числа. Наприклад, "23" зустрічається в "213", але не в "312".

Для розв'язання цієї задачі ми використаємо динамічне програмування.

Для того, щоб зрозуміти, що покласти в кожну комірку таблиці нам треба згадати, що в десятковій системі $abcd = ((a \times 10 + b) \times 10 + c) \times 10 + d. $ І ще, що для будь-якого додатного $N $ $abcd \% N = ((a \times 1000) \% N + (b \times 100) \% N + (c \times 10) \% N + d \% N) \% N. $

Зауважимо, що якщо $a \% N$ рівне нулю, то $a$ кратне $N$, тобто, якщо $a$ присутнє в $B,$ то ми вже маємо врахувати це додавши одиничку до відповіді. Власне, саме таку перевірку для всіх цифр окрім нуля ми і робимо коли ініціалізуємо перший рядок таблиці.

Виходить, що кожна комірка $(i, j)$ представляє всі можливі числа, що закінчуються на $B[i],$ зустрічаються в $B[0..i]$ і мають залишок $j$ при діленні на $N.$

const auto Mod = 1000000007ll;

const auto MaxDigitsCount = 5000;
const auto MaxN = 1000;
const auto Base = 10;

long long DP[MaxDigitsCount + 1][MaxN + 1];
long long Next[MaxDigitsCount + 1][Base];

class Solution {

    static const auto kNoPosistion = -1;

    // Швидкий пошук необхідної цифри починаючи із заданої позиції
    void InitTransitions(string b)
    {
        const auto len = (int)b.length();
        
        for (auto d = 0; d < Base; ++d)
        {
            Next[len][d] = kNoPosistion;
        }

        for (auto i = len - 1; i >= 0; --i)
        {
            for (auto d = 0; d < Base; ++d)
            {
                Next[i][d] = Next[i + 1][d];
            }
            Next[i][b[i] - '0'] = i;
        }
    }

public:
    int count(string B, int N)
    {
        InitTransitions(B);

        auto result = (long long)0;

        const auto len = (int)B.length();
        
        for (auto i = 0; i < len + 1; ++i)
        {
            for (auto j = 0; j < N + 1; ++j)
            {
                DP[i][j] = 0;
            }
        }

        // Ініціалізуємо перший рядок таблиці
        for (auto d = 1; d < Base; ++d)
        {
            const auto nextDigitPosition = Next[0][d];

            if (nextDigitPosition != kNoPosistion)
            {
                const auto reminder = d % N;
                DP[nextDigitPosition][reminder] += 1;

                if (reminder == 0)
                {
                    result += 1;
                }
            }
        }

        for (auto position = 0; position < len - 1; ++position)
        {
            for (auto reminder = 0; reminder < N; ++reminder)
            {
                for (auto d = 0; d < Base; ++d)
                {
                    const auto nextDigitPosition = Next[position + 1][d];
                    
                    if (nextDigitPosition != kNoPosistion)
                    {
                        const auto appended = reminder * 10 + d;
                        const auto newReminder = appended % N;
                        const auto upToD = DP[position][reminder];

                        (DP[nextDigitPosition][newReminder] += upToD) %= Mod;

                        if (newReminder == 0)
                        {
                            (result += upToD) %= Mod;
                        }
                    }
                }
            }
        }

        return result;
    }
};

#pragma comment(linker, "/STACK:200000000")

int main()
{
    // В дужках правильна відповідь.
    cout << Solution().count("12", 2) << " (2)" << endl;
    cout << Solution().count("14", 7) << " (1)" << endl;
    cout << Solution().count("10", 2) << " (1)" << endl;
    cout << Solution().count("12345678", 2) << " (170)" << endl;
    cout << Solution().count("1111111111", 7) << " (1)" << endl;
    cout << Solution().count("1357913579135791357913579", 2) << " (0)" << endl;
    cout << Solution().count("1122334455", 6) << " (20)" << endl;
    cout << Solution().count("1020402", 24) << " (6)" << endl;
    cout << Solution().count("123456789012345678901234567890", 1) << " (62224120)" << endl;
    
}

суботу, 19 вересня 2020 р.

Жорданова нормальна форма і мінімальний многочлен.

Мінімальний многочлен

Нехай задана дійсна матриця $A$ розміром $n \times n.$

Лема М1
Існує дійсний многочлен $p$ такий, що $p(A) = 0.$
Простір $M_{n\times n}(\mathbb R) - $ це $n^2$-вимірний векторний простір. Візьмімо $\mathbf I, \mathbf A, \mathbf A^2, \dots, \mathbf A^{n^2}.$ Всього маємо $n^2 + 1$ елементів, отже вони лінійно залежні: $\mu_0 \mathbf I + \mu_1 \mathbf A + \dots + \mu_{n^2}\mathbf A^{n^2} = 0.$ Тоді можна записати многочлен $\ p(t) = \sum \mu t^i$ такий, що $p(A) = 0.$

Нехай серед усіх таких многочленів такий, щоб мав найменший степінь і був нормованим, тобто перший коефіцієнт був $1.$ Називатимемо такий многочлен *мінімальним многочленом* матриці $\mathbf A.$

Лема М2
Мінімальний многочлен - унікальний.
Якби ми мали два таких многочлени, то остача від їхнього ділення була б многочленом меншого степеня, що давав би нуль для $\mathbf A.$
Лема М3
Якщо $p$ це деякий многочлен, для якого $p(A) = 0,$ тоді $m$ ділить $p.$
Зрозуміло, що $\deg(p) \ge \deg(m),$ тоді можна записати таке: $p = mq + r$ для певних многочленів $q, r$ таких, що $\deg(r) < \deg(m).$ Тоді $r(\mathbf A) = p(\mathbf A) - m(\mathbf A)q(\mathbf A) = 0,$ але ж $r$ має меншу степінь ніж $m,$ це протирічить мінімальності $m.$
Лема М4
Якщо $\lambda $ це в-значення $\mathbf A,$ тоді воно корінь $m.$
Для певного $v \ne 0$ маємо, що $\mathbf Av = \lambda v.$ Також $\mathbf A^k v^k = \mathbf A^{k-1} \lambda v = \cdots = \lambda^k v,$ для будь-кого цілого $k \ge 0.$ Отже, для будь-якого многочлена $p(A)v = p(\lambda)v,$ в тому числі і для мінімального. З того, що $m(\mathbf A) = 0$ і $v \ne 0$ випливає, що $m(\lambda) = 0.$

Теорема Гамільтона — Келі: Нехай $\chi(t) = \det(\mathbf A - t\mathbf I)$ буде характеристичним многочленом $\mathbf A.$ Тоді $\chi(\mathbf A) = 0.$

Лема М5
Якщо $\chi(t) = \prod_{i=1}^r (t-\lambda_i)^{a_i},$ то що ми можемо сказати про $m(t)?$
Згідно з Гамільтоном — Келі, М3 і М4 маємо, що $m$ може мати лише такий вигляд $\prod_{i=1}^r(t-\lambda_i)^{c_i},$ де $1 \le c_i \le a_i.$
Лема М6
$\chi(t)$ і $m(t) - $ інваріанти щодо зміни базису.
Нехай $\mathbf B = \mathbf P^{-1} \mathbf A \mathbf P.$ Зауважимо, що для довільного многочлена $p(\mathbf A) = \mathbf P^{-1}p(\mathbf B) \mathbf P,$ бо $\mathbf B^k = \mathbf P^{-1} \mathbf A^k \mathbf P.$ Тому $p(\mathbf A) = 0$ тоді і тільки тоді, коли $p(\mathbf B) = 0.$ Отже мінімальні многочлени для $\mathbf A$ і $\mathbf B$ ділять один одного і, відповідно, рівні.

Жорданові клітини і жорданова нормальна форма

$\mathbf J = \begin{pmatrix} a & 1 & & \\ & a & 1 \\ & & a & 1 \\ & & & a \\ \end{pmatrix} - $ жорданова клітина, $\begin{pmatrix} \mathbf J_1 & & & \\ & \mathbf J_2 & & \\ & & \ddots & \\ & & & \mathbf J_k \\ \end{pmatrix} - $ жорданова нормальна форма.

Лема Ж1: Нехай $\mathbf J - $ жорданова клітина з $\lambda$ на діагоналі. Тоді $\lambda - $ єдине в-значення з власним простором розмірності $1.$

Лема Ж2: Нехай $\mathbf J$ це жорданова клітина розміру $k \times k,$ тоді $(\mathbf J - \lambda \mathbf I)^n = 0$ для $n = k$ і не рівна нулю для $n < k.$

Лема Ж3: Якщо $A, B$ це блочні матриці з однаковими розмірами блоків, тоді

$\begin{pmatrix} \mathbf A_1 & & \\ & \ddots & \\ & & \mathbf A_n \end{pmatrix} \begin{pmatrix} \mathbf B_1 & & \\ & \ddots & \\ & & \mathbf B_n \end{pmatrix} = \begin{pmatrix} \mathbf A_1\mathbf B_1 & & \\ & \ddots & \\ & & \mathbf A_n\mathbf B_n \end{pmatrix}$
Теорема
Припустимо, що задана $\mathbf A$ в нормальній жордановій формі. Припустимо, що відмінні діагональні елементи це $\lambda_1, \dots, \lambda_k,$ і $\lambda_i$ з'являється $a_i$ раз. Припустимо, що кількість жорданових клітин з $\lambda_i$ на діагоналі становить $g_i.$ І припустимо, що найбільша жорданова клітина з $\lambda_i$ на діагоналі має розмір $c_i.$ Тоді
  1. $\chi(t) = \prod_{i=1}^k(t-\lambda_i)^{a_i}$
  2. $m(t) = \prod_{i=1}^k(t-\lambda_i)^{с_i}$ (степінь $(t-\lambda_i)$ в $m(t)$ це розмір найбільшої $\lambda_i$-клітини)
  3. $\dim\ker(\mathbf A - \lambda_i \mathbf I) = g_i$ (вимірність в-простору $\lambda_i$ це кількість $\lambda_i$-клітин)
  1. Раз $\mathbf A$ в ЖНФ, то вона верхньотрикутна, отже її в-значення у неї на головній діагоналі.
  2. Згідно з першим пунктом і М4 маємо, що $m(t) = \prod_{i=1}^k(t-\lambda_i)^{k_i}$ для деяких $k_i.$ Тепер, для того, щоб _вбити_ найбільшу клітину для $j$-го в-значення потрібно $\mathbf A - \lambda_j \mathbf I$ піднести до $c_j$ степеня. Отже, якщо $k_i = c_i,$ то в $(A-\lambda_i)^{k_i}$ всі $\lambda_i$-клітини нульові. Отже, використовуючи Ж3, щоб отримати нулі для $\lambda_p$ і $\lambda_q$ нам треба взяти $(\mathbf A - \lambda_p \mathbf I)^{c_p}(\mathbf A - \lambda_q \mathbf I)^{c_q}.$ Роблячи це для всіх в-значень ми отримуємо нульову матрицю.
  3. Згідно з Ж1 маємо по одному виміру в-простору для кожної жорданової клітини.

вівторок, 8 січня 2019 р.

Довідничок для перетворення Фур'є

Період 2$\pi$

Фур'є
Припустімо, що $f(t)$ має період $2\pi$, тоді \begin{equation} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left(a_n\cos(nt) + b_n \sin(nt)\right)\label{eq:FourierPi} \end{equation} де коефіцієнти $a_0, a_1, \dots, b_1, \dots$ обчислюються так: \begin{align} a_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\cos(nt)dt\label{eq:anpi}\\ b_n &= \frac{1}{\pi}\int_{-\pi}^{\pi}f(t)\sin(nt)dt\label{eq:bnpi}\\ \end{align} Знак рівності може не виконуватись в місцях розривів функції.
Означення
В \eqref{eq:FourierPi} частота $n=1$ називається фундаментальною, частоти з $n > 1$ називаються гармоніками.

Ортогональність базису

$$ \int_{-\pi}^{\pi}\cos(nt)\cos(nt)dt = \frac{1}{2} \int_{-\pi}^{\pi}\cos^2(nt)+\sin^2(nt)dt = \frac{1}{2}n2\frac{\pi}{n}=\pi $$ Через це $\pi$ ми і ділимо в \eqref{eq:anpi} і \eqref{eq:bnpi} на $\pi$.

Тепер розглянемо з $n\neq m$: \begin{align} \int_{-\pi}^{\pi}\cos(nt)\cos(mt)dt &= \frac{1}{2}\int_{-\pi}^{\pi}\left(\cos((n+m)t) + \cos((n-m)t)\right)dt &= 0\\ \int_{-\pi}^{\pi}\sin(nt)\cos(mt)dt &= \frac{1}{2}\int_{-\pi}^{\pi}\left(\sin((n+m)t) + \sin((n-m)t)\right)dt &= 0\\ \end{align}

Період 2L

Фур'є
Тепер припустимо, що $f(t)$ має період $P = 2L$. Тоді ряд Фур'є набуває вигляду \begin{equation} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\left(a_n\cos(n\frac{\pi}{L}t) + b_n\sin(n\frac{\pi}{L}t)\right). \end{equation} А коефіцієнти такого: \begin{align} a_n &= \frac{L}{\pi}\int_{-L}^{L}f(t)\cos(n\frac{\pi}{L}t)dt\\ b_n &= \frac{L}{\pi}\int_{-L}^{L}f(t)\sin(n\frac{\pi}{L}t)dt\\ \end{align}

Дискретне перетворення Фур'є (ДПФ)

Базис в $\mathbb{C}^N$

Розглянемо вектори $\mathbf{w}^{(k)} = [\dots, e^{j\frac{2\pi}{N}ik}, \dots]'$. Фундаментальною частотою тут буде $\frac{2\pi}{N}k$. Можна побачити, що складові вектора $\mathbf{w}^{(1)}$ роблять один оборот навколо одиничного кола. Найшвидше ми рухаємось по колу коли $k = N/2$ далі, зі збільшенням $k$ швидкість падає, а ознакою того що фундаметальна частота перетнула $\pi$ буде протилежний знак уявної частини.
Базис в $\mathbb{C}^N$
$N$ сигналів $\mathbf{w}^{(k)}_n = e^{j\frac{2\pi}{N}nk}$, де $n,k \in [0, N)$ - це ортогональний базис в $\mathbb{C}^N$.
Для $k\ne h$ \begin{align*} \langle\mathbf{w}^{(k)}\mathbf{w}^{(h)}\rangle &= \sum_{n=0}^{N-1}e^{-j\frac{2\pi}{N}nk}e^{j\frac{2\pi}{N}nh}\\ &= \begin{cases} N,& k = h\\ \frac{1 - e^{j\frac{2\pi}{N}N(h-k)}}{1 - e^{j\frac{2\pi}{N}(h-k)}} = 0,& k \ne h\\ \end{cases} \end{align*} Тут ми скористались тим, що $(1+x+\cdots+x^{N-1})(1-x) = 1-x^N$. А також фактом, що $e^{j2\pi(h-k)} = 1$ в результаті чого числівник стає нулем.
Отже, $\mathbf{w}^{(k)}$ формують ортогональний базис для $\mathbb{C}^N$. Але цей базис не ортонормальний. Коефіцієнтом ортонормалізації буде $\frac{1}{\sqrt{N}}$.

Частоти базису

Якщо розглянути як складові базисних векторів рухаються уздовж одиничного кола, то помітно, що спочатку зі сростанням індексу базисного вектора маркер, що позначає складову, рухається все швидше, досягаючі найбільшої швидкості коли $k = N/2 - 1$ і $k = N/2$ для парного $N$. Після цього ми спостерігаємо ефект колеса з фільмів, тобто коли колесо крутиться в зворотній бік, це пов'язано з частотою кадрів. Так і тут, коли маркер робить більше ніж півоберта за крок, то нам здається, що він обертається в зворотньому напрямку з меншою швидкістю.
Python-код генерації анімації.
import numpy as np
import matplotlib.pylab as plt
from matplotlib.animation import FuncAnimation

N = 32
t1 = np.linspace(0,N-1,N)
t10 = np.linspace(0,N-1,N*10)

fig, (ax1, ax2) = plt.subplots(2,1)
fig.set_tight_layout(True)

ax1.set_xlim(-1,N)
ax2.set_xlim(-1,N)

def update(k):
    ax1.cla()
    ax2.cla()
    ax1.set_ylabel("Re")
    ax2.set_ylabel("Im")
    if (k < N):
        if k < N/2:
            ax1.plot(t10, np.cos(2*np.pi/N*(k)*t10))
            ax2.plot(t10, np.sin(2*np.pi/N*(k)*t10))
        else:
            ax1.plot(t10, np.cos(2*np.pi/N*(N - k)*t10))
            ax2.plot(t10, -np.sin(2*np.pi/N*(N - k)*t10))
        ax1.set_title('Базисний вектор #' + str(k) + ' в С^' + str(N))
        markerline, stemlines, baseline = ax1.stem(t1, np.cos(2*np.pi/N*k*t1))
        plt.setp(baseline, color='black', linewidth=2)
        plt.setp(markerline, color='orange')
        plt.setp(stemlines, color='orange', linewidth=2)
        markerline, stemlines, baseline = ax2.stem(t1, np.sin(2*np.pi/N*k*t1), 'violet')
        plt.setp(baseline, color='black', linewidth=2)
        plt.setp(markerline, color='orange')
        plt.setp(stemlines, color='orange', linewidth=2)

anim = FuncAnimation(fig, update, frames=np.arange(0, N-1 + 5), interval=1000)
plt.show()
Зауважте, що уявна частина змінила знак в другій половині. Бо частота стала від'ємною.

До коефіцієнтів і назад

\begin{equation} X_k = \langle\mathbf{w}^{(k)},\mathbf{x}\rangle\label{eq:analysis}. \end{equation} \begin{equation} \mathbf{x} = \frac{1}{N}\sum_{k=0}^{N-1}X_k\mathbf{w}^{(k)}\label{eq:synthesis}. \end{equation}

Зміна базису в матричній формі

Ми можемо записати зміну базису в матричній формі. Якщо позначити $W=e^{-j\frac{2\pi}{N}}$, то \begin{equation} \mathbf{W} = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & W & W^2 & \dots & W^{N-1}\\ 1 & W^2 & W^4 & \dots & W^{2(N-1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & W^{N-1} & W^{2(N-1)} & \dots & W^{(N-1)^2}\\ \end{bmatrix} \end{equation} Фактично, в кожному рядку маємо спряжений до базового вектора. Тепер формули \eqref{eq:analysis} і \eqref{eq:synthesis} можна записати в матричній формі: \begin{equation} \mathbf{X} = \mathbf{W}\mathbf{x}. \end{equation} \begin{equation} \mathbf{x} = \frac{1}{N}\mathbf{W}^H\mathbf{X}. \end{equation}

Перетворення Фур'є дискретного часу (ПФДЧ)

Якщо збільшувати кількість точок зростає, то фундаментальна частота $2\pi/N$ стає дуже малою і, чисто інтуїтивно, ми можемо казати, що $\lim_{N\to\infty}(2\pi/N)k = \omega$.
Формально
  • $x[n] \in l_2(\mathbb{Z})$ (сумовна з квадратом послідовність)
  • означимо функцію \begin{equation} F(\omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\label{eq:dtft} \end{equation}
  • обернена (коли $F(\omega)$ існує): \begin{equation} x[n] = \frac{1}{2\pi}\int_{-\pi}^{\pi}F(\omega)e^{j\omega n}d\omega,\quad n\in\mathbb{Z}\label{eq:idtft} \end{equation}
Зазвичай замість $F(\omega)$ пишуть $X(e^{j\omega})$. З поміж іншого, при такому записі краще видно $2\pi$ періодичність. Прийнято представляти $X(e^{j\omega})$ на $[-\pi, \pi]$.

Властивості

  • лінійність: \begin{equation}\mbox{ПФДЧ}\{\alpha x[n] + \beta y[n]\} = \alpha X(e^{j\omega})+\beta Y(e^{j\omega})\end{equation}
  • часовий зсув: \begin{equation}\mbox{ПФДЧ}\{x[n-M]\} = e^{-j\omega M}X(e^{j\omega})\label{eq:dtft_timeshift}\end{equation}
  • модуляція: \begin{equation}\mbox{ПФДЧ}\{e^{j\omega_0 n} x[n]\} = X(e^{j(\omega - \omega_0)})\label{eq:dtft_modulation}\end{equation}
  • зворотній час: \begin{equation}\mbox{ПФДЧ}\{x[-n]\} = X(e^{-j\omega})\end{equation}
  • спряження: \begin{equation}\mbox{ПФДЧ}\{x^*[n]\} = X^*(e^{-j\omega})\end{equation}
Властивості \eqref{eq:dtft_timeshift} і \eqref{eq:dtft_modulation} двоїсті одна до одної.

ПФДЧ одиниці

Очевидно, що стала $1$ несумовна з квадратом. Але $\mbox{ДПФ}\{1\} = N\delta[n]$. Було б добре могти обчислити її ПФДЧ. Для цього нам знадобиться дельта-функція Дірака. В ній нас цікавить така властивість: \begin{equation} (\delta; f) = \int_{-\infty}^{+\infty}\delta(x-a)f(x)dx = f(a), \end{equation} де $f(x)$ - це будь-яка неперервна функція. Таким чином ми можемо видобути значення функції в будь-якій точці.

Нам потрібно застосовувати дельта-функцію Дірака в просторі $2\pi$-періодичних функцій, тому нам потрібна така послідовність імпульсів:

\begin{equation} \tilde{\delta}(\omega) = 2\pi\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k) \end{equation}
Python-код генерації зображення
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.ticker as tckr

f,(ax,ax2)=plt.subplots(2,1)
ax2.axis('off')
x=np.linspace(-4*np.pi,4*np.pi,5)
y=[1,1,1,1,1]
ax.stem(x/np.pi,y, markerfmt='^')
ax.xaxis.set_major_formatter(tckr.FormatStrFormatter('%g $\pi$'))
ax.xaxis.set_major_locator(tckr.MultipleLocator(base=1.0))
plt.show()
Добре, тепер давайте розглянемо обернене перетворення Фур'є дискретного часу цієї послідовності імпульсів. Згідно з \eqref{eq:idtft}: \begin{align} \mbox{ОПФДЧ}\{\tilde{\delta}(\omega)\} &= \frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{\delta}(\omega)e^{j\omega n}d\omega\nonumber\\ &= \int_{-\pi}^{\pi}\delta(\omega)e^{j\omega n}d\omega\nonumber\\ &= e^{j\omega n}|_{\omega=0}\nonumber\\ & = 1\label{eq:idtft_one} \end{align} Давайте перевіримо це за допомогою чисельного аналізу. Розглянемо частинних сум $S_k = \sum_{n=-k}^k e^{-j\omega n}$.
Python-код для відображення частинних сум
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(1,1)

N = 500
x = np.linspace(-np.pi, np.pi, N)

def update(k):
    re = np.zeros(N)
    im = np.zeros(N)
    for n in range(-k, k+1):
        re += np.cos(x*n)
        im += np.sin(x*n)
    y = np.multiply(re, re) + np.multiply(im, im)
    y = np.sqrt(y)
    ax.cla()
    ax.set_title('Частинні суми k = ' + str(k))
    ax.set_xlim(-np.pi, np.pi)
    ax.set_ylim(0,100)
    ax.plot(x, y);

anim = FuncAnimation(fig, update, frames=np.arange(0, 50), interval=100)
plt.show()
З \eqref{eq:idtft_one} за допомогою властивості \eqref{eq:dtft_modulation} маємо: \begin{equation} \mbox{ПФДЧ}\{e^{j\omega_0 n}\} = \tilde{\delta}(\omega - \omega_0). \end{equation} А використовуючи формулу Ейлера можна знайти ще дві: $$\mbox{ПФДЧ}\{\cos \omega_0 n\} = [\tilde\delta(\omega-\omega_0) + \tilde\delta(\omega+\omega_0)]/2$$ $$\mbox{ПФДЧ}\{\sin \omega_0 n\} = -j[\tilde\delta(\omega-\omega_0) - \tilde\delta(\omega+\omega_0)]/2$$

Вбудовування сигналу скінченної довжини

Природне спектральне представлення сигналу скінченної довжини це ДПФ ($X[k]$). Існує два способи вбудувати такий сигнал в нескіченну послідовність:
  • періодичне продовження: $\tilde{x}[n] = n[n\mod N]$
  • продовження зі скінченним носієм: $\bar{x}[n] = \begin{cases} x[n] & 0\le n \le N-1\\ 0 & \mbox{інакше}\\ \end{cases}$

Зараз ми розглянемо зв'язок між ДПФ скінченного сигналу $X[n]$ і ПФДЧ продовженого сигналу $\tilde{X}(e^{j\omega})$.

Випадок періодичного сигналу

Починаючи з формули \eqref{eq:dtft} маємо: \begin{align*} \tilde{X}(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}\\ &= \sum_{n=-\infty}^{\infty}\left(\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N} nk}\right)e^{-j\omega n}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]\sum_{n=-\infty}^{\infty}e^{j\frac{2\pi}{N} nk}e^{-j\omega n}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]~\mbox{ПФДЧ}\{e^{j\frac{2\pi}{N} k}\}\\ &=\frac{1}{N}\sum_{k=0}^{N-1}X[k]\tilde{\delta}(w - \frac{2\pi}{N} k)\\ \end{align*}

Випадок сигналу зі скінченним носієм

Спочатку розглянемо сигнал $$\bar{r}[n] = \begin{cases}1&0\le n \le N-1\\0& \mbox{інакше}\end{cases}$$ Його ПФДЧ дорівнює $\sum_{n=0}^{N-1}e^{-j\omega n}$, яке ми позначатимемо $\bar{R}(e^{j\omega})$. \begin{align*} \bar{X}(e^{j\omega}) &= \sum_{n=-\infty}^{\infty}\bar{x}[n]e^{-j\omega n}\\ &= \sum_{n=0}^{N-1}x[n]e^{-j\omega n}\\ &= \sum_{n=0}^{N-1}\left(\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N} nk}\right)e^{-j\omega n}\\ &= \frac{1}{N}\sum_{k=0}^{N-1}X[k]\sum_{n=0}^{N-1}e^{j\frac{2\pi}{N} nk}e^{-j\omega n}\\ &= \frac{1}{N}\sum_{k=0}^{N-1}X[k]\bar{R}(e^{\omega - \frac{2\pi}{N} k})\\ \end{align*} Зараз докладніше розглянемо $\bar{R}(e^{j\omega})$: \begin{align*} \sum_{n=0}^{N-1}e^{-j\omega n} &= \frac{1 - e^{-j\omega N}}{1 - e^{-j\omega}}\\ &=\frac{e^{-j\omega N/2}\left(e^{j\omega N/2}-e^{-j\omega N/2}\right)}{e^{-j\omega N}\left(e^{j\omega/2}-e^{-j\omega/2}\right)}\\ &=\frac{\sin(\omega N/2)}{\sin(\omega / 2)}e^{-j\omega (N-1)/2}\\ \end{align*}

вівторок, 25 грудня 2018 р.

Хвильки (вейвлети) Хаара в 1D

Хвильки (вейвлети) - це математичне приладдя для ієрархіного розкладення зображень. У цій статті я розповім про найпростішу форму хвильок - базис Хаара в одновимірному просторі.

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

$$[9\ 7\ 3\ 5]$$

Ми можемо зберегти це зображення залишивши лише середні значення для кожної двійки пікселів:

$$[8\ 4]$$

Очевидно, що частина інформації загубилась. Щоб отримати початкові чотири пікселі нам також треба зберегати коефіцієнти деталізації. Тут такими коефіцієнтами є $1$ для першої пари і $-1$ для другої пари значень.

РоздільністьСереднє значенняКоеф. деталізації
$4$$[\ 9\ 7\ 3\ 5\ ]$
$2$$[\ 8\ 4\ ]$$[\ 1\ {-1}\ ]$
$1$$[\ 6\ ]$$[\ 2\ ]$
Отже, для одновимірного базиса Хаара хвильковие перетворення нашого зображення виглядає так: $$[\ 6\ 2\ 1\ -1\ ]$$ Зауважте, що інформація не була ані втрачена, ані набута. Початкове зображення мало 4 коефіцієнти і перетворене також. Початкове зображення можна легко відновити з його перетворення.

Базис хвилькового перетворення

Розглянемо дві ортогональні функції:
  1. Маштабувальна - $\psi(x)$.
  2. Хвилькова (материнська) - $\phi(x)$.
Які задаються так: \begin{equation} \phi(s) := \begin{cases} 1 & 0 \le x < 1\\ 0 & \mbox{інакше.} \end{cases} \end{equation} \begin{equation} \psi(s) := \begin{cases} 1 & 0 \le x < \frac{1}{2}\\ -1 & \frac{1}{2} \le x < 1\\ 0 & \mbox{інакше.} \end{cases} \end{equation} Фактично, $[\ 8\ 4\ 1\ -1\ ]$ кодує $[\ 9\ 7\ 3\ 5\ ]$ у такому вигляді: \begin{equation} 8\phi(x) + \psi(x) + 4\phi(x-1) - \psi(x-1).\label{dec9735pre} \end{equation} УВАГА! Зараз ми зробимо невеличкий виверт. Ми розглядатимемо зображення як кусково-неперервні функції на проміжку $[0, 1)$. Наприклад, одновимірне зображення $[\ 9\ 7\ 3\ 5\ ]$ виглядатиме так:
Отже, тепер ми вже не можемо сказати, що \eqref{dec9735pre} предсталяє $[\ 9\ 7\ 3\ 5\ ]$. Тепер нам потрібно відмаштабувати наші функції. У загальному випдаку для зображення розміру $2^n$ на першому кроці треба використовувати $\phi_n^i(x) = \phi(2^nx-i)$ і $\psi_n^i(x) = \psi(2^nx-i)$.

Базис Хаара має чудову властивість ортогональності, тобто всі базисні функції - $\phi_0^0, \psi_0^0, \phi_1^1, \psi_1^1, \dots$ ортогональні. Нам також хотілось би нормалізувати базис.

Означення
Базисна функція $u(x)$ нормалізована, якщо $$ \langle u|u\rangle = \int_0^1 u(x)u(x)dx = 1.$$
Базис Хаара можна нормалізувати підправивши наявне визначення базисних функцій так: \begin{align} \phi_i^j(x):=2^{j/2}\phi(2^jx-i),\\ \psi_i^j(x):=2^{j/2}\psi(2^jx-i). \end{align} В такому базисі ненормалізовані коефіцієнти $[\ 6\ 2\ 1\ -1\ ]$ стають нормалізованими $$[\quad 6\quad 2\quad \frac{1}{\sqrt{2}}\quad -\frac{1}{\sqrt{2}}\quad ].$$ Ніж очислювати ненормалізовані коейіцієнти і потім нормалізовувати їх, можна включити нормалізацію в алгоритм.

Реалізація

Наступний код створює масив 512 символів завдовжки і заповнює його значеннями функції Веєрштраса на проміжку $[-1, 1]$:
#include <algorithm>
#include <cassert>
#include <vector>
#include <fstream>

double Weierstrass(double x) {
    const int summandsCount = 100;
    const double a = 0.7;
    const double b = 11;
    double aRaised = 1;
    double bRaised = 1;
    double res = 0;
    for (int i = 0; i < summandsCount; ++i) {
        aRaised *= a;
        bRaised *= b;
        res += aRaised * cos(bRaised * atan(1) * 4 * x);
    }
    return res;
}

const auto kLevelsCount = 9;
const auto kN = 1 << 9;
const auto kSqrt2 = sqrt(2);

void save(const std::string& name, const std::vector<double>& data) {
    std::ofstream out(name);
    for (int i = 0; i < (int)data.size(); ++i) out << data[i] << std::endl;
}

void DecompositionStep(std::vector<double>& data, std::vector<double> storage, int size) {
    for (int i = 0; i < size; i = i + 2) {
        storage[i / 2] = (data[i] + data[i + 1]) / kSqrt2;
        storage[size / 2 + i / 2] = (data[i] - data[i + 1]) / kSqrt2;
    }
    std::copy_n(storage.cbegin(), size, data.begin());
}

void Decomposition(std::vector<double>& data, std::vector<double>& storage) {
    assert(data.size() == storage.size() && "At the end lod levels take 1 position and differences (data.size() - 1) positions");
    std::transform(data.cbegin(), data.cend(), data.begin(), [s = data.size()](double y) { return y / sqrt(s); });
    auto currSize = (int)data.size();
    while (currSize > 1) {
        DecompositionStep(data, storage, currSize);
        currSize /= 2;
    }
}

int main()
{
    std::vector<double> data(kN);
    const auto interval = std::make_pair(-0.5, 0.5);
    const auto length = interval.second - interval.first;
    for (int i = 0; i < kN; ++i) {
        double x = interval.first + length * i / (kN - 1);
        data[i] = Weierstrass(x);
    }
    
    std::vector<double> storage(data.size());
    Decomposition(data, storage);
    save("compressed.dat", data);
    return 0;
}
А цей код MATLAB'у зчитує утворений файл і відновлює початкові дані починаючи з найгрубішого рівня. По мірі відновлення він записує зображення кожного рівня в файл анімації.
h = figure;
axis tight manual
filename = 'haarWeierstrass.gif';
hold on

xStart = -1;
xEnd = 1;

load 'C:\Yola\VC\Test17\Test17\compressed.dat' compressed -ascii
L = length(compressed);
lodLevelsCount = log(L)/log(2) + 1;
storage = zeros(size(compressed));
for i = 1:lodLevelsCount
    lodStart = 1;
    detailStart = 2^(i-1) + 1;
    localY = compressed(1:detailStart - 1);
    x = linspace(xStart, xEnd, length(localY));
    plot(x, localY);
    title(['Рівень деталізації ' num2str(i-1) ', ' num2str(2^(i-1)) ' точок'])
    drawnow
    frame = getframe(h);
    im = frame2im(frame); 
    [imind,cm] = rgb2ind(im,256); 
    if i == 1 
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf); 
    else 
        imwrite(imind,cm,filename,'gif','WriteMode','append'); 
    end 
    
    if (i < lodLevelsCount)
        % відновлення наступного рівня
        m = 2^((i-1)/2);
        for j = 1:detailStart - 1
            storage(2*j - 1) = compressed(j) + m * compressed(detailStart - 1 + j);
            storage(2*j) = compressed(j) - m * compressed(detailStart - 1 + j);
        end
        compressed(1:2*(detailStart-1)) = storage(1:2*(detailStart-1));
        pause(1);
        clf();
    end
end

четвер, 24 травня 2018 р.

Розв'язки до вправ з "An Introduction to the Analysis of Algorithms" Роберта Седжвіка. Розділ II.

Вправа 2.2
Зовнішній цикл робить $N_{max}$ ітерацій. Внутрішній $\frac{N_{max}(N_{max}-1)}{2}$. На кожній ітерації внутрішнього циклу виконується два додавання і одне ділення. В підсумку $3\frac{N_{max}^2(N_{max}-1)}{2}$.
Вправа 2.3
int Cn(int n) {
    if (n == 1) return 2;
    return (1 + 1. / n) * Cn(n - 1) + 2;
}
На кожному рівні рекурсії ми маємо дві операції з $C_n$ - множення і додавання. Всього таких рівнів $N-1$, отже $2(n-1)$.
Вправа 2.6
\[ a_n = q \operatorname{Fib}(n-1) + p \operatorname{Fib}(n-2) \]
Вправа 2.7
\[ a_n = q \operatorname{Fib}(n-1) + p \operatorname{Fib}(n-2) + r \operatorname{Fib}(n-2) \] p, q, p + q + r, p + 2q + r, 2p + 3q + 2r, 3p + 5q + 3r, 5p + 8q + 5r.
Вправа 2.8
\[ f^*(a_{n-1}, a_{n-2}) = f(a_{n-1}, a_{n-2}) - f(0,0) \] Маємо, що $f^*$ однорідна. Нехай \begin{align*} u_n &= f^*(u_{n-1}, u_{n-2}),\quad \text{для}\ n > 1, u_0 = 1, u_1 = 0,\\ v_n &= f^*(v_{n-1}, v_{n-2}),\quad \text{для}\ n > 1, v_0 = 0, v_1 = 1. \end{align*} Тоді \[ a_n = a_0 u_n + a_1 v_n + f(0,0)\operatorname{Fib}(n-2). \]
Вправа 2.9
\[ a_n = \frac{n}{n+2}a_{n-1} = \prod_{k=1}^n \frac{k}{k+2} = \frac{1\cdot 2\cdot \cdots n}{3\cdot 4 \cdots (n+2)} = \frac{1\cdot 2}{(n+1)(n+2)}. \]
Вправа 2.10
\[ a_n = 1 + \underbrace{(-1) + 2}_1 + \underbrace{(-3) + 4}_1 \cdots, \] \[ a_n= \begin{cases} 1 + \frac{n}{2} & n\ \text{парне}\\ 1 - \frac{n}{2} & n\ \text{непарне} \end{cases} \]
Вправа 2.11
\begin{align*} n(n-1)a_n &= (n-1)(n-2)a_{n-1} + 2(n-1)\\ (n-1)(n-2)a_{n-1} &= (n-2)(n-3)a_{n-4} + 2(n-3)\\ \dots&\\ 2\cdot 1\cdot a_2 &= 0\cdot(-1)\cdot 1 + 2\cdot(n-(n-2)-1) \end{align*} Тепер складемо всі рівності. \[ n(n-1)a_n = 2\sum_{k=1}^{n-1}k = n(n-1) \] Отримуємо: \[ a_n = 1 \]
Вправа 2.12
\begin{align*} 2^{-n}a_n &= 2^{-(n-1)}a_{n-1} + 2^{-n}\\ \dots &\\ 2^{-2}a_2 &= 2^{-1}\cdot 1 + 2^{-2}\\ \end{align*} Тепер складемо всі рівності. \[ 2^{-n}a_n = 2^{-1}\cdot 1 + \sum_{k=2}^{n}2^{-k} \] Отримуємо: \[ a_n = 2^{n-1} + \sum_{k=2}^n2^{n-k} = 2^{n-1} + \sum_{k=0}^{n-2}2^k = 2^n - 1. \]
Вправа 2.13
Тут нам треба звернути увагу на те, що на відміну від теореми 2.1, нулю дорівнює не $a_1$, а $a_0$. Тому сума починається з $j = 0$: \begin{align*} a_n &= 1 + \sum_{j=0}^{n-1}\frac{j+1}{j+2}\cdots\frac{n}{n+1}\\ &= 1 + \sum_{j=0}^{n-1}\frac{j+1}{n+1} = 1 + \frac{1}{n+1}\sum_{j=0}^{n-1} (j+1)\\ &= \frac{1}{n+1}\sum_{j=1}^{n+1}j = \frac{n+2}{2}. \end{align*}
Вправа 2.14
Тут сума починається з $t+1$, також присутній член $a_t$, бо не сказано, що він дорівнює нулю. \begin{equation} a_n = y_n + \sum_{j=t+1}^{n-1} y_{j} x_{j+1} \dots x_n + a_tx_{t+1}\dots x_n.\label{eq:rr_at} \end{equation}
Вправа 2.15
\begin{align*} a_n &= 2 \sum_{j=1}^{n-1} \frac{j+2}{j+1} \dots \frac{n+1}{n} + 2\\ &= 2 \sum_{j=1}^{n-1}\frac{n+1}{j+1} + 2\\ &= 2(n+1)\left(\sum_{j=1}^{n-1}\frac{1}{j+1}+\frac{1}{n+1}\right)\\ &= 2(n+1)(H_{n+1}-1). \end{align*} Під час перевірки можна скористатись тим, що $nH_n - n = \sum_{k=1}^{n-1}H_k$. \begin{align*} 2n\sum_{k=1}^nH_k &= 2(n+1)\sum_{k=1}^{n-1}H_k + 2n\\ nH_n &= \sum_{k=1}^{n-1}H_k + n \end{align*}
Вправа 2.16
Скориставшись теоремою 2.1: \begin{align*} a_n = 12 H_n + 12 \sum_{j=5}^{n-1} H_j \frac{(j-3)(j-2)(j-1)j}{n(n-1)(n-2)(n-3)}, \end{align*} або ми можемо використати сумувальний множник $(n-1)(n-2)(n-3)$: \begin{align*} n(n-1)(n-2)(n-3)a_n &= (n-1)(n-2)(n-3)(n-4)a_{n-1}\\ &\phantom{=}+ 12 n(n-1)(n-2)(n-3) H_n,\\ (5\cdot 4\cdot 3\cdot 2)a_5 &= 12 (5\cdot 4\cdot 3\cdot 2) H_5. \end{align*} Сумуючи отримуємо той самий вислід.
Вправа 2.17
Зауважимо, що кожен 2-вузол містить одну вершину, а кожен 3-вузол містить дві вершини, тобто розпадається на два 2-вузли. \begin{align*} A_N &= A_{N-1}\left(\frac{N-6}{N}\right) +2 \\ &= 2 + 2 \sum_{j=1}^{N-1}\left(\frac{j-5}{j+1}\right)\dots\left(\frac{N-7}{N-1}\right)\\ \end{align*} Тут і у попередній вправі $x_n$ дуже схожі. І виконуючи теорему 2.1 нам може захотітись дати таку відповідь: \[ A_N = 2 + 2 \sum_{j=1}^{N-1} \frac{(j-5)(j-4)(j-3)(j-2)(j-1)j}{N(N-1)(N-2)(N-3)(N-4)(N-5)}, \] але це не правильно, бо ми не можемо ділити на нуль. Ми мусимо використати \eqref{eq:rr_at}. Отже, відповіддю буде \begin{align*} A_N &= 2 + 2 \sum_{j=6}^{N-1} \frac{(j-5)(j-4)(j-3)(j-2)(j-1)j}{N(N-1)(N-2)(N-3)(N-4)(N-5)}\\ &\phantom{=} + A_5\frac{0\cdot1\cdot2\cdot3\cdot4\cdot5}{N(N-1)(N-2)(N-3)(N-4)(N-5)},\quad \text{для }\ N > 5. \end{align*} Нам навіть не потрібно вручну обчислювати $A_5$, бо його потрібно множити на $0$.
Щоб краще уявити собі кількість 2-вузлів і перевірити формулу, можна запустити таку програму:
float tree23(int n) {
    if (n == 0) return 0;
    float prev = tree23(n - 1);
    return prev - 2 * prev / n + 2 * (1 - 2 * prev / n);
}

float tree23_sum(int n) {
    float sum = 0.;
    for (int j = 6; j < n; ++j) {
        float num = 1.f, denom = 1.f;
        for (int k = 0; k < 6; ++k) {
            num *= j - k;
            denom *= n - k;
        }
        sum += num / denom;
    }
    return 2 + 2 * sum;
}

int main() {
    for (int i = 6; i < 30; ++i) {
        cout << setw(6) << setprecision(3) << fixed << tree23(i) << " " << tree23_sum(i) << endl;
    }
    return 0;
}

 2.000 2.000
 2.286 2.286
 2.571 2.571
 2.857 2.857
 3.143 3.143
 3.429 3.429
 3.714 3.714
 4.000 4.000
 4.286 4.286
 4.571 4.571
 4.857 4.857
 5.143 5.143
 5.429 5.429
 5.714 5.714
 6.000 6.000
 6.286 6.286
 6.571 6.571
 6.857 6.857
 7.143 7.143
 7.429 7.429
 7.714 7.714
 8.000 8.000
 8.286 8.286
 8.571 8.571
Вправа 2.18
\begin{align*} a_1 &= \frac{1 + 0\cdot a_0}{1+1\cdot a_0} = \frac{\operatorname{Fib}_1 + \operatorname{Fib}_0\cdot a_0}{\operatorname{Fib}_2 + \operatorname{Fib}_1\cdot a_0}\\ a_2 &= \frac{1 + 1\cdot a_0}{2+1\cdot a_0} = \frac{\operatorname{Fib}_2 + \operatorname{Fib}_1\cdot a_0}{\operatorname{Fib}_3 + \operatorname{Fib}_2\cdot a_0}\\ \dots&\\ a_n &= \frac{\operatorname{Fib}_n + \operatorname{Fib}_{n-1}\cdot a_0}{\operatorname{Fib}_{n+1} + \operatorname{Fib}_n\cdot a_0}\\ \end{align*} Зі збільшенням $n$ перші доданки починають домінувати і ми отримуємо $\frac{\operatorname{Fib}_{n}}{\operatorname{Fib}_{n+1}}$, тоді як ми знаємо, що $\lim_{n\to \infty}\frac{\operatorname{Fib}_{n+1}}{\operatorname{Fib}_{n}} = \varphi$ - золотий перетин, а $\frac{1}{\varphi}$ і є потрібним результатом, з яким $\lim_{n\to \infty} b_n = 0$.
Вправа 2.19
На рисунку праворуч показани графіки функцій $ y = \cos(x)$ і $x = \cos(y)$ або інакше $y = \arccos(x)$. Точка перетину має абсцису $0.73908\dots$.
Вправа 2.20
Щоб зрозуміти чому формула працює можна розглянути такий її варіант: $$2a_na_{n-1} = a_{n-1}^2+b$$ Якщо $b$ від'ємне, то нема причин очікуватись збіжність $a_n$ до $b$.
Вправа 2.21
  1. Спершу доведемо, що $a_n \in \Theta(1/n)$. \begin{align*} \frac{1}{a_n} &= \frac{1}{a_{n-1}}\frac{1}{1-a_{n-1}}&\\ \frac{1}{a_n} &< \frac{1}{a_{n-1}}\frac{n-1}{n-2}&\text{використали}\ a_n < \frac{1}{n}\\ a_n &> \frac{1}{n-1}a_3&\text{де}\ a_3 = 39/256. \end{align*} Отже, для достатньо великих $n$, $a_n > \frac{1}{n-1}$.
  2. Тепер оцінимо $c$:
        float a = 0.5f;
        for (int i = 1; i < 1000; ++i) {
            a = a * (1 - a);
            cout << setw(15) << a << " " << a * i << endl;
        }
    
    Така програма видає нам:
    ...
        0.000996334 0.991353
        0.000995342 0.99136
        0.000994351 0.991368
        0.000993362 0.991376
        0.000992376 0.991383
    
    Отже, схоже що $c$ дорівнює $1$.
  3. З пункту 1 і з книги ми знаємо, що $$\frac{1}{n-1} < a_n < \frac{1}{n}$$ отже $na_n \to 1$.
Вправа 2.22
  1. На початку $\sin(x)$ зростає повільніше ніж $x$, тобто $x > \sin(x)$ окрім як в $0$. Отже послідовність $a_n$ обмежена спадна функція, тобто вона має границю. $\sin(x)$ - неперервна, отже ця границя в стаціонарній точці, єдина така точка в $[0,1]$ це $0$.
  2. Тепер доведемо, що $a_n = O(1/\sqrt{n})$.
Вправа 2.23
Розглянемо $$a_n = \alpha + \varepsilon_{n} = f(a_{n-1}) = f(\alpha + \varepsilon_{n-1}).$$ Використаємо розклад в ряд Тейлора: $$f(\alpha + \varepsilon_{n-1}) = f(\alpha) + f'(\alpha) * \varepsilon_{n-1} + f''(\alpha) * \varepsilon_{n-1} / 2!.$$ Зауважимо, що $f(\alpha) = \alpha$. \begin{equation} \varepsilon_{n} = f'(\alpha) * \varepsilon_{n-1} + f''(\alpha) * \varepsilon_{n-1} / 2!\label{eq:cnvrg_rate_drv} \end{equation} Отже, якщо $f'(\alpha) > 1$, то ряд розбігається.
Вправа 2.24
Використовуючи \eqref{eq:cnvrg_rate_drv} маємо квадратну збіжність якщо $f'(\alpha)$ рівне $0$ з коефіцієнтом $f''(\alpha)/2$. Аналогічно можна отримати коефіцієнти для $f'(\alpha) \ne 0$.
Вправа 2.25
Для того, щоб збагнути як утворити необхідне рекурентне відношення розглянемо таке рекурентне відношення: \begin{equation} a_n - 3*a_{n-1} = 4*(a_{n-1} - 3*a_{n-2}),\quad a_0 = 0, a_1 = 1.\label{eq:rr_25_43} \end{equation} Якщо пройти рекурсивно до базового випадку, маємо: \begin{align*} a_n - 3 a_{n-1} &= 4^{n-1}\\ 3^1(a_{n-1} - 3 a_{n-2}) &= 3^1 4^{n-2}\\ \dots &\\ 3^i(a_{n-i} - 3 a_{n-i-1}) &= 3^i 4^{n-i-1}\\ \dots &\\ 3^{n-1}(a_1 - 3 a_0) &= 3^{n-1} \end{align*} Складаємо і скорочуємо: \begin{align*} a_n = \sum_{i=1}^n 4^{n-i}3^{i-1} \end{align*} Зауважимо, що $4^n = 1\dot 4^{n-1} + 3\cdot 4^{n-1}$ або ж $3^{n-1}\dot 4^0 + 3^n\cdot 4^0 = 4\cdot 3^n$. Застосовуючи останню формулу маємо для $n=5$: $$ \underbrace{4^4\cdot 3^0 + \underbrace{4^3\cdot3 + \underbrace{4^2\cdot3^2+ \underbrace{4\cdot3^3 + \underbrace{3^4\cdot 4^0 + 3^5\cdot 4^0 }_{3^4\cdot 4} }_{3^3\cdot 4^2} }_{3^2\cdot 4^3} }_{3^1\cdot 4^4} }_{4^5} $$ Зверніть увагу, що тут ми додали $3^5$, отже формула така: $$\sum_{i=1}^n 4^{n-i}3^{i-1} = 4^n - 3^n.$$ Тобто розв'язок до \refeq{eq:rr_25_43} такий $a_n = 4^n - 3^n$. Гм.. це маже те, що нам треба! А ось і наша відповідь: $$a_n - 3*a_{n-1} = 4*(a_{n-1} - 3*a_{n-2})+2^{n-1}.$$
long a(int n) { 
    switch (n) {
    case 0: return 1;
    case 1: return 3;
    default: return powl(4, n) - powl(3, n) + powl(2,n);
    }
}

    for (int n = 2; n < 25; ++n) {
        long lhs = a(n) - 3 * a(n - 1);
        long rhs = 4 * (a(n - 1) - 3 * a(n - 2)) + powl(2,n-1);
        cout << setw(2) << n << ": " << setw(10) << lhs << " " << rhs << endl;
    }

 2:          2 2
 3:         12 12
 4:         56 56
 5:        240 240
 6:        992 992
 7:       4032 4032
...
Вправа 2.26
  1. Ввести $b_n = a_n - a_{n-1}$, це дасть нам однорідне рекурентне відношення.
  2. Розписати в стовпчик усі $a_i - a_{i-1}$ і скласти, завдяки телескопуванню, ми зможемо отримати формулу для $a_n$, спрощення цієї формули може бути досить складним.
Вправа 2.27
Характеристичний многочлен такий: $$\beta^2 - 5\beta + 6 = 0.$$ Коренями є $2, 3$. Значить маємо розв'язки у вигляді $c_12^n + c_23^n$. Звідси бачимо, що для розв'язку $2^n$ $c_1 = 1, c_2 = 0$, а розв'язок $2^n - 1$ отримати не можливо.
Вправа 2.28
Корені $2, -i, i$.
  1. $a_0 = a_1 = a_2 = 0$.
  2. $a_0 = a_1 = a_2 = 1$.
  3. $a_0 = i, a_1 = -1, a_2 = i$.
Вправа 2.29
Коренями є $1+\sqrt{5}, 1-\sqrt{5}$. Відповідь $\frac{(1+\sqrt{5})^n - (1-\sqrt{5})^n}{2\sqrt{5}}$.
Вправа 2.30
Перепишемо рекурентне відношення у вигляді $$a_n - a_{n-1} = a_{n-1} - a_{n-2}.$$ Ми бачимо, що зі зростанням $n$ на одиницю $a$ зростає на сталу. Звідси нам легко знайти відповіді.
  1. $a_n = n$.
  2. $a_n = 1$.
Вправа 2.31
Тут корені це $\frac{1\pm i\sqrt{3}}{2}$. Отже, маємо два рівняння: \begin{cases} c_1 + c_2 &= 0\\ c_1\frac{1 + i\sqrt{3}}{2} + c_2\frac{1 - i\sqrt{3}}{2} &= 1 \end{cases} Що дає $c_1 = -c_2 = \frac{1}{i\sqrt{3}}$. І, власне, розв'язок: $$a_n = \frac{1}{i\sqrt{3}}\left(\frac{1 + i\sqrt{3}}{2}\right)^n -\frac{1}{i\sqrt{3}}\left(\frac{1 - i\sqrt{3}}{2}\right)^n$$
Вправа 2.32
Корені тут $\frac{1}{2}, \frac{1\pm\sqrt{3}}{2}$. \begin{cases} c_1+c_2+c_3 &= 0\\ c_1\frac{1}{2}+c_2 \frac{1+\sqrt{3}}{2}+c_3 \frac{1-\sqrt{3}}{2} &= 1\\ c_1\left(\frac{1}{2}\right)^2+c_2 \left(\frac{1+\sqrt{3}}{2}\right)^2+c_3 \left(\frac{1-\sqrt{3}}{2}\right)^2 &= 2 \end{cases}
Вправа 2.33
$$a_n = 4(a_{n-2}-a_{n-4}).$$ При $a_0 = 1, a_1 = -1, a_2 = 2, a_3 = -2$. Для парних маємо $1,2,4,8,\dots$ і для непарних те саме тільки з мінусом.
Вправа 2.34
>> p = [1 -1 -1 -1];
>> r = roots(p)
r =

   1.83929 + 0.00000i
  -0.41964 + 0.60629i
  -0.41964 - 0.60629i

>> A = [r.^0'; r.^1'; r.^2'];
>> b = [0 0 1]';
>> c = A\b
c =

   0.182804 + 0.000000i
  -0.091402 - 0.340547i
  -0.091402 + 0.340547i

>> F20_exact = r.^19' * c
F20_exact = 1.9513e+004 + 6.5683e-013i
>> F20_approx = r(1)^19' * c(1)
F20_approx = 1.9513e+004 + 6.5683e-013i
Вправа 2.35
Домножити обидва боки на $(n-2)!$. Маємо Фібоначчі: $$n!a_n = (n-1)!a_{n-1} + (n-2)!a_{n-2}.$$ Отже, $a_n = F_n / n!$.
Вправа 2.36
Входження $t_i$ означає відсутність $a_i$ і $a_{i+1}$. Всі одночлени - унікальні.
Вправа 2.37
Характеристичний многочлен для послідовності $b_n$ такий: $x^2 - \frac{1}{2}x - \frac{1}{2} = 0$. Корені - $1, -\frac{1}{2}$. \begin{cases} c_1+c_2 &= 0\\ c_1 + c_2\left(-\frac{1}{2}\right) & = 1 \end{cases} Звідси, $c_1 = -c_2 = \frac{2}{3}$. І відповідь $$b_n = \frac{2}{3} - \frac{2}{3}\left(-\frac{1}{2}\right)^n.$$ $$a_n = 2^{\frac{2}{3} - \frac{2}{3}\left(-\frac{1}{2}\right)^n}.$$
Вправа 2.38
Введемо заміну змінних $b_n = a_n^2$. Маємо: $$b_n = b_{n-1} + 1.$$ Це відношення має розв'язок $b_n = n$. Відповідно, $a_n = \sqrt{n}$.
Вправа 2.39
Для отримання формули в для $a_0 = 3, a_0 = 4$ необхідно підставити конкретні значення в формулу наведену в книзі. У випадку коли $a_0 = \frac{3}{2}$ послідовність "хаотично" рухається в $[-2, 2 ]$.
Вправа 2.40
Для великих $\epsilon$: \begin{align*} a_n &\approx \left(\frac{1}{2}(2+\epsilon + \sqrt{(2+\epsilon)^2 - 4}\right)^{2^n}\\ &= \left(\frac{1}{2}(2+\epsilon+\sqrt{4+4\epsilon+\epsilon^2 -4})\right)^{2^n}&\quad\text{в підкорінневому виразі домінує}\ \epsilon^2\\ &\approx \left(1 + \epsilon\right)^{2^n}\\ &\approx \epsilon^{2^n}. \end{align*}
Вправа 2.41
Якщо $b_n = f\cdot a_n + g$, тоді $a_n = \frac{b_n-g}{f}$. \begin{align*} \frac{b_n - g}{f} &= \alpha\left(\frac{b_{n-1} - g}{f}\right)^2 + \beta\frac{b_{n-1} - g}{f} + \gamma\\ b_n &= \frac{\alpha}{f}(b_{n-1} - 2b_{n-1}g + g^2)+\beta b_{n-1}-\beta g + \gamma f + g\\ b_n &= \underbrace{\frac{\alpha}{f}}_{=1}b_{n-1}^2 - \underbrace{\left(\frac{2\alpha g}{f}-\beta\right)}_{=0}b_{n-1} + \underbrace{\frac{\alpha g^2}{f} - \beta g + \gamma f + g}_{=-2} \end{align*} З цього ми миожемо зробити такі висновки: $f = \alpha$, $g = \frac{\beta}{2}$. А те, що $\frac{\alpha g^2}{f} - \beta g + \gamma f + g = -2$ дозволяє нам обрати значення для $\alpha, \beta, \gamma$, ми маємо два ступені свободи тут.
Наприклад, $\beta = 2, \alpha = 1, \gamma = -2$. $$a_n = a_{n-1}^2 + 2 a_{n-1} - 2$$ Отже, $f = 1$, $g = 1$ і $b_n = a_b + 1$. Перевіримо: \begin{align*} a_n + 1 &= (a_{n-1} + 1)^2 - 2\\ & = a_{n-1}^2 + 2 a_{n-1} + 1 - 2 \end{align*} Що дає $$a_n = a_{n-1}^2 + 2 a_{n-1} -2.$$
Вправа 2.42
\begin{align*} a_n &= 2 a_{n-1} \sqrt{1-a_{n-1}^2}&\quad n>0, a_0=\frac{1}{2}\\ a_n^2 &= 4 a_{n-1}^2 (1-a_{n-1}^2)\\ b_n &= 4b_{n-1}-4b_{n-1}^2&\quad b_n=a_n^2 \end{align*} Тут ми можемо використати результат попередньої вправи з $\alpha = -4$, $\beta = -4$, $\gamma = 0$. Тоді $c_n = -4b_n+2$. Маємо $$c_n = c_{n-1}^2 - 2.$$ А це вже алгоритм розподілення регістрів з книги.
Вправа 2.43
Ми можемо переписати нашу рекуренцію так: $$a_n = \underbrace{\frac{\alpha a_{n-1}}{\gamma a_{n-1} + \delta}}_{a'_{n}} + \underbrace{\frac{\beta}{\gamma a_{n-1} + \delta}}_{a''_{n}}.$$ Тепер ми можемо працювати з двома рекуренціями і скласти їх результати. \begin{align*} a'_n &= \frac{\alpha a'_{n-1}}{\gamma a'_{n-1} + \delta}\\ \frac{1}{a'_n} &= \frac{\gamma}{\alpha} + \frac{\delta}{\alpha}\frac{1}{a'_{n-1}}&\quad b'_n = \frac{1}{a'_n}\\ b'_n &= \frac{\gamma}{\alpha} + \frac{\delta}{\alpha}b'_{n-1}\\ \end{align*} Тепер подібно до теореми 2.1: $$b'_n = \frac{\gamma}{\alpha}\sum_{i=0}^{n-1}\left(\frac{\delta}{\alpha}\right)^i + \left(\frac{\delta}{\alpha}\right)^nb'_0.$$ Друга рекуренція майже як в прикладі з книги: \begin{align*} a''_n &= \frac{\beta}{\gamma a''_{n-1} + \delta}&\quad a''_n = \frac{b''_{n-1}}{b''_n}\\ b''_n &= \frac{\delta}{\beta}b''_{n-1} + \frac{\gamma}{\beta}b''_{n-2}\\ \end{align*} Тут ми знайдемо два корені і нам буде потрібно вибрати значення коефіцієнтів залежно початкових умов. Тобто, початкова умова $a_0 = a'_0 + a''_0 = 1$ має виконуватись.
Вправа 2.44
Розпишемо перші члени послідовності $(a_n)$: \begin{align*} a_0 &= 1\\ a_1 &= \frac{1}{s_1+t_1}\\ a_2 &= \frac{s_1+t_1}{s_2s_1+s_2t_1 + t_2}\\ a_3 &= \frac{s_2s_1+s_2t_1 + t_2}{s_3s_2s_1+s_3s_2t_1 + s_3t_2 + t3s_1 + t_3t_1}\\ a_4 &= \frac{s_3s_2s_1+s_3s_2t_1 + s_3t_2 + t_3s_1 + t_3t_1}{s_4s_3s_2s_1+s_4s_3s_2t_1 + s_4s_3t_2 + s_3t_2s_1 + s_3t_2t_1 + t_4s_2s_1+t_4s_2t_1 + t_4t_2}\\ \end{align*} Тут ми вже бачимо взірець за яким можна будувати наступні члени. Ми можемо записати послідовність трошки інакше: $$a_n = \frac{s_{n-1}\operatorname{den}{a_{n-2}} + t_{n-1}\operatorname{nom}{a_{n-2}}}{s_{n}\operatorname{den}{a_{n-1}} + t_{n}\operatorname{nom}{a_{n-1}}}$$ Зверніть увагу на послідовність опрацьовану в вправі 2.36, вона дуже подібна. Якщо $a_n = \frac{b_{n}}{b_{n+1}}$, то $$b_n = s_{n-1}b_{n-1}+t_{n-1}b_{n-2},\quad n>1,\ b_0 = b_1 = 1.$$
Вправа 2.45
Найперше ми маємо знайти форму $a_n$, яка дасть нам входження $n^3$ в $a_n - \left(2\sum_{j=0}^{n-1}\right)/n$. Для цього ми можемо вибрати $n(n-1)(n-2)$. \begin{align*} \sum_{k=0}^n k(k-1)(k-2) & = \sum_{k=0}^n (k^2 - 3k^2 + 2k)\\ &= \frac{n^2(n+1)^2}{4} - \frac{n(n+1)(2n+1)}{2} + n(n+1)\\ &= \frac{1}{4}(n^4 - 2n^3 - n^2 + 2n). \end{align*} Перевіримо: $$a_3 = (3\cdot 2 \cdot 1) = 6 = \frac{1}{4}(81 - 54 - 9 + 6).\quad\text{збіг!}$$ Отже, \begin{align*} a_n - \left(2\sum_{j=0}^{n-1}\right)/n & = n(n-1)(n-2)\\ &\phantom{=}- \frac{1}{2}\big((n-1)^4 - 2(n-1)^3 - (n-1)^2 + 2(n-1)\big)/n \end{align*} Тепер маючи компонент з $n^3$ ми можемо за допомогою нескладної, але нудної алгебри зкомпонувати його з уже відомими розв'язками, щоб отримати $f(n) = n^3$.
Вправа 2.46
Розв'яжіть рекурентне відношення для швидкого сортування медіани трьох використовуючи репертуарний метод. $$c_n = n+1+\sum_{k=1}^n\frac{(n-k)(k-1)}{\binom{n}{3}}(c_{k-1} + c_{n-k}),\quad n>3.$$ Перепишемо: \begin{align*} c_n &= n+1+\frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(n-(k+1))((k+1)-1))c_{k}\\ &=n+1+\frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk- k^2 - k)c_{k} \end{align*} Заповнимо репертуарну таблицю:
  • Для $c_n = 1$:
\begin{align*} &\phantom{=\ }1 - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk - k^2 - k)\\ &= 1 - \frac{2}{\binom{n}{3}}\left(\frac{n^2(n-1)}{2} - \frac{(n-1)n(2n-1)}{6} + \frac{n(n-1)}{2}\right)\\ &= 1 - 2\frac{3\cdot 2}{n(n-1)(n-2)}\left(\frac{1}{6}(n^3-3n^2+2n)\right)\\ &= -1. \end{align*}
  • Для $c_n = n$:
\begin{align*} &\phantom{=\ }n - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk^2 - k^3 - k^2)\\ &= n - \frac{2}{\binom{n}{3}}\left((n-1)\frac{(n-1)n(2n-1)}{6} - \frac{(n-1)^2n^2}{4}\right)\\ &= n - 2\frac{3\cdot 2}{n(n-1)(n-2)}\frac{1}{12}n(n-1)^2(n-2)\\ &= 1. \end{align*}
  • Для $c_n = (n-1)^{\underline t}$:
\begin{align*} &\phantom{=\ }(n-1)^{\underline t} - \frac{2}{\binom{2}{3}}\sum_{1\lt k\lt n}(k-1)(n-k)(k-2)^{\underline t}\\ &= (n-1)^{\underline t} - \frac{2}{\binom{n}{3}}\sum_{1\lt k\lt n}(n-k)(k-1)^{\underline {t+1}}\\ &= (n-1)^{\underline t} - \frac{2(t+1)!}{\binom{n}{3}}\sum_{1\lt k\lt n}(n-k)\binom{k-1}{t+1}\\ &= (n-1)^{\underline t} - \frac{2(t+1)!}{\binom{n}{3}}\sum_{1\lt k\lt n}\binom{n}{t+3}\\ &= (n-1)^{\underline t} - \frac{2(t+1)!}{\binom{n}{3}}\binom{n}{t+3}\\ &= (n-1)^{\underline t} - \frac{12}{(t+2)(t+3)}(n-3)^{\underline{t}}\\ \end{align*} Бррр... Лінійний час видобути не так просто. І не дивно, ми ж очікуємо, що ця рекурсія виконується за $n \lg n$.
  • Для $c_n = \sum_{k=1}^n H_n$:
\begin{align*} &\phantom{=\ }n - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}((n-1)k - k^2)\sum_{i=1}^k H_i\\ \end{align*} Не завершено.
Вправа 2.47
Ми бачимо, що $$\frac{2}{n+1} < a_n < \frac{2}{n}.$$ Тепер підставимо значення для $a_{n-1}$ в рекуренцію: \begin{align*} \frac{2}{n+\frac{2}{n}} &< a_n < \frac{2}{n + \frac{2}{n-1}}\\ \frac{2n}{n^2+2} &< a_n < \frac{2(n-1)}{n^2 - n + 2}\\ \end{align*} Але $\frac{2}{n+1} < \frac{2n}{n^2+2}$ і $\frac{2(n-1)}{n^2 - n + 2} < \frac{2}{n}$, тобто завдяки бутстрепу вдалось звузити діапазон можливих начень $a_n$.
Вправи з 48 і далі нерозв'язані.