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

ضرب مونتگومری

بسیاری از الگوریتم‌ها در نظریه اعداد، مانند آزمایش اول بودن یا تجزیه اعداد، و در رمزنگاری، مانند RSA، به عملیات‌های زیادی به پیمانه یک عدد بزرگ نیاز دارند. محاسبه یک ضرب مانند $x y \bmod{n}$ با الگوریتم‌های معمول بسیار کند است، زیرا برای دانستن اینکه چند بار باید $n$ از حاصل‌ضرب کم شود، به یک عملیات تقسیم نیاز دارد. و تقسیم، به خصوص با اعداد بزرگ، یک عملیات بسیار پرهزینه است.

ضرب (پیمانه‌ای) مونتگومری روشی است که امکان محاسبه سریع‌تر چنین ضرب‌هایی را فراهم می‌کند. به جای تقسیم حاصل‌ضرب و کم کردن چندباره $n$، این روش مضارب $n$ را اضافه می‌کند تا بیت‌های کم‌ارزش را حذف کند و سپس به سادگی بیت‌های کم‌ارزش را نادیده می‌گیرد.

نمایش مونتگومری

با این حال، ضرب مونتگومری رایگان نیست. این الگوریتم تنها در فضای مونتگومری کار می‌کند. و ما باید قبل از شروع ضرب، اعداد خود را به آن فضا تبدیل کنیم.

برای این فضا، ما به یک عدد صحیح مثبت $r \ge n$ نیاز داریم که نسبت به $n$ اول باشد، یعنی $\gcd(n, r) = 1$. در عمل، ما همیشه $r$ را به صورت $2^m$ برای یک عدد صحیح مثبت $m$ انتخاب می‌کنیم، زیرا در این صورت عملیات ضرب، تقسیم و پیمانه $r$ را می‌توان با استفاده از شیفت‌ها و دیگر عملیات بیتی به طور کارآمد پیاده‌سازی کرد. $n$ تقریباً در تمام کاربردها یک عدد فرد خواهد بود، زیرا تجزیه یک عدد زوج کار دشواری نیست. بنابراین هر توانی از $2$ نسبت به $n$ اول خواهد بود.

نمایانگر $\bar{x}$ از عدد $x$ در فضای مونتگومری به صورت زیر تعریف می‌شود:

$$\bar{x} := x \cdot r \bmod n$$

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

درون فضای مونتگومری، شما هنوز هم می‌توانید بیشتر عملیات را به صورت معمول انجام دهید. شما می‌توانید دو عنصر را با هم جمع کنید ($x \cdot r + y \cdot r \equiv (x + y) \cdot r \bmod n$)، تفریق کنید، برابری را بررسی کنید و حتی بزرگترین مقسوم‌علیه مشترک یک عدد با $n$ را محاسبه کنید (زیرا $\gcd(n, r) = 1$). همه اینها با الگوریتم‌های معمول انجام می‌شود.

اما این موضوع در مورد ضرب صدق نمی‌کند.

ما انتظار داریم نتیجه به این صورت باشد:

$$\bar{x} * \bar{y} = \overline{x \cdot y} = (x \cdot y) \cdot r \bmod n.$$

اما ضرب معمولی به ما این نتیجه را می‌دهد:

$$\bar{x} \cdot \bar{y} = (x \cdot y) \cdot r \cdot r \bmod n.$$

بنابراین، ضرب در فضای مونتگومری به صورت زیر تعریف می‌شود:

$$\bar{x} * \bar{y} := \bar{x} \cdot \bar{y} \cdot r^{-1} \bmod n.$$

کاهش مونتگومری

ضرب دو عدد در فضای مونتگومری نیازمند محاسبه کارآمد $x \cdot r^{-1} \bmod n$ است. این عملیات کاهش مونتگومری نامیده می‌شود و با نام الگوریتم REDC نیز شناخته می‌شود.

از آنجا که $\gcd(n, r) = 1$، می‌دانیم که دو عدد $r^{-1}$ و $n^{\prime}$ با شرط $0 < r^{-1}, n^{\prime} < n$ وجود دارند که:

$$r \cdot r^{-1} + n \cdot n^{\prime} = 1.$$

هر دو عدد $r^{-1}$ و $n^{\prime}$ را می‌توان با استفاده از الگوریتم اقلیدسی تعمیم‌یافته محاسبه کرد.

با استفاده از این اتحاد، می‌توانیم $x \cdot r^{-1}$ را به صورت زیر بنویسیم:

$$\begin{aligned} x \cdot r^{-1} &= x \cdot r \cdot r^{-1} / r = x \cdot (-n \cdot n^{\prime} + 1) / r \\ &= (-x \cdot n \cdot n^{\prime} + x) / r \equiv (-x \cdot n \cdot n^{\prime} + l \cdot r \cdot n + x) / r \bmod n\\ &\equiv ((-x \cdot n^{\prime} + l \cdot r) \cdot n + x) / r \bmod n \end{aligned}$$

این هم‌نهشتی‌ها برای هر عدد صحیح دلخواه $l$ برقرار است. این بدان معناست که ما می‌توانیم یک مضرب دلخواه از $r$ را به $x \cdot n^{\prime}$ اضافه یا از آن کم کنیم، یا به عبارت دیگر، می‌توانیم $q := x \cdot n^{\prime}$ را به پیمانه $r$ محاسبه کنیم.

این موضوع الگوریتم زیر را برای محاسبه $x \cdot r^{-1} \bmod n$ به ما می‌دهد:

function reduce(x):
    q = (x mod r) * n' mod r
    a = (x - q * n) / r
    if a < 0:
        a += n
    return a

از آنجا که $x < n \cdot n < r \cdot n$ (حتی اگر $x$ حاصل یک ضرب باشد) و $q \cdot n < r \cdot n$، می‌دانیم که $-n < (x - q \cdot n) / r < n$. بنابراین عملیات پیمانه نهایی با استفاده از یک بررسی و یک جمع پیاده‌سازی می‌شود.

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

کاربرد دوم کاهش مونتگومری، بازگرداندن یک عدد از فضای مونتگومری به فضای عادی است.

ترفند معکوس سریع

برای محاسبه کارآمد معکوس $n^{\prime} := n^{-1} \bmod r$، می‌توانیم از ترفند زیر استفاده کنیم (که از روش نیوتن الهام گرفته شده است):

$$a \cdot x \equiv 1 \bmod 2^k \Longrightarrow a \cdot x \cdot (2 - a \cdot x) \equiv 1 \bmod 2^{2k}$$

این را می‌توان به راحتی اثبات کرد. اگر داشته باشیم $a \cdot x = 1 + m \cdot 2^k$، آنگاه داریم:

$$\begin{aligned} a \cdot x \cdot (2 - a \cdot x) &= 2 \cdot a \cdot x - (a \cdot x)^2 \\ &= 2 \cdot (1 + m \cdot 2^k) - (1 + m \cdot 2^k)^2 \\ &= 2 + 2 \cdot m \cdot 2^k - 1 - 2 \cdot m \cdot 2^k - m^2 \cdot 2^{2k} \\ &= 1 - m^2 \cdot 2^{2k} \\ &\equiv 1 \bmod 2^{2k}. \end{aligned}$$

این بدان معناست که می‌توانیم با $x=1$ به عنوان معکوس $a$ به پیمانه $2^1$ شروع کنیم، این ترفند را چند بار اعمال کنیم و در هر تکرار، تعداد بیت‌های صحیح $x$ را دو برابر کنیم.

پیاده‌سازی

با استفاده از کامپایلر GCC، وقتی هر سه عدد، اعداد صحیح ۶۴ بیتی باشند، هنوز هم می‌توانیم $x \cdot y \bmod n$ را به طور کارآمد محاسبه کنیم، زیرا این کامپایلر از اعداد صحیح ۱۲۸ بیتی با نوع‌های __int128 و __uint128 پشتیبانی می‌کند.

long long result = (__int128)x * y % n;

اما نوعی برای اعداد صحیح ۲۵۶ بیتی وجود ندارد. بنابراین در اینجا یک پیاده‌سازی برای ضرب ۱۲۸ بیتی نشان خواهیم داد.

using u64 = uint64_t;
using u128 = __uint128_t;
using i128 = __int128_t;

struct u256 {
    u128 high, low;

    static u256 mult(u128 x, u128 y) {
        u64 a = x >> 64, b = x;
        u64 c = y >> 64, d = y;
        // (a*2^64 + b) * (c*2^64 + d) =
        // (a*c) * 2^128 + (a*d + b*c)*2^64 + (b*d)
        u128 ac = (u128)a * c;
        u128 ad = (u128)a * d;
        u128 bc = (u128)b * c;
        u128 bd = (u128)b * d;
        u128 carry = (u128)(u64)ad + (u128)(u64)bc + (bd >> 64u);
        u128 high = ac + (ad >> 64u) + (bc >> 64u) + (carry >> 64u);
        u128 low = (ad << 64u) + (bc << 64u) + bd;
        return {high, low};
    }
};

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;
    }

    u128 init(u128 x) {
        x %= mod;
        for (int i = 0; i < 128; i++) {
            x <<= 1;
            if (x >= mod)
                x -= mod;
        }
        return x;
    }

    u128 reduce(u256 x) {
        u128 q = x.low * inv;
        i128 a = x.high - u256::mult(q, mod).high;
        if (a < 0)
            a += mod;
        return a;
    }

    u128 mult(u128 a, u128 b) {
        return reduce(u256::mult(a, b));
    }

    u128 mod, inv;
};

تبدیل سریع

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

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

$$\bar{x} := x \cdot r \bmod n = x \cdot r^2 / r = x * r^2$$

تبدیل یک عدد به این فضا، صرفاً یک ضربِ آن عدد در $r^2$ داخل همان فضا است. بنابراین می‌توانیم $r^2 \bmod n$ را از قبل محاسبه کنیم و به جای ۱۲۸ بار شیفت دادن عدد، فقط یک ضرب انجام دهیم.

در کد زیر، ما r2 را با -n % n مقداردهی اولیه می‌کنیم که معادل $r - n \equiv r \bmod n$ است، سپس آن را ۴ بار شیفت می‌دهیم تا به $r \cdot 2^4 \bmod n$ برسیم. این عدد را می‌توان به عنوان $2^4$ در فضای مونتگومری تفسیر کرد. اگر آن را ۵ بار به توان دو برسانیم، به $(2^4)^{2^5} = (2^4)^{32} = 2^{128} = r$ در فضای مونتگومری می‌رسیم که دقیقاً همان $r^2 \bmod n$ است.

struct Montgomery {
    Montgomery(u128 n) : mod(n), inv(1), r2(-n % n) {
        for (int i = 0; i < 7; i++)
            inv *= 2 - n * inv;

        for (int i = 0; i < 4; i++) {
            r2 <<= 1;
            if (r2 >= mod)
                r2 -= mod;
        }
        for (int i = 0; i < 5; i++)
            r2 = mul(r2, r2);
    }

    u128 init(u128 x) {
        return mult(x, r2);
    }

    u128 mod, inv, r2;
};