Category Archives: Maths

Prime Counting Function and Chebyshev Bounds

English, Maths, , ,

The distribution of primes plays a central role in number theory. The famous mathematician Gauss had conjectured that the number of primes between \(1\) and \(n\) is roughly \(n/\log n\). This estimation gets more and more accurate as \(n\to \infty\). We use \(\pi(n)\) to denote the number of primes between \(1\) and \(n\). So, mathematically, Gauss’s conjecture is equivalent to the claim

\[\lim_{n\to\infty}\frac{\pi(n)}{n/\log n}=1\]

This conjecture (currently known as the Prime Number Theorem) is notoriously difficult to prove. The first proof appeared almost a century after it was stated.

However, there are much simpler methods to bound \(\pi(n)\) above and below. Today I shall show such a bound that gives the correct order of magnitude for \(\pi(n)\) up to some constants. In particular, I shall prove that

\[(1+2\log 4)\frac{n}{\log n}\ge \pi(n)\ge \frac{\log 2}{2}\frac{n}{\log n}\]

for \(n\ge 12\). This proof technique is originally due to nineteenth century mathematician Pafnuty Chebyshev. He used a more complicated version of this technique to show that

\[1.1\frac{n}{\log n}\ge \pi(n)\ge 0.92\frac{n}{\log n}\]

for sufficiently large \(n\)s.

The Lower Bound

Suppose \(n\) is even, ie. \(n=2k\). Let us consider the prime factorization of \(\binom{2k}{k}\). All of its prime factors are between \(1\) and \(2k\). Therefore,

\[\binom{2k}{k}=p_1^{a_1}\cdots p_{\pi(2k)}^{a_{\pi(2k)}}\]

Claim: For every prime \(p_i\), \(p_i^{a_i}\leq 2k\)

Proof: We shall use this very simple identity \[\lfloor a+b\rfloor-\lfloor a\rfloor-\lfloor b\rfloor\leq1\]

It is well-known that the largest exponent of a prime dividing \(r!\) is

\[\left\lfloor\frac{r}{p_i}\right\rfloor+\left\lfloor\frac{r}{p_i^2}\right\rfloor+\cdots\]

Since \(\binom{2k}{k}=\frac{(2k)!}{(k!)^2}\), we have the following relationship

\begin{align}
a_i=&\left\lfloor\frac{2k}{p_i}\right\rfloor-2\left\lfloor\frac{k}{p_i}\right\rfloor+\left\lfloor\frac{2k}{p_i^2}\right\rfloor-2\left\lfloor\frac{k}{p_i^2}\right\rfloor+\cdots \\
&\leq 1+1+\cdots
\end{align}

Now there are as many \(1\)s as there are powers of \(p_i\leq 2k\). Hence, \(p_i^{a_i}\leq 2k\).

(Q.E.D)

So, we can conclude,

\[\binom{2k}{k}=p_1^{a_1}\cdots p_{\pi(2k)}^{a_{\pi(2k)}}\leq (2k)^{\pi(2k)}\]

It is easy to prove \(\binom{2k}{k}\ge 2^k\), hence

\[(2k)^{\pi(2k)}\ge 2^k\implies \pi(2k)\ge \frac{\log 2}{2}\cdot\frac{2k}{\log 2k}\]

The same bound can be derived for \(n=2k+1\) using \(\binom{2k+1}{k+1}\).  Therefore, the lower bound is

\[\boxed{\pi(n)\ge \frac{\log 2}{2}\cdot\frac{n}{\log n}}\]

The Upper Bound

Claim: For all \(n\ge 2\), \[\prod_{p\leq n\\ p\ \text{prime}}p\leq 4^n\]

Proof: We induct on \(n\). The base case trivially holds. \(n=2k-1\to 2k\) also holds because the product on the left hand side does not increment. So, we only need to prove \(n=2k\to 2k+1\)
\begin{align}
\prod_{p \leq 2k+1\\ p\ \text{prime}}p &= \prod_{p \leq k+1\\ p\ \text{prime}}p \cdot\prod_{k+2\leq p \leq 2k+1\\ p\ \text{prime}}p\\
&\leq 4^{k+1}\cdot \binom{2k+1}{k}\quad (1)
\end{align}
The first part of the product is \(\leq 4^{k+1}\) because of the inductive hypothesis. The second part is just the product of all the primes in between \(k+2\) and \(2k+1\). It’s easy to see that they all divide \(\binom{2k+1}{k}\), hence their product is \(\leq \binom{2k+1}{k}\).

It’s easy to show \(\binom{2k+1}{k}<4^k\). Substituting this bound back to (1) gives us

\[\prod_{p \leq 2k+1\\ p\ \text{prime}}p \leq 4^{k+1}\cdot 4^k=4^{2k+1}\quad (\text{Q.E.D.})\]

Now take a number \(2\leq m\leq n\). Consider the primes that are between \(m\) and \(n\). There are \(\pi(n)-\pi(m)\) of them, each of which is \(\ge m\). Consequently,

\[m^{\pi(n)-\pi(m)}\leq \prod_{m\leq p \leq n\\ p\ \text{prime}}p\leq \prod_{1\leq p \leq n\\ p\ \text{prime}}p\leq 4^n\]

Now take log of both sides

\[(\pi(n)-\pi(m))\log m\leq n\log 4\]

Rewrite this as

\[\pi(n)\leq \pi(m)+\frac{n\log 4}{\log m}\leq m+\frac{n\log 4}{\log m}\]

Set \(m=\frac{n}{\log n}\) and simplify the expression

\[\pi(n)\leq \frac{n}{\log n}\left(1+\frac{\log 4\cdot\log n}{\log(n/\log n)}\right)\]

Easy to show \(\frac{\log n}{\log(n/\log n)}\leq 2\) for \(n\ge 12\). (In fact you can show much tighter bound for sufficiently large \(n\).) Hence, the upper bound is

\[\boxed{\pi(n)\leq \left(1+2\log 4\right)\frac{n}{\log n}}\]

Multiplying Two Polynomials with Fast Fourier Transform

Computer Science, English, Maths,

Polynomial multiplication is one of the most important problems in mathematics and computer science. It can be formally defined as follows:

You are given two polynomials roughly of equal size\[A(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1}\]\[B(x)=b_0+b_1x+\dots+b_{n-1}x^{n-1}\]find a polynomial \[C(x)=c_0+c_1x+\dots+c_{2n-2}x^{2n-2}\] such that \(A(x)\cdot B(x)=C(x)\). At first, this may look like an easy problem, since you already know how to multiply polynomials. You just sum the product of every two terms, right?

def multiply(A, B):
	"""
	A = list of coefficients of A(x)
	B = list of coefficients of B(x)
	""" 
	n = len(A) - 1
	# C can have 2*n-2 degree at most
	C = [0 for i in range(2*n-1)]
	for i in range(n):
		for j in range(n):
	    	# adding a_i*b_j to c_{i+j}
	    	C[i+j] += A[i]*B[j]

It definitely works, but has asymptotic runtime of \(O(n^2)\). That’s quite bad. Can we do it faster?

Turns out we can! We can in fact multiply two polynomials in \(O(n\log n)\) time. How is that possible? It all comes down to looking at the problem from another perspective.

Change of Perspective

So far we only saw polynomials as a list of coefficients. However, a polynomial of degree \(n\) can also be viewed as a curve in \(n+1\) dimensional space. It is quite possible define it in terms of the points it passes through. For instance, take the line \(P(x)=2x+1\) in 2d space. We only need to know two points on this line, say \((0, 1)\) and \((1, 3)\), to identify it uniquely. To determine a polynomial with degree \(n\), It has been proven that you need to know its values at \(n+1\) distinct points. This is called the Point Value Representation of Polynomials.

A polynomial \(P(x)\) in this form can be saved like \(\{0 : P(0), \dots, n : P(n)\}\), although you can evaluate \(P(x)\) at any complex point. (This will come handy later)

Note that multiplication is really cheap in this form. Given
\[A(x)=\{0:A(0),\dots (2n-2): A(2n-2)\}\] and \[B(x)=\{0:B(0),\dots (2n-2): B(2n-2)\}\] we can compute \(C(i)=A(i)\cdot B(i)\) at every point \(i\) from \(0\) to \(2n-2\). For instance, \(C(0)=A(0)\cdot B(0)\).This costs \(O(n)\) time. So, we are happy.

Wait, are we? We do not know the coefficients of \(C(x)\)! Evaluation at arbitrary points is very costly in point value form. So, we could not, for instance, raise \(C(x)\) to some arbitrary power. What if we could find a way to go back and forth between coefficient representation and point value representation? We could get the best of both worlds.

This is what Fast Fourier Transform allows us to do.

Before I move on, please note a couple of things.

  • Prior to multiplication, we must evaluate both \(A(x)\) and \(B(x)\) at at least \(2n-1\) points even though each of them have degree \(n\). This is because \(C(x)\) has degree \(2n-2\).
  • It is not necessary to evaluate the polynomials exactly at \(1,2,\dots,2n-2\). There is nothing special about these numbers. However, it is necessary that we evaluate both \(A(x)\) and \(B(x)\) at the same set of points.

Coefficients to Point Value Representation

The only way to convert a polynomial with degree \(n\) to point value representation is to evaluate it at \(n+1\) distinct points. However, polynomial evaluation at a single point stills costs \(O(n)\) in coefficient representation. For instance, you can evaluate \(A(p)\) at any point \(p\) like this: \[A(p)=a_0+p(a_1+p(a_2+\cdots+pa_{n-1}))\] Then for \(n+1\) points, we would end up with \(O(n^2)\) runtime. However, you can save a lot of this cost by algebraic manipulation. Let me give you an example.

Suppose we want to compute \(P(x)=2+3x-4x^2+2x^3+5x^4\) at \(1, 2, -1, -2\). We can evaluate the odd powers and even powers separately since even powers of negative numbers is same as that of positive numbers. To be more precise, define \(Q(x)=2-4x+5x^2\) and \(R(x)=3x+2x^2\).

Note that we only kept the order of the coefficients from \(P(x)\), but not the associated powers \(x\). This gives an algebraic relationship as follows:
\[P(x)=Q(x^2)+xR(x^2)\] Now, note that \(Q((-1)^2)=Q(1^2), Q((-2)^2)=Q(2^2)\) and so on. Therefore, we will only need to compute \(Q(1), Q(4), R(1), R(4)\) in order to evaluate \(P(x)\) at four points. Both \(Q(x)\) and \(R(x)\) has half the degree as \(P(x)\), hence we’re saving almost half the work.

In order for trick to work, we must select our list of points very carefully. If the polynomial were much bigger, we would have to do that odd-even split multiple times. Every time we would like to have numbers with opposite signs, so that when squared, they give the same value. You can easily see that no initial set of real numbers can satisfy such requirements. This is why we choose roots of unity. In particular, we choose \(m=2^k\)th roots of unity such that \(2n\leq m\leq 4n\). (This is essentially the smallest power of two that’s greater than \(2n\)). They have this magic property that if \(w\) is a \(2^k\)th root, then so is \(-w\) regardless of which \(k\) you choose.

Point Value to Coefficient Representation

If you know some linear algebra, you might notice that evaluating a polynomial at \(n+1\) point essentially defines a linear transformation, which can be represented by a matrix. Due to special nature of roots of unity, the inverse matrix consists of conjugates of previously used roots scaled by a factor of \(1/n\). I won’t show it here in details, but if you are interested you should check out this blogpost.

In English, the equations tell us that inverse of Fast Fourier Transform is another Fast Fourier Transform where every point of polynomial evaluation must be replaced by its conjugate, and the final result has to be scaled by a factor of \(1/n\).

Time Complexity

Fast Fourier Transform is essentially a divide-and-conquer algorithm. The divide step corresponds to odd-even coefficient split. Each of the resulting subproblems has cost \(T(n/2)\). The combining step has cost \(O(n)\). (There are \(2n-2\) points and you have to compute the algebraic identity at each point.) So, the recurrence relationship should be
\[T(n)=2T\left(\frac{n}{2}\right)+O(n)\] Using master theorem, \(T(n)=O(n\log(n))\). Since FFT is the most expensive step in this approach, runtime of the entire process is dominated by \(O(n\log n)\). Therefore, polynomial multiplication also costs \(O(n\log n)\).

Summary

  • Given \(A(x)\) and \(B(x\) in coefficient representation, convert them to point value representation with Fast Fourier Transform. (\(O(n\log n)\)))
  • Compute \(C(x)=A(x)\cdot B(x)\). (\(O(n)\))
  • Go back to coefficients representation with Inverse Fast Fourier Transform. (\(O(n\log n)\)))

Generating All Subsets of Size k in Python

Computer Science, English, Maths, , ,

Suppose you are given a set \(S\) with \(n\) elements and you need to generate every subset of size \(k\) from it. For instance, if \(S=\{3,4,5\}\) and \(k=2\), then the answer would be \(\{3,4\}\), \(\{3,5\}\), \(\{4,5\}\). So, how would you do that in Python?

First of all, this is a really simple exercise. Stuff like this often comes up when someone writes a moderately large script. I wanna talk about this one in particular because it has a combinatorial solution.

Before I show it to you, I highly recommend you think about it on your own.

Solution

Notice that if you convert the set to a list, every element is assigned to a unique index between \(0\) and \(n-1\). So, it suffices to generate all the subsets of size \(k\) from \(\{0,1,\ldots,n-1\}\).

Now, every such subset falls into one of two categories: either it contains \(n-1\) or it doesn’t.

If \(n-1\) belongs to the subset, the rest of its elements form a subset of size \((k-1)\) in \(\{0,1,\ldots,n-2\}\).

If \(n-1\) doesn’t belong to the subset, then that subset is also a subset of \(\{0,1,\ldots,n-2\}\).

This is essentially a recursion. So, we are done.

Pretty neat, eh?

def gen_subset(n, k):
    if n >= 0 and n >= k >= 0:
        if n == 0 or k == 0: 
            yield set()
        elif n == k:
            # returns {0, 1, ..., n-1}
            yield set(range(n))
        else:
            # ksubsets without n-1
            yield from gen_subset(n-1, k)
            # ksubsets with n-1
            for subset in gen_subset(n-1, k-1):
                subset.add(n-1)
                yield subset
    else: 
        raise ValueError("Illegal Parameters")

A Special Case of Zsigmondy’s Theorem

English, Maths, , , , ,

Theorem (Zsigsmondy): 

  1. For two coprime positive integers \(a\) and \(b\) and for any positive integer \(n\), \(a^n-b^n\) has a prime divisor that does not divide \(a^k-b^k\) for all positive integers \(k< n\) with the following exceptions:
    • \(n = 2\) and \(a+b\) is a power of two.
    • \(n=6, a=2, b=1\).
  2. For two coprime positive integers \(a\) and \(b\) and for any positive integer \(n\), \(a^n+b^n\) has a prime divisor that does not divide \(a^k+b^k\) for all positive integers \(k< n\) with the following exception:
    • \(n=3, a=2, b=1\).

This theorem is very helpful in solving many olympiad number theory problems. However, its use is often frowned upon, as this theorem is quite hard to prove using only elementary mathematics. (I am aware of the proof using Cyclotomic Polynomials, but that’s not “elementary enough” in my opinion.)

In olympiad mathematics, one can get away in most of the cases by just using this theorem for a fixed pair of \(n\) and \(k\). This modification allows us to provide a very simple proof using LTE. Here I’ll restate this particular case, and then prove the first part. The second part can be proven analogously.

Theorem (Zsigsmondy Special Case): 

  1. For two coprime positive integers \(a\) and \(b\) and any two positive integers \(n\) and \(k\) with \(k<n\), \(a^n-b^n\) has a prime divisor that does not divide \(a^k-b^k\) with the following exception:
    • \(n = 2, k=1\) and \(a+b\) is a power of two.
  2. For two coprime positive integers \(a\) and \(b\) and any two positive integers \(n\) and \(k\) with \(k<n\), \(a^n+b^n\) has a prime divisor that does not divide \(a^k+b^k\) with the following exception:
    • \(n=3, k=1,a=2, b=1\).

Proof of 1:

Suppose \(a^n-b^n\) and \(a^k-b^k\) share same set of prime divisors. This implies \(a^n-b^n\) and \(a^{\gcd(n,k)}-b^{\gcd(n,k)}\) also share same set of prime divisors. Assuming \(A=a^{\gcd(n,k)}\) and \(B=b^{\gcd(n,k)}\), we could follow the rest of this argument and get to a contradiction. So, without loss of generality, we may assume \(k=1\).

Now we consider two cases:

Case 1: \(n\) is a power of \(2\).
If \(n=2\) and \(a+b\) is a power of two, we arrive at one of the listed exceptions. Now note that,
\(\begin{align*}&\gcd(a+b, a^2+b^2)\\
=&\gcd(a+b, (a+b)^2-a^2-b^2)\\
=&\gcd(a+b,2ab)\\
=&\gcd(a+b,2)\end{align*}\)
This implies both \(a+b\) and \(a^2+b^2\) can’t be powers of two unless \(a+b=2\implies a=b=1\). Hence at least one of them must have an odd prime divisor. Furthermore, this prime divisor does not divide \(a-b\) since \(\gcd(a-b, a^2+b^2)=\gcd(a-b,2)\) (following the previous steps). However, if \(n>2\), then \(a^2+b^2|a^n-b^n\), implying that odd prime will divide \(a^n-b^n\). This is a contradiction.

Case 2: \(n=2^md\) with \(d>1\) being odd.
Without loss of generality, we may assume \(n>1\) is odd, since \(a^d-b^d|a^n-b^n\) and it is sufficient to show that \(a^d-b^d\) has a prime divisor that does not divide \(a-b\).  From LTE, \(v_p(a-b)+v_p(n)=v_p(a^n-b^n)\) for each odd prime \(p|a-b\). Furthermore, \[\frac{a^n-b^n}{a-b}\equiv na^{n-1}\equiv 1\pmod 2\] implying \(v_2(a-b)=v_2(a^n-b^n)\). Therefore, we can conclude \[\frac{a^n-b^n}{a-b}=\prod_{p|\gcd(n,a-b)}p^{v_p(n)}\leq n\] which is impossible for \(n>1\). This raises another contradiction and we are done.

ইন্টিজার ফ্যাক্টোরাইজেশন: পর্ব ১ পোলার্ডের রো (rho) এলগোরিদম

Bangla, Computer Science, Maths, , , , , , , , ,

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

গত পর্বে আমরা গায়ের জোরে কিছু সংখ্যার উৎপাদক বের করার চেষ্টা করেছিলাম। দুঃখজনকভাবে আমাদের বাচ্চা প্রোগ্রামটি \(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}\) ধারাটি হয়

\(i\) \(0\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\) \(7\) \(8\) \(9\) \(10\) \(11\)
\(a_i\) \(1\) \(2\) \(5\) \(26\) \(10\) \(14\) \(23\) \(8\) \(7\) \(21\) \(\color{red}7\) \(\color{red}{21}\)

দেখতেই পাচ্ছ, \(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}}\).

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

Bangla, Computer Science, Maths

ইন্টিজার ফ্যাক্টোরাইজেশন: পর্ব ১ পোলার্ডের রো (rho) এলগোরিদম

আমরা কেন ইন্টিজার ফ্যাক্টোরাইজেশন করি

ইন্টেজার ফ্যাক্টোরাইজেশন অর্থ হচ্ছে একটা পূর্ণসংখ্যাকে উৎপাদকে বিশ্লেষণ করা। যেমন \(36 = 4\times 9\). এখানে \(4\) আর \(36\) উভয়েই \(36\)-এর উৎপাদক।

ছোটবেলায় আমাদের সবাই ইন্টেজার ফ্যাক্টোরাইজেশন করতাম। কেন জানো? পরীক্ষায় নম্বর পাওয়ার জন্য! প্রশ্নে একটা সংখ্যা দেওয়া থাকত আর আমরা দুই থেকে শুরু করে একটা একটা সংখ্যা দিয়ে প্রশ্নে দেওয়া সংখ্যাটাকে ভাগ দিতাম। যখন ভাগ যেত, তখন বোঝা যেত যে একটা উৎপাদক (বা গুণনীয়ক) পেয়ে গেছি। খুবই রুটিনমাফিক একটা কাজ, তাই না?

তবে বাচ্চাদের পরীক্ষার প্রশ্নে আসা ছাড়াও ইন্টিজার ফ্যাক্টোরাইজেশনের অন্যরকম একটা গুরুত্ব আছে। যেমন- অনলাইনে তোমার, আমার, সকলের ব্যক্তিগত তথ্য পাজি লোকদের হাত থেকে সুরক্ষিত রাখতে গোপন সংকেতের আড়ালে রাখা হয়। এরকম অনেক গোপন সংকেত বানাতে ও ভাঙতে দুইশ, তিনশ ডিজিটের বিশাল বিশাল সংখ্যাকে উৎপাদকে বিশ্লেষণ করতে হয়, কিংবা প্রাইম খুঁজে বের করতে হয়। এসব হিসেব এত বড় হয় যে খাতা কলমে করতে গেলে আমাদের চুল পেকে যাবে, পাতার পর পাতা ফুরিয়ে যাবে। সেজন্য, আমরা তখন কম্পিউটারের সাহায্য নেই।

কিন্তু সমস্যা কি জানো? কম্পিউটার খুব বোকা। যদি আমি বলি, “কম্পিউটার, আমার একটা হিসেবে \(731\)-এর উৎপাদকে বিশ্লেষণ লাগবে। তুমি চট করে সেটা বের করে দাও তো দেখি!”, সে হা করে বসে থাকবে। কারণ সে তো আর জানে না কীভাবে উৎপাদকে বিশ্লেষণ করতে হয়!

কম্পিউটারকে দিয়ে উৎপাদকে বিশ্লেষণ করানোর জন্য তাকে তোমার প্রাইমারির বাচ্চাদের মত একটু একটু করে বুঝিয়ে বলতে হবে, “কম্পিউটার, \(2\) থেকে \(730\) পর্যন্ত প্রতিটি সংখ্যা দিয়ে \(731\)-কে ভাগ দাও। যদি কোনটা দিয়ে নিঃশেষে ভাগ যায়, তবে সেই সংখ্যা এবং ভাগফল আমাকে দেখাও।” 

def factor(N):
    for d in range(2, abs(N) + 1):
        if N % d == 0: return d

N = int(raw_input())
d = 1 if(2 > N > -2) else factor(N)
print d, '*', N / d

কম্পিউটার যেহেতু বাংলা পারে না, তাই সে বোঝে এমন একটি ভাষা পাইথনে আমরা আমাদের কথাগুলো লিখেছি। এই কোডে একটি পূর্ণসংখ্যা \(N\)-কে ইনপুট হিসেবে দেওয়া হবে। আর কম্পিউটার সাথে সাথে \(N\)-কে উৎপাদকে বিশ্লেষণ করে খুশিমনে আমাদের জানিয়ে দেবে। যত বড় সংখ্যাই হোক না কেন, আমাদের ছেলেবেলায় শেখা নিয়ম দিয়ে তার উৎপাদক বের করা সম্ভব।

উপরের প্রোগ্রামে \(731\) ইনপুট দিলে কম্পিউটার আমাদের বলে দেবে যে \(731 = 17\times 43\).

এবার একটা কাজ কর। উপরের প্রোগ্রামটি চালিয়ে \(N = 18446744073709551619\) ইনপুট দাও। কী হল?

প্রোগ্রামটি আর থামছে না।

এর কারণ হচ্ছে উপরের প্রোগ্রামে কম্পিউটারকে \(2\) থেকে \(N-1\) পর্যন্ত সকল সংখ্যা দিয়েই \(N\)-কে ভাগ (প্রকৃতপক্ষে মডুলাস নির্ণয়) করতে হচ্ছে। কম্পিউটারের প্রতিটা ভাগ করতে কিছুটা সময় লাগে। তাই \(N\) যত বড় হবে, কম্পিউটারকে তত বেশি সংখ্যা দিয়ে ভাগ করতে হবে, এবং প্রোগ্রামের রান করতে সময় (রানটাইম) বেশি লাগবে। রানটাইম \(N\)-এর সাথে পাল্লা দিয়ে কত দ্রুত বাড়বে সেটার একটা ধারণাও আমরা পেতে পারি।

ধর \(N\) সংখ্যাটিতে \(k\) সংখ্যক ডিজিট আছে। আমরা আসন্নভাবে \(N\approx 10^k\) ধরতে পারি। যদি \(N\) একটি মৌলিক সংখ্যা হয়, তবে প্রোগ্রাম ১-এ \(N\)-কে \(2\) থেকে \(N-1\) পর্যন্ত সকল সংখ্যা দিয়েই ভাগ করা হবে। যদি প্রতিবার ভাগ করতে সর্বোচ্চ \(c\) পরিমাণ সময় লাগে, তবে উপরের প্রোগ্রামটি রান করতে সর্বোচ্চ \(Nc=10^{k}c\) সময় লাগবে। এর অর্থ হচ্ছে, ডিজিটের সংখ্যা বাড়ার সাথে সাথে উপরের প্রোগ্রামের রানটাইমও সূচকীয়ভাবে বৃদ্ধি পায়। [বা, বিগ O নোটেশনে উপরের প্রোগ্রামের রানটাইম \(O(10^{k})\)]

প্রোগ্রামিং-এ রানটাইম সূচকীয়ভাবে বাড়া খুবই দুঃখের একটি কথা। কারণ ইনপুটের একটুখানি পরিবর্তনে প্রোগ্রামের রানটাইমে রাতদিন পরিবর্তন হয়ে যেতে পারে।

মনে কর তোমার কাছে যদি কোন মহাক্ষমতাধর সুপার কম্পিউটার আছে যে \(200\) ডিজিটের একটা সংখ্যার উৎপাদকে বিশ্লেষণ প্রোগ্রাম ১ দিয়ে \(1\) সেকেন্ডে বের করতে পারে। তাহলে \(210\) ডিজিটের একটা সংখ্যার উৎপাদকে বিশ্লেষণ করতে কম্পিউটারটির সময় লাগবে \(10^{210-200}=10^{10}\) সেকেন্ড \(\approx 317\) বছর!

আমদের জীবনটা খুব ছোট, জানো? একটা প্রোগ্রামের আউটপুট দেখে যাওয়ার জন্য আমাদের \(317\) বছর বেঁচে থাকা সম্ভব না। তাই খুব ফাস্ট কোন এলগোরিদম দরকার।

প্রোগ্রাম ১-কেই অনেকখানি ফাস্ট বানানো যায় যদি আমরা \(N-1\) এর পরিবর্তে \(\sqrt{N}\) পর্যন্ত প্রতিটা সংখ্যা দিয়ে ভাগ দেই। কেননা প্রতিটি কৃত্রিম সংখ্যার \(\sqrt N\)-এর মাঝেই অন্তত একটা উৎপাদক থাকে। আর যদি উৎপাদক না পাওয়া যায়, তবে সংখ্যাটি প্রাইম। এই পরিবর্তনের ফলে আমাদের এলগোরিদমের রানটাইম হবে \(10^{\frac k 2}c\) [বা \(O\left(10^{\frac k 2}\right)\)]. কিন্তু এটাও যে সূচকীয় হারে বাড়ে!

এর থেকে ফাস্ট কি কোন এলগোরিদম নেই? উত্তরটা হচ্ছে, হ্যাঁ আছে। তবে সেটি জানতে হলে তোমার মডুলার এরিথমেটিক সম্পর্কে হালকা পাতলা ধারণা লাগবে।

যদি \(a\) এবং \(b\) কে \(n\) দিয়ে ভাগ করলে একই ভাগশেষ থাকে, তাহলে আমরা লিখি \(a \equiv b\pmod{n}\). এটাকে পড়া হয় \(a\) is congruent to \(b\) modulo \(n\). যেমন: \(57\equiv 112\pmod{55}\).

লক্ষ্য কর যে, যদি \(n\) দিয়ে \(a\)-কে ভাগ করলে \(r\) ভাগশেষ থাকে, তবে \(a\equiv r\pmod{n}\) লিখা যায়। যেমন: \(57\equiv 2\pmod{55}\).

যেকোন পূর্ণসংখ্যাকে যদি \(55\) দিয়ে ভাগ কর, তবে তুমি \(0\) থেকে \(54\) এর মাঝে কোন একটি ভাগশেষ পাবে। তাহলে, তুমি যদি অনেকগুলো র‍্যান্ডম পূর্ণসংখ্যার নাও, এদের মাঝে \(\pmod {55}\)-এ কনগ্রুয়েন্ট, অর্থাৎ \(55\) দিয়ে ভাগ করলে একই ভাগশেষ থাকে, এমন দুটি পূর্ণসংখ্যা পাওয়ার কিছুটা সম্ভাবনা থাকে, তাই না? যত বেশি সংখ্যা নেবে, এই সম্ভাবনাও তত বাড়বে।

আবার, একই সংখ্যাগুলির মাঝে \(\pmod {5}\)-এ কনগ্রুয়েন্ট দুটি সংখ্যা পাওয়ার সম্ভাবনাও থাকে। কিন্তু এই সম্ভাবনা \(\pmod{55}\)-এ কনগ্রুয়েন্ট হওয়ার সম্ভাবনার চেয়ে অনেক বেশি। কেননা, \(\pmod{5}\)-এ তো খালি \(0\) থেকে \(4\) পর্যন্ত ভাগশেষ আসতে পারে।

এই ধারণাকে ব্যবহার করে উৎপাদকে বিশ্লেষণের একটি চমৎকার এলগোরিদম রয়েছে, যার নাম পোলার্ডের রো(\(\rho\)) এলগোরিদম। পরবর্তী পর্বে এটি নিয়ে আলোচনা করব।

Guidelines to Open a Successful Math Club

Maths, Uncategorized,

[I run Mymensingh Parallel Math School (MPMS), a free-for-all math club whose members (including me) have consistently won many prizes  in Regional, National and International Mathematical Olympiads. Naturally, many math enthusiasts often ask me for tips on starting math clubs. So, 4 years ago, I wrote the following short note. I have decided to republish it here for future references.]

Unlike seeking for philosopher’s stone, starting a math club is very easy. You see, it requires only two things: a few classrooms and some members. Even a ‘teacher’ is not needed because the classes can find their way all by themselves! Although first a few month they need an instructor to teach them how to walk all-alone.

Through this tender age you have to be a bit careful. Try to drive the class in a relaxed mood so that everybody feels free to join in the discussion. It’s a very, very important thing for encouraging everyone to share his comments and ideas with others. In fact continuing this practice one day enables them to take the duster and the chalks in own hands. Moreover, we advise you to look for the true dedicated persons amongst the members. They will run it and never leave it. In this process MPMS has got its perfect structure.

May be the mostly asked questions are all about selecting class-topic. And yeah, this is too difficult to answer. Actually you can choose anything of your choice. There is neither any limitation nor any boundary. Math may appear boring sometimes. So why you don’t look at physics or literature then? In MPMS, physics, chemistry, economics, literature, biology-even philosophy take place frequently in the chaos. Also it is wise to discuss attractive topics like combinatorics at the early stage. Later you may go through harder things day by day depending on everybody’s capability to cope with it.

So…turn off your computer and come out to the field. See you at the BdMO National… 🙂

An interesting application of AM-GM to prove the upper bound for the harmonic numbers

Maths, Uncategorized, , ,

Harmonic series is a very well-known divergent series defined by
\[\sum_{k=1}^\infty \frac 1 k =  1+\frac 1 2 +\frac 1 3+\frac 1 4\cdots\]

The partial sum of first \(n\) terms of this series is called the \(n\)th harmonic number. In other words, \[H_n = \sum_{k=1}^n\frac 1 k\]

Although there is no nice and simple formula for computing \(H_n\), we can, however, do some approximation. For example, Euler proved \[H_n = \ln(n)+\gamma+\varepsilon_n < \ln(n)+1\]

where \(\gamma\) is the Euler–Mascheroni constant and \(\varepsilon_n\approx \frac 1 {2n}\). I shall prove this bound here with AM-GM Inequality, which says, for positive real \(a_i\)s, that
\[\frac 1 n \left(\sum_{i=1}^n a_i\right)\ge \sqrt[n]{\prod_{i=1}^n a_i}\]

So we can write \[ 1+\sum_{k=2}^n \frac{k-1}{k}> n\sqrt[n]{\prod_{k=2}^n \frac{k-1}{k}}=n\sqrt[n]{\frac{1}{n}}\]
The inequality is strict since the \(\frac {k-1}k\)s are not equal. Also note that

\[\sum_{k=2}^n \frac{k-1}{k}=\sum_{k=1}^n 1-\frac{1}{k}=n-H_n\]

Plugging it back to the main inequality, we derive that

\[1+n-H_n>n\sqrt[n]{\frac{1}{n}}\]

And by rewriting

\[1+n\left(1-\sqrt[n]{\frac{1}{n}}\right)>H_n\]

Now, all we need to prove is

\[\ln(n)\ge n\left(1-\sqrt[n]{\frac{1}{n}}\right)\]

Assume \(e^m = n\) and \(\dfrac {m} {e^m} = x\). The previous inequality can be rewritten as

\[m \ge e^m(1-e^{-x})\]

It can be further simplified to

\[x+e^{-x}\ge 1\]

Now, the derivative of \(f(x)=x+e^{-x}\) is \(1-e^{-x}\), which is non-negative in \([0, \infty )\). And we also have \(f(0)=1\). Therefore, we definitely have \(\displaystyle x+e^{-x}\ge 1\) and we are done.

Remarks:

  1. \[\ln(n)\ge n\left(1-\sqrt[n]{\frac{1}{n}}\right)\]
    A weird inequality. You can even see from the previous proof that
    \[\lim_{n\to\infty}n\left(1-\sqrt[n]{\frac{1}{n}}\right) = \ln (n)\]
    Perhaps, with some luck, it can be used to derive some more interesting results. Gotta think later.
  2. In my proof, the \(1\) beside \(\ln(n)\) just sat there to initiate the AM-GM. If we replace it by some function \(\lambda(n)\) and do the calculation, we might be able to optimize our inequality down to the \(\gamma\) bound. But I’m done for today. Anyone interested might give it a shot.