CODEWITHSUNDEEP
matlab
user1210230
user1210230

Reputation: 233

Is this a bug relating to sub2ind / ind2sub?

I was trying to do sub2ind on single-precision variables that I came across the following strange behavior. For example, when I try:

[a b] = ind2sub([50000 50000], sub2ind([50000 50000], single(1000), single(1000)))

I get:

a = 1001
b = 1000

Is this a bug or I am missing something? I know this can be because of an overflow somewhere in matlab's code but it shouldn't happen, right?

I get the same wrong behavior from 64-bit (glnxa64) R2012a, R2011a, R2010b, R2010a, but correct results from 32-bit (glnx86) R2010b.

Upvotes: 2

Views: 982

Answers (2)

Chris
Chris

Reputation: 1531

The reason this happens is the following.

On line 35 of ind2sub.m is 

vi = rem(ndx-1, k(i)) + 1;

where ndx is the index that is passed in. Thus ndx is 49951000 via the call from sub2ind. Now, when you pass in a single precision value, it forces matlab to evaluate all math in single precision. Thus, compare the difference in what happens on l.35.

K>> 49951000-1
ans =
        49950999

K>> single(49951000)-1
ans =
        49951000

This subtraction of a small number from a large number is the issue. So no, this isn't a bug, it's a limitation to the precision of single-precision floating point. This might help some.

edit: As pointed out by Rasman, a nice way of showing this is with eps

eps(single(49951000))=4

so any value added or subtracted from single(49951000) that is on the range (-4,4) will result in 495100 being returned due to the accuracy of single-precision.

Upvotes: 6

Rasman
Rasman

Reputation: 5359

Been thinking some more, especially looking a gnovice's comment which sparked everything.

so here's the thing: a single is a 32-bit signed floating point number. Thus, there is a limitation onto how precise you can get your data

eps(single(49951000)) = 4, while eps(single(50000^2)) = 256

Thus your range is off, and you will get quantization especially as you increase the values of the indices: 

for i = 1000:1010
    sub2ind([50000 50000], single(i), single(1000))
end

If you really want a 32 bit number, I suggest you use uint32, as you will be able to represent all data points (a square matrix of order 50000 has 2.5e+09 data points > 2^32 (=4.3e9)) .

for i = 1000:1010
    sub2ind([50000 50000], uint32(i), uint32(1000))
end

Upvotes: 1

Related Questions