function [k,l,mindupl] = HRtaxiNr(p,M)
% HRTaxiNr Hardy's "uninteresting" taxinumber 1729
%   Ramanujan: It is the smallest numbers that can be expressed
%   as the sum of cubes of two numbers.
%   [k,l,mini]=HRtaxinr(3)  returns the solution.
%   This can be used to solve the same prpblem with other powers as well
%   as long as computing resources allow.
%   Works instantly also if p=4, which was Hardy's next question.
%    Ramanujan:"It must be a very large number.
%   Can be genaralized to variuos directions. In addition to higher powers
%   there have been some investigations for more than two summands.
%   Exercise: Extend this to the case of 3 summands, then more ... 

if nargin==1
    switch p
        case 2
         M=9;
        case 3
         M=16;
        case 4    
         M=200;
        otherwise    
         warning(['No default M-value for p = ' num2str(p)])    
    end
end   

i=1:M;
j=1:M;
[I,J]=meshgrid(i,j);
%
% Compute a matrix of sums of power p for each mesh-point:
%
powersums=I.^p + J.^p;

%% Remove upper diagonal to avoid "symmetry-matches"
%
pickLD=logical(tril(ones(M))); % Lower diagonal logical matrix.
powV=powersums(:); % powersums-matrix columwise into long vector.            % 
powLD=powV(pickLD(:)); % Pick lower diagonal entries.
%
sortedS=sort(powLD);
differences=diff(sortedS);
% Note that difference-vector has 0, at first member of duplicate.
pick=differences==0;
duplicates=sortedS(pick(1:end-1)); 
%%
mindupl=min(duplicates);
[k,l]=find(mindupl==powersums);
 
end

