Matlab: как представить действительное число как двоичное

Вы, конечно, можете вычислить это в реальном числовом пространстве, но вы рискуете столкнуться с проблемами точности (в зависимости от начальной точки). Если вы заинтересованы в изучении орбит, вы можете предпочесть работать с представлением рациональных дробей. Есть более эффективные способы сделать это, но следующий код иллюстрирует один из способов вычисления ряда, полученного из этой карты. Вы увидите определение периода-n на странице 2 вашей ссылки 2. Из этого кода вы должны увидеть, как вы можете легко работать в реальном числовом пространстве в качестве альтернативы (в этом случае функция Matlab rat восстановит рациональное приближение от вашего реального числа).

[EDIT] Теперь с явной двоичной последовательностью!

% start at some point on period-n orbit
period = 6;
num = 3;
den = 2^period-1;

% compute for this many steps of the sequence
num_steps = 20;

% for each step
for n = 1:num_steps

    % * 2
    num = num * 2;

    % mod 1
    if num >= den
        num = num - den;
    end

    % simplify rational fraction
    g = gcd(num, den);
    if g > 1
        num = num / g;
        den = den / g;
    end

    % recover 8-bit binary representation
    bits = 8;
    q = 2^bits;
    x = num / den * q;
    b = dec2bin(x, bits);

    % display
    fprintf('%4i / %4i  ==  0.%s\n', num, den, b);

end

Ах... для полноты, вот версия с действительным значением. Чистые математики должны сейчас отвести взгляд.

% start at some point on period-n orbit
period = 6;
num = 3;
den = 2^period-1;

% use floating point approximation
x = num / den;

% compute for this many steps of the sequence
num_steps = 20;

% for each step
for n = 1:num_steps

    % apply map
    x = mod(x*2, 1);

    % display
    [num, den] = rat(x);
    fprintf('%i / %i\n', num, den);

end

И, в довершение всего, почему эта реализация быстрая, но глупая? (СОВЕТ: попробуйте установить num_steps на 50)...

% matlab vectorised version
period = 6;
num = 3;
den = 2^period-1;
x = zeros(1, num_steps);
x(1) = num / den;
y = filter(1, [1 -2], x);
[a, b] = rat(mod(y, 1));
disp([a' b']);

Хорошо, это должен быть ответ, а не вопрос, так что давайте ответим на мои вопросы...

Это быстро, потому что он использует встроенную (и высокооптимизированную) функцию filter Matlab для обработки итерации (то есть на практике итерация выполняется в C, а не в M-скрипте). Всегда стоит помнить filter в Matlab, я постоянно удивляюсь тому, как его можно использовать с пользой для приложений, которые не выглядят как проблемы с фильтрацией. Однако filter не может выполнять условную обработку и не поддерживает арифметику по модулю, так что же нам сойдет с рук? Просто потому, что эта карта имеет свойство, согласно которому целые периоды на входе сопоставляются с целыми периодами на выходе (поскольку операция карты умножается на целое число).

Это глупо, потому что очень быстро решает вышеупомянутые проблемы с точностью. Установите num_steps на 50 и посмотрите, как он начнет давать неправильные ответы. Происходит следующее: число внутри операции фильтра становится настолько большим (порядка 10^14), что бит, о котором мы действительно заботимся (дробная часть), больше не может быть представлен в одной и той же переменной двойной точности.

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

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

[N,D] = rat(x0) ;

Когда у вас есть числитель N и знаменатель D, очень легко вычислить ряд x(k+1)=mod(2*x(k), 1) , и вам даже не понадобится цикл.

для части 2*x(k) это означает, что все Numerator(k) будут умножены на последовательную степень 2, что можно сделать путем матричного умножения (или bsxfun для любителя функции):
поэтому 2*x(k) => в Matlab N.*(2.^(0:n-1)) (N равно скаляр, числитель x0, n — это количество членов, которые вы хотите вычислить).

Операцию Mod1 также легко перевести в рациональное число: mod(x,1)=mod(Nx,Dx)/Dx (где Nx и Dx являются числителем и знаменателем x.

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

xn = mod( N.*(2.^(0:n-1).'),D) ;

но для визуального комфорта иногда лучше упростить, поэтому рассмотрим следующую функцию:

function y = dyadic_rat(x0,n)

   [N,D] = rat(x0) ;                   %// get Numerator and Denominator of first element
   xn = mod( N.*(2.^(0:n-1).'),D) ;    %'// calculate all Numerators
   G = gcd( xn , D ) ;                 %// list all "Greatest common divisor"
   y = [xn./G D./G].' ;                %'// output simplified Numerators and Denominators

Если я начну с примера, приведенного в вашей вики-ссылке (x0=11/24), я получу:

>> y = dyadic_rat(11/24,8)
y =
    11    11     5     2     1     2     1     2
    24    12     6     3     3     3     3     3

Если я начну с примера, приведенного Rattus Ex Machina (x0=3/(2^6-1)), я также получу тот же результат:

>> y = dyadic_rat(3/63,8)
y =
     1     2     4     8    16    11     1     2
    21    21    21    21    21    21    21    21