% The nth roots of unity are solutions to the
% equation
%
%      n
%     z  = 1;
%
% that is, they are all complex numbers that
% equal 1 when raised to the power n.
%
% Equivalently, they are the roots of the
% polynomial
%
%      n
%     z  - 1.
%
% This is a polynomial of degree n, so the
% Fundamental Theorem of Algebra tells us that
% it has n complex roots, counting multiplicity.
%
% We will see that these n roots are all
% different. Thus, there are exactly n distinct
% nth roots of unity.


% This assigns 'n' the value 10.
%
% You may choose 'n' to be any positive integer.

n = 10


% This is how in MATLAB and Octave we can check
% that:
%
%   1. 'n' has a finite value;
%   2. 'n' has an integer value; and
%   3. 'n' is positive.
%
% The function 'assert(...)' causes an error if
% the condition passed to it is false.

assert(isfinite(n) && (n == round(n)))
assert(n > 0)


% Recall that the complex number 1 can be
% represented in polar form as
%
%     1 phase 0.
%
% However, arguments are not unique. Rotating
% through 2*pi radians, or 360 degrees, returns
% us to the same point.
%
% Therefore, all of the following represent the
% complex number 1:
%
%     1 phase 0
%     1 phase 2*pi
%     1 phase 4*pi
%     1 phase 6*pi
%           ...
%
% More generally,
%
%     1 = exp(1j*2*k*pi)
%
% for any integer k.


% Suppose that z is an nth root of unity. Write
% z in polar form as
%
%     z = r exp(1j*phi).
%
% Raising this to the power n gives
%
%      n    n
%     z  = r  exp(1j*n*phi).
%
% We require this to equal 1.
%
% Comparing magnitudes gives
%
%      n
%     r  = 1.
%
% Since a magnitude is nonnegative, this means
%
%     r = 1.
%
% Thus, every root of unity lies on the unit
% circle.


% Comparing arguments, we require
%
%     n*phi = 2*k*pi
%
% for some integer k.
%
% Dividing by n gives
%
%           2*k*pi
%     phi = -------
%              n
%
% or, in degrees,
%
%           360*k
%     phi = ----- degrees.
%             n
%
% Thus, the nth roots of unity are
%
%                    2*k*pi
%     exp( 1j * ------------- )
%                       n
%
% for
%
%     k = 0, 1, 2, ..., n - 1.


% Define
%
%             2*pi
%     theta = ----.
%               n
%
% This is the angle between consecutive nth
% roots of unity.

format long

theta = 2*pi/n


% In degrees, this angle is

theta_degrees = 360/n


% Thus, the roots are separated from one another
% by exactly
%
%     360/n degrees.
%
% For n = 10, for example, the roots are
% separated by 36 degrees.


% We can now calculate all n roots.
%
% This loop is executed n times. On successive
% executions, k takes the values
%
%     0, 1, 2, ..., n - 1.

fprintf('\nThe %d roots of unity:\n', n)

for k = 0:(n - 1)

    % The kth root is
    %
    %     cos(k*theta) + j*sin(k*theta).
    %
    % Equivalently, using the complex
    % exponential, it is
    %
    %     exp(1j*k*theta).

    z = cos(k*theta) + 1j*sin(k*theta)

    % Since z is an nth root of unity, raising
    % it to the power n should give 1.

    z^n

end


% There is an especially useful way to describe
% all of these roots.
%
% Define the complex number
%
%                    2*pi
%     omega = exp(1j*-----).
%                      n
%
% In degrees, this is
%
%                360
%     1 phase    --- degrees.
%                 n
%
% Thus, omega is the point on the unit circle
% obtained by rotating counterclockwise through
% 360/n degrees.

omega = exp(1j*2*pi/n)


% Since
%
%                   2*pi
%     omega = exp(1j-----),
%                     n
%
% raising omega to the kth power gives
%
%      k               2*k*pi
%     omega = exp(1j * -------).
%                         n
%
% Consequently, every nth root of unity is a
% power of omega.
%
% The n roots are therefore
%
%          0       1       2              n-1
%     omega,  omega,  omega,  ...,  omega.
%
% Since
%
%          0
%     omega = 1,
%
% we can also write them as
%
%                         2
%     1, omega, omega, ..., omega^(n - 1).


fprintf('\nThe roots as powers of omega:\n')

for k = 0:(n - 1)
    omega^k
end


% In our phase notation, we can summarize this
% particularly neatly.
%
% If
%
%                  360 degrees
%     omega = 1 phase -----------,
%                       n
%
% then the nth roots of unity are
%
%          0       1       2              n-1
%     omega,  omega,  omega,  ...,  omega.
%
% That is, they are
%
%                                k
%             /        360      \
%             | 1 phase --- deg |
%             \         n       /
%
% for
%
%     k = 0, 1, 2, ..., n - 1.


% Why do we stop at n - 1?
%
% Consider the next power:
%
%      n
%     omega.
%
% Its argument is
%
%         360
%     n * --- = 360 degrees,
%          n
%
% so
%
%      n
%     omega = 1.
%
% But this is exactly the root we started with:
%
%      n       0
%     omega = omega = 1.
%
% After this, the roots simply repeat:
%
%      n+1       1
%     omega   = omega,
%
%      n+2       2
%     omega   = omega,
%
% and so on.


% This also shows why the n roots
%
%          0       1              n-1
%     omega,  omega,  ...,  omega
%
% are all different.
%
% Their arguments are
%
%       360     2*360              (n - 1)*360
%     0, ---,   -----,  ...,       -----------
%        n        n                       n
%
% degrees.
%
% All of these angles are at least 0 degrees
% and less than 360 degrees, and no two are
% equal.
%
% Thus, we have found n distinct roots. By the
% Fundamental Theorem of Algebra, there cannot
% be any others.


% There is a beautiful geometric interpretation.
%
% Every nth root of unity has magnitude 1, so
% every root lies on the unit circle.
%
% Consecutive roots are separated by an angle
% of
%
%     360/n degrees.
%
% Therefore, the n roots are equally spaced
% around the unit circle.
%
% Connecting consecutive roots produces a
% regular n-sided polygon centered at the
% origin.
%
% For example:
%
%     n = 2    opposite points
%     n = 3    equilateral triangle
%     n = 4    square
%     n = 5    regular pentagon
%     n = 6    regular hexagon
%
% and so on.


% For example, the fourth roots of unity are
%
%     1 phase   0 degrees =  1
%     1 phase  90 degrees =  j
%     1 phase 180 degrees = -1
%     1 phase 270 degrees = -j.
%
% Thus,
%
%      4
%     z  = 1
%
% has the four roots
%
%     1, j, -1, -j.


% The nth roots of unity also have considerable
% symmetry.
%
% Since the roots are equally spaced around the
% unit circle, multiplying any root by omega
% rotates it counterclockwise to the next root:
%
%          k          k+1
%     omega * omega = omega.
%
% Multiplication by omega therefore cycles
% through all of the roots.


% Also recall that, on the unit circle, the
% multiplicative inverse and complex conjugate
% of a number are equal:
%
%     1/z = conj(z).
%
% Therefore, the conjugate of every root of
% unity is also a root of unity.
%
% In terms of powers of omega,
%
%                k        n-k
%     conj(omega ) = omega
%
% for k = 1, 2, ..., n - 1.


% Finally, because these are precisely the roots
% of
%
%      n
%     z  - 1,
%
% the polynomial can be factored completely as
%
%      n
%     z  - 1
%
%       = (z - 1)
%         (z - omega)
%                    2
%         (z - omega )
%             ...
%                        n-1
%         (z - omega     ).
%
% This is the Fundamental Theorem of Algebra
% made particularly visible: the degree-n
% polynomial has exactly n linear factors over
% the complex numbers.


%%%%%%%%%%%%%%%%%%% IMPORTANT %%%%%%%%%%%%%%%%%%%
% The claim is that these are n different roots
% of unity, so each of these n values, when
% raised to the power n, should equal 1.
%
% When you look at the output, however, a value
% may not look exactly like
%
%     1 + 0j.
%
% Instead, it may look like, for example,
%
%     9.999999999999999e-01
%       - 1.276756478318930e-15j.
%
% This initially may not look like 1 + 0j, but
% look closely at its components.
%
% The real component is
%
%     0.9999999999999999,
%
% which is extremely close to 1.
%
% The imaginary component is
%
%     -1.276756478318930e-15,
%
% or
%
%                              -15
%     -1.276756478318930 x 10,
%
% which is extremely close to zero.
%
% The reason is that floating-point
% representations and floating-point arithmetic
% are not perfectly precise.
%
% A floating-point number stores only a finite
% number of significant binary digits. Thus,
% even constants such as pi can only be stored
% approximately.
%
% The exact value of pi begins
%
%     3.1415926535897932384626433832795...
%
% while the stored floating-point approximation
% is approximately
%
%     3.141592653589793.
%
% Consequently, calculations involving pi,
% sine, cosine, multiplication, and powers all
% introduce small rounding errors.
%
% Thus, when checking numerical calculations,
% we often ask whether a result is sufficiently
% close to the expected value rather than
% requiring exact equality.
%%%%%%%%%%%%%%%%%%% IMPORTANT %%%%%%%%%%%%%%%%%%%