Essentially, this is more about finite Fourier transforms (although you could argue that this counts as representation theory of finite abelian groups).

Our matrix fulfils the identity A^k = 1. Therefore all eigenvalues must fulfill this equation as well.

We conclude the eigenvalues are of the form exp(2pi*I*t/k) where t goes from 0 to k-1 and I is the imaginary unit. Say the eigenvalue exp(2pi*I*t/k) has multiplicity a_t.

The trace of a matrix is the sum of all its eigenvalues, and taking powers of a matrix has the effect that the eigenvalues are taken to this power as well.

So, A^j has the eigenvalue exp(2pi*I*j*t/k) with multiplicity a_t and its trace is:

tr(A^j) = sum (t = 0 to k-1) a_t * exp(2pi*I*j*t/k)

This is actually just the discrete Fourier transform of the function taking t to a_t. Therefore, we want the inverse Fourier transform of the function taking j to tr(A^j). Call the latter f_j.

The inverse transform is given by:

a_t = 1/k * sum (t = 0 to k-1) f_j * exp(-2pi*I*j*t/k).

I know that a_t is an integer, so I just calculate the real part of this sum and do some rounding with this result. Furthermore, I modded out k from j*t - the result of exp(...) is the same but errors due to big arguments are less probable to occur.

Then I print the eigenvalue exp(2pi*I*t/k) simply a_j times (my code calculates the occurrences of exp(-2pi*I*t/k) but this just reverses the order of the eigenvalues).

By the way - the number n is not needed for the calculations.

Here is my code:

cout.precision(7);

    double x ;
    
    vector<double> real ;
    vector<double> imag ;
    
    double PI = 3.14159265358979 ;
    
    int k, n ;
    
    cin >> n >> k  ;
    
    double a,b ;
    
    for(int i = 0 ; i < k ; i++){
        cin >> a >> b ;
        real.push_back(a) ;
        imag.push_back(b) ;
    }
    
    for(int i = 0 ; i < k ; i++){
        double sumreal = 0 ;
        
        for(int j = 0 ; j < k ; j++){
            int modded = (i*j) % k ; //
            sumreal = sumreal + real[j]*cos(2*modded*PI/k) - imag[j]*sin(2*modded*PI/k) ;
        }
        
        sumreal = sumreal/k + 0.5 ;
        
        int howmany ;
        
        howmany = sumreal ;
        
        double x = cos(2*i*PI/k) ;
        double y = -sin(2*i*PI/k) ;
        
        for(int l = 0 ; l < howmany ; l++){
            cout << fixed << x << " " << y << endl ;
        }
        
    }