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

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

بیا این مسئله رو در نظر بگیریم: می‌خوایم $ax^y \pmod{2^d}$ رو حساب کنیم. فرض کن که $a$، $x$، $y$ و $d$ عدد صحیح هستن، $d \geq 3$ و مهم‌تر از همه، $x$ فرده.

الگوریتمی که می‌خوایم در موردش حرف بزنیم، این مسئله رو با $O(d)$ تا عمل جمع و عملیات بیتی (دودویی) و فقط یه دونه ضرب در $y$ حل می‌کنه.

ببین، به خاطر ساختار خاص گروه ضربی به پیمانه $2^d$، هر عدد فرد $x$ که شرط $x \equiv 1 \pmod 4$ رو داشته باشه، میشه اینجوری نمایشش داد:

$$ x \equiv b^{L(x)} \pmod{2^d}, $$

اینجا $b$ یه عدده که $b \equiv 5 \pmod 8$. برای اینکه کارمون راحت‌تر بشه و کلیت مسئله هم حفظ بشه، فرض می‌کنیم که $x \equiv 1 \pmod 4$. چرا؟ چون اگه $x \equiv 3 \pmod 4$ بود، خیلی راحت می‌تونیم با جایگزین کردن $x$ با $-x$ و $a$ با $(-1)^{y} a$، مسئله رو به همون حالت $x \equiv 1 \pmod 4$ تبدیل کنیم. پس چیزی رو از دست نمیدیم. با این حساب، $ax^y$ این شکلی میشه:

$$ a x^y \equiv a b^{yL(x)} \pmod{2^d}. $$

ایده اصلی این الگوریتم اینه که از این حقیقت که داریم به پیمانه $2^d$ کار می‌کنیم، استفاده کنیم تا محاسبه $L(x)$ و $b^{y L(x)}$ رو ساده‌تر کنیم. حالا یه نکته‌ای هست: به دلایلی که جلوتر می‌بینی، ما به جای خود $L(x)$، با $4L(x)$ کار می‌کنیم. البته این مقدار به پیمانه $2^d$ گرفته میشه، نه $2^{d-2}$.

توی این مقاله، پیاده‌سازی رو برای عددهای صحیح ۳۲ بیتی با هم می‌بینیم. فرض کن:

  • mbin_log_32(r, x) تابعی باشه که $r+4L(x) \pmod{2^d}$ رو حساب می‌کنه؛
  • mbin_exp_32(r, x) تابعی باشه که $r b^{\frac{x}{4}} \pmod{2^d}$ رو حساب می‌کنه؛
  • mbin_power_odd_32(a, x, y) تابعی باشه که $ax^y \pmod{2^d}$ رو حساب می‌کنه.

با این فرض‌ها، تابع mbin_power_odd_32 اینجوری پیاده‌سازی میشه:

uint32_t mbin_power_odd_32(uint32_t rem, uint32_t base, uint32_t exp) {
    if (base & 2) {
        /* پایه رو منفی در نظر می‌گیریم */
        base = -base;
        /* چک می‌کنیم ببینیم نتیجه باید منفی بشه یا نه */
        if (exp & 1) {
            rem = -rem;
        }
    }
    return (mbin_exp_32(rem, mbin_log_32(0, base) * exp));
}

محاسبه 4L(x) از روی x

خب، فرض کن $x$ یه عدد فرده که $x \equiv 1 \pmod 4$. میشه اون رو اینجوری نمایش داد:

$$ x \equiv (2^{a_1}+1)\dots(2^{a_k}+1) \pmod{2^d}, $$

که اینجا $1 < a_1 < \dots < a_k < d$. تابع $L(\cdot)$ برای هر کدوم از این ضریب‌ها خوش‌تعریفه، چون همه‌شون به پیمانه ۴ برابر با ۱ هستن. در نتیجه:

$$ 4L(x) \equiv 4L(2^{a_1}+1)+\dots+4L(2^{a_k}+1) \pmod{2^{d}}. $$

پس، اگه ما مقدار $t_n = 4L(2^n+1)$ رو برای همه $n$ هایی که بین $1$ و $d$ هستن از قبل حساب کرده باشیم، می‌تونیم $4L(x)$ رو برای هر $x$ دلخواهی حساب کنیم.

برای عددهای صحیح ۳۲ بیتی، می‌تونیم از این جدول از پیش محاسبه شده استفاده کنیم:

const uint32_t mbin_log_32_table[32] = {
    0x00000000, 0x00000000, 0xd3cfd984, 0x9ee62e18,
    0xe83d9070, 0xb59e81e0, 0xa17407c0, 0xce601f80,
    0xf4807f00, 0xe701fe00, 0xbe07fc00, 0xfc1ff800,
    0xf87ff000, 0xf1ffe000, 0xe7ffc000, 0xdfff8000,
    0xffff0000, 0xfffe0000, 0xfffc0000, 0xfff80000,
    0xfff00000, 0xffe00000, 0xffc00000, 0xff800000,
    0xff000000, 0xfe000000, 0xfc000000, 0xf8000000,
    0xf0000000, 0xe0000000, 0xc0000000, 0x80000000,
};

در عمل، ما یه رویکرد یه کم متفاوت‌تر از چیزی که بالا گفتم رو به کار می‌بریم. به جای اینکه تجزیه خود $x$ رو پیدا کنیم، میایم $x$ رو پشت سر هم توی $2^n+1$ ضرب می‌کنیم تا وقتی که به پیمانه $2^d$ برابر با $1$ بشه. اینجوری، در واقع ما نمایش $x^{-1}$ رو پیدا می‌کنیم، یعنی:

$$ x (2^{a_1}+1)\dots(2^{a_k}+1) \equiv 1 \pmod {2^d}. $$

برای این کار، روی $n$ از ۲ تا $d-1$ یه حلقه می‌زنیم. اگه بیت $n$-امِ $x$ فعلی یک بود، $x$ رو در $2^n+1$ ضرب می‌کنیم. این کار توی C++ خیلی راحت با x = x + (x << n) انجام میشه. این عملیات بیت‌های کمتر از $n$ رو دست‌نخورده باقی میذاره، ولی بیت $n$-ام رو صفر می‌کنه (چون $x$ فرده).

با در نظر گرفتن همه اینا، تابع mbin_log_32(r, x) اینجوری پیاده‌سازی میشه:

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }

    return r;
}

یادت باشه که $4L(x) = -4L(x^{-1})$، برای همین به جای اینکه $4L(2^n+1)$ رو جمع کنیم، اون رو از $r$ (که اولش صفره) کم می‌کنیم.

محاسبه x از روی 4L(x)

دقت کن که برای $k \geq 1$ این رابطه برقراره:

$$ (a 2^{k}+1)^2 = a^2 2^{2k} +a 2^{k+1}+1 = b2^{k+1}+1, $$

که از این رابطه (با چند بار به توان رسوندن) میشه نتیجه گرفت که:

$$ (2^a+1)^{2^b} \equiv 1 \pmod{2^{a+b}}. $$

با استفاده از این نتیجه، میشه فهمید که مرتبه ضربی $2^n+1$ یکی از مقسوم‌علیه‌های $2^{d-n}$ هست.

این به نوبه خودش یعنی $L(2^n+1)$ باید بر $2^{n}$ بخش‌پذیر باشه. چرا؟ چون مرتبه $b$ برابر $2^{d-2}$ هست و مرتبه $b^y$ میشه $2^{d-2-v}$ (که اینجا $2^v$ بزرگترین توان ۲ هست که $y$ رو عاد می‌کنه). پس باید این شرط برقرار باشه:

$$ 2^{d-k} \equiv 0 \pmod{2^{d-2-v}}, $$

که یعنی $v$ باید بزرگتر یا مساوی $k-2$ باشه. این قضیه یه کم کار رو سخت می‌کنه. برای همین بود که از اول گفتیم $L(x)$ رو در $4$ ضرب می‌کنیم تا از این مشکل خلاص بشیم. حالا اگه $4L(x)$ رو داشته باشیم، می‌تونیم با چک کردن بیت‌هاش، اون رو به صورت یکتا به مجموعی از جمله‌های $4L(2^n+1)$ تجزیه کنیم. اگه بیت $n$-ام یک بود، نتیجه رو در $2^n+1$ ضرب می‌کنیم و مقدار $4L(x)$ فعلی رو به اندازه $4L(2^n+1)$ کم می‌کنیم.

بنابراین، mbin_exp_32 این شکلی پیاده‌سازی میشه:

uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n < 32; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }

    return r;
}

بهینه‌سازی‌های بیشتر

اگه دقت کنی که $4L(2^{d-1}+1)=2^{d-1}$ و برای $2n \geq d$ این رابطه برقراره:

$$ (2^n+1)^2 \equiv 2^{2n} + 2^{n+1}+1 \equiv 2^{n+1}+1 \pmod{2^d}, $$

که بهمون اجازه میده نتیجه بگیریم برای $2n \geq d$ رابطه $4L(2^n+1)=2^n$ برقراره، میشه تعداد تکرارهای حلقه رو نصف کرد. یعنی می‌تونی الگوریتم رو با اجرای حلقه فقط تا $\frac{d}{2}$ ساده کنی و بعدش با استفاده از همین نکته، بقیه کار رو با چندتا عملیات بیتی ساده انجام بدی:

uint32_t mbin_log_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            x = x + (x << n);
            r -= mbin_log_32_table[n];
        }
    }

    r -= (x & 0xFFFF0000);

    return r;
}

uint32_t mbin_exp_32(uint32_t r, uint32_t x) {
    uint8_t n;

    for (n = 2; n != 16; n++) {
        if (x & (1 << n)) {
            r = r + (r << n);
            x -= mbin_log_32_table[n];
        }
    }

    r *= 1 - (x & 0xFFFF0000);

    return r;
}

چطوری جدول لگاریتم رو بسازیم؟

برای ساختن اون جدول لگاریتم، میشه الگوریتم Pohlig–Hellman رو برای حالتی که پیمانه توانی از ۲ هست، یه کم تغییر داد.

کار اصلی ما اینجا اینه که $x$ رو طوری حساب کنیم که $g^x \equiv y \pmod{2^d}$، که اینجا $g=5$ و $y$ یه عددی از جنس $2^n+1$ هست.

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

$$ g^{2^k x} \equiv y^{2^k} \pmod{2^d}. $$

یادت باشه که مرتبه $g$ از $2^d$ بزرگتر نیست (در واقع، از $2^{d-2}$ هم بزرگتر نیست، ولی برای سادگی با $2^d$ کار می‌کنیم). پس اگه $k=d-1$ رو بذاریم، سمت چپ یا $g^1$ میشه یا $g^0$. این به ما اجازه میده کم‌ارزش‌ترین بیت $x$ رو با مقایسه $y^{2^k}$ با $g$ پیدا کنیم. حالا فرض کن $x=x_0 + 2^k x_1$، که $x_0$ بخش معلوم و $x_1$ بخش نامعلوم ماست. اونوقت:

$$ g^{x_0+2^k x_1} \equiv y \pmod{2^d}. $$

اگه دو طرف رو در $g^{-x_0}$ ضرب کنیم، داریم:

$$ g^{2^k x_1} \equiv (g^{-x_0} y) \pmod{2^d}. $$

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

منابع