پرش به محتویات

کسرهای مسلسل (Continued fraction)

کسر مسلسل (Continued fraction) یه راه برای نمایش یه عدد حقیقی به صورت یه دنباله‌ی خاص و همگرا از عددهای گویاست. این کسرها توی برنامه‌نویسی رقابتی خیلی به درد می‌خورن، چون حساب کتاب کردنشون راحته و می‌تونی خیلی موثر ازشون استفاده کنی تا بهترین تقریب گویای ممکن برای یه عدد حقیقی رو پیدا کنی (یعنی از بین همه‌ی کسرهایی که مخرجشون از یه حدی کوچیکتره، بهترینش رو پیدا کنی).

تازه، کسرهای مسلسل یه ربط خیلی نزدیکی هم به الگوریتم اقلیدس دارن که باعث میشه توی یه عالمه مسئله‌ی نظریه اعداد به کار بیان.

نمایش کسر مسلسل

تعریف

فرض کن $a_0, a_1, \dots, a_k \in \mathbb Z$ باشن و $a_1, a_2, \dots, a_k \geq 1$. اون‌وقت، به عبارتِ:

$$r=a_0 + \frac{1}{a_1 + \frac{1}{\dots + \frac{1}{a_k}}},$$

میگن نمایش کسر مسلسل برای عدد گویای $r$ و خلاصه‌اش رو اینجوری می‌نویسن: $r=[a_0;a_1,a_2,\dots,a_k]$.

مثال

فرض کن $r = \frac{5}{3}$ رو داریم. دو تا راه برای نشون دادنش به صورت کسر مسلسل وجود داره:

$$ \begin{align} r = [1;1,1,1] &= 1+\frac{1}{1+\frac{1}{1+\frac{1}{1}}},\\ r = [1;1,2] &= 1+\frac{1}{1+\frac{1}{2}}. \end{align} $$

میشه ثابت کرد که هر عدد گویایی رو میشه دقیقاً به ۲ روش به صورت کسر مسلسل نشون داد:

$$r = [a_0;a_1,\dots,a_k,1] = [a_0;a_1,\dots,a_k+1].$$

یه نکته‌ی دیگه هم اینه که طول $k$ این کسر مسلسل برای $r=\frac{p}{q}$ از مرتبه‌ی $k = O(\log \min(p, q))$ هست.

دلیل این موضوع وقتی که ساختار کسر مسلسل رو با جزئیات بیشتری بررسی کنیم، روشن‌تر میشه.

تعریف

حالا فرض کن $a_0,a_1,a_2, \dots$ یه دنباله از اعداد صحیح باشه، طوری که $a_1, a_2, \dots \geq 1$. و فرض کن $r_k = [a_0; a_1, \dots, a_k]$. اون‌وقت به این عبارت:

$$r = a_0 + \frac{1}{a_1 + \frac{1}{a_2+\dots}} = \lim\limits_{k \to \infty} r_k.$$

میگن نمایش کسر مسلسل برای عدد گنگ $r$ و خلاصه‌اش رو اینجوری می‌نویسن: $r = [a_0;a_1,a_2,\dots]$.

یادت باشه که برای $r=[a_0;a_1,\dots]$ و یه عدد صحیح $k$، این رابطه برقراره: $r+k = [a_0+k; a_1, \dots]$.

یه مشاهده‌ی مهم دیگه اینه که اگه $a_0 > 0$ باشه، $\frac{1}{r}=[0;a_0, a_1, \dots]$ و اگه $a_0 = 0$ باشه، $\frac{1}{r} = [a_1; a_2, \dots]$.

تعریف

توی تعریف بالا، به عددهای گویای $r_0, r_1, r_2, \dots$ میگن همگراها (convergents)ی $r$.

به همین ترتیب، هر $r_k = [a_0; a_1, \dots, a_k] = \frac{p_k}{q_k}$، $k$-امین همگرای $r$ نامیده میشه.

مثال

بیا $r = [1; 1, 1, 1, \dots]$ رو در نظر بگیریم. میشه با استقرا ثابت کرد که $r_k = \frac{F_{k+2}}{F_{k+1}}$، که اینجا $F_k$ دنباله فیبوناچیه که اینجوری تعریف میشه: $F_0 = 0$, $F_1 = 1$ و $F_{k} = F_{k-1} + F_{k-2}$. از فرمول Binet هم می‌دونیم که:

$$r_k = \frac{\phi^{k+2} - \psi^{k+2}}{\phi^{k+1} - \psi^{k+1}},$$

که اینجا $\phi = \frac{1+\sqrt{5}}{2} \approx 1.618$ همون نسبت طلاییه و $\psi = \frac{1-\sqrt{5}}{2} = -\frac{1}{\phi} \approx -0.618$. بنابراین،

$$r = 1+\frac{1}{1+\frac{1}{1+\dots}}=\lim\limits_{k \to \infty} r_k = \phi = \frac{1+\sqrt{5}}{2}.$$

جالبه بدونی که توی این مورد خاص، یه راه دیگه برای پیدا کردن $r$ اینه که این معادله رو حل کنی:

$$r = 1+\frac{1}{r} \implies r^2 = r + 1. $$

تعریف

فرض کن $r_k = [a_0; a_1, \dots, a_{k-1}, a_k]$. به عددهای $[a_0; a_1, \dots, a_{k-1}, t]$ برای $1 \leq t \leq a_k$ میگن نیم‌همگرا (semiconvergents).

ما معمولاً به (نیم)همگراهایی که از $r$ بزرگترن میگیم (نیم)همگراهای بالایی و به اونایی که از $r$ کوچیکترن میگیم (نیم)همگراهای پایینی.

تعریف

برای تکمیل بحث همگراها، خارج‌قسمت‌های کامل رو هم اینجوری تعریف می‌کنیم: $s_k = [a_k; a_{k+1}, a_{k+2}, \dots]$.

به همین ترتیب، به $s_k$ میگیم $k$-امین خارج‌قسمت کامل $r$.

از تعریف‌های بالا، میشه نتیجه گرفت که $s_k \geq 1$ برای $k \geq 1$.

اگه $[a_0; a_1, \dots, a_k]$ رو به عنوان یه عبارت جبری در نظر بگیریم و به جای $a_i$ ها اعداد حقیقی دلخواه بذاریم، به این رابطه می‌رسیم:

$$r = [a_0; a_1, \dots, a_{k-1}, s_k].$$

به‌طور خاص، $r = [s_0] = s_0$. از طرف دیگه، می‌تونیم $s_k$ رو اینجوری بنویسیم:

$$s_k = [a_k; s_{k+1}] = a_k + \frac{1}{s_{k+1}},$$

این یعنی می‌تونیم $a_k = \lfloor s_k \rfloor$ و $s_{k+1} = (s_k - a_k)^{-1}$ رو از روی $s_k$ حساب کنیم.

دنباله‌ی $a_0, a_1, \dots$ همیشه به خوبی تعریف شده، مگر اینکه $s_k=a_k$ بشه که این فقط وقتی اتفاق میفته که $r$ یه عدد گویا باشه.

بنابراین، نمایش کسر مسلسل برای هر عدد گنگ $r$ به طور یکتا تعریف میشه.

پیاده‌سازی

توی کدهایی که می‌نویسیم، عمدتاً فرض می‌کنیم که با کسرهای مسلسل متناهی سر و کار داریم.

از روی $s_k$، برای رسیدن به $s_{k+1}$ اینجوری عمل می‌کنیم:

$$s_k =\left\lfloor s_k \right\rfloor + \frac{1}{s_{k+1}}.$$

از این عبارت، خارج‌قسمت کامل بعدی یعنی $s_{k+1}$ اینجوری به دست میاد:

$$s_{k+1} = \left(s_k-\left\lfloor s_k\right\rfloor\right)^{-1}.$$

برای $s_k=\frac{p}{q}$ این یعنی:

$$ s_{k+1} = \left(\frac{p}{q}-\left\lfloor \frac{p}{q} \right\rfloor\right)^{-1} = \frac{q}{p-q\cdot \lfloor \frac{p}{q} \rfloor} = \frac{q}{p \bmod q}. $$

بنابراین، محاسبه‌ی نمایش کسر مسلسل برای $r=\frac{p}{q}$ دقیقاً مراحل الگوریتم اقلیدس برای $p$ و $q$ رو دنبال می‌کنه.

از این موضوع این هم نتیجه میشه که $\gcd(p_k, q_k) = 1$ برای $\frac{p_k}{q_k} = [a_0; a_1, \dots, a_k]$. در نتیجه، همگراها همیشه کسرهای ساده‌نشدنی هستن.

auto fraction(int p, int q) {
    vector<int> a;
    while(q) {
        a.push_back(p / q);
        tie(p, q) = make_pair(q, p % q);
    }
    return a;
}
def fraction(p, q):
    a = []
    while q:
        a.append(p // q)
        p, q = q, p % q
    return a

نتایج کلیدی

برای اینکه بیشتر انگیزه بگیری که کسرهای مسلسل رو یاد بگیری، اینجا چند تا نکته‌ی کلیدی و باحال رو برات میگم.

رابطه بازگشتی

برای همگراهای $r_k = \frac{p_k}{q_k}$، این رابطه‌ی بازگشتی باحال برقراره که باعث میشه بتونی خیلی سریع حسابشون کنی:

$$\frac{p_k}{q_k}=\frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}},$$

که اینجا $\frac{p_{-1}}{q_{-1}}=\frac{1}{0}$ و $\frac{p_{-2}}{q_{-2}}=\frac{0}{1}$ رو به عنوان مقادیر اولیه در نظر می‌گیریم.

انحراف‌ها

انحراف $r_k = \frac{p_k}{q_k}$ از $r$ رو میشه به طور کلی اینجوری تخمین زد:

$$\left|\frac{p_k}{q_k}-r\right| \leq \frac{1}{q_k q_{k+1}} \leq \frac{1}{q_k^2}.$$

اگه دو طرف رو در $q_k$ ضرب کنیم، یه تخمین دیگه به دست میاد:

$$|p_k - q_k r| \leq \frac{1}{q_{k+1}}.$$

از رابطه‌ی بازگشتی بالا میشه نتیجه گرفت که $q_k$ حداقل با سرعتی برابر با اعداد فیبوناچی رشد می‌کنه.

توی عکس زیر می‌تونی ببینی که همگراهای $r_k$ چجوری به $r=\frac{1+\sqrt 5}{2}$ نزدیک میشن:

$r=\frac{1+\sqrt 5}{2}$ با خط‌چین آبی نشون داده شده. همگراهای فرد از بالا بهش نزدیک میشن و همگراهای زوج از پایین.
پوش‌های مشبکه

پوش‌های محدب (convex hulls) نقاط بالا و پایین خط $y=rx$ رو در نظر بگیر.

همگراهای فرد $(q_k;p_k)$ رأس‌های پوش بالایی هستن، در حالی که همگراهای زوج $(q_k;p_k)$ رأس‌های پوش پایینی هستن.

تمام رأس‌های صحیح روی این پوش‌ها به شکل $(q;p)$ هستن، طوری که:

$$\frac{p}{q} = \frac{tp_{k-1} + p_{k-2}}{tq_{k-1} + q_{k-2}}$$

برای یه عدد صحیح $0 \leq t \leq a_k$. به عبارت دیگه، مجموعه‌ی نقاط شبکه‌ای روی پوش‌ها با مجموعه‌ی نیم‌همگراها یکیه.

توی عکس زیر، می‌تونی همگراها و نیم‌همگراهای (نقاط خاکستری وسط) $r=\frac{9}{7}$ رو ببینی.

بهترین تقریب‌ها

فرض کن $\frac{p}{q}$ کسریه که مقدار $\left|r-\frac{p}{q}\right|$ رو با شرط $q \leq x$ برای یه $x$ مشخص، کمینه می‌کنه.

اون‌وقت $\frac{p}{q}$ حتماً یه نیم‌همگرای $r$ هست.

این نکته‌ی آخر به ما اجازه میده تا بهترین تقریب‌های گویای $r$ رو با بررسی نیم‌همگراهاش پیدا کنیم.

در ادامه توضیحات بیشتر و کمی شهود و تفسیر برای این نکات رو پیدا می‌کنی.

همگراها

بیا یه نگاه دقیق‌تر به همگراهایی که قبلاً تعریف کردیم بندازیم. برای $r=[a_0, a_1, a_2, \dots]$، همگراهاش اینا هستن:

$$\begin{gather} r_0=[a_0],\\r_1=[a_0, a_1],\\ \dots,\\ r_k=[a_0, a_1, \dots, a_k]. \end{gather}$$

همگراها مفهوم اصلی کسرهای مسلسل هستن، پس مطالعه‌ی ویژگی‌هاشون خیلی مهمه.

برای یه عدد $r$، $k$-امین همگرای اون یعنی $r_k = \frac{p_k}{q_k}$ رو میشه اینجوری حساب کرد:

$$r_k = \frac{P_k(a_0,a_1,\dots,a_k)}{P_{k-1}(a_1,\dots,a_k)} = \frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}},$$

که اینجا $P_k(a_0,\dots,a_k)$ یه ادامه‌دهنده (continuant) هست، که یه چندجمله‌ای چندمتغیره است و اینجوری تعریف میشه:

$$P_k(x_0,x_1,\dots,x_k) = \det \begin{bmatrix} x_k & 1 & 0 & \dots & 0 \\ -1 & x_{k-1} & 1 & \dots & 0 \\ 0 & -1 & x_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & 1 \\ 0 & 0 & \dots & -1 & x_0 \end{bmatrix}_{\textstyle .}$$

بنابراین، $r_k$ یه جور میانه (mediant) وزنی از $r_{k-1}$ و $r_{k-2}$ هست.

برای اینکه فرمول‌ها درست کار کنن، دو تا همگرای اضافی هم تعریف می‌کنیم: $r_{-1} = \frac{1}{0}$ و $r_{-2} = \frac{0}{1}$.

توضیح کامل

صورت و مخرج $r_k$ رو میشه به عنوان چندجمله‌ای‌های چندمتغیره از $a_0, a_1, \dots, a_k$ در نظر گرفت:

$$r_k = \frac{P_k(a_0, a_1, \dots, a_k)}{Q_k(a_0,a_1, \dots, a_k)}.$$

از تعریف همگراها داریم:

$$r_k = a_0 + \frac{1}{[a_1;a_2,\dots, a_k]}= a_0 + \frac{Q_{k-1}(a_1, \dots, a_k)}{P_{k-1}(a_1, \dots, a_k)} = \frac{a_0 P_{k-1}(a_1, \dots, a_k) + Q_{k-1}(a_1, \dots, a_k)}{P_{k-1}(a_1, \dots, a_k)}.$$

از این رابطه نتیجه میشه که $Q_k(a_0, \dots, a_k) = P_{k-1}(a_1, \dots, a_k)$. این ما رو به رابطه‌ی زیر می‌رسونه:

$$P_k(a_0, \dots, a_k) = a_0 P_{k-1}(a_1, \dots, a_k) + P_{k-2}(a_2, \dots, a_k).$$

در ابتدا، $r_0 = \frac{a_0}{1}$ و $r_1 = \frac{a_0 a_1 + 1}{a_1}$، بنابراین:

$$\begin{align}P_0(a_0)&=a_0,\\ P_1(a_0, a_1) &= a_0 a_1 + 1.\end{align}$$

برای اینکه فرمول‌ها سازگار باشن، راحت‌تره که $P_{-1} = 1$ و $P_{-2}=0$ رو تعریف کنیم و به طور صوری بگیم $r_{-1} = \frac{1}{0}$ و $r_{-2}=\frac{0}{1}$.

از آنالیز عددی، می‌دونیم که دترمینان یه ماتریس سه‌قطری دلخواه

$$T_k = \det \begin{bmatrix} a_0 & b_0 & 0 & \dots & 0 \\ c_0 & a_1 & b_1 & \dots & 0 \\ 0 & c_1 & a_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & c_{k-1} \\ 0 & 0 & \dots & b_{k-1} & a_k \end{bmatrix}$$

میتونه به صورت بازگشتی به شکل $T_k = a_k T_{k-1} - b_{k-1} c_{k-1} T_{k-2}$ محاسبه بشه. با مقایسه‌ی این با $P_k$، یه عبارت مستقیم به دست میاریم:

$$P_k = \det \begin{bmatrix} x_k & 1 & 0 & \dots & 0 \\ -1 & x_{k-1} & 1 & \dots & 0 \\ 0 & -1 & x_2 & . & \vdots \\ \vdots & \vdots & . & \ddots & 1 \\ 0 & 0 & \dots & -1 & x_0 \end{bmatrix}_{\textstyle .}$$

این چندجمله‌ای به خاطر ارتباط نزدیکش با کسرهای مسلسل، به عنوان ادامه‌دهنده (continuant) هم شناخته میشه. ادامه‌دهنده اگه دنباله روی قطر اصلی برعکس بشه، تغییر نمی‌کنه. این یه فرمول جایگزین برای محاسبه‌اش به ما میده:

$$P_k(a_0, \dots, a_k) = a_k P_{k-1}(a_0, \dots, a_{k-1}) + P_{k-2}(a_0, \dots, a_{k-2}).$$

پیاده‌سازی

ما همگراها رو به صورت یه جفت دنباله‌ی $p_{-2}, p_{-1}, p_0, p_1, \dots, p_k$ و $q_{-2}, q_{-1}, q_0, q_1, \dots, q_k$ حساب می‌کنیم:

auto convergents(vector<int> a) {
    vector<long long> p = {0, 1};
    vector<long long> q = {1, 0};
    for(auto it: a) {
        p.push_back(p.back() * it + p[p.size() - 2]);
        q.push_back(q.back() * it + q[q.size() - 2]);
    }
    return make_pair(p, q);
}
def convergents(a):
    p = [0, 1]
    q = [1, 0]
    for it in a:
        p.append(p[-1]*it + p[-2])
        q.append(q[-1]*it + q[-2])
    return p, q

درخت‌های کسرهای مسلسل

دو تا راه اصلی وجود داره که میشه همه‌ی کسرهای مسلسل ممکن رو توی ساختارهای درختی مفید جا داد.

درخت اشترن-بروکوت (Stern-Brocot)

درخت اشترن-بروکوت یه درخت جستجوی دودوییه که شامل همه‌ی اعداد گویای مثبت و متمایزه.

شکل کلی درخت اینجوریه:

تصویر از Aaron Rotenberg تحت مجوز CC BY-SA 3.0

کسرهای $\frac{0}{1}$ و $\frac{1}{0}$ به طور "مجازی" به ترتیب در سمت چپ و راست درخت قرار دارن.

بعد، کسر توی هر گره، میانه (mediant) $\frac{a+c}{b+d}$ دو کسر $\frac{a}{b}$ و $\frac{c}{d}$ بالای اون گره هست.

رابطه‌ی بازگشتی $\frac{p_k}{q_k}=\frac{a_k p_{k-1} + p_{k-2}}{a_k q_{k-1} + q_{k-2}}$ یعنی نمایش کسر مسلسل، مسیر رسیدن به $\frac{p_k}{q_k}$ رو توی درخت کدگذاری می‌کنه. برای پیدا کردن $[a_0; a_1, \dots, a_{k}, 1]$، باید $a_0$ بار به راست بری، $a_1$ بار به چپ، $a_2$ بار به راست و همینطور تا $a_k$ ادامه بدی.

والدِ $[a_0; a_1, \dots, a_k,1]$ کسریه که با برداشتن یه قدم به عقب در آخرین جهتی که رفتی به دست میاد.

به عبارت دیگه، وقتی $a_k > 1$ باشه، والدش میشه $[a_0; a_1, \dots, a_k-1,1]$ و وقتی $a_k = 1$ باشه، والدش میشه $[a_0; a_1, \dots, a_{k-1}, 1]$.

بنابراین، فرزندان $[a_0; a_1, \dots, a_k, 1]$ اینا هستن: $[a_0; a_1, \dots, a_k+1, 1]$ و $[a_0; a_1, \dots, a_k, 1, 1]$.

بیا درخت اشترن-بروکوت رو شماره‌گذاری کنیم. به رأس ریشه شماره ۱ رو میدیم. بعد، برای یه رأس $v$، شماره‌ی فرزند چپش با تغییر بیت پیشروی $v$ از ۱ به ۱۰ و برای فرزند راست، با تغییر بیت پیشرو از ۱ به ۱۱ به دست میاد:

توی این شماره‌گذاری، نمایش کسر مسلسل یه عدد گویا، کدگذاری طول اجرا (run-length encoding) شماره‌ی دودویی اون رو مشخص می‌کنه.

برای $\frac{5}{2} = [2;2] = [2;1,1]$، شماره‌اش $1011_2$ هست و کدگذاری طول اجرای اون، با در نظر گرفتن بیت‌ها به ترتیب صعودی، میشه $[2;1,1]$.

یه مثال دیگه $\frac{2}{5} = [0;2,2]=[0;2,1,1]$ هست که شماره‌اش $1100_2$ هست و کدگذاری طول اجرای اون در واقع $[0;2,2]$ میشه.

جالبه بدونی که درخت اشترن-بروکوت، در واقع یه treap هست. یعنی، از نظر مقدار $\frac{p}{q}$ یه درخت جستجوی دودوییه، اما از نظر $p$ و $q$ هر دو یه هیپ (heap) هست.

مقایسه کسرهای مسلسل

فرض کن دو تا کسر $A=[a_0; a_1, \dots, a_n]$ و $B=[b_0; b_1, \dots, b_m]$ داری. چطوری بفهمیم کدومشون کوچیکتره؟

راه‌حل

فعلاً فرض کن $A$ و $B$ گنگ هستن و نمایش کسرهای مسلسلشون یه مسیر بی‌نهایت رو توی درخت اشترن-بروکوت نشون میده.

همونطور که قبلاً گفتیم، توی این نمایش $a_0$ تعداد حرکت‌های به راست، $a_1$ تعداد حرکت‌های متوالی به چپ و غیره رو نشون میده. بنابراین، وقتی $a_k$ و $b_k$ رو مقایسه می‌کنیم، اگه $a_k = b_k$ باشه باید بریم سراغ مقایسه‌ی $a_{k+1}$ و $b_{k+1}$. در غیر این صورت، اگه توی مسیرهای راست هستیم، باید ببینیم آیا $a_k < b_k$ و اگه توی مسیرهای چپ هستیم، باید ببینیم آیا $a_k > b_k$ تا بگیم آیا $A < B$.

به عبارت دیگه، برای $A$ و $B$ گنگ، $A < B$ خواهد بود اگر و تنها اگر $(a_0, -a_1, a_2, -a_3, \dots) < (b_0, -b_1, b_2, -b_3, \dots)$ با مقایسه‌ی لغت‌نامه‌ای (lexicographical comparison).

حالا، با استفاده‌ی صوری از $\infty$ به عنوان یه عنصر از نمایش کسر مسلسل، می‌تونیم اعداد گنگ $A-\varepsilon$ و $A+\varepsilon$ رو شبیه‌سازی کنیم، یعنی عناصری که کوچیکتر (بزرگتر) از $A$ هستن، اما بزرگتر (کوچیکتر) از هر عدد حقیقی دیگه‌ای هستن. به طور خاص، برای $A=[a_0; a_1, \dots, a_n]$، یکی از این دو عنصر رو میشه به صورت $[a_0; a_1, \dots, a_n, \infty]$ و اون یکی رو به صورت $[a_0; a_1, \dots, a_n - 1, 1, \infty]$ شبیه‌سازی کرد.

اینکه کدوم یکی به $A-\varepsilon$ و کدوم به $A+\varepsilon$ مربوط میشه، با زوجیت $n$ یا با مقایسه‌ی اون‌ها به عنوان اعداد گنگ مشخص میشه.

# بررسی می‌کنه که آیا a < b با فرض اینکه a[-1] = b[-1] = infty و a != b
def less(a, b):
    a = [(-1)**i*a[i] for i in range(len(a))]
    b = [(-1)**i*b[i] for i in range(len(b))]
    return a < b

# [a0; a1, ..., ak] -> [a0, a1, ..., ak-1, 1]
def expand(a):
    if a: # a خالی = inf
        a[-1] -= 1
        a.append(1)
    return a

# a-eps, a+eps رو برمی‌گردونه
def pm_eps(a):
    b = expand(a.copy())
    a.append(float('inf'))
    b.append(float('inf'))
    return (a, b) if less(a, b) else (b, a)

بهترین نقطه درونی

فرض کن $\frac{0}{1} \leq \frac{p_0}{q_0} < \frac{p_1}{q_1} \leq \frac{1}{0}$ رو داری. عدد گویای $\frac{p}{q}$ رو پیدا کن طوری که $(q; p)$ از نظر لغت‌نامه‌ای کوچیکترین باشه و $\frac{p_0}{q_0} < \frac{p}{q} < \frac{p_1}{q_1}$.

راه‌حل

از دیدگاه درخت اشترن-بروکوت، این یعنی ما باید پایین‌ترین جد مشترک (LCA) $\frac{p_0}{q_0}$ و $\frac{p_1}{q_1}$ رو پیدا کنیم. به خاطر ارتباط بین درخت اشترن-بروکوت و کسر مسلسل، این LCA تقریباً با بزرگترین پیشوند مشترک نمایش‌های کسر مسلسل برای $\frac{p_0}{q_0}$ و $\frac{p_1}{q_1}$ یکیه.

بنابراین، اگه $\frac{p_0}{q_0} = [a_0; a_1, \dots, a_{k-1}, a_k, \dots]$ و $\frac{p_1}{q_1} = [a_0; a_1, \dots, a_{k-1}, b_k, \dots]$ اعداد گنگ باشن، LCA برابر میشه با $[a_0; a_1, \dots, \min(a_k, b_k)+1]$.

برای اعداد گویای $r_0$ و $r_1$، یکی از اون‌ها می‌تونه خودش LCA باشه که ما رو مجبور می‌کنه حالت‌های مختلفی رو بررسی کنیم. برای ساده کردن راه‌حل برای $r_0$ و $r_1$ گویا، میشه از نمایش کسر مسلسل $r_0 + \varepsilon$ و $r_1 - \varepsilon$ که توی مسئله‌ی قبلی به دست آوردیم، استفاده کرد.

# کوچیکترین (q, p) از نظر لغت‌نامه‌ای رو پیدا می‌کنه
# طوری که p0/q0 < p/q < p1/q1
def middle(p0, q0, p1, q1):
    a0 = pm_eps(fraction(p0, q0))[1]
    a1 = pm_eps(fraction(p1, q1))[0]
    a = []
    for i in range(min(len(a0), len(a1))):
        a.append(min(a0[i], a1[i]))
        if a0[i] != a1[i]:
            break
    a[-1] += 1
    p, q = convergents(a)
    return p[-1], q[-1]

GCJ 2019, Round 2 - New Elements: Part 2

$N$ تا جفت عدد صحیح مثبت $(C_i, J_i)$ داری. باید یه جفت عدد صحیح مثبت $(x, y)$ پیدا کنی طوری که $C_i x + J_i y$ یه دنباله‌ی اکیداً صعودی باشه.

از بین همه‌ی این جفت‌ها، کوچیکترین جفت از نظر لغت‌نامه‌ای رو پیدا کن.

راه‌حل

اگه صورت مسئله رو بازنویسی کنیم، $A_i x + B_i y$ باید برای همه‌ی $i$ ها مثبت باشه، که اینجا $A_i = C_i - C_{i-1}$ و $B_i = J_i - J_{i-1}$.

بین این معادلات، چهار گروه مهم برای $A_i x + B_i y > 0$ داریم:

  1. $A_i, B_i > 0$ رو میشه نادیده گرفت چون ما دنبال $x, y > 0$ هستیم.
  2. $A_i, B_i \leq 0$ باعث میشه جواب "IMPOSSIBLE" باشه.
  3. $A_i > 0$, $B_i \leq 0$. این قیدها معادل $\frac{y}{x} < \frac{A_i}{-B_i}$ هستن.
  4. $A_i \leq 0$, $B_i > 0$. این قیدها معادل $\frac{y}{x} > \frac{-A_i}{B_i}$ هستن.

فرض کن $\frac{p_0}{q_0}$ بزرگترین $\frac{-A_i}{B_i}$ از گروه چهارم و $\frac{p_1}{q_1}$ کوچیکترین $\frac{A_i}{-B_i}$ از گروه سوم باشه.

مسئله الان اینه که با داشتن $\frac{p_0}{q_0} < \frac{p_1}{q_1}$، یه کسر $\frac{p}{q}$ پیدا کنیم طوری که $(q;p)$ از نظر لغت‌نامه‌ای کوچیکترین باشه و $\frac{p_0}{q_0} < \frac{p}{q} < \frac{p_1}{q_1}$.

    def solve():
    n = int(input())
    C = [0] * n
    J = [0] * n
    # p0/q0 < y/x < p1/q1
    p0, q0 = 0, 1
    p1, q1 = 1, 0
    fail = False
    for i in range(n):
        C[i], J[i] = map(int, input().split())
        if i > 0:
            A = C[i] - C[i-1]
            B = J[i] - J[i-1]
            if A <= 0 and B <= 0:
                fail = True
            elif B > 0 and A < 0: # y/x > (-A)/B if B > 0
                if (-A)*q0 > p0*B:
                    p0, q0 = -A, B
            elif B < 0 and A > 0: # y/x < A/(-B) if B < 0
                if A*q1 < p1*(-B):
                    p1, q1 = A, -B
    if p0*q1 >= p1*q0 or fail:
        return 'IMPOSSIBLE'

    p, q = middle(p0, q0, p1, q1)
    return str(q) + ' ' + str(p)

درخت کالکین-ویلف (Calkin-Wilf)

یه راه یه کم ساده‌تر برای سازماندهی کسرهای مسلسل توی یه درخت دودویی، درخت کالکین-ویلف هست.

شکل کلی درخت اینجوریه:

تصویر از Olli Niemitalo, Proz تحت مجوز CC0 1.0

توی ریشه‌ی درخت، عدد $\frac{1}{1}$ قرار داره. بعد، برای یه رأس با عدد $\frac{p}{q}$، فرزندانش $\frac{p}{p+q}$ و $\frac{p+q}{q}$ هستن.

برخلاف درخت اشترن-بروکوت، درخت کالکین-ویلف یه درخت جستجوی دودویی نیست، پس نمیشه ازش برای جستجوی دودویی روی اعداد گویا استفاده کرد.

توی درخت کالکین-ویلف، والد مستقیم کسر $\frac{p}{q}$، وقتی $p>q$ باشه میشه $\frac{p-q}{q}$ و در غیر این صورت میشه $\frac{p}{q-p}$.

برای درخت اشترن-بروکوت، از رابطه‌ی بازگشتی برای همگراها استفاده کردیم. برای اینکه بین کسر مسلسل و درخت کالکین-ویلف ارتباط برقرار کنیم، باید رابطه‌ی بازگشتی برای خارج‌قسمت‌های کامل رو به یاد بیاریم. اگه $s_k = \frac{p}{q}$ باشه، اون‌وقت $s_{k+1} = \frac{q}{p \mod q} = \frac{q}{p-\lfloor p/q \rfloor \cdot q}$.

از طرف دیگه، اگه به طور مکرر از $s_k = \frac{p}{q}$ به والدش توی درخت کالکین-ویلف بریم (وقتی $p > q$)، به $\frac{p \mod q}{q} = \frac{1}{s_{k+1}}$ می‌رسیم. اگه این کار رو ادامه بدیم، به $s_{k+2}$، بعدش $\frac{1}{s_{k+3}}$ و همینطور الی آخر می‌رسیم. از این میشه نتیجه گرفت که:

  1. وقتی $a_0> 0$، والد مستقیم $[a_0; a_1, \dots, a_k]$ توی درخت کالکین-ویلف میشه $\frac{p-q}{q}=[a_0 - 1; a_1, \dots, a_k]$.
  2. وقتی $a_0 = 0$ و $a_1 > 1$، والد مستقیمش میشه $\frac{p}{q-p} = [0; a_1 - 1, a_2, \dots, a_k]$.
  3. و وقتی $a_0 = 0$ و $a_1 = 1$، والد مستقیمش میشه $\frac{p}{q-p} = [a_2; a_3, \dots, a_k]$.

به همین ترتیب، فرزندان $\frac{p}{q} = [a_0; a_1, \dots, a_k]$ اینا هستن:

  1. $\frac{p+q}{q}=1+\frac{p}{q}$، که برابر با $[a_0+1; a_1, \dots, a_k]$ هست،
  2. $\frac{p}{p+q} = \frac{1}{1+\frac{q}{p}}$، که برای $a_0 > 0$ برابر با $[0, 1, a_0, a_1, \dots, a_k]$ و برای $a_0=0$ برابر با $[0, a_1+1, a_2, \dots, a_k]$ هست.

جالبه بدونی، اگه رأس‌های درخت کالکین-ویلف رو به ترتیب جستجوی سطح-اول (level-order) شماره‌گذاری کنیم (یعنی ریشه شماره ۱ داره و فرزندان رأس $v$ به ترتیب شماره‌های $2v$ و $2v+1$ رو دارن)، شماره‌ی یه عدد گویا توی درخت کالکین-ویلف با شماره‌اش توی درخت اشترن-بروکوت یکی میشه.

بنابراین، اعداد توی سطوح یکسان درخت اشترن-بروکوت و درخت کالکین-ویلف یکی هستن، اما ترتیبشون از طریق جایگشت بیت-معکوس (bit-reversal permutation) با هم فرق داره.

همگرایی

برای یه عدد $r$ و $k$-امین همگرای اون $r_k=\frac{p_k}{q_k}$ این فرمول برقراره:

$$r_k = a_0 + \sum\limits_{i=1}^k \frac{(-1)^{i-1}}{q_i q_{i-1}}.$$

به طور خاص، این یعنی:

$$r_k - r_{k-1} = \frac{(-1)^{k-1}}{q_k q_{k-1}}$$

و

$$p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1}.$$

از این می‌تونیم نتیجه بگیریم که:

$$\left| r-\frac{p_k}{q_k} \right| \leq \frac{1}{q_{k+1}q_k} \leq \frac{1}{q_k^2}.$$

نامساوی دوم به این دلیله که $r_k$ و $r_{k+1}$ معمولاً در دو طرف مختلف $r$ قرار دارن، بنابراین:

$$|r-r_k| = |r_k-r_{k+1}|-|r-r_{k+1}| \leq |r_k - r_{k+1}|.$$
توضیح کامل

برای تخمین زدن $|r-r_k|$، با تخمین تفاوت بین همگراهای مجاور شروع می‌کنیم. طبق تعریف:

$$\frac{p_k}{q_k} - \frac{p_{k-1}}{q_{k-1}} = \frac{p_k q_{k-1} - p_{k-1} q_k}{q_k q_{k-1}}.$$

با جایگزین کردن $p_k$ و $q_k$ توی صورت کسر با روابط بازگشتی‌شون، داریم:

$$\begin{align} p_k q_{k-1} - p_{k-1} q_k &= (a_k p_{k-1} + p_{k-2}) q_{k-1} - p_{k-1} (a_k q_{k-1} + q_{k-2}) \\&= p_{k-2} q_{k-1} - p_{k-1} q_{k-2},\end{align}$$

بنابراین صورت کسر $r_k - r_{k-1}$ همیشه برابر با قرینه‌ی صورت کسر $r_{k-1} - r_{k-2}$ هست. این به نوبه‌ی خودش، برای

$$r_1 - r_0=\left(a_0+\frac{1}{a_1}\right)-a_0=\frac{1}{a_1},$$

برابر با ۱ هست، پس:

$$r_k - r_{k-1} = \frac{(-1)^{k-1}}{q_k q_{k-1}}.$$

این یه نمایش جایگزین از $r_k$ به عنوان مجموع جزئی یه سری بی‌نهایت به ما میده:

$$r_k = (r_k - r_{k-1}) + \dots + (r_1 - r_0) + r_0 = a_0 + \sum\limits_{i=1}^k \frac{(-1)^{i-1}}{q_i q_{i-1}}.$$

از رابطه‌ی بازگشتی نتیجه میشه که $q_k$ به صورت یکنواخت حداقل با سرعتی برابر با اعداد فیبوناچی افزایش پیدا می‌کنه، بنابراین:

$$r = \lim\limits_{k \to \infty} r_k = a_0 + \sum\limits_{i=1}^\infty \frac{(-1)^{i-1}}{q_i q_{i-1}}$$

همیشه به خوبی تعریف شده، چون سری زیربنایی همیشه همگراست. جالبه بدونی، سری باقیمانده

$$r-r_k = \sum\limits_{i=k+1}^\infty \frac{(-1)^{i-1}}{q_i q_{i-1}}$$

به خاطر سرعت کاهش $q_i q_{i-1}$، هم‌علامت با $(-1)^k$ هست. از این رو $r_k$ با اندیس زوج از پایین به $r$ نزدیک میشه در حالی که $r_k$ با اندیس فرد از بالا بهش نزدیک میشه:

همگراهای $r=\phi = \frac{1+\sqrt{5}}{2}=[1;1,1,\dots]$ و فاصله‌شون از $r$.

از این عکس می‌تونیم ببینیم که:

$$|r-r_k| = |r_k - r_{k+1}| - |r-r_{k+1}| \leq |r_k - r_{k+1}|,$$

بنابراین فاصله‌ی بین $r$ و $r_k$ هیچوقت بزرگتر از فاصله‌ی بین $r_k$ و $r_{k+1}$ نیست:

$$\left|r-\frac{p_k}{q_k}\right| \leq \frac{1}{q_k q_{k+1}} \leq \frac{1}{q_k^2}.$$

اقلیدس گسترش‌یافته؟

سه تا عدد صحیح $A, B, C \in \mathbb Z$ داری. اعداد صحیح $x, y \in \mathbb Z$ رو پیدا کن طوری که $Ax + By = C$.

راه‌حل

اگرچه این مسئله معمولاً با الگوریتم اقلیدس گسترش‌یافته حل میشه، اما یه راه‌حل ساده و مستقیم با کسرهای مسلسل وجود داره.

فرض کن $\frac{A}{B}=[a_0; a_1, \dots, a_k]$. بالاتر ثابت کردیم که $p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1}$. با جایگزین کردن $p_k$ و $q_k$ با $A$ و $B$، داریم:

$$Aq_{k-1} - Bp_{k-1} = (-1)^{k-1} g,$$

که اینجا $g = \gcd(A, B)$ هست. اگه $C$ بر $g$ بخش‌پذیر باشه، اون‌وقت راه‌حل میشه $x = (-1)^{k-1}\frac{C}{g} q_{k-1}$ و $y = (-1)^{k}\frac{C}{g} p_{k-1}$.

# (x, y) رو برمی‌گردونه طوری که Ax+By=C
# فرض میشه که چنین (x, y) وجود داره
def dio(A, B, C):
    p, q = convergents(fraction(A, B))
    C //= A // p[-1] # تقسیم بر ب.م.م(A, B)
    t = -1 if len(p) % 2 == 0 else 1 # در متن اصلی (-1)^(k-1) بود. len(p) = k+3. پس k-1 زوج است اگر len(p) فرد باشد
    return t*C*q[-2], -t*C*p[-2]

تبدیلات کسری خطی

یه مفهوم مهم دیگه برای کسرهای مسلسل، به اصطلاح تبدیلات کسری خطی هست.

تعریف

یه تبدیل کسری خطی یه تابعی مثل $f : \mathbb R \to \mathbb R$ هست طوری که $f(x) = \frac{ax+b}{cx+d}$ برای مقادیری از $a,b,c,d \in \mathbb R$.

ترکیب $(L_0 \circ L_1)(x) = L_0(L_1(x))$ از تبدیلات کسری خطی $L_0(x)=\frac{a_0 x + b_0}{c_0 x + d_0}$ و $L_1(x)=\frac{a_1 x + b_1}{c_1 x + d_1}$ خودش یه تبدیل کسری خطی هست:

$$\frac{a_0\frac{a_1 x + b_1}{c_1 x + d_1} + b_0}{c_0 \frac{a_1 x + b_1}{c_1 x + d_1} + d_0} = \frac{a_0(a_1 x + b_1) + b_0 (c_1 x + d_1)}{c_0 (a_1 x + b_1) + d_0 (c_1 x + d_1)} = \frac{(a_0 a_1 + b_0 c_1) x + (a_0 b_1 + b_0 d_1)}{(c_0 a_1 + d_0 c_1) x + (c_0 b_1 + d_0 d_1)}.$$

معکوس یه تبدیل کسری خطی، اونم یه تبدیل کسری خطی هست:

$$y = \frac{ax+b}{cx+d} \iff y(cx+d) = ax + b \iff x = -\frac{dy-b}{cy-a}.$$

DMOPC '19 Contest 7 P4 - Bob and Continued Fractions

یه آرایه از اعداد صحیح مثبت $a_1, \dots, a_n$ داری. باید به $m$ تا پرس‌وجو جواب بدی. هر پرس‌وجو محاسبه‌ی $[a_l; a_{l+1}, \dots, a_r]$ هست.

راه‌حل

اگه بتونیم کسرهای مسلسل رو به هم بچسبونیم، می‌تونیم این مسئله رو با درخت بازه‌ها (segment tree) حل کنیم.

به طور کلی این درسته که $[a_0; a_1, \dots, a_k, b_0, b_1, \dots, b_m] = [a_0; a_1, \dots, a_k, [b_0; b_1, \dots, b_m]]$. اینجا $b_0$ باید جایگزین $x$ بشه.

بیا $L_{k}(x) = [a_k; x] = a_k + \frac{1}{x} = \frac{a_k\cdot x+1}{1\cdot x + 0}$ رو نشون بدیم. توجه کن که $L_k(\infty) = a_k$. با این نمادگذاری، داریم:

$$[a_0; a_1, \dots, a_k, x] = [a_0; [a_1; [\dots; [a_k; x]]]] = (L_0 \circ L_1 \circ \dots \circ L_k)(x) = \frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}.$$

بنابراین، مسئله به محاسبه‌ی این خلاصه میشه:

$$(L_l \circ L_{l+1} \circ \dots \circ L_r)(\infty).$$

ترکیب تبدیلات شرکت‌پذیره، پس می‌تونیم توی هر گره از یه درخت بازه‌ها، ترکیب تبدیلات زیردرختش رو حساب کنیم.

تبدیل کسری خطی یک کسر مسلسل

فرض کن $L(x) = \frac{ax+b}{cx+d}$. نمایش کسر مسلسل $[b_0; b_1, \dots, b_m]$ از $L(A)$ رو برای $A=[a_0; a_1, \dots, a_n]$ حساب کن.

این کار به ما امکان محاسبه‌ی $A + \frac{p}{q} = \frac{qA + p}{q}$ و $A \cdot \frac{p}{q} = \frac{p A}{q}$ رو برای هر $\frac{p}{q}$ میده.

راه‌حل

همونطور که بالاتر اشاره کردیم، $[a_0; a_1, \dots, a_k] = (L_{a_0} \circ L_{a_1} \circ \dots \circ L_{a_k})(\infty)$، از این رو $L([a_0; a_1, \dots, a_k]) = (L \circ L_{a_0} \circ L_{a_1} \circ \dots L_{a_k})(\infty)$.

بنابراین، با اضافه کردن متوالی $L_{a_0}$، $L_{a_1}$ و غیره، می‌تونیم این رو حساب کنیم:

$$(L \circ L_{a_0} \circ \dots \circ L_{a_k})(x) = L\left(\frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}\right)=\frac{a_k' x + b_k'}{c_k' x + d_k'}.$$

از اونجایی که $L(x)$ معکوس‌پذیره، توی $x$ هم یکنواست. بنابراین، برای هر $x \geq 0$ داریم که $L(\frac{p_k x + p_{k-1}}{q_k x + q_{k-1}})$ بین $L(\frac{p_k}{q_k}) = \frac{a_k'}{c_k'}$ و $L(\frac{p_{k-1}}{q_{k-1}}) = \frac{b_k'}{d_k'}$ قرار داره.

علاوه بر این، برای $x=[a_{k+1}; \dots, a_n]$ این مقدار برابر با $L(A)$ هست. از این رو، $b_0 = \lfloor L(A) \rfloor$ بین $\lfloor L(\frac{p_k}{q_k}) \rfloor$ و $\lfloor L(\frac{p_{k-1}}{q_{k-1}}) \rfloor$ هست. وقتی این دو تا برابر باشن، اون‌ها هم برابر با $b_0$ هستن.

توجه کن که $L(A) = (L_{b_0} \circ L_{b_1} \circ \dots \circ L_{b_m})(\infty)$. با دونستن $b_0$، می‌تونیم $L_{b_0}^{-1}$ رو با تبدیل فعلی ترکیب کنیم و به اضافه کردن $L_{a_{k+1}}$، $L_{a_{k+2}}$ و غیره ادامه بدیم، و دنبال توافق کف‌های جدید بگردیم، که از اون می‌تونیم $b_1$ و غیره رو استنتاج کنیم تا زمانی که همه‌ی مقادیر $[b_0; b_1, \dots, b_m]$ رو به دست بیاریم.

محاسبات روی کسرهای مسلسل

فرض کن $A=[a_0; a_1, \dots, a_n]$ و $B=[b_0; b_1, \dots, b_m]$. نمایش‌های کسر مسلسل $A+B$ و $A \cdot B$ رو حساب کن.

راه‌حل

ایده‌ی اینجا شبیه به مسئله‌ی قبلیه، اما به جای $L(x) = \frac{ax+b}{cx+d}$ باید تبدیل کسری دوخطی $L(x, y) = \frac{axy+bx+cy+d}{exy+fx+gy+h}$ رو در نظر بگیری.

به جای $L(x) \mapsto L(L_{a_k}(x))$، تبدیل فعلی خودت رو به صورت $L(x, y) \mapsto L(L_{a_k}(x), y)$ یا $L(x, y) \mapsto L(x, L_{b_k}(y))$ تغییر میدی.

بعد، بررسی می‌کنی که آیا $\lfloor \frac{a}{e} \rfloor = \lfloor \frac{b}{f} \rfloor = \lfloor \frac{c}{g} \rfloor = \lfloor \frac{d}{h} \rfloor$ و اگه همه‌شون یکی بودن، از این مقدار به عنوان $c_k$ توی کسر حاصل استفاده می‌کنی و تبدیل رو اینجوری تغییر میدی:

$$L(x, y) \mapsto \frac{1}{L(x, y) - c_k}.$$

تعریف

یه کسر مسلسل $x = [a_0; a_1, \dots]$ رو متناوب (periodic) میگن اگه $x = [a_0; a_1, \dots, a_k, x]$ برای یه $k$ مشخص.

یه کسر مسلسل $x = [a_0; a_1, \dots]$ رو در نهایت متناوب (eventually periodic) میگن اگه $x = [a_0; a_1, \dots, a_k, y]$، که اینجا $y$ متناوبه.

برای $x = [1; 1, 1, \dots]$ داریم $x = 1 + \frac{1}{x}$، پس $x^2 = x + 1$. یه ارتباط کلی بین کسرهای مسلسل متناوب و معادلات درجه دو وجود داره. این معادله رو در نظر بگیر:

$$ x = [a_0; a_1, \dots, a_k, x].$$

از یه طرف، این معادله یعنی نمایش کسر مسلسل $x$ با دوره‌ی تناوب $k+1$ متناوبه.

از طرف دیگه، با استفاده از فرمول همگراها، این معادله یعنی:

$$x = \frac{p_k x + p_{k-1}}{q_k x + q_{k-1}}.$$

یعنی، $x$ یه تبدیل کسری خطی از خودشه. از این معادله نتیجه میشه که $x$ یه ریشه‌ی معادله‌ی درجه دومه:

$$q_k x^2 + (q_{k-1}-p_k)x - p_{k-1} = 0.$$

یه استدلال مشابه برای کسرهای مسلسلی که در نهایت متناوب هستن هم برقراره، یعنی $x = [a_0; a_1, \dots, a_k, y]$ برای $y=[b_0; b_1, \dots, b_m, y]$. در واقع، از معادله‌ی اول نتیجه می‌گیریم که $x = L_0(y)$ و از معادله‌ی دوم که $y = L_1(y)$، که اینجا $L_0$ و $L_1$ تبدیلات کسری خطی هستن. بنابراین،

$$x = (L_0 \circ L_1)(y) = (L_0 \circ L_1 \circ L_0^{-1})(x).$$

میشه حتی ثابت کرد (و اولین بار لاگرانژ این کار رو کرد) که برای هر معادله‌ی درجه دوم دلخواه $ax^2+bx+c=0$ با ضرایب صحیح، ریشه‌اش $x$ یه کسر مسلسل در نهایت متناوبه.

گنگی درجه دو

کسر مسلسل $\alpha = \frac{x+y\sqrt{n}}{z}$ رو پیدا کن که اینجا $x, y, z, n \in \mathbb Z$ و $n > 0$ یه مربع کامل نیست.

راه‌حل

برای $k$-امین خارج‌قسمت کامل $s_k$ از یه عدد، به طور کلی داریم که:

$$\alpha = [a_0; a_1, \dots, a_{k-1}, s_k] = \frac{s_k p_{k-1} + p_{k-2}}{s_k q_{k-1} + q_{k-2}}.$$

بنابراین،

$$s_k = -\frac{\alpha q_{k-1} - p_{k-1}}{\alpha q_k - p_k} = -\frac{q_{k-1} y \sqrt n + (x q_{k-1} - z p_{k-1})}{q_k y \sqrt n + (xq_k-zp_k)}.$$

با ضرب صورت و مخرج در $(xq_k - zp_k) - q_k y \sqrt n$، از شر $\sqrt n$ توی مخرج خلاص میشیم، پس خارج‌قسمت‌های کامل به این شکل در میان:

$$s_k = \frac{x_k + y_k \sqrt n}{z_k}.$$

بیا $s_{k+1}$ رو با فرض اینکه $s_k$ رو می‌دونیم، پیدا کنیم.

اول از همه، $a_k = \lfloor s_k \rfloor = \left\lfloor \frac{x_k + y_k \lfloor \sqrt n \rfloor}{z_k} \right\rfloor$ اشتباهه. باید $\lfloor \frac{x_k + y_k \sqrt n}{z_k} \rfloor$ باشه. اما چون $y_k \sqrt n$ گنگه، محاسبه‌اش سخته. توی کد از $\lfloor \frac{x_k + y_k \lfloor \sqrt n \rfloor}{z_k} \rfloor$ استفاده شده که برای $\alpha = \sqrt n$ کار می‌کنه.

$$s_{k+1} = \frac{1}{s_k-a_k} = \frac{z_k}{(x_k - z_k a_k) + y_k \sqrt n} = \frac{z_k (x_k - z_k a_k) - y_k z_k \sqrt n}{(x_k - z_k a_k)^2 - y_k^2 n}.$$

بنابراین، اگه $t_k = x_k - z_k a_k$ رو نشون بدیم، خواهیم داشت:

$$\begin{align}x_{k+1} &=& z_k t_k, \\ y_{k+1} &=& -y_k z_k, \\ z_{k+1} &=& t_k^2 - y_k^2 n.\end{align}$$

نکته‌ی خوب در مورد این نمایش اینه که اگه $x_{k+1}, y_{k+1}, z_{k+1}$ رو بر بزرگترین مقسوم‌علیه مشترکشون ساده کنیم، نتیجه منحصر به فرد میشه. بنابراین، می‌تونیم ازش برای بررسی اینکه آیا حالت فعلی قبلاً تکرار شده و همچنین برای بررسی اینکه اندیس قبلی که این حالت رو داشته، کجا بوده، استفاده کنیم.

توی کد زیر، $a_k$ به صورت (x * n0 + y) // z محاسبه میشه که برای $\sqrt n$ درسته چون توی اون حالت $x, y, z$ ضرایب $\frac{x\sqrt n + y}{z}$ هستن. کد زیر برای محاسبه‌ی نمایش کسر مسلسل $\alpha = \sqrt n$ هست:

import math
# محاسبه کسر مسلسل sqrt(n)
def sqrt_cf(n):
    if int(math.sqrt(n))**2 == n:
        return [int(math.sqrt(n))]

    m = 0
    d = 1
    a0 = int(math.sqrt(n))
    a = [a0]

    seen = {}

    while (m, d) not in seen:
        seen[(m,d)] = len(a)
        m = d * a[-1] - m
        d = (n - m**2) // d
        a.append(int((a0 + m) / d))

    start_period = seen[(m, d)]
    return a[:start_period], a[start_period:] # پیش‌دوره، دوره

با استفاده از یه تابع step مشابه اما با $x$, $y$ و $z$ اولیه‌ی متفاوت، میشه اون رو برای $\frac{x+y \sqrt{n}}{z}$ دلخواه هم حساب کرد.

Tavrida NU Akai Contest - Continued Fraction

$x$ و $k$ رو داری، $x$ یه مربع کامل نیست. فرض کن $\sqrt x = [a_0; a_1, \dots]$، باید $\frac{p_k}{q_k}=[a_0; a_1, \dots, a_k]$ رو برای $0 \leq k \leq 10^9$ پیدا کنی.

راه‌حل

بعد از محاسبه‌ی دوره‌ی تناوب $\sqrt x$، میشه $a_k$ رو با استفاده از توان‌رسانی دودویی (binary exponentiation) روی تبدیل کسری خطی ناشی از نمایش کسر مسلسل حساب کرد. برای پیدا کردن تبدیل حاصل، دوره‌ی تناوب به اندازه‌ی $T$ رو به یه تبدیل واحد فشرده می‌کنی و اون رو $\lfloor \frac{k-1}{T}\rfloor$ بار تکرار می‌کنی، بعدش به صورت دستی اون رو با تبدیلات باقیمانده ترکیب می‌کنی.

# کد با فرض وجود تابع sqrt_cf از مثال قبل
x, k = map(int, input().split())

mod = 10**9+7

# ترکیب (A[0]*x + A[1]) / (A[2]*x + A[3]) و (B[0]*x + B[1]) / (B[2]*x + B[3])
def combine(A, B):
    # A, B ماتریس‌های 2x2 هستند
    # (a b) (e f) = (ae+bg af+bh)
    # (c d) (g h)   (ce+dg cf+dh)
    return [ (A[0]*B[0]+A[1]*B[2]) % mod, (A[0]*B[1]+A[1]*B[3]) % mod, 
             (A[2]*B[0]+A[3]*B[2]) % mod, (A[2]*B[1]+A[3]*B[3]) % mod ]

def bpow(A, n):
    res = [1, 0, 0, 1] # ماتریس همانی
    base = A
    while n > 0:
        if n % 2 == 1:
            res = combine(res, base)
        base = combine(base, base)
        n //= 2
    return res

pre_period, period = sqrt_cf(x)

# تبدیل نهایی: p_k / q_k
# [a_0; a_1, ..., a_k]
# p_k = a_k p_{k-1} + p_{k-2}
# q_k = a_k q_{k-1} + q_{k-2}
# (p_k) = (a_k 1) (p_{k-1})
# (q_k)   (1  0) (q_{k-1})

# ماتریس برای کل تبدیل
final_matrix = [1, 0, 0, 1]

# بخش پیش‌دوره
k_rem = k
for i in range(len(pre_period)):
    if k_rem < 0: break
    mat_a = [pre_period[i], 1, 1, 0]
    final_matrix = combine(final_matrix, mat_a)
    k_rem -= 1

if k_rem >= 0:
    # ماتریس برای یک دوره کامل
    period_matrix = [1, 0, 0, 1]
    for val in period:
        mat_a = [val, 1, 1, 0]
        period_matrix = combine(period_matrix, mat_a)

    num_periods = (k_rem + 1) // len(period)
    final_matrix = combine(final_matrix, bpow(period_matrix, num_periods))
    k_rem -= num_periods * len(period)

    # بخش باقیمانده دوره
    for i in range(k_rem + 1):
        mat_a = [period[i], 1, 1, 0]
        final_matrix = combine(final_matrix, mat_a)

# [p_k, p_{k-1}], [q_k, q_{k-1}]
# p_{-2}=0, p_{-1}=1
# q_{-2}=1, q_{-1}=0
pk = final_matrix[0] * 1 + final_matrix[1] * 0
qk = final_matrix[2] * 1 + final_matrix[3] * 0

print(str(pk % mod) + '/' + str(qk % mod))

تعبیر هندسی

بیا فرض کنیم $\vec r_k = (q_k;p_k)$ بردار مربوط به همگرای $r_k = \frac{p_k}{q_k}$ باشه. اونوقت، این رابطه‌ی بازگشتی رو داریم:

$$\vec r_k = a_k \vec r_{k-1} + \vec r_{k-2}.$$

حالا فرض کن $\vec r = (1;r)$. با این حساب، هر بردار $(x;y)$ با شیبش یعنی $\frac{y}{x}$ مطابقت داره.

با استفاده از مفهوم ضرب شبه‌اسکالر (cross product) $(x_1;y_1) \times (x_2;y_2) = x_1 y_2 - x_2 y_1$، میشه نشون داد (توضیحش رو پایین‌تر می‌بینی) که:

$$s_k = -\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r} = \left|\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r}\right|.$$

معادله‌ی آخر به این دلیله که $r_{k-1}$ و $r_{k-2}$ در دو طرف مختلف $r$ قرار دارن، پس ضرب‌های شبه‌اسکالر $\vec r_{k-1}$ و $\vec r_{k-2}$ با $\vec r$ علامت‌های متفاوتی دارن. با در نظر داشتن $a_k = \lfloor s_k \rfloor$، فرمول برای $\vec r_k$ الان اینجوری میشه:

$$\vec r_k = \vec r_{k-2} + \left\lfloor \left| \frac{\vec r \times \vec r_{k-2}}{\vec r \times \vec r_{k-1}}\right|\right\rfloor \vec r_{k-1}.$$

توجه کن که $\vec r_k \times r = (q_k;p_k) \times (1;r) = q_k r - p_k$، بنابراین:

$$a_k = \left\lfloor \left| \frac{q_{k-1}r-p_{k-1}}{q_{k-2}r-p_{k-2}} \right| \right\rfloor.$$
توضیح

همونطور که قبلاً اشاره کردیم، $a_k = \lfloor s_k \rfloor$، که اینجا $s_k = [a_k; a_{k+1}, a_{k+2}, \dots]$. از طرف دیگه، از رابطه‌ی بازگشتی همگراها نتیجه می‌گیریم که:

$$r = [a_0; a_1, \dots, a_{k-1}, s_k] = \frac{s_k p_{k-1} + p_{k-2}}{s_k q_{k-1} + q_{k-2}}.$$

اگه این رو به شکل برداری بنویسیم، اینجوری میشه:

$$\vec r \parallel s_k \vec r_{k-1} + \vec r_{k-2},$$

این یعنی $\vec r$ و $s_k \vec r_{k-1} + \vec r_{k-2}$ هم‌خط هستن (یعنی شیب یکسانی دارن). با گرفتن ضرب شبه‌اسکالر هر دو طرف با $\vec r$، داریم:

$$0 = s_k (\vec r_{k-1} \times \vec r) + (\vec r_{k-2} \times \vec r),$$

که فرمول نهایی رو به ما میده:

$$s_k = -\frac{\vec r_{k-2} \times \vec r}{\vec r_{k-1} \times \vec r}.$$

الگوریتم کشیدن بینی

هر بار که $\vec r_{k-1}$ رو به بردار $\vec p$ اضافه می‌کنی، مقدار $\vec p \times \vec r$ به اندازه‌ی $\vec r_{k-1} \times \vec r$ زیاد میشه.

بنابراین، $a_k=\lfloor s_k \rfloor$ حداکثر تعداد صحیح بردارهای $\vec r_{k-1}$ هست که میشه به $\vec r_{k-2}$ اضافه کرد بدون اینکه علامت ضرب خارجی با $\vec r$ تغییر کنه.

به عبارت دیگه، $a_k$ حداکثر تعداد دفعاتیه که می‌تونی $\vec r_{k-1}$ رو به $\vec r_{k-2}$ اضافه کنی بدون اینکه از خط تعریف شده توسط $\vec r$ رد بشی:

همگراهای $r=\frac{7}{9}=[0;1,3,2]$. نیم‌همگراها با نقاط میانی بین فلش‌های خاکستری مطابقت دارن.

توی عکس بالا، $\vec r_2 = (3;4)$ با اضافه کردن مکرر $\vec r_1 = (1;1)$ به $\vec r_0 = (1;0)$ به دست میاد.

وقتی دیگه نشه $\vec r_1$ رو به $\vec r_0$ اضافه کرد بدون اینکه از خط $y=rx$ رد بشیم، میریم اون طرف خط و به طور مکرر $\vec r_2$ رو به $\vec r_1$ اضافه می‌کنیم تا $\vec r_3 = (7;9)$ رو به دست بیاریم.

این روش بردارهایی با طول نمایی تولید می‌کنه که به خط نزدیک میشن.

به خاطر این ویژگی، روش تولید بردارهای همگرای متوالی توسط بوریس دلونی الگوریتم کشیدن بینی (Algorithm of pulling noses) نامیده شده.

اگه به مثلثی که روی نقاط $\vec r_{k-2}$، $\vec r_{k}$ و $\vec 0$ رسم شده نگاه کنیم، می‌بینیم که دو برابر مساحتش برابره با:

$$|\vec r_{k-2} \times \vec r_k| = |\vec r_{k-2} \times (\vec r_{k-2} + a_k \vec r_{k-1})| = a_k |\vec r_{k-2} \times \vec r_{k-1}| = a_k.$$

در ترکیب با قضیه پیک، این یعنی هیچ نقطه‌ی شبکه‌ای (lattice point) دقیقاً داخل مثلث وجود نداره و تنها نقاط شبکه‌ای روی مرزش $\vec 0$ و $\vec r_{k-2} + t \cdot \vec r_{k-1}$ برای همه‌ی اعداد صحیح $t$ طوری که $0 \leq t \leq a_k$ هستن. وقتی این‌ها رو برای همه‌ی $k$های ممکن به هم وصل کنیم، یعنی هیچ نقطه‌ی صحیحی توی فضای بین چندضلعی‌های تشکیل شده توسط بردارهای همگرای با اندیس زوج و فرد وجود نداره.

این به نوبه‌ی خودش یعنی $\vec r_k$ با ضرایب فرد یه پوش محدب (convex hull) از نقاط شبکه‌ای با $x \geq 0$ بالای خط $y=rx$ رو تشکیل میده، در حالی که $\vec r_k$ با ضرایب زوج یه پوش محدب از نقاط شبکه‌ای با $x > 0$ زیر خط $y=rx$ رو تشکیل میده.

تعریف

این چندضلعی‌ها به عنوان چندضلعی‌های کلاین هم شناخته میشن، به نام فلیکس کلاین که اولین بار این تعبیر هندسی رو برای کسرهای مسلسل پیشنهاد کرد.

مثال‌های مسئله

حالا که مهم‌ترین نکات و مفاهیم رو یاد گرفتیم، وقتشه که بریم سراغ چند تا مثال خاص از مسائل.

پوش محدب زیر خط

پوش محدب نقاط شبکه‌ای $(x;y)$ رو طوری پیدا کن که $0 \leq x \leq N$ و $0 \leq y \leq rx$ برای $r=[a_0;a_1,\dots,a_k]=\frac{p_k}{q_k}$.

راه‌حل

اگه مجموعه‌ی نامحدود $0 \leq x$ رو در نظر می‌گرفتیم، پوش محدب بالایی توسط خود خط $y=rx$ داده میشد.

اما، با قید اضافی $x \leq N$ باید در نهایت از خط منحرف بشیم تا پوش محدب مناسبی رو حفظ کنیم.

فرض کن $t = \lfloor \frac{N}{q_k}\rfloor$، اون‌وقت $t$ نقطه‌ی شبکه‌ای اول روی پوش بعد از $(0;0)$ برابر با $\alpha \cdot (q_k; p_k)$ برای عدد صحیح $1 \leq \alpha \leq t$ هستن.

اما $(t+1)(q_k; p_k)$ نمی‌تونه نقطه‌ی شبکه‌ای بعدی باشه چون $(t+1)q_k$ از $N$ بزرگتره.

برای رسیدن به نقاط شبکه‌ای بعدی توی پوش، باید به نقطه‌ی $(x;y)$ بریم که با کمترین حاشیه از $y=rx$ منحرف میشه، در حالی که $x \leq N$ حفظ بشه.

پوش محدب نقاط شبکه‌ای زیر $y=\frac{4}{7}x$ برای $0 \leq x \leq 19$ از نقاط $(0;0), (7;4), (14;8), (16;9), (18;10), (19;10)$ تشکیل شده.

فرض کن $(x; y)$ آخرین نقطه‌ی فعلی توی پوش محدب باشه. اون‌وقت نقطه‌ی بعدی $(x'; y')$ طوریه که $x' \leq N$ و $(x'; y') - (x; y) = (\Delta x; \Delta y)$ تا حد امکان به خط $y=rx$ نزدیک باشه. به عبارت دیگه، $(\Delta x; \Delta y)$ مقدار $r \Delta x - \Delta y$ رو با شرط $\Delta x \leq N - x$ و $\Delta y \leq r \Delta x$ بیشینه می‌کنه.

نقاطی مثل این روی پوش محدب نقاط شبکه‌ای زیر $y=rx$ قرار دارن. به عبارت دیگه، $(\Delta x; \Delta y)$ باید یه نیم‌همگرای پایینی $r$ باشه.

به این ترتیب، $(\Delta x; \Delta y)$ به شکل $(q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)$ برای یه عدد فرد $i$ و $0 \leq t < a_i$ هست.

برای پیدا کردن چنین $i$، می‌تونیم همه‌ی $i$های ممکن رو از بزرگترین شروع کنیم و از $t = \lfloor \frac{N-x-q_{i-1}}{q_i} \rfloor$ برای $i$ طوری که $N-x-q_{i-1} \geq 0$ استفاده کنیم.

با $(\Delta x; \Delta y) = (q_{i-1}; p_{i-1}) + t \cdot (q_i; p_i)$، شرط $\Delta y \leq r \Delta x$ با خواص نیم‌همگرا حفظ میشه.

و $t < a_i$ برقرار خواهد بود چون ما قبلاً نیم‌همگراهای به دست اومده از $i+2$ رو تموم کردیم، از این رو $x + q_{i-1} + a_i q_i = x+q_{i+1}$ بزرگتر از $N$ هست.

حالا که می‌تونیم $(\Delta x; \Delta y)$ رو به $(x;y)$ برای $k = \lfloor \frac{N-x}{\Delta x} \rfloor$ بار اضافه کنیم قبل از اینکه از $N$ رد بشیم، بعدش نیم‌همگرای بعدی رو امتحان خواهیم کرد.

// [ah, ph, qh] را برمی‌گرداند به طوری که نقاط r[i]=(qh[i], ph[i]) پوش محدب بالایی را تشکیل می‌دهند
// از نقاط مشبکه روی 0 <= x <= N و 0 <= y <= r * x، که r = [a0; a1, a2, ...]
// و ah[i] نقطه صحیح روی پاره‌خط بین r[i-1] و r[i] وجود دارد.
auto hull(auto a, long long N) {
    auto [p, q] = convergents(a);
    long long t = N / q.back();
    vector<long long> ah = {t};
    vector<long long> ph = {0, t * p.back()};
    vector<long long> qh = {0, t * q.back()};

    for(int i = q.size() - 2; i >= 0; i--) {
        if(i % 2) {
            while(qh.back() + q[i - 1] <= N) {
                t = (N - qh.back() - q[i - 1]) / q[i];
                long long dp = p[i - 1] + t * p[i];
                long long dq = q[i - 1] + t * q[i];
                long long k = (N - qh.back()) / dq;
                ah.push_back(k);
                ph.push_back(ph.back() + k * dp);
                qh.push_back(qh.back() + k * dq);
            }
        }
    }
    return make_tuple(ah, ph, qh);
}
# [ah, ph, qh] رو برمی‌گردونه طوری که نقاط r[i]=(qh[i], ph[i]) پوش محدب بالایی رو تشکیل میدن
# از نقاط شبکه‌ای روی 0 <= x <= N و 0 <= y <= r * x، که r = [a0; a1, a2, ...]
# و ah[i] تعداد نقاط صحیح روی پاره‌خط بین r[i-1] و r[i] هست.
def hull(a, N):
    p, q = convergents(a)
    t = N // q[-1]
    ah = [t]
    ph = [0, t*p[-1]]
    qh = [0, t*q[-1]]
    # p,q شامل p_(-2), p_(-1) هستن، پس اندیس همگراها با اندیس آرایه متفاوته
    for i in reversed(range(2, len(q))):
        # همگراهای پایینی (که r_k < r) اندیس k فرد دارن
        # اندیس در آرایه: k+2
        if (i-2) % 2 == 1:
            # q[i-2] = q_{i-4}, q[i-1] = q_{i-3}
            # ما به q_{j} و q_{j-1} نیاز داریم که j فرد است. پس q[i], q[i-1]
            while qh[-1] + q[i-1] <= N:
                t = (N - qh[-1] - q[i-1]) // q[i]
                dp = p[i-1] + t*p[i]
                dq = q[i-1] + t*q[i]
                k = (N - qh[-1]) // dq
                ah.append(k)
                ph.append(ph[-1] + k * dp)
                qh.append(qh[-1] + k * dq)
    return ah, ph, qh

Timus - Crime and Punishment

اعداد صحیح $A$, $B$ و $N$ رو داری. $x \geq 0$ و $y \geq 0$ رو پیدا کن طوری که $Ax + By \leq N$ و $Ax + By$ بیشترین مقدار ممکن باشه.

راه‌حل

توی این مسئله $1 \leq A, B, N \leq 2 \cdot 10^9$ هست، پس میشه اون رو توی $O(\sqrt N)$ حل کرد. اما، یه راه‌حل $O(\log N)$ با کسرهای مسلسل وجود داره.

برای راحتی، جهت $x$ رو با جایگزینی $x \mapsto \lfloor \frac{N}{A}\rfloor - x$ برعکس می‌کنیم، طوری که الان باید نقطه‌ی $(x; y)$ رو پیدا کنیم که $0 \leq x \leq \lfloor \frac{N}{A} \rfloor$، $By - Ax \leq N \pmod A$ و $By - Ax$ بیشترین مقدار ممکن باشه. $y$ بهینه برای هر $x$ مقداری برابر با $\lfloor \frac{Ax + (N \bmod A)}{B} \rfloor$ داره.

برای اینکه به طور کلی‌تر باهاش برخورد کنیم، یه تابع می‌نویسیم که بهترین نقطه رو روی $0 \leq x \leq N$ و $y = \lfloor \frac{Ax+B}{C} \rfloor$ پیدا می‌کنه.

ایده‌ی اصلی راه‌حل توی این مسئله اساساً مسئله‌ی قبلی رو تکرار می‌کنه، اما به جای استفاده از نیم‌همگراهای پایینی برای انحراف از خط، از نیم‌همگراهای بالایی برای نزدیک‌تر شدن به خط بدون عبور از اون و بدون نقض $x \leq N$ استفاده می‌کنی. متأسفانه، بر خلاف مسئله‌ی قبلی، باید مطمئن بشی که وقتی به خط $y=\frac{Ax+B}{C}$ نزدیک میشی ازش رد نمیشی، پس باید موقع محاسبه‌ی ضریب $t$ نیم‌همگرا این رو در نظر داشته باشی.

# (x, y) طوری که y = (A*x+B) // C,
# Cy - Ax بیشینه و 0 <= x <= N.
def closest(A, B, C, N):
    # y <= (A*x + B)/C <=> diff(x, y) <= B
    def diff(x, y):
        return C*y-A*x
    a = fraction(A, C)
    p, q = convergents(a)
    # نقطه شروع (0, B//C)
    best_x, best_y = 0, B//C
    curr_x, curr_y = 0, B//C

    # همگراهای بالایی (k زوج)
    for i in range(2, len(q)):
        if (i-2) % 2 == 0:
            # امتحان q[i] و q[i-1]
            dq, dp = q[i], p[i]
            if curr_x + dq > N:
                # امتحان نیم‌همگراهای بین r_{i-2} و r_i
                dq_s, dp_s = q[i-1], p[i-1]
                if curr_x + dq_s <= N:
                    # حداکثر ضریب t
                    t = (N - curr_x - dq_s) // dq

                    # ضریب t برای اینکه از خط عبور نکنیم
                    if diff(dq, dp) < 0:
                        t = min(t, (B - diff(curr_x+dq_s, curr_y+dp_s)) // abs(diff(dq, dp)))

                    curr_x += dq_s + t*dq
                    curr_y += dp_s + t*dp
                    if diff(curr_x, curr_y) > diff(best_x, best_y):
                        best_x, best_y = curr_x, curr_y
                break

            # ضریب k برای اینکه از خط عبور نکنیم
            k = (N-curr_x) // dq
            if diff(dq, dp) < 0:
                k = min(k, (B - diff(curr_x, curr_y)) // abs(diff(dq, dp)))

            curr_x += k*dq
            curr_y += k*dp
            if diff(curr_x, curr_y) > diff(best_x, best_y):
                best_x, best_y = curr_x, curr_y
            if k < (N-curr_x) // dq : # نتونستیم کامل پیش بریم
                break
    return best_x, best_y

def solve_crime(A, B, N):
    if N < 0: return 0, 0
    # max (Ax+By) s.t. Ax+By <= N
    # بدون از دست دادن کلیت A <= B
    if A > B: A, B = B, A
    if A == 0: return 0, N//B

    # x_new = N/A - x, y_new = y
    # A(N/A - x) + By <= N -> N-Ax+By <= N -> By - Ax <= 0
    # می‌خواهیم By - Ax را بیشینه کنیم
    x, y = closest(B, 0, A, N // B)

    return (N - B*y) // A, y

June Challenge 2017 - Euler Sum

$\sum\limits_{x=1}^N \lfloor ex \rfloor$ رو حساب کن، که اینجا $e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, \dots, 1, 2n, 1, \dots]$ عدد اویلر هست و $N \leq 10^{4000}$.

راه‌حل

این مجموع برابر با تعداد نقاط شبکه‌ای $(x;y)$ هست طوری که $1 \leq x \leq N$ و $1 \leq y \leq ex$.

بعد از ساختن پوش محدب نقاط زیر $y=ex$، این تعداد رو میشه با استفاده از قضیه پیک حساب کرد:

// مجموع floor(r * x) برای x از 1 تا N و r = [a0; a1, a2, ...]
long long sum_floor(auto a, long long N) {
    N++;
    auto [ah, ph, qh] = hull(a, N);

    // تعداد نقاط مشبکه درون یک ذوزنقه قائم‌الزاویه عمودی
    // روی نقاط (0; 0) - (0; y1) - (dx; y2) - (dx; 0) که
    // a+1 نقطه صحیح روی پاره‌خط (0; y1) - (dx; y2) دارد.
    auto picks = [](long long y1, long long y2, long long dx, long long a) {
        long long b = y1 + y2 + a + dx; // نقاط مرزی
        long long A = (y1 + y2) * dx; // دو برابر مساحت
        return (A - b + 2) / 2 + b - (y2 + 1); // نقاط داخلی + نقاط روی دو ضلع
    };

    long long ans = 0;
    for(size_t i = 1; i < qh.size(); i++) {
        ans += picks(ph[i - 1], ph[i], qh[i] - qh[i - 1], ah[i - 1]);
    }
    return ans - N;
}
# مجموع floor(r * x) برای x از 1 تا N و r = [a0; a1, a2, ...]
def sum_floor_func(a, N):
    # N شامل خود نقطه پایانی نیست
    ah, ph, qh = hull(a, N)

    # تعداد نقاط مشبکه درون یک ذوزنقه قائم‌الزاویه عمودی
    # روی نقاط (0; 0) - (0; y1) - (dx; y2) - (dx; 0) که
    # a+1 نقطه صحیح روی پاره‌خط (0; y1) - (dx; y2) دارد.
    def picks(y1, y2, dx, a):
        b = y1 + y2 + a + dx
        A = (y1 + y2) * dx
        return (A - b + 2) // 2 + b - y2

    ans = 0
    for i in range(1, len(qh)):
        # مجموع floor(r*x) برای qh[i-1] < x <= qh[i]
        y_start = ph[i-1]
        y_end = ph[i]
        x_start = qh[i-1]
        x_end = qh[i]
        dx = x_end - x_start
        # تعداد نقاط روی پاره‌خط از (x_start, y_start) تا (x_end, y_end) برابر است با
        # gcd(dx, dy) + 1. در اینجا ah[i-1] تعداد پله هاست.
        ans += picks(y_start, y_end, dx, ah[i-1])

    return ans - (N + 1)

NAIPC 2019 - It's a Mod, Mod, Mod, Mod World

با داشتن $p$, $q$ و $n$, $\sum\limits_{i=1}^n (p \cdot i \bmod q)$ رو حساب کن.

راه‌حل

این مسئله به مسئله‌ی قبلی تبدیل میشه اگه یادت باشه که $a \bmod b = a - \lfloor \frac{a}{b} \rfloor b$. با این نکته، مجموع به این شکل در میاد:

$$\sum\limits_{i=1}^n \left(p \cdot i - \left\lfloor \frac{p \cdot i}{q} \right\rfloor q\right) = \frac{pn(n+1)}{2}-q\sum\limits_{i=1}^n \left\lfloor \frac{p \cdot i}{q}\right\rfloor.$$

اما، جمع کردن $\lfloor rx \rfloor$ برای $x$ از $1$ تا $N$ چیزیه که ما از مسئله‌ی قبلی بلدیم چجوری انجامش بدیم.

void solve(long long p, long long q, long long N) {
    cout << p * N * (N + 1) / 2 - q * sum_floor(fraction(p, q), N) << "\n";
}
def solve_mod_sum(p, q, N):
    return p * N * (N + 1) // 2 - q * sum_floor_func(fraction(p, q), N)

Library Checker - Sum of Floor of Linear

با داشتن $N$, $M$, $A$ و $B$, $\sum\limits_{i=0}^{N-1} \lfloor \frac{A \cdot i + B}{M} \rfloor$ رو حساب کن.

راه‌حل

این مسئله از نظر فنی دردسرسازترین مسئله تا اینجاست.

میشه از همون رویکرد استفاده کرد و پوش محدب کامل نقاط زیر خط $y = \frac{Ax+B}{M}$ رو ساخت.

ما قبلاً می‌دونیم چجوری این رو برای $B = 0$ حل کنیم. علاوه بر این، قبلاً یاد گرفتیم چجوری این پوش محدب رو تا نزدیکترین نقطه‌ی شبکه‌ای به این خط توی بازه‌ی $[0, N-1]$ بسازیم (این کار توی مسئله‌ی "Crime and Punishment" بالا انجام شده).

حالا باید توجه کنیم که وقتی به نزدیکترین نقطه به خط رسیدیم، می‌تونیم فرض کنیم که خط در واقع از اون نزدیکترین نقطه عبور می‌کنه، چون هیچ نقطه‌ی شبکه‌ای دیگه‌ای توی $[0, N-1]$ بین خط واقعی و خطی که یه کم به پایین منتقل شده تا از نزدیکترین نقطه عبور کنه، وجود نداره.

به این ترتیب، برای ساختن پوش محدب کامل زیر خط $y=\frac{Ax+B}{M}$ توی $[0, N-1]$، می‌تونیم اون رو تا نزدیکترین نقطه به خط توی $[0, N-1]$ بسازیم و بعدش ادامه بدیم انگار که خط از این نقطه عبور می‌کنه، و از الگوریتم ساخت پوش محدب با $B=0$ دوباره استفاده کنیم:

# پوش محدب نقاط شبکه‌ای (x, y) طوری که M*y <= A*x+B
def general_hull(A, B, M, N):
    a = fraction(A, M)
    p, q = convergents(a)
    ah = []
    ph = [B // M]
    qh = [0]

    # پیدا کردن نزدیکترین نقطه
    curr_x, curr_y = 0, B // M
    best_x, best_y = 0, B // M
    max_diff = M * curr_y - A * curr_x

    for i in range(2, len(q)):
        if (i-2) % 2 == 0: # همگراهای بالایی
            dq, dp = q[i], p[i]
            if curr_x + dq >= N:
                # نیم‌همگراها
                dq_s, dp_s = q[i-1], p[i-1]
                if curr_x + dq_s < N:
                    t = (N - 1 - curr_x - dq_s) // dq
                    # بررسی عبور از خط
                    if M*dp - A*dq < 0:
                        t = min(t, (B - (M* (curr_y+dp_s) - A*(curr_x+dq_s))) // abs(M*dp-A*dq))
                    curr_x += dq_s + t*dq
                    curr_y += dp_s + t*dp
                    if M * curr_y - A * curr_x > max_diff:
                        best_x, best_y = curr_x, curr_y
                        max_diff = M * curr_y - A * curr_x
                break

            k = (N-1 - curr_x) // dq
            if M*dp - A*dq < 0:
                k = min(k, (B - (M*curr_y - A*curr_x)) // abs(M*dp - A*dq))
            curr_x += k*dq
            curr_y += k*dp
            if M*curr_y - A*curr_x > max_diff:
                best_x, best_y = curr_x, curr_y
                max_diff = M*curr_y - A*curr_x
            if k < (N-1-curr_x)//dq: break

    # محاسبه مجموع تا بهترین نقطه
    # ...
    # سپس مجموع از بهترین نقطه تا N
    # ... این راه حل پیچیده است. راه حل استاندارد برای این مسئله ساده تر است.
    # راه حل استاندارد:
    if A == 0: return (B // M) * N
    if A >= M: return (A // M) * N * (N - 1) // 2 + (B // M) * N + solve_floor_sum(N, M, A % M, B % M)
    if B >= M: return N * (B // M) + solve_floor_sum(N, M, A, B % M)
    m = (A * (N-1) + B) // M
    return (N - 1) * m - solve_floor_sum(m, A, M, -B - 1 + A)

def solve_floor_sum(N, M, A, B):
    # ... (کد کامل برای تابع بازگشتی)
    pass

OKC 2 - From Modular to Rational

یه عدد گویای $\frac{p}{q}$ وجود داره طوری که $1 \leq p, q \leq 10^9$. تو می‌تونی مقدار $p q^{-1}$ به پیمانه‌ی $m \sim 10^9$ رو برای چند تا عدد اول $m$ بپرسی. $\frac{p}{q}$ رو پیدا کن.

فرمول‌بندی معادل: $x$ رو پیدا کن که $Ax \pmod M$ رو برای $1 \leq x \leq N$ کمینه کنه.

راه‌حل

به خاطر قضیه‌ی باقیمانده‌ی چینی، پرسیدن نتیجه به پیمانه‌ی چند تا عدد اول مثل اینه که اون رو به پیمانه‌ی حاصلضربشون بپرسی. به همین دلیل، بدون اینکه کلیت مسئله از بین بره، فرض می‌کنیم که باقیمانده رو به پیمانه‌ی یه عدد به اندازه‌ی کافی بزرگ $m$ می‌دونیم.

ممکنه چند تا راه‌حل ممکن $(p, q)$ برای $p \equiv qr \pmod m$ برای یه باقیمانده‌ی $r$ داده شده وجود داشته باشه. اما، اگه $(p_1, q_1)$ و $(p_2, q_2)$ هر دو راه‌حل باشن، اون‌وقت $p_1 q_2 \equiv p_2 q_1 \pmod m$ هم برقراره. با فرض اینکه $\frac{p_1}{q_1} \neq \frac{p_2}{q_2}$، این یعنی $|p_1 q_2 - p_2 q_1|$ حداقل $m$ هست.

توی صورت مسئله به ما گفته شده که $1 \leq p, q \leq 10^9$، پس اگه هم $p_1, q_1$ و هم $p_2, q_2$ حداکثر $10^9$ باشن، اون‌وقت تفاوت حداکثر $2 \cdot 10^{18}$ هست. برای $m > 2 \cdot 10^{18}$ این یعنی راه‌حل $\frac{p}{q}$ با $1 \leq p, q \leq 10^9$ به عنوان یه عدد گویا منحصر به فرده.

بنابراین، مسئله به این خلاصه میشه که با داشتن $r$ به پیمانه‌ی $m$، هر $q$ رو پیدا کنیم طوری که $1 \leq q \leq 10^9$ و $qr \pmod m \leq 10^9$.

این عملاً مثل اینه که $q$ رو پیدا کنیم که کمترین مقدار ممکن $qr \bmod m$ رو برای $1 \leq q \leq 10^9$ به دست بده.

برای $qr = km + b$ این یعنی باید یه جفت $(q, k)$ پیدا کنیم طوری که $1 \leq q \leq 10^9$ و $qr - km \geq 0$ کمترین مقدار ممکن باشه.

از اونجایی که $m$ ثابته، می‌تونیم بر اون تقسیم کنیم و مسئله رو اینجوری بازنویسی کنیم: $q$ رو پیدا کن طوری که $1 \leq q \leq 10^9$ و $\frac{r}{m} q - k \geq 0$ کمترین مقدار ممکن باشه.

از دیدگاه کسرهای مسلسل، این یعنی $\frac{k}{q}$ بهترین تقریب دیوفانتی به $\frac{r}{m}$ هست و کافیه فقط نیم‌همگراهای پایینی $\frac{r}{m}$ رو بررسی کنیم.

# پیدا کردن Q که Q*r mod m رو برای 1 <= Q <= n < m کمینه می‌کنه
def mod_min(r, n, m):
    a = fraction(r, m)
    p, q = convergents(a)
    min_val = m
    best_q = 1
    # بررسی همگراها
    for i in range(2, len(q)):
        if q[i] > n: break
        val = (q[i] * r) % m
        if val < min_val:
            min_val = val
            best_q = q[i]

    # بررسی نیم‌همگراهای آخرین همگرای معتبر
    # i-1 آخرین همگرایی است که q[i-1] <= n
    last_valid_i = i - 1
    if last_valid_i >= 2:
        # نیم‌همگراهای بین r_{i-3} و r_{i-2}
        # ما به همگراهای با اندیس زوج نیاز داریم (بالایی ها)
        # چون r - p/q < 0 -> rq - p < 0 (mod m)
        # ولی ما دنبال rq - mp > 0 هستیم پس همگراهای پایینی (اندیس فرد)
        if (last_valid_i - 2) % 2 == 1:
            prev_q, prev_p = q[last_valid_i-1], p[last_valid_i-1]
            curr_q, curr_p = q[last_valid_i], p[last_valid_i]

            t = (n - prev_q) // curr_q
            test_q = prev_q + t*curr_q
            val = (test_q*r)%m
            if val < min_val:
                min_val = val
                best_q = test_q

    return best_q, min_val

مسائل تمرینی