Как я получаю реальные целочисленные переполнения в MATLAB/октаве?

Хм, да ...

На самом деле, я смог решить проблему с помощью моей специальной подпрограммы "переполнение" ... Теперь она работает мучительно медленно, но без неожиданного поведения! Моя ошибка заключалась в отсутствии раунда (), поскольку Matlab / Octave будет вносить небольшие ошибки.

Но если кто-то знает более быстрое решение, я был бы рад попробовать!

function ret = overflow_sg(arg,bw)

    % remove possible rounding errors, and prepare returnvalue (if number is inside boundaries, nothing will happen)
    ret = round(arg);

    argsize = size(ret);

    for i = 1:argsize(1)
        for j = 1:argsize(2)
            ret(i,j) = flow_sg(ret(i,j),bw);
        end
    end

end%function

%---

function ret = flow_sg(arg,bw)
    ret = arg;
    while (ret < (-2^(bw-1)))
        ret = ret + 2^bw;
    end

    % Check for overflows:
    while (ret > (2^(bw-1)-1))
        ret = ret - 2^bw;
    end
end%function

Я выполнил следующий фрагмент кода

test = int32(2^31-12);
for i = 1:24
    test = test + int32(1)
end

с неожиданными результатами. Похоже, что для Matlab intmax ('int32') + 1 == intmax ('int32') . Я использую 2010a на 64-битной Mac OS X.

Не уверен, что это ответ, еще одно подтверждение того, что Matlab ведет себя нелогично. Однако в документации к функции intmax () говорится:

Любое значение, превышающее значение, возвращаемое intmax, насыщается до значения intmax при преобразовании в 32-битное целое число.

Я полагаю, что Matlab ведет себя, как описано в документации.

Если 64 бит достаточно, чтобы не переполняться, и вам нужно много таких индексов, возможно, сделайте так:

function ret = overflow_sg(arg,bw)
  mask = int64(0);
  for i=1:round(bw)
    mask = bitset(mask,i);
  end
  topbit = bitshift(int64(1),round(bw-1));
  subfrom = double(bitshift(topbit,1))


  ret = bitand( int64(arg) , mask );
  i = (ret >= topbit);
  ret(i) = int64(double(ret(i))-subfrom);
  if (bw<=32)
    ret = int32(ret);
  end
end

Почти все делается как вычисление матрицы, и многое делается с битами, и все делается за один шаг (без циклов while), так что это должно быть довольно быстро. Если вы собираетесь заполнить его rand, вычтите 0.5, поскольку предполагается, что он должен округлять до целочисленных значений (а не усекать).

Если вы хотите получить числовые операции в стиле C, вы можете использовать функцию MEX для прямого вызова операторов C, и по определению они будут работать как типы данных C.

Этот метод намного работает больше, чем переопределения gnovice, но он должен лучше интегрироваться в большую кодовую базу и безопаснее, чем изменение определения встроенных типов, поэтому я думаю, что его следует упомянуть для полнота.

Вот файл MEX, который выполняет операцию C "+" над массивом Matlab. Сделайте одно из них для каждого оператора, для которого требуется поведение в стиле C.

/* c_plus.c - MEX function: C-style (not Matlab-style) "+" operation */

#include "mex.h"
#include "matrix.h"
#include <stdio.h>

void mexFunction(
                 int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]
                 )
{
    mxArray     *out;
    /* In production code, input/output type and bounds checks would go here. */
    const mxArray     *a = prhs[0];
    const mxArray     *b = prhs[1];
    int         i, n;
    int *a_int32, *b_int32, *out_int32;
    short *a_int16, *b_int16, *out_int16;

    mxClassID datatype = mxGetClassID(a);
    int n_a = mxGetNumberOfElements(a);
    int n_b = mxGetNumberOfElements(b);
    int         a_is_scalar = n_a == 1;
    int         b_is_scalar = n_b == 1;
    n = n_a >= n_b ? n_a : n_b;
    out = mxCreateNumericArray(mxGetNumberOfDimensions(a), mxGetDimensions(a),
            datatype, mxIsComplex(a));

    switch (datatype) {
        case mxINT32_CLASS:
            a_int32 = (int*) mxGetData(a);
            b_int32 = (int*) mxGetData(b);
            out_int32 = (int*) mxGetData(out);
            for (i=0; i<n; i++) {
                if (a_is_scalar) {
                    out_int32[i] = a_int32[i] + b_int32[i];
                } else if (b_is_scalar) {
                    out_int32[i] = a_int32[i] + b_int32[0];
                } else {
                    out_int32[i] = a_int32[i] + b_int32[i];
                }
            }
            break;
        case mxINT16_CLASS:
            a_int16 = (short*) mxGetData(a);
            b_int16 = (short*) mxGetData(b);
            out_int16 = (short*) mxGetData(out);
            for (i=0; i<n; i++) {
                if (a_is_scalar) {
                    out_int16[i] = a_int16[0] + b_int16[i];
                } else if (b_is_scalar) {
                    out_int16[i] = a_int16[i] + b_int16[0];
                } else {
                    out_int16[i] = a_int16[i] + b_int16[i];
                }
            }
            break;
        /* Yes, you'd have to add a separate case for every numeric mxClassID... */
        /* In C++ you could do it with a template. */
        default:
            mexErrMsgTxt("Unsupported array type");
            break;
    }

    plhs[0] = out;
}

Затем вам нужно выяснить, как вызвать его из кода Matlab. Если вы пишете весь код, вы можете просто везде вызывать «c_plus (a, b)» вместо «a + b». В качестве альтернативы вы можете создать свой собственный числовой класс-оболочку, например. @cnumeric, который содержит числовой массив Matlab в своем поле и определяет plus () и другие операции, которые вызывают соответствующую функцию MEX стиля C.

classdef cnumeric
    properties
        x % the underlying Matlab numeric array
    end
    methods
        function obj = cnumeric(x)
            obj.x = x;
        end

        function out = plus(a,b)
            [a,b] = promote(a, b); % for convenience, and to mimic Matlab implicit promotion
            if ~isequal(class(a.x), class(b.x))
                error('inputs must have same wrapped type');
            end
            out_x = c_plus(a.x, b.x);
            out = cnumeric(out_x);
        end

        % You'd have to define the math operations that you want normal
        % Matlab behavior on, too
        function out = minus(a,b)
            [a,b] = promote(a, b);
            out = cnumeric(a.x - b.x);
        end

        function display(obj)
            fprintf('%s = \ncnumeric: %s\n', inputname(1), num2str(obj.x));
        end

        function [a,b] = promote(a,b)
        %PROMOTE Implicit promotion of numeric to cnumeric and doubles to int
            if isnumeric(a); a = cnumeric(a); end
            if isnumeric(b); b = cnumeric(b); end
            if isinteger(a.x) && isa(b.x, 'double')
                b.x = cast(b.x, class(a.x));
            end
            if isinteger(b.x) && isa(a.x, 'double')
                a.x = cast(a.x, class(b.x));
            end
        end
    end

end

Затем оберните свои числа в @cnumeric там, где вы хотите поведение int в стиле C, и произведите с ними математические расчеты.

>> cnumeric(int32(intmax))
ans = 
cnumeric: 2147483647
>> cnumeric(int32(intmax)) - 1
ans = 
cnumeric: 2147483646
>> cnumeric(int32(intmax)) + 1
ans = 
cnumeric: -2147483648
>> cnumeric(int16(intmax('int16')))
ans = 
cnumeric: 32767
>> cnumeric(int16(intmax('int16'))) + 1
ans = 
cnumeric: -32768

Вот ваше поведение переполнения в стиле C, изолированное от взлома примитивного типа @ int32. Кроме того, вы можете передать объект @cnumeric другим функциям, которые ожидают регулярные числа, и он будет «работать» до тех пор, пока они будут обрабатывать свои входные данные полиморфно.

Предупреждение о производительности: поскольку это объект, + будет иметь меньшую скорость отправки метода вместо встроенного. Если у вас мало обращений к большим массивам, это будет быстро, потому что фактические числовые операции выполняются в C. Многие вызовы небольших массивов могут замедлить работу, потому что вы много платите накладные расходы на вызов каждого метода.