이전 글까지는 실수 실수행렬과 벡터에 대한 선형대수 연산을 중심으로 다뤘습니다. 이제 복소수(Complex Number)를 지원함으로써 신호처리, 스펙트럼 분석, 물리 시뮬레이션 등 광범위한 분야에서 필수적인 푸리에 변환(FFT, Fast Fourier Transform)을 구현하고, 그 과정에서 다항식 연산(Polynomial Arithmetic) 문제도 함께 다뤄보겠습니다.
이번 글의 핵심 주제:
- 복소수 지원: num_complex 크레이트를 활용하거나 직접 Complex<f64> 타입을 정의해 복소수를 다루고, 벡터/행렬/다항식에서 복소수를 처리할 수 있도록 확장합니다.
- 다항식 연산(Polynomial Arithmetic): 다항식 덧셈, 곱셈, 분할 등을 벡터 형태로 구현합니다. 다항식 곱셈을 빠르게 처리하기 위해 FFT를 활용할 수 있습니다.
- FFT(고속 푸리에 변환): FFT는 O(n log n) 시간에 다항식 곱셈, 신호 변환, 이미지 처리 등에 널리 사용됩니다. Cooley-Tukey 알고리즘을 간단히 구현하고, 이를 다항식 곱셈 가속에 이용해봅니다.
복소수 지원
러스트에는 num_complex 크레이트가 있어 복소수 타입을 쉽게 사용할 수 있습니다.
Cargo.toml에 추가:
[dependencies]
num-complex = "0.4"
// src/complex_support.rs
use num_complex::Complex;
pub type C64 = Complex<f64>;
// 예: C64를 이용한 간단한 연산
// 벡터/행렬을 복소수 타입으로 선언하면 동일한 연산을 지원할 수 있음
이미 구현한 Vector<T>나 Matrix<T>를 제네릭으로 작성했다면, T = C64를 사용해 복소수 벡터/행렬을 다룰 수 있습니다. 다만 LU, QR 등 선형대수 연산을 복소수에 적용하려면 알맞은 수학적 처리(Conjugate Transpose, 헥시단서적 안정성 등)가 필요합니다. 여기서는 FFT에 초점을 맞추어 복소 벡터를 다루는 예를 보여주겠습니다.
다항식 연산(Polynomial)
다항식을 [a0, a1, a2, ...] 형태의 벡터로 표현하면 p(x) = a0 + a1*x + a2*x² + ... 라고 할 수 있습니다.
// src/polynomial.rs
pub struct Polynomial<T> {
coeffs: Vec<T>
}
impl<T: Copy> Polynomial<T> {
pub fn new(coeffs: Vec<T>) -> Self {
Self { coeffs }
}
pub fn degree(&self) -> usize {
if self.coeffs.is_empty() {
0
} else {
self.coeffs.len() - 1
}
}
pub fn evaluate(&self, x: T, mul: fn(T,T)->T, add: fn(T,T)->T) -> T {
// 호너(Horner) 방식으로 다항식 평가
let mut result = self.coeffs[self.degree()];
for &c in self.coeffs[..self.degree()].iter().rev() {
result = add(mul(result,x), c);
}
result
}
}
// 실수 다항식을 예로 들면 mul = |a,b| a*b, add = |a,b| a+b
// 복소수나 다른 타입에도 mul/add 함수를 바꿔 적용 가능
다항식 덧셈 & 곱셈
impl<T: Copy + Default> Polynomial<T> {
pub fn add(&self, other: &Self, add: fn(T,T)->T) -> Self {
let len = self.coeffs.len().max(other.coeffs.len());
let mut result = vec![T::default(); len];
for i in 0..len {
let a = if i < self.coeffs.len() { self.coeffs[i] } else { T::default() };
let b = if i < other.coeffs.len() { other.coeffs[i] } else { T::default() };
result[i] = add(a,b);
}
Self::new(result)
}
}
// 다항식 곱셈을 직관적으로 하면 O(n²) 시간. FFT를 통한 O(n log n) 개선을 목표로 한다면
// 다음 섹션에서 FFT를 구현한 뒤 FFT 기반 곱셈을 보여주겠다.
FFT 구현
Cooley-Tukey FFT 알고리즘을 간단히 구현해 봅시다. FFT는 복소수 벡터를 받아 주파수 영역으로 변환합니다. 여기서는 길이가 2의 거듭제곱인 벡터를 처리하는 단순 FFT를 예로 듭니다.
// src/fft.rs
use crate::complex_support::C64;
use std::f64::consts::PI;
pub fn fft(a: &mut [C64], invert: bool) {
let n = a.len();
// 비트 반전(bit-reversal) 순서 재배치
let mut j = 0;
for i in 1..(n-1) {
let mut bit = n >> 1;
while j & bit != 0 {
j ^= bit;
bit >>= 1;
}
j ^= bit;
if i < j {
a.swap(i,j);
}
}
let mut len = 2;
while len <= n {
let ang = 2.0*PI/(len as f64)*(if invert {-1.0}else{1.0});
let wlen = C64::new(ang.cos(), ang.sin());
let half = len/2;
let mut i = 0;
while i < n {
let mut w = C64::new(1.0,0.0);
for j in 0..half {
let u = a[i+j];
let v = a[i+j+half]*w;
a[i+j] = u+v;
a[i+j+half] = u-v;
w = w*wlen;
}
i += len;
}
len <<= 1;
}
if invert {
for x in a {
*x /= n as f64;
}
}
}
이 FFT 구현은 기본형이고, 수치안정성이나 최적화를 위해 더 많은 작업이 필요할 수 있습니다.
FFT 기반 다항식 곱셈
두 다항식 P, Q를 곱할 때 차수가 N이라 하면 직접 곱셈은 O(N²)입니다. FFT를 이용하면 O(N log N)에 할 수 있습니다.
- P와 Q를 복소수 벡터로 변환(계수만 복소수 실수부에 담음).
- FFT를 각각 적용해 주파수 영역으로 변환
- 주파수 영역에서 요소별 곱(element-wise multiplication)
- IFFT(FFT with invert=true)를 적용해 다시 시간 영역(계수)로 복귀
// src/polynomial.rs (추가)
use crate::fft::fft;
use crate::complex_support::C64;
impl Polynomial<f64> {
pub fn mul_fft(&self, other: &Self) -> Self {
let mut n = 1;
let len = self.coeffs.len() + other.coeffs.len() - 1;
while n < len { n <<= 1; }
let mut fa = vec![C64::new(0.0,0.0); n];
let mut fb = vec![C64::new(0.0,0.0); n];
for i in 0..self.coeffs.len() {
fa[i] = C64::new(self.coeffs[i],0.0);
}
for i in 0..other.coeffs.len() {
fb[i] = C64::new(other.coeffs[i],0.0);
}
fft(&mut fa, false);
fft(&mut fb, false);
for i in 0..n {
fa[i] *= fb[i];
}
fft(&mut fa, true);
let mut result = fa[..len].iter().map(|x| x.re).collect::<Vec<f64>>();
// 반올림 처리 (수치 오차)
for r in &mut result {
*r = (*r).round();
}
Polynomial::new(result)
}
}
예제
use mathlib::polynomial::Polynomial;
fn main() {
let p = Polynomial::new(vec![1.0,2.0,3.0]); // 1+2x+3x²
let q = Polynomial::new(vec![4.0,5.0]); // 4+5x
let pq = p.mul_fft(&q);
// (1+2x+3x²)(4+5x) = 4 +13x +22x² +15x³
assert_eq!(pq.coeffs, vec![4.0,13.0,22.0,15.0]);
}
앞으로의 학습 방향
4편에서는 복소수 지원을 통해 FFT를 구현하고, 이를 이용한 다항식 곱셈 가속화를 예제로 다뤄보았습니다. 앞으로 다음과 같은 영역을 탐색할 수 있습니다.
- 고급 최적화: SIMD, 병렬화(rayong), GPU 가속 등을 통한 성능 향상
- 확률/통계 함수, 최적화 알고리즘: MCMC, 회귀 분석, 최소자승법 구현
- 미분/적분, ODE/PDE 솔버: 수치해석 전반으로 범위 확장
- 자동미분(AD) 및 머신러닝 커스텀 연산: 심볼릭 연산/AD 지원
결론
이번 글에서 복소수 지원, 다항식 연산, FFT를 통한 빠른 다항식 곱셈까지 구현해보며 러스트로 만들어가는 수학/과학 계산 라이브러리가 점점 더 범용적이고 강력해질 수 있음을 확인했습니다. FFT는 신호처리, 이미지처리, 수치해석, 머신러닝까지 광범위하게 응용 가능한 핵심 알고리즘이며, 러스트의 안전성과 성능 잠재력을 활용하면 고성능 FFT 라이브러리 구현도 가능합니다.
앞으로도 다양한 분야의 수치 알고리즘을 러스트로 구현하여, 견고하고 효율적인 수학 라이브러리를 완성해나가 봅시다.
유용한 링크와 리소스
- num-complex 문서: https://docs.rs/num-complex
- FFT 알고리즘 설명: "Numerical Recipes", "Understanding FFT" 관련 자료
- Rust FFT 크레이트 예시: rustfft, ndrustfft 등