问题1 最接近点对问题
(1) 1维
① n最接近点对问题,q一维最近点对问题,直线上的n个点采用排序+扫描方法,时间复杂度为O(nlogn)。
②蛮力算法(不推荐):两两比较, T(n)=(n-1)(n-2)/2=O(n2)。
③分治
时间复杂度为O(nlogn)
主要思路:
将一排仅有x坐标的点划分为两部分(快排中的partition),划分的数m=(m1+m2)/2,m1和m2是数组中的最大值和最小值,类似于中位数。
然后再划分的两边进行查找点对,这时可能会碰到一种情况是跨过中心点m,一边一个的情况,解决办法是找左边数组最大值p,右边数组最小值q,将q-p的值与左右两边找出的两个最小距离相比较。
伪代码:
c++代码:
// 一维
#include <iostream>
#include <stdio.h>
#include <cstdlib>
#include <ctime>
using namespace std;
#include <vector>
const int L = 100;
//点对结构体
struct Pair
{
float d;//点对距离
float d1, d2;//点对坐标
};
float Random();
int input(vector<float>& s);//构造S
float Max(vector<float>& s, int p, int q);
float Min(vector<float>& s, int p, int q);
template <class Type>
void Swap(Type& x, Type& y);
template <class Type>
int Partition(Type s[], Type x, int l, int r);
Pair Cpair(vector<float>& s, int l, int r);
// 一维 1d
//Ques20200927
int main()
{
srand((unsigned)time(NULL));
int m;
vector<float> s(L);
Pair d;
m = input(s); // m是length
d = Cpair(s, 0, m - 1);
cout << endl << "最近点对坐标为: (d1:" << d.d1 << ",d2:" << d.d2 << ")";
cout << endl << "这两点距离为: " << d.d << endl;
return 0;
}
float Random()
{
float result = rand() % 10000;
return result * 0.01;
}
int input(vector<float>& s)
{
int length;
cout << "输入点的数目: ";
cin >> length;
cout << "点集在X轴上坐标为:";
for (int i = 0; i < length; i++)
{
s[i] = Random();//产生随机数组元素
cout << s[i] << " ";
}
return length;
}
float Max(vector<float>& s, int l, int r)//返回s[]中的最大值
{
float s_max = s[l];
for (int i = l + 1; i <= r; i++)
if (s_max < s[i])
s_max = s[i];
return s_max;
}
float Min(vector<float>& s, int l, int r)//返回s[]中的最小值
{
float s_min = s[l];
for (int i = l + 1; i <= r; i++)
if (s_min > s[i])
s_min = s[i];
return s_min;
}
template <class Type>
void Swap(Type& x, Type& y)
{
Type temp = x;
x = y;
y = temp;
}
int Partition(vector<float>& s, float x, int l, int r)//快速排序中的一次划分
{
int i = l - 1, j = r + 1;
while (true)
{
while (s[++i] < x && i < r);
while (s[--j] > x);
if (i >= j)
{
break;
}
Swap(s[i], s[j]);
}
return j;
}
//返回s[]中的具有最近距离的点对及其距离
Pair Cpair(vector<float>& s, int l, int r)
{
Pair min_d = { 99999,0,0 };//最短距离
if (r - l < 1) return min_d;
float m1 = Max(s, l, r), m2 = Min(s, l, r);
float m = (m1 + m2) / 2;//找出点集中的中位数
//将点集中的各元素按与m的大小关系分组
int j = Partition(s, m, l, r);
Pair d1 = Cpair(s, l, j), d2 = Cpair(s, j + 1, r);//递归,求左边最小距离,右边最小距离
float p = Max(s, l, j), q = Min(s, j + 1, r);//左边最大值,右边最小值
//返回s[]中的具有最近距离的点对及其距离
if (d1.d < d2.d) // d1比较小
{
if ((q - p) < d1.d) // 横跨的q - p比d1还要小
{
min_d.d = (q - p);
min_d.d1 = q;
min_d.d2 = p;
return min_d;
}
else return d1;
}
else
{
if ((q - p) < d2.d)
{
min_d.d = (q - p);
min_d.d1 = q;
min_d.d2 = p;
return min_d;
}
else return d2;
}
}
(2)二维点对
二维复杂些,同样用分治,求S1、S2最近点对用x-坐标的中位数去划分点集
令d1、d2为S1、S2中最接近点对的距离,
但是并不需要像一维那样每个都去遍历, 那样复杂度就是O(n^2)了。
压缩最近点对搜索范围
分割直线L:x=m,m为S中各点x坐标的中位数
d=min(d1,d2)
当是一边一个的情况时
真正的距离<d,则一边一个,假设p,q点,那么p和q距直线L的距离均小于d,三角形3边关系,很容易看出来。
p q点关系
由d的意义可知,S2中任何2个S中的点的距离都不小于d。可以得到如下图红框所示的矩形Sp,由此可以推出矩形Sp中最多只有6个S中的点。
事实上,我们可以将矩形Sp的长为2d的边3等分,将它的长为d的边2等分,由此导出6个(d/2)×(2d/3)的矩形
将Sp按纵向三等分、横向2等分划分为6个方格,每个方格中至多有一个S2中的点 ,搜索范围中至多含有S2中的6个点。
对于S1中任一点p,P2中最多只有6个点与它构成最接近点对的候选者。
分治法的合并步骤中,我们最多只需要检查6×n/2=3n对候选者,而不是n^2/4对候选者。
伪代码:
C++代码:
分为按x坐标增序的X数组,y坐标增序的Y数组,以及用来最终存储全部点的Z数组(按区域划分的)。
递归开始,出口是只有两个点或三个点时。
递归体开始,m=(left+right)/2
两个指针分别指向left和m+1
for循环一遍left到right,Z:记录左右两个区中有哪些点。左Z[l:m],右Z[m+1:r],按y坐标从小到大,若它的x坐标>m,分到右半个区.若它的x坐标<m,分到左半个区。
递归求左边的最近的两个点a和b及距离d
递归求右边的最近的两个点ar和br及距离dr
比较两边得到的最小距离dr < d,更新a,b,d。
Merge(Z, Y, l, m, r);//重构数组Y,按点的y坐标大小,记录到Y数组
扫描Y数组中每一个点,将矩形内的点置于Z中,
再搜索Z[l:k-1],查找是否有距离更小的q点。
//最接近点对
//二维
//2d10-2 二维最邻近点对问题
#include <iostream>
#include <stdio.h>
#include <cstdlib>
#include <ctime>
#include<cmath>
using namespace std;
const int M = 50;
//用类PointX和PointY表示依x坐标和y坐标排好序的点
class PointX {
public:
int operator<=(PointX a)const
{
return (x <= a.x);
}
int ID; //点编号
float x, y; //点坐标
};
class PointY {
public:
int operator<=(PointY a)const
{
return(y <= a.y);
}
int p; //同一点在数组x中的坐标
float x, y; //点坐标
};
float Random();
template <class Type>
float dis(const Type& u, const Type& v);
bool Cpair2(PointX X[], int n, PointX& a, PointX& b, float& d);
void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX& a, PointX& b, float& d);
template <typename Type>
void Copy(Type a[], Type b[], int left, int right);
template <class Type>
void Merge(Type c[], Type d[], int l, int m, int r);
template <class Type>
void MergeSort(Type a[], Type b[], int left, int right);
// 二维 2d
int main()
{
srand((unsigned)time(NULL));
int length;
cout << "请输入点对数:";
cin >> length;
PointX X[M];
cout << "随机生成的二维点对为:" << endl;
for (int i = 0; i < length; i++)
{
X[i].ID = i;//点的编号从0开始到length-1
X[i].x = Random();
X[i].y = Random();
cout << "(" << X[i].x << "," << X[i].y << ") ";
}
PointX a;
PointX b;
float d;
Cpair2(X, length, a, b, d);
cout << endl;
cout << "最邻近点对为:(" << a.x << "," << a.y << ")和(" << b.x << "," << b.y << ") " << endl;
cout << "最邻近距离为: " << d << endl;
return 0;
}
//平面上任意两点u和v之间的距离可计算如下
template <class Type>
inline float dis(const Type& u, const Type& v)
{
float dx = u.x - v.x;
float dy = u.y - v.y;
return sqrt(dx * dx + dy * dy);
}
bool Cpair2(PointX X[], int n, PointX& a, PointX& b, float& d)
{//共n个点
if (n < 2) return false;
PointX* tmpX = new PointX[n];
MergeSort(X, tmpX, 0, n - 1);//按点的x坐标从小到大排序
PointY* Y = new PointY[n];
for (int i = 0; i < n; i++) //将数组X中的点复制到数组Y中
{
Y[i].p = i;
Y[i].x = X[i].x;
Y[i].y = X[i].y;
}
PointY* tmpY = new PointY[n];
MergeSort(Y, tmpY, 0, n - 1);//按点的y坐标排序
PointY* Z = new PointY[n];
closest(X, Y, Z, 0, n - 1, a, b, d);
delete[]Y; //
delete[]Z;
delete[]tmpX;
delete[]tmpY;
return true;
}
void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX& a, PointX& b, float& d)
{
if (r - l == 1) //两点的情形
{
a = X[l];
b = X[r];
d = dis(X[l], X[r]);
return;
}
if (r - l == 2) //3点的情形
{
float d1 = dis(X[l], X[l + 1]);
float d2 = dis(X[l + 1], X[r]);
float d3 = dis(X[l], X[r]);
if (d1 <= d2 && d1 <= d3)
{
a = X[l];
b = X[l + 1];
d = d1;
return;
}
if (d2 <= d3)
{
a = X[l + 1];
b = X[r];
d = d2;
}
else {
a = X[l];
b = X[r];
d = d3;
}
return;
}
//多于3点的情形,用分治法
int m = (l + r) / 2;
int f = l, g = m + 1;
//在算法预处理阶段,将数组X中的点依x坐标排序,将数组Y中的点依y坐标排序
//算法分割阶段,将子数组X[l:r]均匀划分成两个不相交的子集,取m=(l+r)/2
//X[l:m]和X[m+1:r]就是满足要求的分割。
for (int i = l; i <= r; i++)
{//Z:记录左右两个区中有哪些点。左Z[l:m],右Z[m+1:r]
if (Y[i].p > m) Z[g++] = Y[i];//按y坐标从小到大,若它的x坐标>m,分到右半个区
else Z[f++] = Y[i];//若它的x坐标<m,分到左半个区
}
closest(X, Z, Y, l, m, a, b, d);//递归求左边的最近的两个点a和b及距离d
float dr;
PointX ar, br;
closest(X, Z, Y, m + 1, r, ar, br, dr);//递归求右边的最近的两个点ar和br及距离dr
if (dr < d)
{
a = ar;
b = br;
d = dr;
}
Merge(Z, Y, l, m, r);//重构数组Y,按点的y坐标大小,记录到Y数组
//d矩形条内的点置于Z中
int k = l;
for (int i = l; i <= r; i++)//扫描Y数组中每一个点
{
if (fabs(X[m].x - Y[i].x) < d)//2d宽
{
Z[k++] = Y[i];
}
}
//搜索Z[l:k-1]
for (int i = l; i < k; i++)
{
for (int j = i + 1; j < k && Z[j].y - Z[i].y < d; j++)//高2d矩形框
{
float dp = dis(Z[i], Z[j]);
if (dp < d)
{
d = dp;
a = X[Z[i].p];
b = X[Z[j].p];
}
}
}
}
template <class Type>
void Merge(Type c[], Type d[], int l, int m, int r)
{
int i = l, j = m + 1, k = l;
while ((i <= m) && (j <= r))
{
if (c[i] <= c[j])
{
d[k++] = c[i++];
}
else
{
d[k++] = c[j++];
}
}
if (i > m)
{
for (int q = j; q <= r; q++)
{
d[k++] = c[q];
}
}
else
{
for (int q = i; q <= m; q++)
{
d[k++] = c[q];
}
}
}
template <class Type>
void MergeSort(Type a[], Type b[], int left, int right)
{
if (left < right)
{
int i = (left + right) / 2;
MergeSort(a, b, left, i);
MergeSort(a, b, i + 1, right);
Merge(a, b, left, i, right);//合并到数组b
Copy(a, b, left, right);//复制回数组a
}
}
template <typename Type>
void Copy(Type a[], Type b[], int left, int right)
{
for (int i = left; i <= right; i++)
a[i] = b[i];
}
时间复杂度分析
来源:国科大L老师PPT
问题2 快速Fourier变换
一些定义
(1)蛮力法
复杂度
计算复杂度:5句做i次乘,4句循环共1+2+…+(2n-1)=2n(2n-1)/2次乘,2句循环2n次,所以T(n)=O(n3)
(2)改进Fourier变换
避免的重复计算,考虑如下代数变换:
对每一个 顺序计算即可。
复杂度
(3)分治法
定义如下两个N次多项式:分别由a(x)偶、奇系数组成
复杂度
根据主定理2,或直接推导,T(n)=O(NlogN)
伪代码
C++代码
// 快速傅里叶变换
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<complex>
using namespace std;
#define cp complex<double>
#define ll long long
#define PI acos(-1.0)
#define MAXN 4000010
cp a[MAXN], b[MAXN], c[MAXN];
int n, m, lim;
inline ll read() {
ll s = 0, w = 1;
char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') w = -1;
for (; isdigit(c); c = getchar()) s = (s << 1) + (s << 3) + (c ^ 48);
return s * w;
}
cp omega(int n, int k) {
return cp{ cos(2 * PI * k / n), sin(2 * PI * k / n) };
}
void fft(cp* a, int n, bool inv) {
if (n == 1) return;
static cp buf[MAXN];
int m = n / 2;
for (register int i = 0; i < m; i++) {
buf[i] = a[2 * i];
buf[i + m] = a[2 * i + 1];
}
for (register int i = 0; i < n; i++)
a[i] = buf[i];
fft(a, m, inv);
fft(a + m, m, inv);
for (register int i = 0; i < m; i++) {
cp x = omega(n, i);
if (inv) x = conj(x);
buf[i] = a[i] + x * a[i + m];
buf[i + m] = a[i] - x * a[i + m];
}
for (register int i = 0; i < n; i++)
a[i] = buf[i];
}
// Ques20201010
int main() {
n = read(), m = read();
// n是 m是
for (register int i = 0; i <= n; i++)
a[i] = { (double)read(), 0 };
for (register int i = 0; i <= m; i++)
b[i] = { (double)read(), 0 };
int lim = 1;
while (lim <= n + m) lim *= 2;
for (int i = n + 1; i <= lim; i++) a[i] = { 0, 0 };
for (int i = m + 1; i <= lim; i++) b[i] = { 0, 0 };
fft(a, lim, true), fft(b, lim, true);
for (register int i = 0; i <= lim; i++)
c[i] = a[i] * b[i];
fft(c, lim, false);
for (register int i = 0; i <= n + m; i++)
printf("%d ", (int)((c[i].real() / lim) + 0.5));
return 0;
}