▶ 书中第六章部分程序,加上自己补充的代码,包括快速傅里叶变换,多项式表示
● 快速傅里叶变换,需要递归
package package01; import edu.princeton.cs.algs4.StdOut;
import edu.princeton.cs.algs4.StdRandom;
import edu.princeton.cs.algs4.Complex; public class class01
{
private static final Complex ZERO = new Complex(0, 0);
private static final double EPSILON = 1e-8;
private class01() { } public static Complex[] fft(Complex[] x)// 正向傅里叶变换
{
int n = x.length;
if (n == 1)
return new Complex[] {x[0]};
if (n % 2 != 0)
throw new IllegalArgumentException("n != 2^k"); Complex[] temp = new Complex[n/2], result = new Complex[n];
for (int k = 0; k < n/2; k++) // 偶数项和奇数分别递归
temp[k] = x[2*k];
Complex[] q = fft(temp);
for (int k = 0; k < n/2; k++)
temp[k] = x[2*k + 1];
Complex[] r = fft(temp);
for (int k = 0; k < n/2; k++) // 合并
{
double kth = -2 * k * Math.PI / n;
Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
result[k] = q[k].plus(wk.times(r[k]));
result[k + n/2] = q[k].minus(wk.times(r[k]));
}
return result;
} public static Complex[] ifft(Complex[] x)// 傅里叶反变换,共轭、正变换、再共轭,最后除以项数
{
int n = x.length;
Complex[] result = new Complex[n], temp;
for (int i = 0; i < n; i++)
result[i] = x[i].conjugate();
result = fft(result);
for (int i = 0; i < n; i++)
result[i] = result[i].conjugate();
for (int i = 0; i < n; i++)
result[i] = result[i].scale(1.0 / n);
return result;
} public static Complex[] cconvolve(Complex[] x, Complex[] y)// 环形卷积
{
if (x.length != y.length)
throw new IllegalArgumentException("x.length != y.length");
int n = x.length;
Complex[] a = fft(x), b = fft(y), c = new Complex[n];
for (int i = 0; i < n; i++)
c[i] = a[i].times(b[i]);
return ifft(c);
} public static Complex[] convolve(Complex[] x, Complex[] y)// 线性卷积,后一半垫 0
{
if (x.length != y.length)
throw new IllegalArgumentException("x.length != y.length");
Complex[] a = new Complex[2 * x.length], b = new Complex[2 * y.length];
for (int i = 0; i < x.length; i++)
a[i] = x[i];
for (int i = x.length; i < 2*x.length; i++)
a[i] = ZERO;
for (int i = 0; i < y.length; i++)
b[i] = y[i];
for (int i = y.length; i < 2*y.length; i++)
b[i] = ZERO;
return cconvolve(a, b);
} private static void show(Complex[] x, String title)
{
StdOut.printf("%s\n-------------------\n",title);
for (int i = 0; i < x.length; i++)
StdOut.println(x[i]);
StdOut.println();
} public static void main(String[] args)
{
int n = 16;
Complex[] x = new Complex[n];
for (int i = 0; i < n; i++)
x[i] = new Complex(StdRandom.uniform(-1.0, 1.0), 0); show(x, "x");
Complex[] y = fft(x);
show(y, "y = fft(x)");
Complex[] z = ifft(y);
show(z, "z = ifft(y)");
Complex[] c = cconvolve(x, x);
show(c, "c = cconvolve(x, x)");
Complex[] d = convolve(x, x);
show(d, "d = convolve(x, x)");
}
}
● 多项式表示
package package01; import edu.princeton.cs.algs4.StdOut; public class class01
{
private int[] coef; // 系数列表
private int degree; // 次数 public class01(int a, int b) // 建立单项式 a*x^b
{
coef = new int[b + 1];
coef[b] = a;
reduce();
} private void reduce() // 记录多项式的次数
{
degree = -1;
for (int i = coef.length - 1; i >= 0; i--)
{
if (coef[i] != 0)
{
degree = i;
return;
}
}
} public int degree()
{
return degree;
} public class01 plus(class01 that) // p(x) + q(x),把系数从低次项向高次项各加一遍
{
class01 poly = new class01(0, Math.max(degree, that.degree));
for (int i = 0; i <= degree; i++)
poly.coef[i] = coef[i];
for (int i = 0; i <= that.degree; i++)
poly.coef[i] += that.coef[i];
poly.reduce();
return poly;
} public class01 minus(class01 that)
{
class01 poly = new class01(0, Math.max(degree, that.degree));
for (int i = 0; i <= degree; i++)
poly.coef[i] = coef[i];
for (int i = 0; i <= that.degree; i++)
poly.coef[i] -= that.coef[i];
poly.reduce();
return poly;
} public class01 times(class01 that)
{
class01 poly = new class01(0, degree + that.degree);
for (int i = 0; i <= degree; i++)
{
for (int j = 0; j <= that.degree; j++)
poly.coef[i + j] += (coef[i] * that.coef[j]);
}
poly.reduce();
return poly;
} public class01 compose(class01 that)// p(q(x)),用了秦九韶
{
class01 poly = new class01(0, 0);
for (int i = degree; i >= 0; i--)
{
class01 term = new class01(coef[i], 0);
poly = term.plus(that.times(poly));
}
return poly;
} public class01 differentiate() // p'(x)
{
if (degree == 0)
return new class01(0, 0);
class01 poly = new class01(0, degree - 1);
poly.degree = degree - 1;
for (int i = 0; i < degree; i++)
poly.coef[i] = (i + 1) * coef[i + 1];
return poly;
} public int evaluate(int x) // p(x) /. {x->a}
{
int p = 0;
for (int i = degree; i >= 0; i--)
p = coef[i] + (x * p);
return p;
} public int compareTo(class01 that)
{
if (degree < that.degree)
return -1;
if (degree > that.degree)
return +1;
for (int i = degree; i >= 0; i--)
{
if (coef[i] < that.coef[i])
return -1;
if (coef[i] > that.coef[i])
return +1;
}
return 0;
} @Override
public String toString()
{
if (degree == -1)
return "0";
if (degree == 0)
return "" + coef[0];
if (degree == 1)
return coef[1] + "x + " + coef[0];
String s = coef[degree] + "x^" + degree;
for (int i = degree - 1; i >= 0; i--)
{
if (coef[i] == 0)
continue;
if (coef[i] > 0)
s += " + " + (coef[i]);
if (coef[i] < 0)
s += " - " + (-coef[i]);
s += (i == 1) ? "x" : ("x^" + i);
}
return s;
} public static void main(String[] args)
{
class01 zero = new class01(0, 0);
class01 p1 = new class01(4, 3);
class01 p2 = new class01(3, 2);
class01 p3 = new class01(1, 0);
class01 p4 = new class01(2, 1);
class01 p = p1.plus(p2).plus(p3).plus(p4); // 4x^3 + 3x^2 + 1 class01 q1 = new class01(3, 2);
class01 q2 = new class01(5, 0);
class01 q = q1.plus(q2); // 3x^2 + 5 StdOut.println("zero(x) = " + zero);
StdOut.println("p(x) = " + p);
StdOut.println("q(x) = " + q);
StdOut.println("p(x) + q(x) = " + p.plus(q));
StdOut.println("p(x) * q(x) = " + p.times(q));
StdOut.println("p(q(x)) = " + p.compose(q));
StdOut.println("p(x) - p(x) = " + p.minus(p));
StdOut.println("0 - p(x) = " + zero.minus(p));
StdOut.println("p(3) = " + p.evaluate(3));
StdOut.println("p'(x) = " + p.differentiate());
StdOut.println("p''(x) = " + p.differentiate().differentiate());
}
}