% function g = generator(n)
%
% Find a generator of the cyclic group U_n = { a in Z_n : gcd(a, n) = 1 } with
% n=2 or n=4 or n of the form p^k or 2*p^k, with p an odd prime.
% If n is not of the above form then U_n is not cyclic and the group is not
% generated by a single element.
%
% (C) 2004 <dirk.nuyens@cs.kuleuven.ac.be>
function g = generator(n)

f = factor(n);
uf = unique(f); % the unique factors of $n$, in increasing order
nf = zeros(size(uf)); % in \texttt{nf} we will put the powers of each unique factor
for i=1:length(nf), nf(i) = length(find(f == uf(i))); end;
% we now have that $n$ = \verb+prod(uf .^ nf)+

if (length(uf) == 2) && (uf(1) == 2) && ...
   (nf(1) == 1) && (mod(uf(2), 2) == 1)
    % $n$ is of the form $2 p^k$, with $p$ an odd prime
    p = uf(2); k = nf(2);
    g = generatorp(p);
    if mod(g, 2) == 0
        g = g + p^k;
    end;
elseif (length(uf) == 1) && (mod(uf(1), 2) == 1)
    % $n$ is of the form $p^k$, with $p$ an odd prime
    p = uf(1); k = nf(1);
    g = generatorp(p);
    if k > 1
        if powmod(g, p-1, p^2) == 1
            % This happens first for the prime 40487,
            % that is for $n = 40487^2 = 1639197169$.
            % The \texttt{powmod.m} routine will then incorrectly
            %  calculate $5^{40487-1} \bmod{1639197169}$ unless using
            % the arbitrary precision integers from Java.
            g = g + p;
        end;
    end;
elseif n == 2
    g = 1;
elseif n == 4
    g = 3;
else
    error(['n must be 2 or 4 or of the form p^k or 2*p^k, ' ...
            'with p an odd prime']);
end;