function [y,idx] = digitrevorder(x, radixbase)
%DIGITREVORDER Permute gpuArray input into digit-reversed order.
%   Y = DIGITREVORDER(X, R)
%   [Y,I] = DIGITREVORDER(X, R)
%
%   See also DIGITREVORDER, BITREVORDER, gpuArray.

% Copyright 2020-2023 The MathWorks, Inc.

% Make sure only X is on GPU, and if not use the standard function
radixbase = gather(radixbase);
if ~isgpuarray(x)
    [y,idx] = digitrevorder(x, radixbase);
    return;
end

% R must be a scalar, real, integer between 2 and 36:
if ~isnumeric(radixbase) ...
        || ((length(radixbase) > 1) ...
        || (~isreal(radixbase)) ...
        || (fix(radixbase) ~= radixbase) ...
        || ((radixbase < 2) || (radixbase > 36)))
    error(message('signal:digitrevorder:MustBePosInteger'));
end

% At this point, if radixbase is 2, use the slightly faster BITREVORDER.
if isequal(radixbase,2)
    if nargout>1
        [y,idx] = bitrevorder(x);
    else
        y = bitrevorder(x);
    end
    return;
end

% Since X is a gpuArray it is guaranteed to be numeric or logical. Just
% check it is 2D.
if ~ismatrix(x)
    error(message('signal:digitrevorder:InvalidDimensions'));
end

if isvector(x)
    N = length(x);
else
    N = size(x,1);
end

% Check to make sure N is an integer power of radixbase:
radixpow = round(log10(N) / log10(radixbase));
if ((radixbase ^ radixpow) ~= N)
    error(message('signal:digitrevorder:MustBeInteger', string(radixbase)));
end

% Max index size for GPU is 2^32-1 so use uint32 for bit manipulation
idx = uint32(gpuArray.colon(1, N))';

idx = arrayfun(@iReverseIdxBaseN, idx, radixbase, radixpow);
if isvector(x)
    y = x(idx);
else
    y = x(idx,:);
end

% Make sure the index output is double if requested
if nargout>1
    idx = double(idx);
end

end

function y = iReverseIdxBaseN(idx, radixbase, radixpow)
y = zeros("like",idx);
outvalue = cast(radixbase.^(radixpow-1),"like",idx);
idx = idx - 1; % 1-based index to 0-based

% Note that every GPU thread does the same number of loop steps so this
% loop is actually pretty efficient.
for bit=1:radixpow
    r = rem(idx,radixbase);
    y = y + r*outvalue;
    idx = (idx-r)/radixbase; % shift right
    outvalue = outvalue/radixbase; % shift right
end
y = y + 1; % 0-based index to 1-based
end
