e = cos θ + i sin θ
The Circle as an Exponential

A single equation uniting five fundamental constants, bridging exponential growth, trigonometry, and the geometry of rotation — often called the most beautiful formula in mathematics.

▼   scroll to explore   ▼

Section 01

Building Blocks

Before we can appreciate Euler's formula, we need to understand each ingredient: the mysterious number e, the imaginary unit i, and the geometry of the complex plane.

Euler's Number e ≈ 2.71828…

The number e is the unique base for which the exponential function equals its own derivative. It arises naturally from compound interest taken to its continuous limit:

e = limn→∞ (1 + 1/n)n

Equivalently, via its Taylor series (which we'll use heavily later):

ex = 1 + x + x²/2! + x³/3! + x⁴/4! + …
Key property: d/dx[ex] = ex. The exponential function is the only function (up to a constant multiple) that equals its own derivative. This makes it the natural language of growth, decay, and oscillation.
The function y = eˣ (cyan) alongside its derivative — they are identical.

Approaching e via Compound Interest

1 2.000
As n increases, (1 + 1/n)ⁿ approaches e = 2.71828…

Key Values

ExpressionValueMeaning
e⁰1Any nonzero base to the 0 is 1
2.71828…Euler's number itself
e−10.36788…Reciprocal of e
e−1Euler's Identity — the crown jewel
e2πi1One full rotation returns to the start

The Imaginary Unit i

We define i as the square root of negative one:

i² = −1   ⟹   i = √(−1)

At first this seems impossible — no real number squares to give a negative. But i is perfectly well-defined as an algebraic object, and introducing it completes the real numbers into the complex numbers.

Powers of i — The Cycle

The powers of i cycle with period 4:

PowerValueGeometric meaning
i⁰1Pointing right (+x)
iPointing up (+y), rotate 90°
−1Pointing left (−x), rotate 180°
−iPointing down (−y), rotate 270°
i⁴1Back to start, rotate 360°
Geometric insight: Multiplying by i rotates a point in the complex plane by 90° counterclockwise. Multiplying by i twice (i.e., by i² = −1) rotates 180°. This is the key to why exponentials can trace circles.
Powers of i rotating around the unit circle.

The Complex Plane

A complex number z = a + bi is plotted on a 2D plane: real part a on the horizontal axis, imaginary part b on the vertical axis.

z = a + bi = r·e = r(cos θ + i sin θ)

Every complex number can be written in either rectangular form (a + bi) or polar form (r e). Here:

  • r = |z| = √(a² + b²) — the magnitude (distance from origin)
  • θ = arg(z) = arctan(b/a) — the argument (angle from positive x-axis)
1.50 1.20
z = 1.50 + 1.20i  |  |z| = 1.92  |  θ = 38.7°

Sine and Cosine on the Unit Circle

On the unit circle, a point at angle θ has coordinates (cos θ, sin θ). These functions describe the projection of circular motion onto the real and imaginary axes:

cos θ = Re(e)      sin θ = Im(e)
Euler's formula says: a point on the unit circle at angle θ is the complex exponential e. The two descriptions (trigonometric and exponential) are identical.
0.80π
cos(θ) = 0.588  |  sin(θ) = 0.809

Taylor Series — Polynomials in Disguise

Any smooth function can be approximated as an infinite polynomial around a point. Around x = 0:

f(x) = f(0) + f'(0)x + f''(0)x²/2! + f'''(0)x³/3! + …

The three key series we'll need are:

FunctionTaylor Series
ex 1 + x + x²/2! + x³/3! + x⁴/4! + x⁵/5! + …
cos x 1 − x²/2! + x⁴/4! − x⁶/6! + …
sin x x − x³/3! + x⁵/5! − x⁷/7! + …
Notice the pattern: cos uses even powers, sin uses odd powers. Euler's formula emerges precisely from substituting into the ex series and sorting the terms by their real/imaginary nature.
3
Showing 3-term Taylor approximation of eˣ

Section 02

The Derivation

We derive Euler's formula by substituting into the Taylor series of ex and then recognizing the resulting pattern.

Start with the Taylor series for eˣ

ex = 1 + x + x²/2! + x³/3! + x⁴/4! + x⁵/5! + x⁶/6! + x⁷/7! + …

This converges for all complex numbers, not just real ones. So we're free to substitute x = iθ.

Substitute x = iθ

e = 1 + (iθ) + (iθ)²/2! + (iθ)³/3! + (iθ)⁴/4! + (iθ)⁵/5! + …

Now we expand each power using the cycling rule i² = −1, i³ = −i, i⁴ = 1, i⁵ = i, …

Expand the powers of (iθ)ⁿ

TermExpand iⁿSimplified
(iθ)⁰i⁰ = 11
(iθ)¹/1!i¹ = i
(iθ)²/2!i² = −1−θ²/2!
(iθ)³/3!i³ = −i−iθ³/3!
(iθ)⁴/4!i⁴ = +1+θ⁴/4!
(iθ)⁵/5!i⁵ = +i+iθ⁵/5!
(iθ)⁶/6!i⁶ = −1−θ⁶/6!
(iθ)⁷/7!i⁷ = −i−iθ⁷/7!

Group real and imaginary terms

e = (1 − θ²/2! + θ⁴/4! − θ⁶/6! + …) + i(θ − θ³/3! + θ⁵/5! − θ⁷/7! + …)

Recognize the Taylor series for cos and sin

cos θ = 1 − θ²/2! + θ⁴/4! − θ⁶/6! + …
sin θ = θ − θ³/3! + θ⁵/5! − θ⁷/7! + …

Conclude: Euler's Formula

e = cos θ + i sin θ

The real part of e is cos θ, and the imaginary part is sin θ. This is not an approximation — it's an exact identity, valid for all real (and complex) θ.


Visual Convergence: Taylor Vectors

Each term in the Taylor series is a complex number (a 2D vector). They chain end-to-end; as we add more terms, the chain spirals in and converges to the point e = cos θ + i sin θ on the unit circle.

1.00π 6
Showing 6-term vector chain converging to e^(iθ)

Section 03

The Unit Circle Explorer

As θ varies from 0 to 2π, the complex number e traces the unit circle exactly once. Watch how the real and imaginary projections generate the cosine and sine waves.

0.32π
e^(i·0.32π) = cos(0.32π) + i·sin(0.32π) = 0.588 + 0.809i
The Real Projection
As the point rotates, its x-coordinate (horizontal shadow) traces out a cosine wave. Cosine starts at 1, dips to −1 at θ = π, and returns.
The Imaginary Projection
The y-coordinate (vertical shadow) traces a sine wave. Sine starts at 0, peaks at 1 when θ = π/2, and oscillates.

Inverse Formulas

By adding and subtracting Euler's formula for +θ and −θ, we recover cosine and sine:

cos θ = (e + e−iθ) / 2
sin θ = (e − e−iθ) / 2i

This is how complex exponentials secretly are trigonometry — they are not separate subjects but two faces of the same mathematical object.


Section 04

Euler's Identity

Setting θ = π in Euler's formula yields the most celebrated equation in mathematics:

e = cos π + i sin π = −1 + i·0 = −1
e + 1 = 0

In this single equation five fundamental constants of mathematics appear exactly once, connected by two operations (addition and exponentiation).

e Euler's number
base of natural log
i Imaginary unit
√(−1)
π Pi
ratio circ./diam.
1 Multiplicative
identity
0 Additive
identity
The point e^(iθ) walks from 1 (θ=0) to −1 (θ=π) along the upper semicircle.
Why is this deep? It is not just that five constants appear — it is that three entirely different mathematical worlds (analysis via e, geometry via π, algebra via i) intersect at the simplest possible relationship. The number i carries the formula into the imaginary axis, π specifies the half-turn, and e provides the exponential scaffold. The result lands exactly at −1, the reflection of 1 through the origin.

Section 05

Rotation & Complex Multiplication

Euler's formula reveals that multiplying by e rotates any complex number by angle φ. This is one of its most powerful applications.

z · e = r·e · e = r·ei(θ+φ)

The magnitudes multiply and the angles add. Rotation is just multiplication by a unit complex number on the circle.

1.20 0.40π
+0.50π

De Moivre's Theorem

A direct consequence of Euler's formula: raising a complex number to an integer power simply multiplies the angle:

(cos θ + i sin θ)n = cos(nθ) + i sin(nθ)
Proof: (e)n = einθ = cos(nθ) + i sin(nθ). ∎ This one-line proof, via Euler's formula, replaces tedious trigonometric induction.

Deriving Trig Identities in Seconds

Using ei(α+β) = e · e and expanding both sides immediately gives the angle addition formulas:

cos(α+β) = cos α cos β − sin α sin β
sin(α+β) = sin α cos β + cos α sin β

No geometric diagrams or induction required — it follows directly from exponential algebra.


Section 06

Roots of Unity

The equation zn = 1 has exactly n solutions (roots) in the complex numbers. Using Euler's formula, they are evenly spaced around the unit circle:

zk = e2πik/n = cos(2πk/n) + i·sin(2πk/n)    for k = 0, 1, …, n−1

These roots form the vertices of a regular n-gon inscribed in the unit circle. This is the foundation of the Discrete Fourier Transform.

6
6th roots of unity — vertices of a regular hexagon

Primitive Root

The first nontrivial root, ω = e2πi/n, is called the primitive nth root of unity. Its powers generate all the others:

1, ω, ω², ω³, …, ωn−1
Sum to zero: The n-th roots of unity always sum to zero — 1 + ω + ω² + … + ωn−1 = 0. Geometrically, the vectors cancel by symmetry. This is the orthogonality property that makes the Fourier transform work.

Roots Listed (n = selected above)

ke^(2πik/n)AngleReIm

Section 07

Phasors & AC Circuits

In electrical engineering and physics, a phasor is a rotating complex exponential representing a sinusoidal signal. Instead of working with A cos(ωt + φ) directly, engineers use the compact representation:

v(t) = Re[V·eiωt]    where V = A·e

Adding two sinusoids of the same frequency becomes simple vector addition of phasors. Differentiation becomes multiplication by — turning differential equations into algebra.

1.00 0.00π
1.0
Single phasor: A=1.00, φ=0, ω=1.0
Impedance: In AC circuit analysis, a capacitor's impedance is ZC = 1/(iωC) and an inductor's is ZL = iωL. The factor of i encodes the 90° phase shift between voltage and current — rotation in the complex plane.

Section 08

Applications Across Science

Euler's formula is not merely beautiful — it is indispensable across mathematics, physics, and engineering.

🌊

Fourier Transform

The DFT and FFT decompose signals into frequencies using e−2πiνt as basis functions. Every spectral analysis tool — audio EQ, MRI, radio — depends on Euler's formula.

AC Circuit Analysis

Complex impedance (Z = R + iX) encodes amplitude and phase simultaneously. Ohm's law, V = IZ, handles reactive components without differential equations.

⚛️

Quantum Mechanics

Wave functions are complex-valued: ψ = Aei(kx−ωt). The Schrödinger equation, propagators, and S-matrices are all expressed via complex exponentials.

🎮

2D/3D Graphics

Multiplying by e rotates a 2D point. Quaternions (the 3D extension) use the same principle to represent 3D rotations without gimbal lock.

📡

Signal Processing

Modulation, filtering, and demodulation use complex exponentials. The analytic signal z(t) = x(t) + iH{x(t)} gives instantaneous amplitude & phase.

🔢

Number Theory

The Riemann ζ function, character sums, and Dirichlet series all use complex exponentials. The Riemann Hypothesis concerns zeros of ζ in the complex plane.

🌀

Fluid Dynamics

Conformal mappings of the complex plane transform complicated geometries into simple ones. Potential flow and airfoil theory are expressed as complex potentials.

🔬

Optics & Diffraction

Fraunhofer diffraction and holography use Fourier optics: the diffraction pattern is the Fourier transform of the aperture function.

🎵

Music & Audio

Every digital audio effect — reverb, pitch-shifting, noise removal — is computed in the frequency domain via the FFT, which is Euler's formula at its core.


Section 09

GNU Octave Examples

Ten complete, runnable GNU Octave scripts exploring every facet of Euler's formula. Each includes comments explaining the mathematics. Copy any example into Octave or a compatible MATLAB environment.

1. Verify Euler's Identity Numerically

Directly compute e and confirm it equals −1 to machine precision.

% ── Euler's Identity: e^(i*pi) + 1 = 0 ──────────────────────────────────
% Compute e^(i*pi) using Octave's built-in exp() and the imaginary unit 1i

result = exp(1i * pi);   % complex exponential: should equal -1

fprintf('e^(i*pi) = %.15f + %.15fi\n', real(result), imag(result));
fprintf('e^(i*pi) + 1 = %.2e  (machine epsilon = %.2e)\n', ...
        abs(result + 1), eps);

% ── Direct component check ────────────────────────────────────────────────
theta = pi;
using_formula = cos(theta) + 1i * sin(theta);   % Euler's formula
using_exp     = exp(1i * theta);                 % exponential form

fprintf('\nEuler formula result:  %+.10f %+.10fi\n', ...
        real(using_formula), imag(using_formula));
fprintf('exp(i*theta) result:   %+.10f %+.10fi\n', ...
        real(using_exp), imag(using_exp));
fprintf('Difference: %.2e\n', abs(using_formula - using_exp));

% ── Explore e^(i*theta) for several values ────────────────────────────────
fprintf('\n  θ/π     |  cos(θ)     |  sin(θ)     |  e^(iθ)\n');
fprintf('%s\n', repmat('-', 1, 55));
for k = 0:8
  t = k * pi / 4;
  z = exp(1i * t);
  fprintf('  %4.2f    |  %+.6f  |  %+.6f  |  %+.4f%+.4fi\n', ...
          k/4, cos(t), sin(t), real(z), imag(z));
end

2. Plot the Unit Circle with e^(iθ)

Visualize how e traces the unit circle, and project the real and imaginary parts to produce cosine and sine waves.

% ── Unit Circle and Wave Projections ─────────────────────────────────────
theta = linspace(0, 2*pi, 360);
z     = exp(1i * theta);          % unit circle as complex exponential
x     = real(z);                   % cos(θ)
y     = imag(z);                   % sin(θ)

figure('Name', 'Euler Formula: Unit Circle', 'Color', 'w');
subplot(1, 3, 1);
  plot(x, y, 'b-', 'LineWidth', 2);  hold on;
  plot(1, 0, 'ro', 0, 1, 'go', -1, 0, 'ko', 0, -1, 'mo', 'MarkerSize', 8);
  axis equal; grid on; axis([-1.4 1.4 -1.4 1.4]);
  title('Unit Circle: e^{i	heta}');
  xlabel('Re(z) = cos	heta'); ylabel('Im(z) = sin	heta');
  legend({'e^{i	heta}','	heta=0','	heta=\pi/2','	heta=\pi','	heta=3\pi/2'}, ...
         'Location', 'SouthEast');

subplot(1, 3, 2);
  plot(theta/pi, x, 'b-', 'LineWidth', 2);
  grid on; xlabel('	heta / \pi'); ylabel('cos	heta = Re(e^{i	heta})');
  title('Real Part: Cosine Wave');
  xticks(0:0.5:2); xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});

subplot(1, 3, 3);
  plot(theta/pi, y, 'r-', 'LineWidth', 2);
  grid on; xlabel('	heta / \pi'); ylabel('sin	heta = Im(e^{i	heta})');
  title('Imaginary Part: Sine Wave');
  xticks(0:0.5:2); xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});

sgtitle('Euler\'s Formula: e^{i	heta} = cos	heta + i sin	heta', 'FontSize', 14);

3. Taylor Series Convergence

Visualize how partial sums of the ex Taylor series converge to the true function, and watch the convergence of e to the unit circle.

% ── Taylor Series Convergence for e^(iθ) ─────────────────────────────────
theta_vals = linspace(0, 2*pi, 500);
N_max = 15;   % max number of terms

% Build partial sums: sum_{k=0}^{N} (iθ)^k / k!
figure('Name', 'Taylor Convergence');
colors = jet(N_max);

subplot(1, 2, 1);
hold on;
for N = [1 2 3 5 8 15]
  partial = zeros(size(theta_vals));
  for k = 0:N
    partial = partial + (1i * theta_vals).^k ./ factorial(k);
  end
  plot(real(partial), imag(partial), 'LineWidth', 1.5, 'Color', colors(N, :));
end
% Exact unit circle
z_exact = exp(1i * theta_vals);
plot(real(z_exact), imag(z_exact), 'k--', 'LineWidth', 2);
axis equal; grid on; axis([-1.5 1.5 -1.5 1.5]);
legend({'N=1','N=2','N=3','N=5','N=8','N=15','Exact'}, 'Location', 'Best');
title('Taylor partial sums trace the unit circle');
xlabel('Re'); ylabel('Im');

subplot(1, 2, 2);
% Error vs number of terms for θ = π
theta_test = pi;
errors = zeros(1, N_max);
partial_sum = 0;
for N = 0:N_max-1
  partial_sum = partial_sum + (1i*theta_test)^N / factorial(N);
  errors(N+1) = abs(partial_sum - exp(1i*theta_test));
end
semilogy(1:N_max, errors, 'bo-', 'LineWidth', 2, 'MarkerFaceColor', 'b');
grid on; xlabel('Number of terms N'); ylabel('|error|');
title('Convergence error at 	heta = \pi');
sgtitle('Taylor Series Convergence of e^{i	heta}');

4. Complex Multiplication as Rotation

Demonstrate that multiplying by e rotates any complex number by angle φ without changing its magnitude.

% ── Complex Multiplication = Rotation ────────────────────────────────────
z = 1.5 + 0.8*1i;          % original complex number
phi_vals = 0 : pi/6 : 2*pi; % rotate in steps of 30°

figure('Name', 'Complex Rotation');
subplot(1, 2, 1);
hold on; axis equal; grid on;
xlim([-2 2]); ylim([-2 2]);
xlabel('Re'); ylabel('Im');
title('z rotated by multiples of 30°');

colors = hsv(length(phi_vals));
for k = 1:length(phi_vals)
  phi = phi_vals(k);
  z_rot = z * exp(1i * phi);   % rotate z by phi
  plot([0 real(z_rot)], [0 imag(z_rot)], '-o', ...
       'Color', colors(k,:), 'LineWidth', 1.5, 'MarkerFaceColor', colors(k,:));
end
plot(real(z), imag(z), 'k*', 'MarkerSize', 12);  % original

% Verify: magnitude is preserved
z_rot_test = z * exp(1i * pi/3);
fprintf('Original z:   %.4f + %.4fi\n', real(z), imag(z));
fprintf('After 60° rot: %.4f + %.4fi\n', real(z_rot_test), imag(z_rot_test));
fprintf('|z| before: %.6f,  |z| after: %.6f  (unchanged)\n', abs(z), abs(z_rot_test));
fprintf('Angle before: %.2f°, angle after: %.2f°  (shifted by 60°)\n', ...
        angle(z)*180/pi, angle(z_rot_test)*180/pi);

subplot(1, 2, 2);
% Demonstrate |z*w| = |z|*|w| and arg(z*w) = arg(z)+arg(w)
n_pts = 200;
z1 = (1.2) * exp(1i * linspace(0, 2*pi, n_pts));  % circle of radius 1.2
z2 = (0.7) * exp(1i * pi/4);                       % fixed rotator
product = z1 .* z2;
plot(real(z1), imag(z1), 'b-', 'LineWidth', 2); hold on;
plot(real(product), imag(product), 'r-', 'LineWidth', 2);
axis equal; grid on;
legend({'z_1 (r=1.2)', 'z_1 	imes e^{i\pi/4} (r=0.84, rotated 45°)'});
title('Multiplying scales AND rotates'); xlabel('Re'); ylabel('Im');

5. De Moivre's Theorem & nth Roots of Unity

% ── De Moivre's Theorem and Roots of Unity ───────────────────────────────
% De Moivre: (cos θ + i sin θ)^n = cos(nθ) + i sin(nθ)

theta = pi / 5;
for n = [1 2 3 5 10]
  via_demoivre   = cos(n*theta) + 1i*sin(n*theta);
  via_direct_pow = (cos(theta)  + 1i*sin(theta))^n;
  via_euler      = exp(1i*n*theta);
  err = max(abs(via_demoivre - via_direct_pow), abs(via_demoivre - via_euler));
  fprintf('n=%2d: DeMoivre=%.6f%+.6fi  |error|=%.2e\n', ...
          n, real(via_demoivre), imag(via_demoivre), err);
end

% ── nth Roots of Unity ────────────────────────────────────────────────────
figure('Name', 'Roots of Unity');
for subplot_n = 1:6
  n = subplot_n + 1;            % n from 2 to 7
  k = 0:n-1;
  roots_n = exp(2i*pi*k/n);    % z^n = 1 solutions: e^(2πi k/n)
  subplot(2, 3, subplot_n);
    circle_t = linspace(0, 2*pi, 200);
    plot(cos(circle_t), sin(circle_t), 'k:'); hold on;
    plot(real(roots_n), imag(roots_n), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
    % Draw polygon
    patch(real(roots_n), imag(roots_n), 'r', 'FaceAlpha', 0.1, 'EdgeColor', 'r');
    for kk = 1:n
      text(real(roots_n(kk))*1.2, imag(roots_n(kk))*1.2, ...
           sprintf('k=%d', kk-1), 'FontSize', 7, 'HorizontalAlignment', 'center');
    end
    axis equal; axis([-1.6 1.6 -1.6 1.6]); grid on;
    title(sprintf('%d-th roots of unity', n));
end
sgtitle('z^n = 1  →  e^{2\pi i k/n}');

% ── Verify sum = 0 ────────────────────────────────────────────────────────
for n = 2:8
  s = sum(exp(2i*pi*(0:n-1)/n));
  fprintf('Sum of %d-th roots of unity = %.2e (should be ~0)\n', n, abs(s));
end

6. Phasor Animation (AC Circuit)

% ── Phasor Animation: AC Circuit with R, L, C ────────────────────────────
% Component values
R   = 100;       % Ohms
L   = 0.1;       % Henry
C   = 1e-4;      % Farad
V0  = 1;         % source amplitude (V)
f_range = logspace(1, 4, 500);   % 10 Hz to 10 kHz

omega = 2*pi*f_range;

% Complex impedances: Z_R, Z_L, Z_C
Z_R = R * ones(size(omega));
Z_L = 1i * omega * L;
Z_C = 1 ./ (1i * omega * C);
Z_total = Z_R + Z_L + Z_C;

% Current phasor I = V/Z  (Ohm's law in complex form)
I = V0 ./ Z_total;

figure('Name', 'RLC Phasor Analysis');

subplot(2,2,1);
  semilogx(f_range, abs(I), 'b-', 'LineWidth', 2); grid on;
  xlabel('Frequency (Hz)'); ylabel('|I| (A)');
  title('Current Amplitude vs Frequency');

subplot(2,2,2);
  semilogx(f_range, angle(I)*180/pi, 'r-', 'LineWidth', 2); grid on;
  xlabel('Frequency (Hz)'); ylabel('Phase (degrees)');
  title('Phase of Current vs Frequency');
  yline(0, 'k--');

% Resonance frequency
f0 = 1 / (2*pi*sqrt(L*C));
fprintf('Resonance frequency: %.2f Hz\n', f0);
fprintf('At resonance: Z_L = %+.2fi, Z_C = %+.2fi  (cancel!)\n', ...
        imag(1i*2*pi*f0*L), imag(1/(1i*2*pi*f0*C)));

% Time-domain phasor at a specific frequency
subplot(2,2,[3,4]);
t = linspace(0, 3/f0, 1000);
omega0 = 2*pi*f0;
V_t = V0 * cos(omega0 * t);
I_phasor = V0 / R;            % at resonance, Z = R
I_t = I_phasor * cos(omega0 * t);
plot(t*1000, V_t, 'b-', 'LineWidth', 2); hold on;
plot(t*1000, I_t*R, 'r--', 'LineWidth', 2);
xlabel('Time (ms)'); ylabel('Amplitude'); grid on;
legend({'V(t) = V_0 cos(\omega_0 t)', 'I(t)	imesR  (in phase at resonance)'});
title(sprintf('Time Domain at Resonance f_0 = %.1f Hz', f0));

7. Trig Identities from Euler's Formula

Derive and verify standard trigonometric identities algebraically using complex exponentials.

% ── Trig Identities via Euler's Formula ──────────────────────────────────
% cos(α+β) = cos α cos β - sin α sin β
% sin(α+β) = sin α cos β + cos α sin β
% Proof: e^{i(α+β)} = e^{iα} · e^{iβ}

alpha = 0.7; beta = 1.2;

lhs = exp(1i*(alpha+beta));
rhs = exp(1i*alpha) * exp(1i*beta);
fprintf('e^{i(a+b)} = %.6f %+.6fi\n', real(lhs), imag(lhs));
fprintf('e^{ia}*e^{ib}= %.6f %+.6fi  |error|=%.2e\n\n', real(rhs), imag(rhs), abs(lhs-rhs));

% Expanding: Real part gives cos addition, imaginary gives sin addition
cos_add_euler  = real(exp(1i*alpha) * exp(1i*beta));
cos_add_trig   = cos(alpha)*cos(beta) - sin(alpha)*sin(beta);
sin_add_euler  = imag(exp(1i*alpha) * exp(1i*beta));
sin_add_trig   = sin(alpha)*cos(beta) + cos(alpha)*sin(beta);
fprintf('cos(a+b): Euler=%.8f  Trig=%.8f\n', cos_add_euler, cos_add_trig);
fprintf('sin(a+b): Euler=%.8f  Trig=%.8f\n\n', sin_add_euler, sin_add_trig);

% ── Double angle formulas ─────────────────────────────────────────────────
% (e^{iθ})^2 = e^{2iθ}  →  expand (cosθ + i sinθ)^2
theta = 0.9;
z2 = (cos(theta) + 1i*sin(theta))^2;
fprintf('cos(2θ): Euler=%.8f  Direct=%.8f\n', real(z2), cos(2*theta));
fprintf('sin(2θ): Euler=%.8f  Direct=%.8f\n\n', imag(z2), sin(2*theta));

% ── Product-to-sum identities ─────────────────────────────────────────────
% cos α cos β = (1/2)[cos(α-β) + cos(α+β)]
% Proof: Re(e^{iα}) * Re(e^{iβ}) = (1/2)Re(e^{i(α-β)} + e^{i(α+β)})
lhs_prod = cos(alpha) * cos(beta);
rhs_prod = 0.5 * (cos(alpha-beta) + cos(alpha+beta));
fprintf('cos(a)cos(b):  Direct=%.8f  Product-to-sum=%.8f\n', lhs_prod, rhs_prod);

% ── Plot several trig identities verified numerically ─────────────────────
t = linspace(0, 4*pi, 1000);
figure('Name', 'Trig Identities via Euler');
subplot(2,1,1);
  euler_expr = real(exp(2i*t));            % Re(e^{2it}) = cos(2t)
  trig_expr  = cos(t).^2 - sin(t).^2;    % classical double angle
  plot(t, euler_expr, 'b-', t, trig_expr, 'r--', 'LineWidth', 2);
  legend({'Re(e^{2it})', 'cos^2(t)-sin^2(t)'}); grid on;
  title('Double angle: cos(2t) = cos^2(t) - sin^2(t)');
subplot(2,1,2);
  euler_expr2 = imag(exp(2i*t));     % Im(e^{2it}) = sin(2t)
  trig_expr2  = 2*cos(t).*sin(t);    % classical form
  plot(t, euler_expr2, 'b-', t, trig_expr2, 'r--', 'LineWidth', 2);
  legend({'Im(e^{2it})', '2cos(t)sin(t)'}); grid on;
  title('Double angle: sin(2t) = 2cos(t)sin(t)');

8. Fourier Series with Complex Exponentials

Decompose a square wave into complex exponential harmonics, demonstrating how eiωt are the building blocks of all periodic signals.

% ── Fourier Series via Complex Exponentials ───────────────────────────────
% Square wave: f(t) = sign(sin(t)), period 2π
% Fourier coefficients: c_n = (1/2π) ∫₀²π f(t) e^{-int} dt
% For square wave: c_n = 0 (even n), c_n = 2/(inπ) (odd n)

T = 2*pi;
t = linspace(0, T, 2000);
f_exact = sign(sin(t));           % exact square wave

% Build partial Fourier sum with N harmonics
N_vals = [1, 3, 7, 15, 31];
figure('Name', 'Fourier Series (Complex Exponentials)');

for idx = 1:length(N_vals)
  N = N_vals(idx);
  f_approx = zeros(size(t));
  for n = -N:N
    if mod(n, 2) ~= 0 && n ~= 0     % only odd harmonics
      c_n = 2 / (1i * n * pi);       % Fourier coefficient
      f_approx = f_approx + real(c_n * exp(1i * n * t));
    end
  end
  subplot(2, 3, idx);
    plot(t, f_exact, 'k:', 'LineWidth', 1); hold on;
    plot(t, f_approx, 'b-', 'LineWidth', 1.5);
    title(sprintf('N = %d harmonics', N));
    xlabel('t'); ylim([-1.5 1.5]); grid on;
    if idx == 1, legend({'exact', 'Fourier sum'}); end
end

% Compute and display Fourier coefficients for first 10 harmonics
subplot(2,3,6);
n_plot = -9:2:9;              % odd harmonics only
c_plot = 2 ./ (1i * n_plot * pi);
stem(n_plot, abs(c_plot), 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b');
xlabel('Harmonic n'); ylabel('|c_n|');
title('Spectrum of square wave'); grid on;
sgtitle('Fourier Series: Square Wave as Sum of e^{int}');

9. 3D Helix: e in the Complex + θ Plane

Visualize the complex exponential in 3D: as θ grows, e traces a helix in the (Re, Im, θ) space.

% ── 3D Helix Visualization of e^{iθ} ─────────────────────────────────────
theta = linspace(0, 6*pi, 2000);
z = exp(1i * theta);
x = real(z);    % cos(θ)
y = imag(z);    % sin(θ)

figure('Name', 'e^{itheta} as a 3D Helix');
subplot(1,2,1);
  plot3(x, y, theta, 'b-', 'LineWidth', 2); hold on;
  % Project onto Re-Im plane (unit circle)
  plot3(x, y, zeros(size(theta)), 'k:', 'LineWidth', 1);
  % Project onto Re-θ plane (cosine wave)
  plot3(x, -ones(size(theta)), theta, 'r:', 'LineWidth', 1);
  % Project onto Im-θ plane (sine wave)
  plot3(ones(size(theta)), y, theta, 'g:', 'LineWidth', 1);
  xlabel('Re(e^{i	heta}) = cos	heta');
  ylabel('Im(e^{i	heta}) = sin	heta');
  zlabel('	heta (radians)');
  title('3D helix: e^{i	heta} with 3 projections');
  legend({'helix','unit circle (Re-Im)','cosine (Re-	heta)','sine (Im-	heta)'});
  grid on; view(35, 25);

subplot(1,2,2);
  % Growing spiral: r(θ) * e^{iθ} where r(θ) = e^{θ/10}
  r = exp(theta/10);           % growing amplitude
  z_spiral = r .* exp(1i*theta);
  plot3(real(z_spiral), imag(z_spiral), theta, 'm-', 'LineWidth', 2);
  xlabel('Re'); ylabel('Im'); zlabel('	heta');
  title('Growing spiral: e^{	heta/10} \cdot e^{i	heta}');
  grid on; view(35, 25);
  
sgtitle('The Complex Exponential in 3D');

10. DFT from Scratch vs. Built-in fft()

Implement the Discrete Fourier Transform using Euler's formula directly, then compare with Octave's optimized FFT.

% ── DFT from Scratch using Euler's Formula ───────────────────────────────
% DFT definition: X[k] = Σ_{n=0}^{N-1} x[n] · e^{-2πi k n / N}
% This is a direct application of complex exponentials as basis functions

function X = my_dft(x)
  N = length(x);
  X = zeros(1, N);
  for k = 0:N-1
    for n = 0:N-1
      X(k+1) = X(k+1) + x(n+1) * exp(-2i*pi*k*n/N);  % Euler here
    end
  end
end

% ── Test signal: sum of two sinusoids ─────────────────────────────────────
N = 64;
n = 0:N-1;
f1 = 5; f2 = 13;               % frequency bins
x = 2*cos(2*pi*f1*n/N) + 0.5*sin(2*pi*f2*n/N);

tic; X_mine  = my_dft(x); t_mine = toc;
tic; X_builtin = fft(x);  t_fft  = toc;

% Verify they match
max_err = max(abs(X_mine - X_builtin));
fprintf('Max error between my_dft and fft(): %.2e\n', max_err);
fprintf('my_dft time: %.4f s,  fft() time: %.6f s\n\n', t_mine, t_fft);

% Plot spectrum
figure('Name', 'DFT via Euler Formula');
k = 0:N/2-1;               % positive frequencies

subplot(2,1,1);
  plot(n, x, 'b-o', 'MarkerSize', 3); grid on;
  xlabel('Sample n'); ylabel('x[n]');
  title('Input signal: 2cos(2\pi 5n/N) + 0.5sin(2\pi 13n/N)');

subplot(2,1,2);
  stem(k, abs(X_builtin(1:N/2))/N*2, 'b', 'LineWidth', 2, ...
       'MarkerFaceColor', 'b'); grid on;
  xlabel('Frequency bin k'); ylabel('|X[k]|/N × 2');
  title(sprintf('DFT Spectrum (peaks at k=%d and k=%d)', f1, f2));
  xline(f1, 'r--'); xline(f2, 'r--');

sgtitle('DFT = Euler\'s formula applied N^2 times');

% ── Show the basis vectors explicitly ─────────────────────────────────────
fprintf('\nDFT basis vector for k=5 (first 8 elements):\n');
k5_basis = exp(-2i*pi*5*n/N);
for j = 0:7
  fprintf('  e^{-2πi·5·%d/%d} = %+.4f %+.4fi\n', j, N, real(k5_basis(j+1)), imag(k5_basis(j+1)));
end