MATLAB 실습

%% DNA 저장 실습 — 1~5단계 누적 인코딩
% 1단계부터 5단계까지 같은 스크립트 아래에 순서대로 이어 붙입니다.
% clear는 맨 처음에만 사용합니다. 각 %% 섹션을 위에서 아래로 실행하세요.
% clear는 이전 실행 변수 제거, clc는 명령창을 비워 이번 단계 로그를 보기 쉽게 합니다.
clear; clc;
%% 0. 설정 — 입력 종류와 파일 이름만 바꿔 실행
% 'text'=한글/문자열, 'image'=사진 파일, 'file'=일반 파일
inputMode = 'text';             % 'text' | 'image' | 'file'
textInput = '나노바이오공학전공';         % inputMode='text'일 때 입력 문자열
inputFile = '';                 % 비우면 파일 선택창 (image / file)
% FASTA는 헤더 포함 형식, TXT는 sequence만 있어 붙여넣기·복호 입력이 간단합니다.
outFasta = 'dnas_out.fasta';
outSequencesTxt = 'dnas_sequences.txt'; % plain TXT: one 148 nt sequence per line
F = 'AGCCTTGTGTCCATCAATCC';     % forward primer (20 nt)
R = 'TGCGCTATGGTTTGGCTAAT';     % reverse primer (20 nt)

%% 1. 입력 → byte (Block 1: 이름은 UTF-8, 사진은 raw byte)
switch lower(inputMode)
    case 'text'
        dataBytes = uint8(unicode2native(textInput, 'UTF-8'));
        sourceName = 'text'; ext = 'txt';
    case {'image','file'}
        if isempty(inputFile)
            [fn, fp] = uigetfile('*.*', '인코딩할 파일 선택');
            if isequal(fn, 0), error('파일 선택을 취소했습니다.'); end
            inputFile = fullfile(fp, fn);
        end
        fid = fopen(inputFile, 'rb');
        if fid < 0, error('파일을 열 수 없습니다: %s', inputFile); end
        dataBytes = fread(fid, Inf, '*uint8').'; fclose(fid);
        [~, sourceName, ext0] = fileparts(inputFile);
        ext = lower(strrep(ext0, '.', ''));
        if isempty(ext), ext = 'bin'; end
    otherwise
        error("inputMode는 'text', 'image', 'file' 중 하나여야 합니다.");
end
dataBytes = reshape(uint8(dataBytes), 1, []);
L = numel(dataBytes);
if L == 0, error('입력 byte가 비어 있습니다.'); end
fprintf('입력: %s | %d byte | 확장자 .%s\n', sourceName, L, ext);
nShow = min(32, L);
fprintf('Block 1 HEX (first %d of %d bytes):', nShow, L); fprintf(' %02X', dataBytes(1:nShow));
if L > nShow, fprintf(' ... (%d more bytes)', L - nShow); end
fprintf('\n');

assert(isa(dataBytes, 'uint8') && L > 0);
fprintf('Block 1 PASS | %d bytes\n', L);

%% 2. 17 byte 행 분할 (Block 2: 빈 자리는 PAD=0x1B)
PAD = uint8(27); ROW_BYTES = 17; NSYM = 8;
% 헤더 없는 포맷이므로 사용자가 decoder에 L과 D를 알려 줍니다.
% 데이터 행 D는 짝수로 올림: 앞 절반과 뒤 절반을 XOR합니다.
% data row는 최소 2개·짝수여야 XOR 쌍을 만들 수 있습니다.
D = max(2, 2 * ceil(ceil(L / ROW_BYTES) / 2));
H = D / 2; N = D + H;
if N > 65535, error('고정 index 한계(65535 strand)를 넘었습니다.'); end
rows = repmat(PAD, N, ROW_BYTES);
for row = 1:D
    first = (row - 1) * ROW_BYTES + 1;
    last = min(row * ROW_BYTES, L);
    if first <= last
        rows(row, 1:last-first+1) = dataBytes(first:last);
    end
end
assert(isequal(reshape(rows(1:D, :)', 1, []), [dataBytes repmat(PAD, 1, 17*D-L)]), 'Block 2 분할 검증 실패');
fprintf('Block 2 PASS | %d data rows × 17 byte (D=%d, H=%d, N=%d)\n', D, D, H, N);

disp('첫 데이터 행 (HEX)'); fprintf('%02X ',rows(1,:)); fprintf('\n');

%% 3. 바깥 XOR 행 추가 (Block 3: 같은 위치 byte끼리 XOR)
% XOR 행의 고정 위치는 D+j (j=1..H)입니다.
for j = 1:H
    rows(D + j, :) = bitxor(rows(j, :), rows(j + H, :));
end
assert(isequal(rows(D+1:end, :), bitxor(rows(1:H, :), rows(H+1:D, :))), 'Block 3 XOR 검증 실패');
fprintf('Block 3 PASS | XOR parity rows = %d\n', H);

disp('첫 XOR 행 (HEX)'); fprintf('%02X ',rows(D+1,:)); fprintf('\n');

%% 4. 고정 index → inner RS(27,19)
% index 1..N을 2 byte big-endian으로 표현합니다.
indexBytesByRow = zeros(N, 2, 'uint8');
messageByRow = zeros(N, 19, 'uint8');
for index = 1:N
    % 높은 byte 먼저 저장: index 3은 [0, 3], decoder가 2 byte를 다시 읽습니다.
    indexBytesByRow(index, :) = uint8([floor(index / 256), mod(index, 256)]);
    messageByRow(index, :) = [rows(index, :) indexBytesByRow(index, :)];
end
assert(isequal(messageByRow(:, 18:19), indexBytesByRow), 'Block 4 index 검증 실패');
fprintf('index 부착: %d rows × 19 bytes\n', N);

% GF(256) 표: primitive polynomial 0x11D, alpha=2.
% 로그 표에서 지수를 더하면 GF 곱셈이 됩니다. 0은 별도로 처리합니다.
powers = zeros(1,512); logs = zeros(1,256); x = 1;
for k = 0:254
    powers(k+1) = x; logs(x+1) = k;
    x = x*2;
    if x >= 256, x = bitxor(x,285); end
end
for k = 255:511, powers(k+1) = powers(k-254); end
% 생성다항식 g(x) = product(x + alpha^r), r=1..8.
generator = 1;
for r = 1:8
    product = zeros(size(generator)); nz = generator ~= 0;
    product(nz) = powers(logs(generator(nz)+1) + logs(powers(r+1)+1) + 1);
    generator = bitxor([generator 0], [0 product]);
end
% 각 19 byte message 뒤에 나눗셈 나머지 8 byte를 붙입니다.
codewords = zeros(N,27,'uint8');
for index = 1:N
    message = double(messageByRow(index,:));
    work = [message zeros(1,8)];
    for k = 1:19
        coefficient = work(k);
        if coefficient ~= 0
            g = generator(2:end); product = zeros(1,8); nz = g ~= 0;
            product(nz) = powers(logs(g(nz)+1)+logs(coefficient+1)+1);
            work(k+(1:8)) = bitxor(work(k+(1:8)),product);
        end
    end
    codewords(index,:) = uint8([message work(20:27)]);
end
% 첫 codeword의 syndrome 8개가 모두 0인지 확인합니다.
syndrome = zeros(1,8);
for r = 1:8
    y = 0;
    for c = double(codewords(1,:))
        if y ~= 0, y = powers(logs(y+1)+r+1); end
        y = bitxor(y,c);
    end
    syndrome(r) = y;
end
assert(isequal(codewords(:,1:19), messageByRow) && ~any(syndrome));
fprintf('Block 4 PASS | %d codewords × 27 bytes, first syndrome = ',N);
fprintf('%d ',syndrome); fprintf('\n');
disp('첫 codeword (HEX)'); fprintf('%02X ',codewords(1,:)); fprintf('\n');

%% 5. byte → DNA → primer → FASTA / sequence TXT 저장
% F와 R은 1단계 설정에서 이미 선언한 프라이머를 그대로 사용합니다.
% F = 'AGCCTTGTGTCCATCAATCC';  % 1단계에서 선언한 forward primer (20 nt)
% R = 'TGCGCTATGGTTTGGCTAAT';  % 1단계에서 선언한 reverse primer (20 nt)
% 위 두 줄은 설명용 주석입니다. 1~4단계를 실행했다면 다시 선언할 필요가 없습니다.
% 아래 [F bodyDNA R]은 F를 108 nt 본문 앞에, R을 뒤에 붙입니다.
% 여기서 R은 가닥 끝에 넣을 서열로 정해져 있어 다시 역상보 변환하지 않습니다.
% 전체 길이: F 20 nt + 본문 108 nt + R 20 nt = 148 nt.
% 00=A, 01=T, 10=G, 11=C. 한 byte를 높은 bit부터 4 nt로 읽습니다.
assert(numel(F)==20 && numel(R)==20 && all(ismember([F R],'ATGC')));
BASES = 'ATGC'; seqs = cell(1,N);
for index = 1:N
    bytes = double(codewords(index,:)); v = zeros(4,27);
    for k = 1:4, v(k,:) = mod(floor(bytes/2^(8-2*k)),4); end
    bodyDNA = BASES(reshape(v,1,[])+1);
    seqs{index} = [F bodyDNA R];
end
assert(all(cellfun(@numel,seqs)==148));
fid = fopen(outFasta,'wt');
if fid < 0, error('FASTA 출력 파일을 열 수 없습니다.'); end
for index = 1:N
    % 복호 정보는 FASTA 이름 줄에 자동 기록합니다. DNA 서열에는 넣지 않습니다.
    if index == 1
        fprintf(fid,'>strand_%05d L=%d D=%d ext=%s\n%s\n',index,L,D,ext,seqs{index});
    else
        fprintf(fid,'>strand_%05d\n%s\n',index,seqs{index});
    end
end
fclose(fid);
fid = fopen(outSequencesTxt,'wt');
if fid < 0, error('TXT 출력 파일을 열 수 없습니다.'); end
for index = 1:N, fprintf(fid,'%s\n',seqs{index}); end
fclose(fid);
fprintf('Block 5 PASS | %d strands × 148 nt\n',N);
fprintf('첫 sequence: %s\n',seqs{1});
fprintf('저장: %s / %s\n',outFasta,outSequencesTxt);
fprintf('디코더 설정: L=%d; D=%d; ext=''%s'';\n',L,D,ext);
% 정확한 byte 비교용 원본을 보관합니다. DNA에 헤더를 넣는 것은 아닙니다.
save('original_bytes.mat','dataBytes','L','D','ext');

%% Decoding.m — 1~5단계 인코더의 파일 복호 예제
% 5단계 인코더의 FASTA / sequence TXT 또는 시퀀싱 FASTQ를 입력합니다.
% 인코더 실행창의 L, D 값을 아래에 그대로 입력해야 원본 파일 길이를 정확히 복원합니다.
clear; clc;

%% 0. 설정 — 인코더 출력값 입력
% 결과 파일은 현재 폴더(MATLAB Drive의 이 스크립트 폴더)에 저장됩니다.
% [바꾸는 곳] 인코더 로그의 L, D, ext를 아래에 옮깁니다.
% 명령창에서 미리 설정해도 clear가 지우므로 이 파일 내부를 수정하세요.
inFile = 'dnas_sequences.txt'; % 인코더 TXT (한 줄에 서열 하나) / FASTA / FASTQ
L = 27;                        % 인코더가 출력한 원본 byte 길이
D = 2;                         % 인코더가 출력한 짝수 data row 수
ext = 'txt';                   % 결과 파일 확장자: txt, jpg, png, ...
outFile = ['decoded.' ext];
F = 'AGCCTTGTGTCCATCAATCC';
R = 'TGCGCTATGGTTTGGCTAAT';
H = D / 2; N = D + H; BASES = 'ATGC';
assert(D >= 2 && mod(D, 2) == 0 && N <= 65535, 'D는 양의 짝수여야 합니다.');
assert(L >= 0 && L <= 17 * D, 'L은 0..17*D 범위여야 합니다.');

%% 1. FASTA / FASTQ / plain sequence 읽기
seqs = read_sequences(inFile);
if isempty(seqs), error('입력 파일에서 서열을 찾을 수 없습니다.'); end
fprintf('읽은 sequence/read: %d\n', numel(seqs));

%% 2. 방향 정리 → primer 확인 → 고정 index로 묶기
groups = repmat({char(zeros(0, 108))}, 1, N);
nRev = 0; nNoP = 0; nBad = 0; nIdx = 0;
for k = 1:numel(seqs)
    s = upper(regexprep(seqs{k}, '\s', '')); s(s == 'U') = 'T';
    % N/갭 등 불명 염기는 byte 변환에서 A로 오인하지 않도록 해당 read를 제외합니다.
    if any(~ismember(s, 'ATGC')), nBad = nBad + 1; continue; end
    if numel(s) < 148 || sum(s(1:20) == F) < 18
        s = reverse_complement(s);
        if numel(s) < 148 || sum(s(1:20) == F) < 18
            nNoP = nNoP + 1; continue;
        end
        nRev = nRev + 1;
    end
    if numel(s) < 148 || sum(s(129:148) == R) < 18
        nBad = nBad + 1; continue;
    end
    body = s(21:128);                         % 108 nt = 27 byte RS codeword
    ib = dna_to_bytes(body(69:76));           % RS message byte 18–19: fixed index
    index = double(ib(1))*256 + double(ib(2));
    if index < 1 || index > N
        nIdx = nIdx + 1; continue;
    end
    groups{index}(end + 1, :) = body;
end

%% 3. 묶음별 consensus → Global RS → read별 Local RS
% 같은 index의 read에서 위치별 최빈 염기를 고른 것이 consensus입니다.
% Global은 consensus 한 가닥에 RS를 적용하는 단계이며 별도 outer RS가 아닙니다.
% 실패하면 각 read를 따로 RS 복호하고 성공한 byte들의 최빈값을 사용합니다.
% RS(27,19)는 8 parity byte를 이용해 최대 4개의 byte 오류를 정정합니다.
rows = zeros(N, 17, 'uint8');
have = false(1, N); how = repmat({'missing'}, 1, N);
for index = 1:N
    G = groups{index};
    if isempty(G), continue; end
    counts = zeros(4, 108);
    for q = 1:4, counts(q, :) = sum(G == BASES(q), 1); end
    [~, best] = max(counts, [], 1);              % 동률은 ATGC 순서
    [msg, ok] = rs_decode(dna_to_bytes(BASES(best)));
    if ok && message_index(msg) == index
        rows(index, :) = msg(1:17); have(index) = true; how{index} = 'global';
        continue;
    end
    localRows = zeros(0, 17, 'uint8');
    for q = 1:size(G, 1)
        [msg, ok] = rs_decode(dna_to_bytes(G(q, :)));
        if ok && message_index(msg) == index
            localRows(end + 1, :) = msg(1:17); %#ok<AGROW>
        end
    end
    if ~isempty(localRows)
        rows(index, :) = uint8(mode(double(localRows), 1));
        have(index) = true; how{index} = 'local';
    end
end

%% 4. 바깥 XOR로 (j, j+H, D+j) 중 하나만 빠진 행 복원
% A XOR B = P이므로 A XOR P = B, B XOR P = A입니다.
% 세 행 중 두 행이 있어야 빠진 하나를 복원할 수 있습니다.
nXor = 0;
for j = 1:H
    triple = [j, j + H, D + j];
    missing = triple(~have(triple));
    if numel(missing) == 1
        known = triple(have(triple));
        rows(missing, :) = bitxor(rows(known(1), :), rows(known(2), :));
        have(missing) = true; how{missing} = 'xor'; nXor = nXor + 1;
    end
end

%% 5. 원래 데이터 행 연결 → 입력 byte 길이 L만큼 파일로 복원
missingData = find(~have(1:D));
if ~isempty(missingData)
    warning('복원되지 않은 data row가 있습니다: %s. 해당 영역은 0 byte로 기록됩니다.', ...
        mat2str(missingData));
end
% XOR 행을 제외한 데이터 D행을 원래 순서로 연결합니다.
% PAD 값과 같은 원본 byte도 있을 수 있으므로 L로 잘라야 정확합니다.
allData = reshape(rows(1:D, :)', 1, []);
data = allData(1:L);
fid = fopen(outFile, 'wb');
if fid < 0, error('출력 파일을 만들 수 없습니다: %s', outFile); end
fwrite(fid, data, 'uint8'); fclose(fid);

%% 6. 검증 보고: 서열 수 · 방향 · primer · RS · XOR · 파일 byte
fprintf('read %d | reverse %d | no_primer %d | bad_length %d | bad_index %d\n', ...
    numel(seqs), nRev, nNoP, nBad, nIdx);
fprintf('strand %d | Global RS %d | Local RS %d | XOR 복원 %d\n', N, ...
    sum(strcmp(how, 'global')), sum(strcmp(how, 'local')), nXor);
if isempty(missingData), fprintf('데이터 행 전부 확보\n');
else, fprintf('누락 data row: %s\n', mat2str(missingData)); end
fprintf('복원 파일: %s (%d byte)\n', outFile, numel(data));
if strcmpi(ext, 'txt')
    fprintf('UTF-8 text: %s\n', native2unicode(data, 'UTF-8'));
elseif ismember(lower(ext), {'jpg','jpeg','png','tif','tiff','bmp'})
    try, figure('Name', outFile); image(imread(outFile)); axis image off;
    catch err, warning('이미지 표시 실패: %s', err.message); end
end

%% 7. 같은 폴더에 원본 byte 기록이 있으면 정확한 일치 확인
% 5단계에서 저장한 original_bytes.mat는 실습 검증용이며 복호 필수 입력은 아닙니다.
% 다른 데이터를 새로 인코딩하면 이 파일도 덮어써집니다.
if isfile('original_bytes.mat')
    original = load('original_bytes.mat','dataBytes');
    if isequal(original.dataBytes,data)
        fprintf('ROUNDTRIP PASS | %d bytes identical\n',numel(data));
    else
        warning('원본 byte와 복원 결과가 다릅니다. 설정값과 누락 행을 확인하세요.');
    end
end

%% 로컬 함수 — FASTA/FASTQ와 RS 복호를 이 파일 하나에서 처리
function seqs = read_sequences(filename)
% FASTA, 4-line FASTQ, plain TXT를 판별해 각 read를 한 줄 문자열로 반환합니다.
t = strrep(fileread(filename), sprintf('\r'), '');
lines = strsplit(t, sprintf('\n')); lines = cellfun(@strtrim, lines, 'UniformOutput', false);
lines = lines(~cellfun(@isempty, lines)); seqs = {};
if isempty(lines), return; end
if lines{1}(1) == '>'                         % FASTA: wrapped sequence 지원
    current = '';
    for k = 1:numel(lines)
        if lines{k}(1) == '>'
            if ~isempty(current), seqs{end + 1} = current; end %#ok<AGROW>
            current = '';
        else
            current = [current lines{k}]; %#ok<AGROW>
        end
    end
    if ~isempty(current), seqs{end + 1} = current; end
elseif lines{1}(1) == '@'                     % FASTQ: wrapped sequence/quality 지원
    k = 1;
    while k <= numel(lines)
        if lines{k}(1) ~= '@', error('FASTQ %d번째 줄이 @ header가 아닙니다.', k); end
        k = k + 1; sequence = '';
        while k <= numel(lines) && lines{k}(1) ~= '+'
            sequence = [sequence lines{k}]; %#ok<AGROW>
            k = k + 1;
        end
        if k > numel(lines), error('FASTQ에서 + 구분자를 찾을 수 없습니다.'); end
        k = k + 1; qualityLength = 0;
        while k <= numel(lines) && qualityLength < numel(sequence)
            qualityLength = qualityLength + numel(lines{k}); k = k + 1;
        end
        if qualityLength ~= numel(sequence), error('FASTQ sequence/quality 길이가 다릅니다.'); end
        seqs{end + 1} = sequence; %#ok<AGROW>
    end
else
    seqs = lines;                               % 한 줄당 sequence 하나
end
end

function s = reverse_complement(s)
% 염기 A↔T, C↔G 치환 후 뒤집어 역상보 서열을 만듭니다.
s = upper(s); s(s == 'A') = 't'; s(s == 'T') = 'a';
s(s == 'G') = 'c'; s(s == 'C') = 'g'; s = upper(fliplr(s));
end

function b = dna_to_bytes(s)
% ATGC를 00/01/10/11로 되돌리고 4 nt마다 1 byte로 묶습니다.
s = upper(char(s));
if mod(numel(s), 4) ~= 0 || any(~ismember(s, 'ATGC'))
    error('DNA-to-byte input must contain only ATGC and have a length divisible by 4.');
end
v = zeros(1, numel(s)); v(s == 'T') = 1; v(s == 'G') = 2; v(s == 'C') = 3;
v = reshape(v, 4, []); b = uint8([64 16 4 1] * v);
end

function index = message_index(msg)
% RS가 복구한 message의 마지막 두 byte(big-endian)를 strand index로 읽습니다.
index = double(msg(18)) * 256 + double(msg(19));
end

function [msg, ok] = rs_decode(codeword)
% RS(27,19), GF(256)/0x11D, roots alpha^1..alpha^8; BM + Chien + Forney.
nsym = 8; msg = uint8([]); ok = false;
r = double(reshape(uint8(codeword), 1, [])); n = numel(r);
if n ~= 27, return; end
S = syndromes(r, nsym);
if ~any(S), msg = uint8(r(1:n - nsym)); ok = true; return; end
C = [1 zeros(1, nsym)]; Bp = C; degree = 0; shift = 1; lastDiscrepancy = 1;
for step = 0:nsym - 1
    discrepancy = S(step + 1);
    for i = 1:degree
        discrepancy = bitxor(discrepancy, gf_mul(C(i + 1), S(step - i + 1)));
    end
    if discrepancy == 0, shift = shift + 1; continue; end
    oldC = C; scale = gf_div(discrepancy, lastDiscrepancy);
    for i = 0:nsym - shift
        C(i + shift + 1) = bitxor(C(i + shift + 1), gf_mul(scale, Bp(i + 1)));
    end
    if 2 * degree <= step
        degree = step + 1 - degree; Bp = oldC;
        lastDiscrepancy = discrepancy; shift = 1;
    else
        shift = shift + 1;
    end
end
locator = C(1:degree + 1);
if degree == 0 || degree > nsym / 2, return; end
[E, ~] = gf_table(); positions = []; roots = [];
for p = 1:n
    x = E(mod(-(n - p), 255) + 1);
    if poly_eval_ascending(locator, x) == 0
        positions(end + 1) = p; roots(end + 1) = x; %#ok<AGROW>
    end
end
if numel(positions) ~= degree, return; end
omega = zeros(1, nsym);
for i = 0:degree
    for j = 0:nsym - i - 1
        omega(i + j + 1) = bitxor(omega(i + j + 1), gf_mul(locator(i + 1), S(j + 1)));
    end
end
derivative = zeros(1, degree);
for i = 1:degree
    if mod(i, 2), derivative(i) = locator(i + 1); end
end
for k = 1:numel(positions)
    denominator = poly_eval_ascending(derivative, roots(k));
    if denominator == 0, return; end
    magnitude = gf_div(poly_eval_ascending(omega, roots(k)), denominator);
    r(positions(k)) = bitxor(r(positions(k)), magnitude);
end
if any(syndromes(r, nsym)), return; end
msg = uint8(r(1:n - nsym)); ok = true;
end

function S = syndromes(r, nsym)
% 수신 codeword를 생성다항식의 근에 대입해 RS syndrome을 계산합니다.
[E, ~] = gf_table(); S = zeros(1, nsym);
for j = 1:nsym
    y = 0; x = E(j + 1);
    for c = r, y = bitxor(gf_mul(y, x), c); end
    S(j) = y;
end
end

function y = poly_eval_ascending(p, x)
% 오름차순 계수 다항식을 GF(256) Horner 방식으로 평가합니다.
y = 0;
for i = numel(p):-1:1, y = bitxor(gf_mul(y, x), p(i)); end
end

function c = gf_mul(a, b)
% GF(256) 곱셈: log/antilog 표로 계산합니다.
if a == 0 || b == 0, c = 0; return; end
[E, Lg] = gf_table(); c = E(Lg(a + 1) + Lg(b + 1) + 1);
end

function c = gf_div(a, b)
% GF(256) 나눗셈: 지수 차이를 255로 순환시킵니다.
if b == 0, error('GF(256): zero divisor'); end
if a == 0, c = 0; return; end
[E, Lg] = gf_table(); c = E(mod(Lg(a + 1) - Lg(b + 1), 255) + 1);
end

function [E, Lg] = gf_table()
% primitive polynomial 0x11D를 이용해 GF(256) lookup table을 한 번 생성합니다.
persistent e l
if isempty(e)
    e = zeros(1, 512); l = zeros(1, 256); x = 1;
    for p = 0:254
        e(p + 1) = x; l(x + 1) = p;
        x = x * 2; if x >= 256, x = bitxor(x, 285); end
    end
    for p = 255:511, e(p + 1) = e(p - 254); end
end
E = e; Lg = l;
end