ইন্টিজার ফ্যাক্টোরাইজেশন পর্ব ০

গত পর্বে আমরা গায়ের জোরে কিছু সংখ্যার উৎপাদক বের করার চেষ্টা করেছিলাম। দুঃখজনকভাবে আমাদের বাচ্চা প্রোগ্রামটি \(20\) ডিজিটের একটি সংখ্যাকে উৎপাদকে বিশ্লেষণ করতে গিয়ে ক্লান্ত হয়ে গিয়েছিল। \(50\) বা \(100\) ডিজিটের কোন সংখ্যাকে উৎপাদকে বিশ্লেষণ করতে হলে সেই প্রোগ্রামের কী অবস্থা হবে ভাব একবার!

মানুষ যখন দেখল বড় বড় সংখ্যার জন্য \(2\) থেকে \(N-1\) [কিংবা \(\sqrt N\)] পর্যন্ত সব সংখ্যা দিয়ে \(N\)-কে এমন ভাগ করা অনেক সময়সাপেক্ষ, তখন তারা ট্রিক খুঁজতে লাগল কী করে আরও কম সংখ্যক সংখ্যা দিয়ে ভাগ করে \(N\)-এর উৎপাদক বের করা যায়। এমন অনেকগুলো ট্রিক রয়েছে, যার মাঝে একটি হল জন্মদিন সমস্যা(পরিশিষ্টে দেখ)। একে কাজে লাগিয়ে জন্ম হয়েছে একটি চালাকচতুর এলগোরিদমের।

পোলার্ডের \(\rho\) এলগোরিদম

মনে কর, আমাদেরকে \(N=55\)-কে উৎপাদকে বিশ্লেষণ করতে হবে। (আমরা যদিও জানি যে \(55\) এর একটি উৎপাদক \(p=5\), তবুও আমরা সেটি এলগোরিদম বের করব।)

আমরা করব কি, র‍্যান্ডম পূর্ণসংখ্যার একটি অসীম ধারা (infinite sequence) \(a_0, a_1, a_2, a_3, a_4,\ldots\) [সংক্ষেপে \(\{a_i\}_{i = 0}^{\infty}\)] থেকে খেয়াল খুশিমত কিছু পদ, যেমন: \(a_0, a_3, a_7, a_{100}\), তুলে নেব। জন্মদিন সমস্যা অনুযায়ী, এই চারটি পদের মাঝে \(\pmod {5}\)-এ কনগ্রুয়েন্ট অন্তত দুটি পদ পাওয়ার সম্ভাবনা\[1-e^{-\frac{4^2}{2\times 5}}\approx 80\%\]এমন দুটি পদ \(a_m\) এবং \(a_n\) যদি সত্যি সত্যি থাকে, তাহলে \[a_m\equiv a_n\pmod{5}\]\[\implies a_m-a_n\equiv 0\pmod{5}\]\[\implies d = \gcd(55, a_m-a_n)\ge 5\] অর্থাৎ, আমরা \(55\)-এর একটি উৎপাদক \(d\) নির্ণয় করে ফেলব। অবশ্য কেউ প্রশ্ন করতেই পারে \(d = 55\) হলে? হ্যাঁ, সেটা হতে পারে। কিন্তু চারটি সংখ্যার মাঝে \(\pmod{55}\)-এ কনগ্রুয়েন্ট দুটি সংখ্যা পাওয়ার সম্ভাবনা মাত্র\[1-e^{-\frac{4^2}{2\times 55}}\approx 13.5\%\]

অতএব, আমরা যদি \(\{a_i\}_{i = 0}^{\infty}\) থেকে প্রতিবার যেকোন দুইটি পদ \(a_i\) এবং \(a_j\) নিয়ে \(\gcd(N, a_i-a_j)\) নির্ণয় করে যেতে থাকি, তাহলে খুবই ভাল সম্ভাবনা আছে যে আমরা \(a_m-a_n\) পেয়ে যাব এবং \(55\)-এর একটি প্রকৃত উৎপাদক \(d\) নির্ণয় করে ফেলব।

তবে ঝামেলা হচ্ছে খাঁটি র‍্যান্ডম সংখ্যার অসীম ধারা ব্যবহার করলে সেটির অনেকগুলো পদ আমাদেরকে মেমরিতে রাখতে হবে, ফলে \(50\) বা \(100\) ডিজিটের কোন সংখ্যার জন্য কোটি কোটি পেটাবাইট মেমরির প্রয়োজন হবে। তাই আমরা এমন কোন সংখ্যার ধারা নেব যার মান প্রকৃতপক্ষে ডিটারমিনিস্টিক, অর্থাৎ সূত্র মেনে চলে, কিন্তু \(\pmod p\)-তে মোটামুটি র‍্যান্ডমলি আসে। (এমন সংখ্যাকে বলে সুডোর‍্যান্ডম সংখ্যা)

উদাহরণস্বরূপ, অসীম ধারা হিসেবে \(a_0=1\), এবং \(a_i=a_{i-1}^2+1\), অর্থাৎ \(1, 2, 5, 26, 677,\ldots\) ধারাটি নেওয়া নেব এবং প্রতিটি পদের জন্য \(\gcd(N, a_i-a_0)\) নির্ণয় করব। লক্ষ্য কর যে, আমাদের \(a_i\)-এর প্রকৃত মানের পরিবর্তে \(a_i\pmod N\) জানলেই চলে।

def rhonaive(N):
    a_0 = d = 1
    a_i = a_0
    while d == 1:
        a_i = (a_i * a_i + 1) % N
        d = gcd(a_i - a_0, N)
    return d

অনেকসময় দেখা যেতে পারে আমাদের নেওয়া সুডোর‍্যান্ডম সংখ্যার ধারাটি “যথেষ্ট র‍্যান্ডম” নয়। তেমন ক্ষেত্রে আমাদের \(a_i = a_{i-1}^2+1\) এর পরিবর্তে আরেকটু জটিল কোন সম্পর্ক ব্যবহার করতে হবে। তাই, সাধারণভাবে(generally) আমরা \(a_i=f(a_{i-1})\) ধরতে পারি যেখানে \(f(x)\) হচ্ছে একটি সুবিধাজনক বহুপদী। যেমন- উপরের প্রোগ্রামে \(f(x)=x^2+1\).

লক্ষ্য কর যে, \(\{a_i\}_{i = 0}^{\infty}\) ধারার প্রতিটি পদ শুধুমাত্র আগের পদের উপর নির্ভরশীল। তাই \(a_m \equiv a_n\pmod p\) একবার পাওয়া গেলে ধারাতে পূর্ববর্তী পদগুলো সাইক্লিকভাবে ঘুরে ফিরে আসতে থাকবে। অর্থাৎ,\(a_{m+i}\equiv a_{n+i}\pmod p\).

তবে এই সাইকেল \(a_0\) থেকে শুরু হতে হবে এমন কোন কথা নেই। যেমন- \(p=29\) এবং \(f(x)=x^2+1\) হলে, \(\pmod p\)-তে \(\{a_i\}_{i = 0}^{\infty}\) ধারাটি হয়

দেখতেই পাচ্ছ, \(a_0\)-র কনগ্রুয়েন্ট কোন পদ দ্বিতীয়বার পাওয়া যায়নি। তাই, প্রোগ্রাম ১-এ সবসময় \(a_0\) থেকে সাইকেল শুরু হবে ধরে নেওয়াটা ভুল ছিল।

এমন ধারাতে সাইকেল (অর্থাৎ \(\pmod p\)-তে কনগ্রুয়েন্ট দুটি পদ) নির্ণয় করতে আরেকটু বুদ্ধি খাটাতে হয়। আমি এখানে দুটি এলগোরিদমের কথা বলে দিচ্ছি। প্রথমটি ফ্লয়েডের এলগোরিদম, দ্বিতীয়টি ব্রেন্টের এলগোরিদম।

ফ্লয়েডের এলগোরিদম

মনে কর \(a_i\) ধারার সাইকেলের দৈর্ঘ্য \(l\). \(n\) যদি যথেষ্ট বড় হয়, তবে,\[a_n\equiv a_{n+l}\equiv a_{n+2l}\equiv\cdots\equiv a_{n+kl}\pmod p\]
আমরা যদি যথেষ্ট বড় কোন \(k\)-র জন্য \(n=kl\) নেই, তবে,\[a_n \equiv a_{n+kl}=a_{2kl}=a_{2n}\pmod p\]
অর্থাৎ, \(a_i\) ধারায় এমন অসীম সংখ্যক পদ \(a_n\) আছে যাদের জন্য \(a_n\equiv a_{2n}\pmod p\). সুতরাং \(\gcd(N, a_{2i}-a_i)\) নির্ণয় করে যেতে থাকলেই একটা সময় \(\pmod p\)-তে কনগ্রুয়েন্ট দুটি পদ পাওয়া যাবে। ফ্লয়েডের এলগোরিদম এই ধারণাটি ব্যবহার করে।

ফ্লয়েডের এলগোরিদমে একটি খরগোশ আর একটি কচ্ছপের মাঝে দৌড় প্রতিযোগিতার চলে। কচ্ছপটি ধারার একটি একটি পদে \(a_0, a_1, a_2,a_3,\ldots\) পর্যায়ক্রমে যেতে থাকে। একইসময়ে, খরগোশটি ধারার \(a_0, a_2, a_4, a_6,\ldots\) পদগুলোতে লাফিয়ে লাফিয়ে যেতে থাকে। আমরা খরগোশের অবস্থান এবং কচ্ছপের অবস্থানের পার্থক্যের সাথে \(N\)-এর গসাগু নির্ণয় করতে থাকি। যখনই গসাগু \(1\) এর চেয়ে বড় হয়ে যায়, এলগোরিদমটি টার্মিনেট করে।

def rhofloyd(N):
    #both hare and tortoise at a_0
    tortoise = hare = d = 1
    while d == 1:
        tortoise = f(tortoise) % N #tortoise steps once
        hare = f(f(hare)) % N #hare steps twice
        d = gcd(hare - tortoise, N)
    return d

ব্রেন্টের এলগোরিদম

ব্রেন্টের এলগোরিদমে অনেকগুলো রাউন্ডে খরগোশ আর কচ্ছপের দৌড় প্রতিযোগিতা চলে। প্রতিযোগিতার \(i\)-তম রাউন্ডে খরগোশ আর কচ্ছপ উভয়েই \(a_{2^i}\)-তম পদে অবস্থান করে। প্রতিযোগিতা শুরুর পর খরগোশ তার জায়গায় আরাম করে ঘুম দেয়, আর কচ্ছপটি টুক টুক করে একটি একটি পদ সামনে এগুতে থাকে। আমরা প্রতিবারই খরগোশের অবস্থান এবং কচ্ছপের অবস্থানের পার্থক্যের সাথে \(N\)-এর গসাগু \(d\) নির্ণয় করি। যখন কচ্ছপ খরগোশের চেয়ে \(2^i\) পদ সামনে এগিয়ে যায়, খরগোশের তড়াক করে জেগে ওঠে এবং সে এক লাফে কচ্ছপের অবস্থানে চলে আসে। এরপর প্রতিযোগিতার \(i+1\)-তম রাউন্ড শুরু হয়।

আগের মতই, \(d> 1\) হয়ে গেলে এলগোরিদমটি টার্মিনেট করবে।

লক্ষ্য কর যে, \(i\) এর মান বাড়ার সাথে সাথে \(2^{i}\) এবং \(2^{i+1}\)-এর মধ্যবর্তী পার্থক্যও বাড়তে থাকে। একটা সময় এই পার্থক্য \(\{a_i\}_{i = 0}^\infty\) ধারার সাইকেলের দৈর্ঘ্য \(l\) থেকে বড় হয়ে যায়। তখন কচ্ছপের কোন এক অবস্থানের জন্য খরগোশ ও কচ্ছপের অবস্থানের পার্থক্য হবে ঠিক \(l\)টি পদ। এই অবস্থাতেই আসলে এলগোরিদমটি টার্মিনেট করবে।

def rhobrent(N):
    #both hare and tortoise are at a_(2^0)
    #power is the nearest 2^i behind tortoise
    tortoise = hare = f(1)
    power, lead, d = 1, 0, 1
    while d == 1:
        if lead == power: #tortoise reaches 2^(i+1)-th term
            hare = tortoise #hare jumps to tortoise's position
            power *= 2 #update power
            lead = 0 #reset lead, since both hare and tortoise 
                     #are in the same position now
        tortoise = f(tortoise) % N #tortoise steps once
        lead += 1 #tortoise leads the race by one more step
        d = gcd(hare - tortoise, N)
    return d

ব্রেন্টের এনালাইসিস অনুসারে, ব্রেন্টের এলগোরিদম ফ্লয়েডের এলগোরিদম থেকে \(24\%\) দ্রুততর।

পূর্ণাঙ্গ প্রোগ্রাম

ব্রেন্টের এলগোরিদম দিয়ে যেহেতু সাইকেল শনাক্তকরণ সবচেয়ে দ্রুত করা যায়, তাই সেটি দিয়েই নিচে একটি পূর্ণাঙ্গ প্রোগ্রাম দেওয়া হল।

from fractions import gcd
from time import time

def f(x):
    return (x**2 + x + 1) % N

def rhobrent(N):
    tortoise = hare = f(1)
    power, lead, d = 1, 0, 1
    while abs(d) == 1:
        if lead == power:
            hare, lead = tortoise, 0
            power *= 2
        tortoise = f(tortoise)
        lead += 1
        d = gcd(tortoise - hare, N)
    return d

N = int(raw_input())
exec_time = time()
d = 1 if(-2 < N < 2) else rhobrent(N)
exec_time = time() - exec_time
print d, '*', N / d
print "Execution time: %.3f s" % exec_time

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <iso646.h> //and

#define MAX 2147483647

long long f(long long x);
long long gcd(long long a, long long b);
long long rhobrent(long long N);

long long N;

int main(void){
    long long d;
    scanf("%lld", &N);
    clock_t exec_time = clock();
    d = (N < 2 and N > -2) ? 1 : rhobrent(N);
    exec_time = clock() - exec_time;
    double t = (double) exec_time / CLOCKS_PER_SEC;
    printf("%lld * %lld\n", d, N/d);
    printf("Execution time: %.3lf s\n", t);
    return 0;
}

long long f(long long x){
    if(x < MAX) 
        return (x * x + x + 1) % N;
    long long y = x, output = x + 1;
    while(y > 0){
        if(y % 2 == 1)
            output += x;
        x *= 2;
        y /= 2;
        output %= N;
        x %= N;
    }
    return output;
}

long long gcd(long long a, long long b){
    long long r;
    while(b != 0){
        r =  a % b;
        a = b;
        b = r;
    }
    return (a > 0) ? a : -a;
}

long long rhobrent(long long N){
    long long tortoise = f(1), hare = f(1), power = 1, lead = 0, d = 1;
    while(d == 1){
        if(lead == power){
            hare = tortoise;
            power *= 2;
            lead = 0;
        }
        tortoise = f(tortoise);
        lead++;
        d = gcd(N, tortoise - hare);
    }
    return d;
}

নোট:

1. উপরের কোডে \(f(x)=x^2+x+1\) এবং ফাংশানের ভেতরেই \(N\) দিয়ে মডুলাসের কাজ শেষ করা হয়েছে।
2. প্রায় সকল আধুনিক কম্পাইলারে long long int-এর ব্যাবধী \([-2^{63}, 2^{63})\). তাই সি প্রোগ্রামে \(f(x)\) ফাংশনটি ঠিকঠাকমত কাজ করার জন্য \[x^2+x+1\leq 2^{63}-1\Longrightarrow x<2^{31.5}-1\] হতে হবে। এজন্য \(x \ge 2^{31}\) -এর জন্য \(y = x\) নামের একটি নতুন ভেরিয়েবল নিয়ে এর প্রতিটি বিটের সাথে আলাদাভাবে \(x\)-এর গুণফল বের করে সেই সমষ্টির মডুলাস নেওয়া হয়েছে।
3. \(\pmod N\) এবং \(\pmod p\)-তে \(N\)-এর সাইকেল একই হলে উপরের প্রোগ্রাম \(N\) কৃত্রিম হওয়া সত্বেও এর উৎপাদক খুঁজে পাবে না। এরকম হলে \(f(x)\) পলিনোমিয়ালটি পরিবর্তন করে আবার চেষ্টা করতে পার।

কমপ্লেক্সিটি

মনে কর, \(\{a_i\}_{i=0}^{\infty}\) একটি সত্যিকারে র‍্যান্ডম ধারা। জন্মদিনের সমস্যা অনুযায়ী, \(\pmod p\)-তে কনগ্রুয়েন্ট দুটি পদ পাওয়ার সম্ভাবনা \(99\%\) হতে হলে \(\sqrt{-2\times\ln(1-0.99)}\sqrt{p}\)-টি পদ নিতে হবে। অতএব, রো এলগোরিদমের কমপ্লেক্সিটি \(O(\sqrt p)\leq O\left(N^{\frac 1 4}\right)\). [অর্থাৎ \(N\) সংখ্যাটিতে \(k\)টি ডিজিট থাকলে রো এলগোরিদমের রানটাইম \(O\left(10^{\frac k 4}\right)\)]
গতপর্বের চেয়ে আরেকটু ভালো বাউন্ড, তাই না? পরবর্তী পর্বে আমরা জন পোলার্ডেরই আরেকটি এলগোরিদম দেখব, যেটির নাম \(p-1\) এলগোরিদম।

পরিশিষ্ট

জন্মদিন সমস্যা

তোমার ফেসবুক ফ্রেন্ডলিস্টে নূন্যতম কতজন থাকলে একথা নিশ্চিত করে বলা যাবে যে অন্তত দুইজন ফ্রেন্ডের জন্মদিন একই তারিখে হওয়ার সম্ভাবনা \(99\%\)?

মনে কর, তোমার ফ্রেন্ডলিস্টে \(m\) জন ফ্রেন্ড আছে, এবং কোন ফ্রেন্ডের জন্মদিন ২৯ ফেব্রুয়ারি নয়। কোন একজন ফ্রেন্ডের জন্মদিন ৬ অক্টোবর হলে আর কোন ফ্রেন্ডের জন্মদিন ৬ অক্টোবর না হওয়ার সম্ভাবনা \(1-\frac {1}{365}\). আরেক ফ্রেন্ডের জন্মদিন মনে কর ৪ আগস্ট। সুতরাং, তৃতীয় কোন ফ্রেন্ডের জন্মদিন ৬ অক্টোবর বা ৪ আগস্ট না হওয়ার সম্ভাবনা \(1-\frac {2}{365}\). অতএব, \(m\) জন ফ্রেন্ডের ভিন্ন ভিন্ন দিনে জন্মদিন পড়ার সম্ভাবনা \[\left(1-\frac {1}{365}\right)\cdots\left(1-\frac {m-1}{365}\right)\]
তাহলে, অন্তত দুজন ফ্রেন্ডের জন্মদিন একই দিনে পড়ার সম্ভাবনা \[1 – \left(1-\frac {1}{365}\right)\cdots\left(1-\frac {m-1}{365}\right)\]
\(x\)-এর মান \(0\)-এর কাছাকাছি হলে আমরা টেইলরের সিরিজ থেকে \(1-x\approx e^{-x}\) লিখতে পারি। এভাবে উপরের সম্ভাবনাটিকে সরল করে লিখা যায় \[1-e^{-\frac {1}{365}}\cdot e^{-\frac {2}{365}}\cdots e^{-\frac {m-1}{365}}\]\[\approx 1-e^{-\frac{m^2}{2\times 365}}\]এখন \(1-e^{-\frac{m^2}{2\times 365}}\ge 0.99\) হতে হলে \[m\ge \sqrt{-2\times 365\ln(1-0.99)}\approx 58\] হতে হবে।

তোমার ফ্রেন্ডলিস্টে মাত্র \(58\) জন মানুষ থাকলেই তুমি \(99\%\) গ্যারান্টি দিয়ে বলতে পারবে যে অন্তত দুইজনের জন্মদিন একই তারিখে! এই বিখ্যাত সমস্যাটির নাম জন্মদিন সমস্যা।

লক্ষ্য কর জন্মদিনের সমস্যার সমাধান থেকে এটা বলা যায় যে, \(0\) থেকে \(N-1\) এর মাঝে \(m\) সংখ্যা নিলে কোন একটি সংখ্যা অন্তত দুইবার নেওয়ার সম্ভাবনা \(1-e^{-\frac{m^2}{2N}}\).