The Internet vs a cycler

I’ve been reading Programming Pearlsfor some time now. It’s a great book, the kind of book that makes you realize what programming really is about.

Among the exercises behind Column 6 in it, there was this strange problem: “At what distances can a courier on a bicycle with a reel of magnetic tape be a more rapid carrier of information than a telephone line that transmits 56,000 bits per second? Than a 1200-bps line?”

The thought that a bicycle can be faster than an electronic communication line seems ridiculous indeed, but it turns out it’s not difficult. For the sake of modernity, let’s replace the ‘magnetic tape’ with a hard disk. 1TB disks are common enough these days that I’ll use that for the calculations.

Now, let’s calculate the time taken for transmitting ‘g’ Gigabytes of data using both the bicycle method and an Internet line, between two places which are ‘d’ kilometres apart. In the bicycle method, we have two kinds of ‘time’s to consider: the time taken to transfer things to and from the hard disk, and the time taken in the bicycle itself. These days, data transfer rates in hard disks lies around 25MB/s from my experience, so I’ll use that as the transfer rate here. Then, 1 gigabyte will take 1GB/25MB ~= 40 seconds. For g GB, it is 40g seconds. This is for one transfer. Since we do this twice (transferring to the hard disk, and later back from it), it’d be 80g seconds in total. In addition, we have the time taken in the bicycle itself. Assuming a minimal speed of 10kmph by the bicycler, the time taken would be d/10 hours, which in seconds would be 360d. So, the total time is 360d + 80g.

For a 56kbps Internet line, the time taken for transferring the same g gigabytes would be  the upload time + download time. From discussions in the Internet, it appears that upload speeds were “from 1/2 to less of the the download speed”. Of course, this is not ‘authoritative’, but is good enough for our purposes. The download speed here is 7KB/s (since 56kbps = 7KB/s). Let’s take the upper limit and take the upload speed to be half of download i.e. 3.5KB/s. Then, time taken = 1000,000g/3.5 + 1000,000g/7 seconds (since 1GB ~= 1000,000KB) = 3000,000g/7. So we’ve got to solve for d in:

360d + 80g < 3000,000g/7

360d < (3000,000 – 560)g/7

d < 2999440g/2520

d < 1200g (approx.)

If g = 1GB, the bicycle is faster for upto a distance of 1200 km.

If we are to use the 1TB hard disk to its full capacity, g ~= 1000GB. Substituting that,

d < 1,200,000 km

which means that upto an astounding distance of about 1200,000 kilometres (that’s about 750,000 miles), a bicycle with a 1TB hard disk is faster than a 56kbps line for transferring 1TB of data!

Now, 56kbps is one of those old dialup lines – the more common internet connection these days is a broadband connection. For purposes of modernity, let’s now assume a 1Mbps connection, which transfers at a rate of 128KB/s. Again, the upload speed is less. And in this case, it appears that upload speeds of 1/4th of download are most common. Hence, let’s assume an uplink of 128/4 = 32KB/s. Then, the inequality is:

360d + 80g < 1000,000g/32 + 1000,000g/128

360d < 4999920g/128

d < 108g

Again, substituting g = 1,

d < 108 km

For 1024 GB (i.e., 1TB data),

d < 110,000 km

Given that the earth’s circumference itself is just 40,000 km, this means that if you need to transfer 1TB of data to anywhere on earth, you’re better off sending a bicycle courier than sending it through today’s internet connections (of course, there’s the issue of crossing the oceans, which we’ll ignore for simplicity ;) )

Now, the same is not true if you need to transfer only 1GB – as we saw above, a bicycle is faster only upto 108 km. So, for what amount of data is a bicycle faster to anywhere on earth? Let’s do one final calculation.

The maximum distance you’d need to travel on earth is half its circumference (can you see why?). This is 20,000 km. So, let’s make d=20,000 km in the inequality and see what g we get.

360*20,000 < 4999920g/128

g > 184 GB

So, the moral of the story is, even if you have a 1Mbps broadband line, if you need to transfer more than 184 GB of data, you’d achieve better speeds to anywhere on earth by carrying a hard disk on a bicycle than by transferring it through the 1Mbps connection!

Finding the Largest prime factor efficiently

Long long ago, I wrote an article on finding the largest prime factor of a number (if you haven't read it, I'd recommend reading it now as much of this article won't make sense without that). It was just a fun little thing I did for solving Project Euler, so I didn't expect it would invite such excellent comments and get posted to Reddit!

And almost every one of the comments pointed out how the algorithm could be made better. I hadn't really thought there would be more optimal solutions than the ones given out by Project Euler, but turns out there are!

The first improvement was one I might have done myself if I had tried to do an iterative version myself. In the recursive version, when I find the number is divisible by (say) 2, I divide out the 2 and call the function recursively. Now, if there's another 2 in there (that is, if the number was divisible by 4), the number would again get divided by 2 and another recursive call would occur. When I'm converting this algorithm to an iterative version too, the most natural way would be:

my $max_factor = 1;
while($n % 2 == 0)
{
$max_factor = 2;
$n /= 2;
}

for(my $i = 3; $i < $n; $i += 2)
{
while($n % $i == 0)
{
$max_factor = $i;
$n /= $i;
}
}

There's a subtle optimization (compared to the recursive version) that has happened here that you probably missed. It's actually a simple one - in the recursive version, when you've divided n by say 2, you pass it back to the same function - which again tries to divide it by 2. This is good as long as n is divisible by 2, but once you've divided out all the 2's, this becomes just a waste of CPU cycles. In the iterative version, this is avoided because, once we know we've divided out all the 2's, we no longer try dividing by that, and instead move on to greater things. That's a Good Thing, as for big numbers, this optimization might save several cycles.

But even this turns out to be a naive implementation. As someone in the math subreddit as well as one of the commenters here pointed out, there's another optimization that can be done here - after 2 and 3, what are the prime numbers you have? 5, 7, 11, 13, 17, 19, 23, 29, etc. Is there some common property these numbers have? Why, yes, they are all either of the form 6*n+1 or 6*n-1. Turns out that it has been proved that all primes after 2 and 3 are of either of these forms (second half of second paragraph). So, instead of trying out all the odd
numbers, we can try only the numbers of these forms. The first code I wrote for this purpose was,

$i = 6;
while($n > ($factor = $i - 1))
{
$n /= $factor while $n % $factor == 0;

$factor = $i + 1;
last if($n == $factor);
$n /= $factor while $n % $factor == 0;

$i += 6;
}

That basically just implements the concept of trying out 6n +/- 1 for each value of n. After that, I remembered a comment on the previous post which had made me realize that this is the same as just adding 2 and 4 in alternate iterations of the check. Just for the heck of it, I tried that too:

$factor = 5;
while($n > $factor)
{
$n /= $factor while $n % $factor == 0;
$factor += 2;
last if($factor == $n);

$n /= $factor while $n % $factor == 0;
$factor += 4;
}

Finally, I put these all in different subroutines in a Perl program, and benchmarked those using Perl's Benchmark.pm module. It's a very easy to use module - my entire benchmarking code consists of

# Upper limit for random number generation
my $rangemax = 1000_000;
cmpthese(-10, {
# Recursive => " largestprimefRecurs((int(rand($rangemax))))",
Naive_Iteration => " largestprimefNaiveIter(int(rand($rangemax)))",
Arcane6nplusminus1 => " largestprimef6n(int(rand($rangemax)))",
Arcane6nadding2and4 => " largestprimef6nadding2and4(int(rand($rangemax)))",
Arcane6nusingCount => "  largestprimef6nusingCount(int(rand($rangemax)))"
});

The -10 there asks Benchmark to run each code segment for at least 10 CPU seconds, so that random fluctuations due to irrelevant factors are smoothed out, and we have reliable results.

You might notice that the call to the recursive version (the original one) is commented out. Due to (what I believe to be) a problem with the Benchmark.pm module, I was not able to benchmark it using this method. Benchmark.pm would report that the recursive function ran thousands of times faster than the newer methods (probably counting the recursive calls too), but other testing revealed it was actually much slower:

sundar@System:~$ time for (( i=0; i<500; i++ )); do perl Recursive.pl; done;
real 0m5.631s
user 0m3.848s
sys 0m1.556s

sundar@System:~$ time for (( i=0; i<500; i++ )); do perl 6nMethod.pl; done;
real 0m3.219s
user 0m1.692s
sys 0m1.344s

We see the recursive version is slower by more than 2 seconds! I'm planning to investigate more about this behaviour of Benchmark.pm (a quick googling didn't turn up much), but for now suffice to know the old algorithm is much slower.

There's another entry in the code there that I haven't talked about: Arcane6nusingCount. [By the way, the 'arcane' naming convention is a tongue-in-cheek reference to this comment where the commenter said "I get the feeling that [the programmer who wrote this post] would consider that working with 6n +/- 1 as too arcane a method..." :) I couldn't start the function name with the digit 6, so I added 'Arcane' before each one. ;)]
Anyway, that Arcane6nusingCount was a version where I tried to implement the 6n +/- 1 algorithm by keeping a counter of the iterations, and incrementing by 2 or 4 based on whether count%2 was 0 or 1. Turns out this extra mod operation kills performance, as you can see in the benchmark results below.

Benchmark Results:
Legend:
Due to Blogger's @&%*&%# formatting, I had to change the names to shorter versions as follows:

Naive_Iteration     : Iter
Arcane6nplusminus1  : 6n+/-1
Arcane6nadding2and4 : Add2,4
Arcane6nusingCount  : Count

The percentages are the amounts by which the algorithm in the horizontal (the row) is better than the one in the vertical (the column). The /s amounts are the number of times this subroutine could have been run in one second, averaged over the given inputs. For this test, the inputs were random numbers from 0 to 1000,000, as seen in the code above.


Trial 1

	Rate	Count	Iter	6n+/-1	Add2,4
Count	45.9/s	--	-24%	-35%	-39%
Iter	60.1/s	31%	--	-15%	-20%
6n+/-1	70.5/s	54%	17%	--	-6%
Add2,4	75.2/s	64%	25%	7%	--

Trial 2

	Rate	Count	Iter	6n+/-1	Add2,4
Count	43.1/s	--	-25%	-41%	-42%
Iter	57.4/s	33%	--	-22%	-23%
6n+/-1	73.6/s	71%	28%	--	-2%
Add2,4	74.9/s	74%	30%	2%	--

We see that the method of adding 2 and 4 directly (Arcane6nadding2and4) works best, but the one using multiples of 6 and dividing by $i-1 and $i+1 (Arcane6nplusminus1) is only slightly worse. And, in my opinion, it is clearer and more elegant. They're both much better than the naive iterative method, and the original recursive solution is not even comparable to these!

So, I guess the moral of the story is, look at ways you can optimize your code even if an 'authoritative source' seems to tell you've reached the optimum. Also, where performance is critical, profile and benchmark your code - they are much more reliable practically than theoretical algorithmic analyses (of course, this is true only when they are done right).

The program along with benchmarking code is available at: http://sundaryourfriend.googlepages.com/largestprimefactorBenchmark.pl
To run it, you need to have the Benchmark.pm Perl module installed (which would probably be available by default).

Have a nice day...

Largest prime factor of a number

Long live the programming overlords! Today is the day I have begun serious programming again... 'serious' is perhaps a bit too much - programming is never actually serious, it's a sort of fun activity that happens to appear boring and difficult and hence obtains good paychecks from the HR types... ;)

Anyway, I happened to come across the Project Euler and registered in the site. The registration was painless and I fell in love with the site immediately. Problems 1 and 2 were trivial (finding the sum of the multiples of 3 or 5 under 1000, finding the sum of even numbers in the Fibonacci series which are less than 4 million).

Problem 3 was a little more involved. I had to find the largest prime factor of 600851475143. It isn't particularly difficult, but it does call for some thinking. The first method I thought of was to

1. divide the number by successive integers upto sqrt(n), and
2. when the remainder is 0, check whether the divisor is prime,
3. if so, store it as answer overwriting the previous value, if any.

An obvious improvement was to come down from sqrt(n) so that the first number to satisfy conditions 2 and 3 would be the answer.
However, by this time, something seemed fishy in this method. Particularly, it occurred to me that a prime number checking was costly and unnecessary here. Instead, you can go about it using a recursive solution:

FindLargestPrimeFactor(n):

1. divide the number by successive integers (each denoted by i) upto sqrt(n),
2. when the remainder is 0, return the maximum among i and FindLargestPrimeFactor(n/i).

It might not be clear at first sight that the algorithm works correctly. There seems to be no check for primality. However, this algorithm exploits the fact that the smallest factor sf of a number n will necessarily be prime. This is because, if it were not prime, its factors would be factors of n and hence sf would not be the smallest factor. Thus, we need only find the smallest factor to have a prime factor. Then, we can divide away this factor and again find the smallest factor. This way, we can have all the prime factors of a number.

In the algorithm, we find a prime factor of n, compare it with the largest prime factor of n/i, and return the larger number. Thus, finally, we'll have the largest prime factor of the number n itself.
Obviously, we need a limit for the recursion. In this case, this limit occurs when i overflows sn. This means that n is prime, and hence we must return n itself as the largest prime number.

This was the logic of my program as I coded and finished it. I entered the answer into the Project Euler website (which is projecteuler.net, by the way). After that, I saw the forums in the website and saw that many had implemented my original algorithm which used prime testing. This more efficient recursive solution had missed the neurons of many... With that boost to my hubris, I turned to the pdf containing the 'overview' of the problem.

Wow, the algorithm described is very close to mine. One improvement is that the recursion has instead become an iteration. This is one area where I stumble every time - converting a recursive program into an iterative one. So I bookmarked this in my mind and moved on.

The next improvement was one that made me feel stupid for not thinking of it earlier. Since I'm going to find only prime factors, I need only use 2 and odd numbers in step 1. So, I can use a separate step for 2 and then start with 3 and increment i by 2. This way, we'll save about half of the iterations...

UPDATE: There are quite a number of other improvements that must be done to this algorithm. For the actual algorithm for this problem (which develops on the algorithm here), see Finding the Largest prime factor efficiently.

The final perl program I wrote was:

use strict;
use warnings;

my $magic = 600851475143;

sub largestprimef($);
sub max($$);

print largestprimef($magic);
sub largestprimef($)
{
my $n = shift;

my $i;
return largestprimef(max(2, $n/2)) if($n % 2 == 0); 

my $sn = int(sqrt($n));

for ($i = 3; $i <= $sn; $i += 2)  { 

  last if($n % $i == 0);

} 

if($i > $sn) #loop ran over, means the number is prime
{
  return $n;
}
else
{
return max($i, largestprimef($n/$i));
}
}

sub max($$)
{
return (sort { $a <=> $b }(@_))[1];
}

A few notes about the program here:

• not 'use'-ing the strict and warning pragmas gave a speed improvement of a few milliseconds
• first I implemented max in the good old comparison way - if i > j, return i else return j. Turns out perl's sort implementation gives faster results than this equivalent speed. I just tried it for the heck of it , and was surprised to see a few milliseconds improvement. Probably the speed improvement came from the removal of 'shift's to get the numbers from @_ to $i and $j.(that was a temporary illusion; in repeated testing both turned out to be equivalent.
• Since I did the program in Linux, the age old newline problem cropped up. When I pasted the code into blogger, it showed the code as one continuous mess. I had to use http://www.hashemall.com/ for converting the newlines. The site dutifully hashed the program, and thankfully returned my original text with MS Windows newlines in place of Unix ones.. (The hash value under sha-1 was b8b7188ad3440084e5db7f44fd70d599f897e33a, in case you are interested... :) )
• The final runtime was around 14 ms.
• UPDATE: The code highlighting is courtesy of Online syntax highlighter 'tohtml' which is an excellent service by the way.

Law of small numbers and such...

The first strong law of small numbers (Gardner 1980, Guy 1988, Guy 1990) states "There aren't enough small numbers to meet the many demands made of them."

The second strong law of small numbers (Guy 1990) states that "When two numbers look equal, it ain't necessarily so." Guy (1988a) gives 35 examples of this statement, and 40 more in Guy (1990). For example, example 35 notes that the first few values of the interpolating polynomial [...] for , 2, ... are 1, 2, 4, 8, 16, .... Thus, the polynomial appears to give the powers of 2, but then continues 31, 57, 99, ... (Sloane's A000127).

I just read this at Wolfram Mathworld; and something flashed within me. What if the differences between this poly and the actual powers of two were significant? Well, being terribly bored and having a lot of time to kill, I set out to explore that. I was stuck at the first step itself - I was in Linux, and didn't know any way to do math explorations like this. Had I been in Windows, I would have used Matlab (damn good and surprisingly easy to use). I'm sure there's some 'popular' Free Software for this kind of job, but I was not in a mood to search for software packages and get into dependency hell just yet. So, I just asked myself, can I do it with what I do know to do in Linux?
Writing C programs. Perl scripts. Shell scripts. For this kind of job, Perl somehow seemed more appropriate. So, I typed out a Perl script to find the differences (here f(n) is the 4th order polynomial given in the above article):

n	f(n)	2^n	diff
6	31	32	1
7	57	64	7
8	99	128	29
9	163	256	93
10	256	512	256
11	386	1024	638
12	562	2048	1486
13	794	4096	3302
14	1093	8192	7099
15	1471	16384	14913
16	1941	32768	30827
17	2517	65536	63019
18	3214	131072	127858
19	4048	262144	258096
20	5036	524288	519252
21	6196	1048576	1042380
22	7547	2097152	2089605
23	9109	4194304	4185195
24	10903	8388608	8377705

No obvious pattern emerged. Then, on second thoughts, there did seem to be some sort of pattern to it. The differences are all increasing. But each one increased by a lesser amount than the previous one. 7/1 is (quite obviously) 7, 29/7 is slightly greater than 4, 93/29 is slightly greater than 3, etc. As the numbers increased, it seemed that they increased by successively smaller factors. The numbers had some kind of geometrical pattern - sort of like a geometric series where you get successive terms by multiplying each term by a constant ratio. However, in this case, the ratio wasn't exactly constant. Instead, it seemed to be tending towards some constant.

To check this out, I then printed the quotient obtained by dividing each difference by its predecessor. Lo and behold, suddenly, there was Light:

n	f(n)	2^n	diff	quotient
6	31	32	1	~
7	57	64	7	7
8	99	128	29	4.14285714285714
9	163	256	93	3.20689655172414
10	256	512	256	2.75268817204301
11	386	1024	638	2.4921875
12	562	2048	1486	2.32915360501567
13	794	4096	3302	2.22207267833109
14	1093	8192	7099	2.14990914597214
15	1471	16384	14913	2.10071841104381
16	1941	32768	30827	2.06712264467243
17	2517	65536	63019	2.04427936549129
18	3214	131072	127858	2.02888017899364
19	4048	262144	258096	2.01861440035039
20	5036	524288	519252	2.01185605356147
21	6196	1048576	1042380	2.00746458367036
22	7547	2097152	2089605	2.00464801703793
23	9109	4194304	4185195	2.00286417767951
24	10903	8388608	8377705	2.00174782775952

So, there you have it: the quotient of dividing each successive difference between 2^n and this f(n) is approaching 2. For n=70 and higher, it was so close to 2 that Perl gave it up and printed it as just 2. So, the quotients are converging towards 2 as far as we can see. Each successive difference, then, is obtained by multiplying the previous difference by a number very close to 2.
That is the kind of beauty math has. You obtain a sequence that looks like powers of 2 for small numbers, you find that it deviates from the actual value for higher powers, then you find that the deviation itself is increasing in powers of two. This thing has some inexplicable beauty to it... :)
So, now that we've found that this thing seems to converge to 2, how do we interpret it? That's the kind of thing for which you need mathematical training. If I had had some good math course on series and sequences, may be I would know what to do with this. May be I'd theoretically prove that the quotients actually converge to 2 as n tends to infinity. May be we can show that this f(n) can indeed be used to find powers of 2.
For now, let's get into the nitty gritty of this number crunching. First, I wrote this cute little program (any claims to its uncuteness are by Agents of Satan; avoid them like plague) :

#!/usr/bin/perl

use warnings;
use strict;

my @range = (6 .. 24); #6 is where the differences start to be non-zero
my $i;
my @diffs;

foreach $i (@range)
{
print $i,"\t\t";
print (( $i**4 - 6*$i**3 + 23*$i**2 - 18*$i + 24 ) / 24, "\t\t");
print 2**($i - 1), "\t\t";
push @diffs, (2**($i - 1) - (( $i**4 - 6*$i**3 + 23*$i**2 - 18*$i + 24 ) / 24));
print $diffs[$#diffs], "\t\t";
$#diffs && print $diffs[$#diffs]/$diffs[$#diffs-1];
print "\n";
}

Then, I ran the output through txt2html with --make_tables option. I got the HTML code for a table, which you can see with naked eye by right clicking here and choosing 'View source' (or whatever option reads close to that). Using the code, I made a blog post named 'Law of small numbers and such...'.

You just finished reading it!

A cool captcha!

Are you human?
This is a question we often get asked when registering for things on the net - the sign up page gives a scrambled image containing some text, and we are asked to type out the characters from that.
This is to make sure that we are indeed human, and not some software created to make fake registrations for illicit usage. The idea is that software programs cannot 'see' the text within the image, while we humans can.
Now, this site has a cool new way of doing it: instead of asking to type out letters from an image, it asks us to solve a math problem. It is usually a simple problem which can be solved by anyone with undergrad level math practice. I got 'what is the least zero of the polynomial x^2 + 3x' (the answer is -3).
While not everyone might be able to solve these problems, the site registration too is not for everyone - it's a registration to download a 'quantum random bit generation service' client; someone downloading that can probably solve these problems hands down...

PS: No, I didn't download the client. I registered only because I heard about this cool method of doing a captcha and wanted to give it a try... :)