《算法》BEYOND 部分程序 part 2

Posted 2021-01-26 cuancuancuanhao

? 书中第六章部分程序，加上自己补充的代码，包括快速傅里叶变换，多项式表示

● 快速傅里叶变换，需要递归

  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 }