拉马努金，是一个在数字上极具天赋的印度数学家。在和朋友的对话中，朋友谈到一个出租车的车牌号是毫无趣味的1729，他立刻反映到说1729可以表示为两组不同的数，并且每组都是两个数的立方和。
出租车数的定义就是：可以用n种不同方式表示为两个立方体之和的最小数。目前已知的出租车数有六个：
T(1) = 2 = 1^3 + 1^3;
T(2) = 1729 = 1^3 + 12^3 = 9^3 + 10^3;
T(3) = 87539319 = 167^3 + 436^3 = 228^3 + 423^3 = 255^3 + 414^3;
T(4) = 6963472309248 = 2421^3 + 19083^3 = 5436^3 + 18948^3 = 10200^3 + 18072^3 = 13322^3 + 16630^3;
T(5) = 48988659276962496 = 38787^3 + 365757^3 = 107839^3 + 362753^3 = 205292^3 + 342952^3 = 221424^3 + 336588^3 = 231518^3 + 331954^3;
T(6) = 24153319581254312065344 = 582162^3 + 28906206^3 = 3064173^3 + 28894803^3 = 85192813^3 + 28657487^3 = 16218068^3 + 27093208^3 = 17492496^3 + 26590452^3 = 18289922^3 + 26224366^3;
如果用程序实现寻找出租车数，首先想到的可能就是用for循环迭代试错，但由上面可见用编程手段寻找出租车数，需要的计算量很大，无脑迭代的运行时间会非常长。MATLAB作为矩阵工作室，内置丰富的函数，现使用一种新的方法，把可能的情况生成到矩阵中，在矩阵中缩小范围，找到满足条件的拉马努金数。
首先我们先编写一个函数寻找满足T(2)的出租车数。编写一个函数TaxicabNum，使[a,b,c,d,M] = TaxicabNum(N)返回满足大于N的最小的出租车数M = a3 + b3 = c3 + d3，即输出M（大于或等于N），以及a,b,c和d的四个值。
主要的思路是：先把可能的情况存成矩阵，再使用matlab函数在矩阵中搜索。
代码如下：
function [a,b,c,d,M] = TaxicabNum(N)
%输入正整数N，找到当M>=N时最小的M，以及a，b，c，d，使a3+b3=c3+d3=M。
limit = ceil(N^(1/3));
for max_no = limit:inf
m = zeros(1,max_no);
n = zeros(1,max_no);
for i = 1:max_no
m(1,i) = i;
end
for j = 1:max_no
n(1,j) = j;
end
V1 = m.^3;
V2 = n.^3;
for x = 1:length(V1)
for y = 1:length(V2)
Space(x,y) = V1(x)+V2(y);
end
end
[row,col] = find(Space>=N);
for h = 1:length(row)
after_retrieving(h) = Space(row(h),col(h));
end
tbl = tabulate(after_retrieving);
Index = find(tbl(:,2)==4);
if Index~=0
[R,C] = find(Space==Index(1));
a=R(1);b=C(1);
c=R(2);d=C(2);
M = a^3+b^3;
return
end
end当输入为36301，在命令行调用[a,b,c,d,M]=TaxicabNum(36031)
运行结果为：a = 34; b = 2; c = 33; d = 15; M = 39312;
代码主要部分都是比较容易看得懂的，其中主要使用的函数是find函数和tabulate函数，感兴趣的可以搜一下这两个函数的用法，尤其是find函数在检索矩阵方面还是比较方便的。其中tabuate函数的返回值是三组列向量，第一列是其括号内矩阵中出现次数大于一的数，第二列是其对应出现的次数，第三列是该数出现的次数所占矩阵内元素总个数的比例。由于所有情况组成的矩阵为一对称矩阵，故满足出租车数的元素应在矩阵中出现了2n次。
根据这个思路，稍微更改一下程序，找到满足T(3)的出租车数。为减少计算量，我将搜索范围限定到87000000到88000000，此处也可手动设置。
代码如下：
clc;
clear all;
T = [87000000,88000000];
limit = ceil(T(2)^(1/3));
m = zeros(1,limit);
n = zeros(1,limit);
for i = 1:limit
m(1,i) = i;
end
for j = 1:limit
n(1,j) = j;
end
V1 = m.^3;
V2 = n.^3;
for x = 1:length(V1)
for y = 1:length(V2)
Space(x,y) = V1(x)+V2(y);
end
end
[row,col] = find(Space>=T(1) & Space<=T(2));
for h = 1:length(row)
after_retrieving(h) = Space(row(h),col(h));
end
tbl = tabulate(after_retrieving);
Index = find(tbl(:,2)==6);
[R,C] = find(Space == Index(1));
a = R(1);b = C(1);
c = R(2);d = C(2);
e = R(3);f = C(3);
M = a^3+b^3;
sprintf('%d %d %d %d %d %d %d',a, b, c, d, e,f, M)%显示a, b, c, d, e,f, M运行结果为：436 167 423 228 414 255 87539319
经过调试验证，在寻找满足T(3)的出租车数时，运行时间为10秒，而使用for循环的运行时间超过了四十分钟，用这种算法寻找拉马努金数发挥了MATLAB矩阵运算的优势，相比使用for循环加条件试错，大大的缩短了运行的时间。