Can I believe the values of these large integers?

% coinStackv4
%
% Geoffrey T Fox
%
% https://projecteuler.net/problem=78
%
% Based on Bob's Excel work
%
% 27 December 2022
%
clear; clc; tic;
mAx = 10000;
daTa = zeros(mAx,mAx);
daTa(:,2) = 2;
daTa(:,1) = 1;
daTa(1,:) = 1;
daTa(3,3) = 3;
for idx = 3:mAx
    daTa(2,idx) = 1 + floor(idx/2);
end
for idx = 3:mAx % mAx
    jMid = 1 + idx;
    for jdx = 3:mAx
        if idx < jdx
            daTa(idx,jdx) = daTa(idx-1,jdx)+ daTa(idx,jdx-idx);
        elseif idx == jdx
            daTa(idx,jdx) = daTa(idx-1,idx) + 1;
        elseif idx > jdx
            daTa(idx,jdx) = daTa(idx-1,jdx);
        end
    end
end
for idx = 1:mAx
    if mod(daTa(idx,idx), 1000000) == 0
        fprintf("%3.0f yields  %15.0f\n",idx,daTa(idx,idx))
        fprintf("above line is divisable  by 1,000,000\n")
    end
end

2301 yields  17022871133751703055227888846952967314604032000000

above line is divisable  by 1,000,000

toc

Elapsed time is 1.925014 seconds.

[Jan]

For which reasons do you hesitate to believe what?

[Geoffrey Fox]

When I run the output to screen shows

2301 yields 17022871133751703055227888846952967314604032000000

but that last integer is much greater than 2^64 which I understand to be the largest integer possibe under 64 bit.

[Stephen23]

"...which I understand to be the largest integer possibe under 64 bit."

The largest possible integer using DOUBLE() class is:

realmax()

ans = 1.7977e+308

https://stackoverflow.com/questions/1848700/biggest-integer-that-can-be-stored-in-a-double

Possibly you are actually thinking of the largest contiguous integer, which is 2^53:

flintmax()

ans = 9.0072e+15

"Can I believe the values of these large integers?"

You should never believe anything that any calculation tells you.

Answer 1

[John D'Errico]

clear; clc; tic;
mAx = 10000;
daTa = zeros(mAx,mAx);
daTa(:,2) = 2;
daTa(:,1) = 1;
daTa(1,:) = 1;
daTa(3,3) = 3;
for idx = 3:mAx
    daTa(2,idx) = 1 + floor(idx/2);
end
for idx = 3:mAx % mAx
    jMid = 1 + idx;
    for jdx = 3:mAx
        if idx < jdx
            daTa(idx,jdx) = daTa(idx-1,jdx)+ daTa(idx,jdx-idx);
        elseif idx == jdx
            daTa(idx,jdx) = daTa(idx-1,idx) + 1;
        elseif idx > jdx
            daTa(idx,jdx) = daTa(idx-1,jdx);
        end
    end
end
for idx = 1:mAx
    if mod(daTa(idx,idx), 1000000) == 0
        fprintf("%3.0f yields  %15.0f\n",idx,daTa(idx,idx))
        fprintf("above line is divisable  by 1,000,000\n")
    end
end

2301 yields  17022871133751703055227888846952967314604032000000

above line is divisable  by 1,000,000

Should you "trust" that value? NO!!!!!!!!!!!

A double that is larger than 2^53-1 is not known exactly as an integer. So while doubles can represent floating point numbers larger than 2^53-1, as integers, that is the absolute limit. In terms of an integer, beyond that point all is garbage when you are working with integers.

I'm sorry, but if you think that number should be an integer, you are kidding yourself if you perform these computations using doubles. The least significant bits of that number will be garbage. Even an int64 or uint64 will fail you, due to overflows.

And sadly, trying to perform that large of a computation with either syms (or my own variable precision tools) will fail, due to memory (AND time) problems. That is a huge array, when each element of the array requires more memory than just a double.

What you actually are trying to do I have no clue, since your code has no comments. It is pretty much unreadable as such. If your goal is to compute moduli. For example, to know if a number is actually divisible by that relatively small integer, then you are doing this the wrong way. You need to compute the modulus at EVERY step. But I really have no clue what you are trying to do in the end.

[Geoffrey Fox]

The code SHOULD be generating only integers since it starts with integers and then only uses addition. However the values do get beyond the integer limit at which point I amnow convinced the vallues are not to be believed.

As for comments - % https://projecteuler.net/problem=78

Thanks to all who responded.

[John D'Errico]

Yes, I just saw that this was a PE problem. And I am absolutely positive that I solved that one, though it has been a few years since I was doing them with any seriousness. I did save the m-files I wrote for many of my solutions. (Ok, I did the sudoku one by hand, solving each puzzle myself. No computer needed.)

The point of this problem is though, that you can simply do EVERYTHING modulo 1e6. If you do that, then the result should just drop out. And it will do so pretty quickly, with no need to resort to high precision arithmetic.

Finally, I just noticed that this appears to be just effectively a computation (mod 1e6) of the number of integer partitions of a number. For that, I recall there is a recursive scheme.

https://mathworld.wolfram.com/PartitionFunctionP.html

http://mathworld.wolfram.com/PentagonalNumber.html

You might want to do a litte reading. But still, do everything mod 1e6.

[Geoffrey Fox]

John, I appreciate the suggestions. I'm an 81 year old trying to keep mentally fit working on various problem.

"Integer Partitions" is a new concept to me soI have som learning to do.

Thanks for your help and encouragement.

[John D'Errico]

I feel like I am approaching that age far more quickly than I want. Anyway, PE problems are great fun, and pretty challenging.

Answer 2

[Image Analyst]

Your data is double. You can make it 64 bit unsigned integer like

daTa = zeros(mAx,mAx, 'uint64'); % Pass in the option to make it 'uint64'.

[Geoffrey Fox]

When I do as suggested, daTa = zeros(mAx,mAx, 'uint64') with mAx = 10000 the code gives

>> daTa(10000,10000)

ans =

uint64

18446744073709551615

But if I run the original way I get

2301 yields 17022871133751703055227888846952967314604032000000

above line is divisable by 1,000,000

Elapsed time is 1.009981 seconds.

>> daTa(10000,10000)

ans =

3.616725132563619e+106

[John D'Errico]

@Geoffrey Fox

The uint64 value you got is the largest integer you can compute using a uint64.

intmax('uint64')

ans = uint64 18446744073709551615

Essentially, you caused an overflow. and uint64 just maxes out at that value when you do. The number you got was actually just 2^64-1.

sym(2)^64

ans = 

18446744073709551616

However, the "number" you computed as a double is meaningless, if you think the lower order digits of that number have any meaning.

[Jan]

Remember, that Matlab stores numbers with about 16 valid digits as doubles in the IEEE754 format. This means, that the trailing zeros of 17022871133751703055227888846952967314604032000000 are not meaningful.

The used approach is not applicable.

Answer 3

[Paul]

What happens if you change this line to work in symbolic world? It will be slow .... (couldn't run here, exceeded the 55 second limit)

%daTa = zeros(mAx,mAx);
daTa = zeros(mAx,mAx,'sym');

[Geoffrey Fox]

Error using zeros

CLASSNAME argument must be a class that supports ZEROS, such as 'double' or 'single'.

Error in coinStackv4 (line 14)

daTa = zeros(mAx,mAx,'sym');

[John D'Errico]

Actually

A = zeros(5,5,'sym')

A = 

is allowed in current releases. However, for large arrays, it WILL be slow.

[Paul]

How about

daTa = sym(zeros(mAx,mAx));

[Geoffrey Fox]

I don't have that toolbox

'sym' requires Symbolic Math Toolbox.