Category Archives: Computer Science

An Online Algorithm to Check for Bipartite Graphs

Bipartite Graphs

A graph \(G(V, E)\) is called bipartite if its vertices can be divided into two groups \(X\) and \(Y\) such that every edge connects one vertex from \(X\) and one vertex from \(Y\). The graph drawn below is bipartite.

Given a graph, can you determine if it is bipartite?

This problem is pretty straightforward because the vertices of a bipartite graph can be colored with two colors such that every edge connects two vertices of different colors (see above). This property is called 2-colorability. The converse is also true, that is if a graph is 2-colorable, it is also bipartite. Hence, it order to check for bipartiteness, we can run a Breadth First Search on the graph and color the vertices as we discover them. It would take \(O(V+E)\) time and storage.

An Online Variant

First of all, what is an online algorithm? An online algorithm is an algorithm that processes the input in a piece by piece in a sequential manner and have access to limited storage. These algorithms are useful in many realistic scenarios.

Now consider an online variant of the previous problem: suppose your algorithm gets a single pass over the edges of a graph and has access to only \(O(V)\) storage. In other words, the edges are fed into the algorithm one at a time, but the algorithm cannot store all of them. Can you determine if this graph is bipartite?

This is a practical problem because a graph with \(|V|\) vertices can have as many as \(|V|^2\) edges. For instance, if a graph has \(10^{10}\) vertices, there can be \(O(10^{20})\) edges. Storing them in the memory might be impossible. The problem above attempts to bypass this issue by getting the edges sequentially, and never storing all of them. I’ll show that it is indeed possible to design such an algorithm with only \(O(E\log V)\) runtime.

We’ll use the 2-colorability property once more. In the beginning, the graph is 2-colorable since all the vertices are disjoint. Let’s color all the vertices to be red. Every time a new edge is added to the graph, it connects two vertices of either different or same colors. In the former case, the graph remains bipartite. In the later case, one of the following two scenarios may occur:

  1. The edge connects two vertices from the same connecting component. Now, if the edge connects two vertices of same color, the graph is no longer 2-colorable. We can output FAIL and stop.
  2. The edge connects two vertices from different connected components. In this case, if the edge connects two edges of the same color, the graph is still 2-colorable. However, we need to flip colors of all the vertices in one of the connected components (see below).

A naive way to do #2 would be to flip color of every vertex in one of the connected components. This may take \(O(V)\) time for each edge, the total runtime being \(O(EV)\). We need a smarter way.

One potential approach is to do the flipping in a lazy manner: if #2 happens, only remember that colors in a connected components need to be flipped, instead of actually flipping them. When the algorithm checks the color of a vertex,

  • go through your records to find the connected components the vertex was previously in.
  • Observe how many of those connected components had their colors flipped.
  • Since all vertices started as red, the number of flips determine current color of that vertex.

This reminds us of a data structure called Union-Find.

A Quick Intro to Union-Find

This data structure is is essentially a list of mutually disjoint sets. It has three methods:

  • make-set: makes a set with a unique id
  • find(x): given element \(x\), finds the id of the set \(x\) belongs to.
  • union(x, y): performs the set-union on the sets \(x\) and \(y\) belong to.

Union-Find is usually implemented as a forest of trees, with root of every tree being its id. Every node contains a pointer to its parent. Roots contain pointers to themselves. find(x) looks for the root of a tree. union(x, y) first does find(x) and find(y), then attaches the root of y as a child to the root of x.

In a naive implementation, the cost of the find and union functions are \(O(n)\). However, those can be reduced to \(O(\log n)\) via a technique called union by rank. In this technique, a shorter tree is always attached to a taller tree. You can easily show (by induction) that such a tree with height \(n\) must contain atleast \(2^n\) nodes. In other words, the time complexity for this data structure is \(O(\log n)\).

Faster Algorithm Using Union Find

In our problem, every connected component can be thought of as a disjoint set. Every node starts in its own set. When two connected components are merged, we perform set-union over their corresponding sets. Now, notice that, the root of every tree can remember if the connected component needs to flip color during the union. When we call find on an element, we’ll pass through all its ancestors, which were roots in some previous connected components. We can check each of them to see if those components had their colors flipped. Since each find costs \(O(\log V)\) time, computing the color of each node requires \(O(\log V)\) time, too. Now, each of the steps #1 and #2 in the previous algorithm takes \(O(\log V)\) time. Hence, the total runtime for our modified algorithm is \(O(E\log V)\).

Appendix

In Union Find, there is another technique, called Path Compression, which (along with Union by Rank) reduces the data structure’s  amortized complexity to \(O(\alpha(n))\). The previous algorithm works with path compression, but needs to be modified. This is left as an exercise to the reader.

A Blind Robot Beside an Infinite Wall

Let’s think about the following problem:

Consider a wall that stretches to infinitely in both directions. There is a robot at position \(0\) and a door at position \(p\in\mathbb Z\) along the wall \((p\neq 0)\). The robot would like to get to the door, but it knows neither \(p\), nor the direction to the door. Furthermore, the robot cannot sense or see the door unless it stands right next to it. Give a deterministic algorithm that minimizes the number of steps the robot needs to take to get to the door.

This problem quite famous in Bangladesh because of Dr. Mohammad Kaykobad, a renowned Bangladeshi computer scientist. He maintains a list of interesting problems like this one to attract high school kids into problem solving.

Dr. Kaykobad’s formulation of this problem is quite ambiguous. He simply asks how the robot should get to the door, without specifying whether the algorithm should be deterministic. It matters as we shall see.

In this article, I shall provide a deterministic algorithm that takes at most \(9|p|\) steps for any \(p\). I shall also prove that any deterministic algorithm will take at least \(8|p|\) steps. I suspect this bound can be improved up to \(9|p|\), but I didn’t have time to think much.

My \(\leq 9|p|\) Deterministic Algorithm

for i = 0, 1, 2, ...:
    walk up to (-2)^i, but don't walk past p
    if you find the door: break loop
    else: get back to 0

Suppose \(p\) is found on \(i\)-th round. Then, \(|p| > 2^{i-2}\), otherwise \(p\) would have been found on \((i-2)\)-th round. Now, it is easy to see the robot takes \(2^1+\dots+2^{i}<2^{i+1}\) steps prior to \(i\)-th round. It takes \(|p|\) more steps to get to \(p\). Therefore, the total number of steps is at most \[|p|+2^{i+1}=|p|+8\cdot2^{i-2}<9|p|\] As \(|p|\to\infty\), total number of steps converges to \(9|p|\).

(Disclaimer: Finding a \(9|p|\) algorithm was a homework problem in my class 6.046. The algorithm and the runtime analysis shown above was part of the official solution. My algorithm was the same, but my runtime analysis was much uglier.)

Other Deterministic Algorithms

In this section, we shall prove that any deterministic algorithm for this problem will take at least \(8|p|\) steps.

In any algorithm, the robot has to walk some steps to the left/right, then turn back and walk some more steps, then turn again, and so on. Without loss of generality, we may assume the optimum algorithm takes first step to the right.

Consider the the optimum deterministic algorithm \(OPT\) with minimum number of turns \(n\). Assume that \(\{a_i\}_{i=1}^n\) denotes the number of steps \(OPT\) takes during the \(i\)-th turn. Define \(A_i=\sum_{j=1}^i(-1)^{j+1}a_j\) to be the position of the robot after \(i\) steps. Also define \(A_0=0\). Suppose \(OPT\) takes fewer than \(8|p|\) steps for all \(p\).

Claim 1: \(a_n>\cdots a_i>\cdots a_1 > 0\) for \(n> i > 1\).

Proof: If the robot takes \(0\) step at any turn, we can simply ignore that turn and merge the steps of the next turn with the steps of the previous turn. Hence, \(a_1,a_i,a_n > 0\).

Suppose \(a_{i-1}-a_{i} \ge 0\) for some \(i\). This means at \(i\)-th step, the robot backtracked to somewhere between \(A_{i-1}\) and \(A_{i-2}\). Then, at the next step, it must go past \(A_{i-1}\), otherwise, \((i+1)\)-th step will revisit some previously visited vertices and it can be removed by walking to \((a_{i-1}-a_i+a_{i+1})\) at \(i\)-th step. This is a contradiction to the minimality of \(OPT\). Hence, we must have \((a_{i-1}-a_i)+a_{i+1} > a_{i-1}\). However, in this case, we can get a better algorithm by taking \(a_{i-1}-a_i+a_{i+1}\) steps at the \((i-1)\)-th turn, eliminating both \(a_i\) and \(a_{i+1}\). This is another contradiction. Hence, \(a_{i-1}-a_{i}< 0\).

\(a_{n}> a_{n-1}\) must be true, otherwise the \(n\)-th step is redundant.

Claim 2: \(a_i=|A_i|+|A_{i-1}|\) for all \(i\ge 1\). Furthermore, \(|A_{i+2}|-|A_{i}|=a_{i+2}-a_{i+1}>0\).

Proof: Trivial to check for odd and even \(i\). This claim also implies \(|A_{i+2}|>0\).

Claim 3: For \(i\ge 2\),
\[\sum_{j=1}^{i}a_j=|A_{i}|+2\sum_{j=1}^{i-1}|A_j|\]

Proof: Apply claim 3 for every \(a_j\).

Claim 4: \(|A_i| \leq 2.5|A_{i-1}|\) for \(i\ge 6\).

Proof: Suppose \(i\) is odd since the other case is analogous. Consider \( p= -|A_{i-1}|- 1 < A_{i-1}\) . The robot does not find \(p\) until at least \((i+1)\)-th turn, since \(|A_{i-1}|>|A_{i-3}|>\cdots>|A_2|\) (claim 2). Thus,
\[
\left(\sum_{j=1}^{i}a_j\right)+a_{i}+1 < 8(|A_{i-1}|+1) \] The left hand side can be thought as follows. First, the robot walks to \(A_{i}\). After \(a_{i}\) more steps in the opposite direction, it reaches \(A_{i-1}\) again. Then it takes 1 more step to reach \(p\). Rewrite the left hand side using claim 2 and 3 to get \[ 1+|A_{i-1}|+2\sum_{j=1}^{i}|A_j|< 8|A_{i-1}|+8 \] It simplifies to \[ 1+2|A_{i}|+2\sum_{j=1}^{i-2}|A_j|< 5|A_{i-1}|+8 \] \[ \implies |A_{i}|< 2.5|A_{i-1}|+3.5-\sum_{j=1}^{i-2}|A_j|<2.5|A_{i-1}| \] The last inequality follows from the fact that \(i\ge 6\) and \(|A_j| > 0\). (Q.E.D)

Notice that claim 4 can be strengthened iteratively.

Claim 5: Suppose we have proven \(|A_i|\leq b|A_{i-1}|\) for some \(b\) and for all sufficiently large \(i\). Let \(b’=(2.5-1/b-1/b^2)\). We claim that \(|A_k|\le b’|A_{k-1}|\) for all sufficiently large \(k\).

Proof: Pick a sufficiently large \(i\) such that \(|A_{i-2}|\leq b|A_{i-3}|\) We proceed like in the previous proof to get
\[
|A_{i}|< 2.5|A_{i-1}|+3.5-\sum_{j=1}^{i-2}|A_j|\qquad(1) \] Notice that \[|A_{i-2}|+|A_{i-3}|\ge \left( \frac{1}{b}+\frac{1}{b^2} \right)|A_{i-1}| \] Since \(i\) is sufficiently large, the right hand side of (1) is at most \[(2.5-1/b-1/b^2)|A_{i-1}|=b'|A_{i-1}|\] Hence, \(|A_{i}| < b' |A_{i-1}|\) . (Q.E.D.)

We start with \(b=2.5\) and repeatedly apply claim 5. We find the following values for \(b\):

\[
\begin{array}{cl}
\text{iter.} & b\\
1 & 1.94\\
2 & 1.71883\\
3 & 1.57973\\
4 & 1.46627\\
5 & 1.35287\\
6 & 1.21445\\
7 & 0.99857
\end{array}
\]

We observe that \(|A_{k}|\leq .999 |A_{k-1}|\) for all sufficiently large \(k\). However, that’s a contradiction, since we know \(|A_{k}| > |A_{k-2}|\). Hence, \(OPT\) has to take at least \(8|p|\) steps in the worst case.

Probabilistic Algorithms

The \(8|p|\) bound does not hold for probabilistic algorithms. In fact, it can be shown that my first algorithm takes \(7|p|\) steps in expectation when the first step is taken in a random direction.

Multiplying Two Polynomials with Fast Fourier Transform

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

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")

Lie Detectors and the Story of Halting Turing Machines

Is it possible to design a machine that detects lies? Sure, people have already built devices called polygraphs that monitor blood pressure, respiration, pulse etc. to determine if a person is giving false information. However, that is not same as “detecting lies” because polygraphs merely measure physical effects of telling lies. On the top of that, they are hella inaccurate. I am talking about REAL lie detectors, ie. machines that can instantly validate the statements themselves. Can we actually make such machines? (*insert new startup idea*)

Well, turns out if we ever could make such machines, they would become incredibly handy tools for scientists and mathematicians. That’s because they could babble all sorts of scientific and mathematical conjectures to such a machine, if it beeped, they would know those conjectures were false. Their conversation might go like this:

Crazy Physicist: Higgs Boson exists.
Machine: *Says Nothing*
Crazy Physicist: Ureka! I was right all along.

Computer Scientist: P=NP #changemymind
Machine: *beep*
Computer Scientist: OMG. P=/=NP. I just became a millionaire.

Clearly, you can see how these lie detectors would trivialize almost every known problem of science and mathematics. Then, a very legit question would be why nobody even tried to build such a machine. Indeed, humankind has spared far more effort after making junks like philosopher stones and elixir of life.

Turns out there is something fundamentally wrong about the design of the lie detectors. Before I get to the details, let me tell you a cool life hack. If you ever meet people that claim to be oracles, do this to them:

You: Will I offer you $10?

Whatever the oracle says, do the opposite.

This clearly shows the oracles cannot predict the future. The same thing is true for our lie detectors. Here is why:

Suppose you brought your friend Joe to show him the lie detector you bought. Now do the following:

You: I’m gonna hand Joe $10.

Whatever the lie detector says, do the opposite.

Therefore, just as above, our precious lie detectors cannot exist. However, they did have a purpose. They helped me to show you a proof technique called Cantor’s Diagonal Argument. It is a really nifty tool to prove many advanced statements in computability theory. Here is the general idea:

  • First you assume a machine with utterly unbelievable functions exists.
  • Then you come up for a statement or input for this machine using its own functions.
  • Finally, show that existence of such an input is contradictory to the existence of the machine itself.

As I’m running out of time, I’ll give you just one example of this technique. First proven by Alan Turing, it’s still one of the most famous proofs in computability theory.

Halting Turing Machine (HTM) is a special type of Turing machine that can instantly determine whether another Turing Machine halts on a specific input. In English, this means HTM is a cool app that can tell you whether one of the apps in your smartphone is gonna become unresponsive forever for a specific environment condition inside your smartphone. HTM would be a handy app, right? If HTM tells you an app is gonna be unresponsive beforehand, you’ll save time by not installing that crappy app. However, like all good things in computability theory, turns out HTM cannot exist either. This is because what if an evil developer created an app like this?

  • Give my app’s code and current environment condition of the smartphone to the HTM
  • Whatever the HTM says about my app, do the opposite inside the app.

Hence HTM is unable to predict what this evil dev’s app will do, and that’s why HTM cannot exist.

That’s it for today. I hope you enjoyed this blogpost. If you want to learn Computability Theory for real, I highly recommend Michael Sipser’s book.

ক্যালকুলেটর কীভাবে log-এর মান বের করে?

ব্যাপারটা একটু ক্যালকুলেটরের মত করে চিন্তা করা যাক। সে খুব সাদাসিধে একটি যন্ত্র। ধরা যাক, সে লগারিদম বোঝে না, কিন্তু বড়ো বড়ো সংখ্যার যোগ, বিয়োগ, গুণ, ভাগ বা বর্গমূল এক নিমিষে বের করে ফেলতে পারে। তাকে লগারিদম না শিখিয়ে আমার \(\log 2017\)-এর মান নির্ণয় করানো দরকার। বলো তো আমি এখন কী করতে পারি?

আমি হয়তো ক্যালকুলেটরকে বলতে পারি যে \(\ln (1+x)\)-এর টেইলর সিরিজ
\[\ln(1+x)=x-\frac{x^2}2+\frac{x^3}3-\cdots\;(|x| < 1)\]

সে ডানপাশের ধারাটির প্রথম কিছু পদের যোগফল নিলে \(\ln (1+x)\)-এর মোটামুটি নিখুঁত মান পাবে। এরপর সহজেই যেকোন ধ্বনাত্নক সংখ্যার \(\log\) নির্ণয় করা সম্ভব, কেননা

\[\log 2017 = \frac{\ln 2017}{\ln 10}=\frac{\ln\left(1-\frac{2016}{2017}\right)}{\ln\left(1-\frac{9}{10}\right)}\]

কিন্তু আমার এভাবে করাটা পছন্দ হল না। কারণ কাজটা বড্ড রসকষহীন।

আরও মজার কিছু ভাবা যাক।

খুব সহজেই দেখা যায় যে \[10^3 < 2017 < 10^4 \Longrightarrow 3 <\log 2017 < 4\] কিন্তু দশমিকের পরের অঙ্কটি কত হবে? \(3.1\) থেকে \(3.9\) পর্যন্ত \(10\) এর পাওয়ার নিয়ে আমি দেখতে পারি কোনটা \(2017\) এর সবচেয়ে কাছাকাছি। এরপর দশমিকের পরের দ্বিতীয় অঙ্কের জন্যও একই কাজ করতে পারি। এভাবে একটি একটি করে অঙ্ক আমি নির্ণয় করে যেতে পারি। কিন্তু প্রতিটি অঙ্কের জন্য এতবার করে \(10\) এর পাওয়ার নেওয়া আসলে অপচয়। একই কাজ বাইনারিতে করলে কেমন হয়? সেক্ষেত্রে প্রতিটি অঙ্ক \(0\) অথবা \(1\).  হিসেব করে দেখা যায় যে \(2017 < 10^{3}\times 10^{\frac{1}{2}}\approx 3162\). অতএব, বাইনারিতে দশমিকের পরের অঙ্কটি  \(0\). এর পরের অঙ্কটি? দেখা যায় যে \(2017 > 10^3\times 10^{\frac{1}{4}} \approx 1778\). অতএব, দশমিকের পরের দ্বিতীয় অঙ্কটি \(1\). বাইনারিতেও আমি একই কাজ পরবর্তী অঙ্কগুলোর জন্যে করে যেতে পারি। এভাবে \(\log\)-এর মান যেকোন ঘর পর্যন্ত নিখুঁতভাবে নির্ণয় করা সম্ভব।

এবার তাহলে ঝটপট কোড লিখে ফেলা যাক!

from math import sqrt
def log(base, num, precision = 40): 
    n, e = float(base), 1.0 
    # Find a power of base greater than num 
    while n <= num: 
        n, e = n**2, e * 2 
    # that's initial guess for log(base, num) 
    approx, log_approx = n, e 
    for i in range(precision): 
        # improve guess 
        n, e = sqrt(n), e / 2 
        if approx / n >= num: 
            approx = approx / n 
            log_approx -= e 
    return log_approx

পুনশ্চ ১: তুমি কি ধরতে পেরেছ উপরের কোডটা যে আসলে বাইনারি সার্চ?

পুনশ্চ ২: সব ক্যালকুলেটর দ্বিতীয় পদ্ধতি অনুসরণ করে না। 

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

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

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

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

ইন্টিজার ফ্যাক্টোরাইজেশন: পর্ব ১ পোলার্ডের রো (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\)) এলগোরিদম। পরবর্তী পর্বে এটি নিয়ে আলোচনা করব।