▶ 书中第六章部分程序，加上自己补充的代码，包括快速傅里叶变换，多项式表示
● 快速傅里叶变换，需要递归
1 package package01; 2 3 import edu.princeton.cs.algs4.StdOut; 4 import edu.princeton.cs.algs4.StdRandom; 5 import edu.princeton.cs.algs4.Complex; 6 7 public class class01 8 { 9 private static final Complex ZERO = new Complex(0, 0); 10 private static final double EPSILON = 1e-8; 11 private class01() { } 12 13 public static Complex[] fft(Complex[] x)// 正向傅里叶变换 14 { 15 int n = x.length; 16 if (n == 1) 17 return new Complex[] {x[0]}; 18 if (n % 2 != 0) 19 throw new IllegalArgumentException("n != 2^k"); 20 21 Complex[] temp = new Complex[n/2], result = new Complex[n]; 22 for (int k = 0; k < n/2; k++) // 偶数项和奇数分别递归 23 temp[k] = x[2*k]; 24 Complex[] q = fft(temp); 25 for (int k = 0; k < n/2; k++) 26 temp[k] = x[2*k + 1]; 27 Complex[] r = fft(temp); 28 for (int k = 0; k < n/2; k++) // 合并 29 { 30 double kth = -2 * k * Math.PI / n; 31 Complex wk = new Complex(Math.cos(kth), Math.sin(kth)); 32 result[k] = q[k].plus(wk.times(r[k])); 33 result[k + n/2] = q[k].minus(wk.times(r[k])); 34 } 35 return result; 36 } 37 38 public static Complex[] ifft(Complex[] x)// 傅里叶反变换，共轭、正变换、再共轭，最后除以项数 39 { 40 int n = x.length; 41 Complex[] result = new Complex[n], temp; 42 for (int i = 0; i < n; i++) 43 result[i] = x[i].conjugate(); 44 result = fft(result); 45 for (int i = 0; i < n; i++) 46 result[i] = result[i].conjugate(); 47 for (int i = 0; i < n; i++) 48 result[i] = result[i].scale(1.0 / n); 49 return result; 50 } 51 52 public static Complex[] cconvolve(Complex[] x, Complex[] y)// 环形卷积 53 { 54 if (x.length != y.length) 55 throw new IllegalArgumentException("x.length != y.length"); 56 int n = x.length; 57 Complex[] a = fft(x), b = fft(y), c = new Complex[n]; 58 for (int i = 0; i < n; i++) 59 c[i] = a[i].times(b[i]); 60 return ifft(c); 61 } 62 63 public static Complex[] convolve(Complex[] x, Complex[] y)// 线性卷积，后一半垫 0 64 { 65 if (x.length != y.length) 66 throw new IllegalArgumentException("x.length != y.length"); 67 Complex[] a = new Complex[2 * x.length], b = new Complex[2 * y.length]; 68 for (int i = 0; i < x.length; i++) 69 a[i] = x[i]; 70 for (int i = x.length; i < 2*x.length; i++) 71 a[i] = ZERO; 72 for (int i = 0; i < y.length; i++) 73 b[i] = y[i]; 74 for (int i = y.length; i < 2*y.length; i++) 75 b[i] = ZERO; 76 return cconvolve(a, b); 77 } 78 79 private static void show(Complex[] x, String title) 80 { 81 StdOut.printf("%s ------------------- ",title); 82 for (int i = 0; i < x.length; i++) 83 StdOut.println(x[i]); 84 StdOut.println(); 85 } 86 87 public static void main(String[] args) 88 { 89 int n = 16; 90 Complex[] x = new Complex[n]; 91 for (int i = 0; i < n; i++) 92 x[i] = new Complex(StdRandom.uniform(-1.0, 1.0), 0); 93 94 show(x, "x"); 95 Complex[] y = fft(x); 96 show(y, "y = fft(x)"); 97 Complex[] z = ifft(y); 98 show(z, "z = ifft(y)"); 99 Complex[] c = cconvolve(x, x); 100 show(c, "c = cconvolve(x, x)"); 101 Complex[] d = convolve(x, x); 102 show(d, "d = convolve(x, x)"); 103 } 104 }
● 多项式表示
1 package package01; 2 3 import edu.princeton.cs.algs4.StdOut; 4 5 public class class01 6 { 7 private int[] coef; // 系数列表 8 private int degree; // 次数 9 10 public class01(int a, int b) // 建立单项式 a*x^b 11 { 12 coef = new int[b + 1]; 13 coef[b] = a; 14 reduce(); 15 } 16 17 private void reduce() // 记录多项式的次数 18 { 19 degree = -1; 20 for (int i = coef.length - 1; i >= 0; i--) 21 { 22 if (coef[i] != 0) 23 { 24 degree = i; 25 return; 26 } 27 } 28 } 29 30 public int degree() 31 { 32 return degree; 33 } 34 35 public class01 plus(class01 that) // p(x) + q(x)，把系数从低次项向高次项各加一遍 36 { 37 class01 poly = new class01(0, Math.max(degree, that.degree)); 38 for (int i = 0; i <= degree; i++) 39 poly.coef[i] = coef[i]; 40 for (int i = 0; i <= that.degree; i++) 41 poly.coef[i] += that.coef[i]; 42 poly.reduce(); 43 return poly; 44 } 45 46 public class01 minus(class01 that) 47 { 48 class01 poly = new class01(0, Math.max(degree, that.degree)); 49 for (int i = 0; i <= degree; i++) 50 poly.coef[i] = coef[i]; 51 for (int i = 0; i <= that.degree; i++) 52 poly.coef[i] -= that.coef[i]; 53 poly.reduce(); 54 return poly; 55 } 56 57 public class01 times(class01 that) 58 { 59 class01 poly = new class01(0, degree + that.degree); 60 for (int i = 0; i <= degree; i++) 61 { 62 for (int j = 0; j <= that.degree; j++) 63 poly.coef[i + j] += (coef[i] * that.coef[j]); 64 } 65 poly.reduce(); 66 return poly; 67 } 68 69 public class01 compose(class01 that)// p(q(x))，用了秦九韶 70 { 71 class01 poly = new class01(0, 0); 72 for (int i = degree; i >= 0; i--) 73 { 74 class01 term = new class01(coef[i], 0); 75 poly = term.plus(that.times(poly)); 76 } 77 return poly; 78 } 79 80 public class01 differentiate() // p'(x) 81 { 82 if (degree == 0) 83 return new class01(0, 0); 84 class01 poly = new class01(0, degree - 1); 85 poly.degree = degree - 1; 86 for (int i = 0; i < degree; i++) 87 poly.coef[i] = (i + 1) * coef[i + 1]; 88 return poly; 89 } 90 91 public int evaluate(int x) // p(x) /. {x->a} 92 { 93 int p = 0; 94 for (int i = degree; i >= 0; i--) 95 p = coef[i] + (x * p); 96 return p; 97 } 98 99 public int compareTo(class01 that) 100 { 101 if (degree < that.degree) 102 return -1; 103 if (degree > that.degree) 104 return +1; 105 for (int i = degree; i >= 0; i--) 106 { 107 if (coef[i] < that.coef[i]) 108 return -1; 109 if (coef[i] > that.coef[i]) 110 return +1; 111 } 112 return 0; 113 } 114 115 @Override 116 public String toString() 117 { 118 if (degree == -1) 119 return "0"; 120 if (degree == 0) 121 return "" + coef[0]; 122 if (degree == 1) 123 return coef[1] + "x + " + coef[0]; 124 String s = coef[degree] + "x^" + degree; 125 for (int i = degree - 1; i >= 0; i--) 126 { 127 if (coef[i] == 0) 128 continue; 129 if (coef[i] > 0) 130 s += " + " + (coef[i]); 131 if (coef[i] < 0) 132 s += " - " + (-coef[i]); 133 s += (i == 1) ? "x" : ("x^" + i); 134 } 135 return s; 136 } 137 138 public static void main(String[] args) 139 { 140 class01 zero = new class01(0, 0); 141 class01 p1 = new class01(4, 3); 142 class01 p2 = new class01(3, 2); 143 class01 p3 = new class01(1, 0); 144 class01 p4 = new class01(2, 1); 145 class01 p = p1.plus(p2).plus(p3).plus(p4); // 4x^3 + 3x^2 + 1 146 147 class01 q1 = new class01(3, 2); 148 class01 q2 = new class01(5, 0); 149 class01 q = q1.plus(q2); // 3x^2 + 5 150 151 StdOut.println("zero(x) = " + zero); 152 StdOut.println("p(x) = " + p); 153 StdOut.println("q(x) = " + q); 154 StdOut.println("p(x) + q(x) = " + p.plus(q)); 155 StdOut.println("p(x) * q(x) = " + p.times(q)); 156 StdOut.println("p(q(x)) = " + p.compose(q)); 157 StdOut.println("p(x) - p(x) = " + p.minus(p)); 158 StdOut.println("0 - p(x) = " + zero.minus(p)); 159 StdOut.println("p(3) = " + p.evaluate(3)); 160 StdOut.println("p'(x) = " + p.differentiate()); 161 StdOut.println("p''(x) = " + p.differentiate().differentiate()); 162 } 163 }