I am trying to generate a plot for the 3x+1 problem for the seed 10^310 but i am unable to get the plot? Why

% Code for 3x+1 Problem
seed = 10^309; % Starting number
x = seed; % Initialize x with the seed
trajectory = x; % Store the trajectory
tic
while x ~= 1
    if mod(x, 2) == 0 % Check if x is even
        x = x / 2;
    else % If x is odd
        x = 3 * x + 1;
    end
    trajectory = [trajectory, x]; % Append x to trajectory
end
toc
% Plot the trajectory
figure;
plot(trajectory, '-o', 'LineWidth', 1, 'MarkerSize', 5,'MarkerFaceColor', 'r');
title('Plot of the Collatz Sequence for positive Seed 10^310');
xlabel('Step');
ylabel('Value');
grid on;

Answer 1

[Matt J]

It's because 1e309 is too large for double float math to handle. The largest integer that can be expressed in double float precision is 2^53. (EDIT - The largest consecutive integer).

mod(1e309, 2)  %Fails

ans = NaN

mod(1 + 2^56 , 2)  %Falsely even number

ans = 0

It's definitely not a problem with plotting, though, as can be seen for smaller seeds below.

% Code for 3x+1 Problem
seed = 1e30; % Starting number
x = seed; % Initialize x with the seed
trajectory = x; % Store the trajectory
tic
while x ~= 1
    if mod(x, 2) == 0 % Check if x is even
        x = x / 2;
    else % If x is odd
        x = 3 * x + 1;
    end
    trajectory = [trajectory, x]; % Append x to trajectory
end
toc

Elapsed time is 0.002391 seconds.

% Plot the trajectory
figure;
plot(trajectory, '-o', 'LineWidth', 1, 'MarkerSize', 5,'MarkerFaceColor', 'r');
title("Plot of the Collatz Sequence for positive Seed " +seed );
xlabel('Step');
ylabel('Value');
grid on;

[Matt J]

One remedy is to use symbolic math:

% Code for 3x+1 Problem
% 3x+1 Problem using Symbolic Math

seed='1e309';
x = sym(seed); % Use symbolic representation of large seed

trajectory=inf(1,2000);
k=1;
trajectory(k)=double(log10(x));

tic
while x ~= 1
    
    k=k+1;
    
    if mod(x, 2) == 0
        x = x / 2;
    else
        x = 3 * x + 1;
    end

    trajectory(k)=double(log10(x));
    
    

end
toc

Elapsed time is 2.866138 seconds.

trajectory(k+1:end)=[];

% Plot the trajectory
figure;
plot(trajectory, '-o', 'LineWidth', 1, 'MarkerSize', 5,'MarkerFaceColor', 'r');
title("Plot of the Collatz Sequence for positive Seed " +seed );
xlabel('Step');
ylabel('log_{10}(Value)');

●

[Steven Lord]

It's because 1e309 is too large for double float math to handle.

Yes, the largest finite double precision value is:

r = realmax

r = 1.7977e+308

How is 1e309 stored in IEEE double precision? It overflows.

q = 1e309

q = Inf

The largest integer that can be expressed in double float precision is 2^53.

Not correct. All exactly repressentable double precision values greater than flintmax are integer values, but not all integer values greater than flintmax are exactly representable in double precision. Let's take flintmax + 1 as an example.

y = flintmax

y = 9.0072e+15

z = y + 1

z = 9.0072e+15

What's the spacing from y to the next largest representable value?

eps(y)

ans = 2

So since 1 is less than 2, flintmax + 1 is not exactly representable as different from flintmax.

y == z

ans = logical
   1

Taking a step of length 3 from flintmax gives the same value as taking a step of length 4.

(y + 3) == (y + 4)

ans = logical
   1

Let's check integer value-ness. Both 2*flintmax and realmax are integer values if we test whether they are equal to their rounded value.

isintegervalue = @(x) x == round(x);
isintegervalue(2*y)

ans = logical
   1

isintegervalue(r)

ans = logical
   1

If you're going to use symbolic numbers, rather than starting off with a piece of text representing the number I'd start off with a value that is exactly representable, convert it into a symbolic variable, then perform symbolic computations with it.

ten = sym(10)

ten = 

10

tenTo309 = ten ^ 309

tenTo309 = 

1.0000e+309

You can even preallocate an array of all 0 values using zeros by saying you want the type to be sym.

symzeros = zeros(1, 10, 'sym') % or

symzeros = 

symzeros2 = zeros(1, 10, 'like', ten)

symzeros2 = 

whos ten tenTo309 symzeros symzeros2

  Name           Size            Bytes  Class    Attributes

  symzeros       1x10                8  sym                
  symzeros2      1x10                8  sym                
  ten            1x1                 8  sym                
  tenTo309       1x1                 8  sym                

[Walter Roberson]

I noticed an oddity in representation that I will be reporting to Mathworks

ten = sym(10)

ten = 

10

tenTo309 = ten ^ 309

tenTo309 = 

1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

When Steve executed, the output representation was 1.0000e+309 . That was due to an error in initializing the symbolic engine for MATLAB Answers.

@Steven Lord please do not re-run and save your post; I need it as proof of bad output.