function [ch,cind,onert,sirt,lj,norder] = charac(multab,G,ncl,...
h,nfirst,classl,inel,steer,prime)
%
% The irreducible characters are calculated according to Dixon's method.
%
% Section 1
%
% The order, norder(I), and the powers, npow(I,K), of the elements in each
% class are found. That is, the index of the class to which groupelement 
% Q^K belongs, if Q belongs to class I. nconj(I) is the index of the class
% of inverse elements of class I. The array nar is used several times for 
% different variables.
%
fid = fopen('Irrepres.m', 'a');
pi2 = pi*2;
sqr = zeros(1,520);
nconj(1) = 1;
if G > 1
  for I = 1:ncl
    K = nfirst(I);
    L = h(I);
    M = K + L - 1;
    for I1 = K:M
      N = classl(I1);
      cind(N) = I;
    end
  end
  fprintf(fid, 'The class index for each group element is:\n');
  fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', cind);
  fprintf(fid, '\n');
%
% cind(N) = I if group element N belongs to class I
%
  for I = 1:ncl
    onert(I) = 0;
    K = nfirst(I);
    K = classl(K);
    J = K;
    norder(I) = 1;
    npow(I,1) = I;
    N = 1;
    if J ~=1
      while J ~=1
        N = N + 1;
        J = multab(J,K);
        if J ~=1
          npow(I,N) = cind(J);
        end
      end
      npow(I,N) = 1;
      norder(I) = N;
      M = inel(K);
      nconj(I) = cind(M);
    end
  end
  fprintf(fid, 'The order of each class is:\n');
  fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', norder);
  fprintf(fid, '\n');
  fprintf(fid, 'The powers of each class is:\n');
  fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', npow);
  fprintf(fid, '\n');
  fprintf(fid, 'The conjugated class of each class is:\n');
  fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', nconj);
  fprintf(fid, '\n');
%
% Section 2
%
% Find the exponent ex, the least common multiple of norder(I), I = 1:ncl.
%
  ex = 1;
  for K = 2:ncl
    N = norder(K);
    M = ex;
%
% Determine the greatest common divisor gcd of N and M. Then the least
% common multiple is N*M/gcd.
%
    ex = N*M/gcd(N,M);
  end
  fprintf(fid, '\n');
  fprintf(fid, 'The exponent of the group is: %3d\n\n', ex);
%
% Section 3
%
% Determine prime P, such that ex divides P - 1.
%
  maxh = fix(2*sqrt(G)) + 1;
  for I = 1:100
    NP1 = prime(I) - 1;
    if (prime(I) >= maxh) & ((fix(NP1/ex)*ex - NP1) == 0)
      P = prime(I);
      break;
    end
  end
  fprintf(fid, 'The prime P, such that ex divides P - 1 is: %3d\n\n', P); 
%
% Section 4
%
% Determine a primitive root Z, such that Z^ex == 1 (modulus P), but 
% Z^I ~= 1 (mod P) for I < ex.
% Determine indices ind(I), such that ind(I) == K if Z^K == I (mod P)
% and powers nup(K), such that nup(K) == I if Z^K == I (mod P).
% Inverses inv(I), such that inv(I)*I == 1 (mod P)
%
  M = P - 1;
  for root = 2:M
    I = 1;
    J = root;
    while J ~= 1
      I = I + 1;
      J = modulus(J*root,P);
    end
    if I == ex
      break;
    end
  end
  if root > M
    fprintf(fid, 'No primitive root\n\n');
    return;
  end
  Zprim = root;
  J = Zprim;
  I = 1;
  while I <= ex
    ind(J) = I;
    nup(I) = J;
    K = modulus(I*I,P);
    if sqr(K) == 0
      sqr(K) = I;
    end
    I = I + 1;
    J = modulus(J*Zprim,P);
  end
  inv(1) = 1;
  for I = 2:M
    for L = 2:M
      if modulus(I*L,P) == 1;
        inv(I) = L;
        break;
      end
    end
  end
  if steer(6)
    fprintf(fid, 'The primitive root of the group is:%3d\n\n', Zprim);
    fprintf(fid, 'The indices of the group are:\n');
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', ind);
    fprintf(fid, '\n');
    fprintf(fid, 'The powers of the group are:\n');
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', nup);
    fprintf(fid, '\n');
    fprintf(fid, 'The inverses modulus P are:\n');
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', inv);
    fprintf(fid, '\n');
    fprintf(fid, 'The square roots modulus P are:\n');
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', sqr);
    fprintf(fid, '\n\n');
  end
%
% Section 5
%
% The eigenspaces of the first classconstant matrix are formed and the
% eigenspaces of dimension higher than 1 are reduced by operation with
% other classconstant matrices.
% The eigenvectors are stored as rows, one should multiply with the
% transposed classconstant matrices.
% The eigenvectors lie in w(I,J).
% If the first row of an eigenspace is row S, then the last row of this
% eigenspace is row nar(S).
% We start with the space of the unit class; all operations modulus P.
%
  w = eye(ncl);
  nar = zeros(1,ncl);
  nar(1) = ncl;
%
% Form classconstant matrix c(LC,I,L) and store it in mat(I,L)
%
  LC = 2;
  while LC <= ncl
    mat = zeros(ncl);
    nflc = nfirst(LC);
    nllc = nflc + h(LC) - 1;
    for I = 1:ncl
      nfi = nfirst(I);
      nli = nfi + h(I) - 1;
      for MJ = nflc:nllc
        NJ = classl(MJ);
        for IH = nfi:nli
          NI = classl(IH);
          L = multab(NJ,NI);
          for K = 1:ncl
            nfk = nfirst(K);
            nfe = classl(nfk);
            if L == nfe
              mat(I,K) = mat(I,K) + 1;
              break;
            end
          end
        end
      end
    end
    S = 1;
    ncl2 = 2*ncl;
    while nar(S) ~= 0
      while nar(S) ~= S
        F = S;
        T = nar(S);
        for I = S:T
          for J = 1:ncl
            u(I,J) = 0;
%
% The result of the operation of the classconstant matrix on the 
% eigenspace is stored in u(I,J).
%
            for K = 1:ncl
              u(I,J) = u(I,J) + w(I,K)*mat(J,K);
%
% Note transpose of mat.
%
            end
            u(I,J) = modulus(u(I,J),P);
          end
        end
        nh = 1;
        while nh <= P
          e = nh - 1;
          for I = S:T
            for J = 1:ncl
              K = J + ncl;
              v(I,J) = w(I,J);
              v(I,K) = modulus(u(I,J) + e*w(I,J),P);
            end
          end
%
% If -e is an eigenvalue of mat(I,J), then the right part of v(I,K) now
% contains the sum of eigenspaces minus the vectors belonging to this
% eigenvalue. Then there are linear combinations of rows, which vanish.
% The corresponding linear combinations of w(I,J) (stored in the left
% part of v(I,K)), will be common eigenvectors of the present
% classconstant matrix and all preceding ones.
%
          J = F;
          for  I = S:T
            K = ncl + 1;
            while K <= ncl2
              if v(I,K) ~= 0
                Y = v(I,K);
                invs = inv(Y);
                for L = 1:ncl2
                  v(I,L) = modulus(v(I,L)*invs,P);
                end
%
% This makes v(I,K) == 1
%
                for row = (I+1):T
                  Y = modulus(P - v(row,K),P);
                  for L = 1:ncl2
                    v(row,L) = modulus(v(row,L) + v(I,L)*Y,P);
                  end
                end
                break;
              end
%
% This makes v(row,K) == 0
%
              K = K + 1;
            end
            if K > ncl2
              z(F,1:ncl) = v(I,1:ncl);
% 
% z(I,K) is a temporary storage of produced common eigenvectors
%
              F = F + 1;
            end
          end
          nar(J) = F - 1;
          if F > T
            break;
          end
          nh = nh + 1;
        end
        if F > T
          break;
        end
        if nh > P
          fprintf(fid, '\n');
          fprintf(fid, 'Not enough eigencolumns\n\n');
          return;
        end
      end
      w(S:T,1:ncl) = z(S:T,1:ncl);
      if nar(S) ~= ncl
        S = nar(S) + 1;
      else
        break;  
      end 
    end
    if nar(S) == 0
      fprintf(fid, '\n\');
      fprintf(fid, 'Error, nar(S) == 0\n');
      return;
    end
    Y = LC;
    for I = 1:ncl
      if nar(I) ~= I
        LC = LC + 1;
        break;
      end
    end
    if Y == LC
      break;
    end
  end 
%
% Section 6
% 
% Calculation of the characters. The dimensions are stored in lj(K).
% The coefficients of the powers of the primitive roots are stored
% in v(I,J).
%
  for I = 1:ncl
    if w(I,1) ~= 0
      L = inv(w(I,1));
      for K = 1:ncl
        w(I,K) = modulus(w(I,K)*L,P);
      end
    else
      fprintf(fid, 'Error, w(I,1) == 0\n');
      return;
    end
  end
  for I = 1:ncl
    J = h(I);
    hinv(I) = inv(J);
    K = norder(I);
    nar(I) = ex/K;
    nordin(I) = inv(K);
  end
  for K = 1:ncl
    norm2 = 0;
    for I = 1:ncl
      J = nconj(I);
      X = modulus(w(K,I)*w(K,J),P);
      X = modulus(X*hinv(I),P);
      norm2 = modulus(norm2 + X,P);
    end
    deg2 = modulus(G*inv(norm2),P);
%
% deg2 is the square of the dimension of the irreducible representation.
%
    lj(K) = sqr(deg2);
    for N = 1:ncl
      u(K,N) = modulus(lj(K)*w(K,N)*hinv(N),P);
    end
%
% u(K,I) contains the irreducible characters modulus P.
%
    for I1 = 1:ncl
      M = norder(I1);
      for J = 1:M
        J1 = J - 1;
        X = J1*nar(I1);
        Y = 0;
        for L2 = 1:M
          L = L2 - 1;
          S = modulus(-L*X,ex);
          S = modulus(S + ex,ex);
          if L == 0
            L = M;
          end
          T = npow(I1,L);
          if S == 0
            Y = modulus(Y + u(K,T),P);
          else 
            Y = modulus(Y + u(K,T)*nup(S),P);
          end
        end
        v(I1,J) = modulus(nordin(I1)*Y,P);
      end
      if steer(2)
%
% A nondegenerate eigenvalue or an eigenvalue with the lowest possible 
% degeneracy, is registered.
%
        if onert(K) ~= 1
          sirtk3 = 1;
          if v(I1,1) ~= 1
            degen = v(I1,1);
            if v(I1,1) == 0
              degen = 10;
            end
            J2 = 2;
            while J2 <= M
              if v(I1,J2) < degen
                if v(I1,J2) ~= 0 
                  degen = v(I1,J2);
                  sirtk3 = J2;
                  if v(I1,J2) == 1
                    break;
                  end
                end
              end
              J2 = J2 + 1;
            end
          end
          if v(I1,1) == 1 | J2 <= M
            onert(K) = 1;
            sirt(K,1) = I1;
            sirt(K,2) = M;
            sirt(K,3) = sirtk3;
          end
          if v(I1,1) ~= 1 & J2 > M
            if I1 == 1 | degen < minimum(K)
              sirt(K,1) = I1;
              sirt(K,2) = M;
              sirt(K,3) = sirtk3;
              minimum(K) = degen;
            end
          end
        end
      end
      if steer(6)
        fprintf(fid,...
'Character (%3d, %3d) in %3d th roots of unity is: %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', K, I1, M, v(I1,1:M));   
        fprintf(fid, '\n\n');
      end
      if I1 ~= 1
        ch(K,I1) = v(I1,1);
%
% The zeroth power of Z is 1, so the first term in the addition is
% 1*v(I1,1)
%
        for L1 = 2:M
          if v(I1,L1) ~= 0
            croot = pi2*(L1 - 1)*i/M;
            croot = exp(croot);
            ch(K,I1) = ch(K,I1) + croot*v(I1,L1);        
          end
        end
      else
        ch(K,1) = lj(K);
      end
      if abs(imag(ch(K,I1))) < 0.0001
        ch(K,I1) = real(ch(K,I1));
      end
      if abs(ch(K,I1)) < 0.0001
        ch(K,I1) = 0;
      end
    end
  end
%
% Section 7
%
% Check the results. Output
%
  L = sum(lj(1:ncl).^2);
  if L ~= G
    fprintf(fid, 'Error, L ~= G\n');
    fprintf(fid, 'L = %3d G = %3d lj = %3d\n', L, G, lj);
    return;
  end 
  if steer(7)
    fprintf(fid, 'lj = \n\n');  
    fprintf(fid, '%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', lj);
    fprintf(fid, '\n\n');
    for J = 1:ncl
      for JJ = 1:ncl
        chord(1,J,JJ) = real(ch(J,JJ));
        chord(2,J,JJ) = imag(ch(J,JJ));
      end
    end
    for J = 1:ncl
      fprintf(fid, 'ch( %3dL ) =  \n', J);
      fprintf(fid, '(%9.5f %9.5f)  (%9.5f %9.5f)  (%9.5f %9.5f)\n', chord(1:2,J,1:ncl));
      fprintf(fid, '\n\n');
    end
    fprintf(fid, '\n');
  end       
else
  lj(1) = 1;
  ch(1,1) = 1;
  cind = 1;
  onert = 1;
  sirt = 1;
  norder = 1;
  fprintf(fid, 'lj = %3d\n\n',lj);
  fprintf(fid, 'ch( 1,1 ) = %7.4f\n', ch(1,1));
end
fclose(fid);
return