UVa 1718 - Tile Cutting

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution
    1. 4.1. FFT
    2. 4.2. NTT

Problem

請問在網格中,若頂點設置在網格邊上,構成平行四邊形的面積介於 $[l, r]$,最多的方法數和面積分別為多少。

Sample Input

1
2
3
2
4 4
2 6

Sample Output

1
2
48
620

Solution

由於是平行四邊形,假設矩形被劃分的兩個對稱三角形的底高 $a, b, c, d$,且 $a, b, c, d > 0$。

得到平行四邊形面積為

$$\begin{align*} \text{Area} &= (a + c)(b + d) - 2 \times 0.5 \times a b - 2 \times 0.5 \times c d \\ &= ad + bc \end{align*}$$

最後找到構成方法數為 ways of area = #(a, b, c, d) sat. ad + bc = Area。分別統計兩數相乘為 $x$ 的方法數 $W_x$,可以構成多項式 $W_1 x + W_2 x^2 + \cdots W_n x^n$,多項式相乘後,得到的係數就是所有的方法數。套用 FFT 的卷積計算即可。

FFT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <assert.h>
#include <algorithm>
#include <complex>
using namespace std;

struct Complex {
double x, y;
Complex(double x = 0, double y = 0): x(x), y(y) {}
Complex operator-(const Complex &v) const {
return Complex(x-v.x,y-v.y);
}
Complex operator+(const Complex &v) const {
return Complex(x+v.x,y+v.y);
}
Complex operator*(const Complex &v) const {
return Complex(x*v.x-y*v.y,x*v.y+y*v.x);
}
Complex operator/(const int &v) const {
return Complex(x/v,y/v);
}
double real() {
return x;
}
};

template<typename T> class TOOL_FFT {
public:
typedef unsigned int UINT32;
#define MAXN (1048576<<1)
Complex p[2][MAXN];
int pre_n;
T PI;
TOOL_FFT() {
pre_n = 0;
PI = acos(-1);
}
int NumberOfBitsNeeded(int PowerOfTwo) {
for (int i = 0; ; ++i) {
if (PowerOfTwo & (1 << i)) {
return i;
}
}
}
inline UINT32 FastReverseBits(UINT32 a, int NumBits) {
a = ( ( a & 0x55555555U ) << 1 ) | ( ( a & 0xAAAAAAAAU ) >> 1 ) ;
a = ( ( a & 0x33333333U ) << 2 ) | ( ( a & 0xCCCCCCCCU ) >> 2 ) ;
a = ( ( a & 0x0F0F0F0FU ) << 4 ) | ( ( a & 0xF0F0F0F0U ) >> 4 ) ;
a = ( ( a & 0x00FF00FFU ) << 8 ) | ( ( a & 0xFF00FF00U ) >> 8 ) ;
a = ( ( a & 0x0000FFFFU ) << 16 ) | ( ( a & 0xFFFF0000U ) >> 16 ) ;
return a >> (32 - NumBits);
}

void FFT(bool InverseTransform, vector<Complex>& In, vector<Complex>& Out) {
// simultaneous data copy and bit-reversal ordering into outputs
int NumSamples = In.size();
int NumBits = NumberOfBitsNeeded(NumSamples);
for (int i = 0; i < NumSamples; ++i) {
Out[FastReverseBits(i, NumBits)] = In[i];
}
// the FFT process
for (register int i = 1; i <= NumBits; i++) {
int BlockSize = 1<<i, BlockEnd = BlockSize>>1, BlockCnt = NumSamples/BlockSize;
for (register int j = 0; j < NumSamples; j += BlockSize) {
Complex *t = p[InverseTransform];
for (register int k = 0; k < BlockEnd; k++, t += BlockCnt) {
Complex a = (*t) * Out[k+j+BlockEnd];
Out[k+j+BlockEnd] = Out[k+j] - a;
Out[k+j] = Out[k+j] + a;
}
}
}
// normalize if inverse transform
if (InverseTransform) {
for (int i = 0; i < NumSamples; ++i) {
Out[i] = Out[i] / NumSamples;
}
}
}
void prework(int n) {
if (pre_n == n)
return ;
pre_n = n;
p[0][0] = Complex(1, 0);
p[1][0] = Complex(1, 0);
for (register int i = 1; i < n; i++) {
p[0][i] = Complex(cos(2*i*PI / n ) , sin(2*i*PI / n ));
p[1][i] = Complex(cos(2*i*PI / n ) , -sin(2*i*PI / n ));
}
}
vector<T> convolution(Complex *a, Complex *b, int n) {
prework(n);
vector<Complex> s(a, a+n), d1(n), d2(n), y(n);
vector<T> ret(n);
FFT(false, s, d1);
s[0] = b[0];
for (int i = 1, j = n-1; i < n; ++i, --j)
s[i] = b[j];
FFT(false, s, d2);
for (int i = 0; i < n; ++i) {
y[i] = d1[i] * d2[i];
}
FFT(true, y, s);
for (int i = 0; i < n; ++i) {
ret[i] = s[i].real();
}
return ret;
}
};
TOOL_FFT<double> tool;

Complex a[MAXN], b[MAXN];
vector<double> c;
long long ret[1048576];
long long pr[1048576] = {};
int main() {
int m;
const int N = 500000;
for (m = 1; m < (N<<1); m <<= 1);
/*
(a + c)(b + d) - 2 * 0.5 * (a b) - 2 * 0.5 * (c * d)
= ad + bc = area
a, b, c, d > 0,
ways of area = #(a, b, c, d) sat. ad + bc = area
*/

for (int i = 1; i <= N; i++) {
for (int j = i; j <= N; j += i)
pr[j]++;
}

memset(a, 0, sizeof(a[0]) * m);
memset(b, 0, sizeof(b[0]) * m);

for (int i = 1; i <= N; i++) {
a[i] = Complex(pr[i], 0);
b[m-i] = Complex(pr[i], 0);
}

c = tool.convolution(a, b, m);

for (int i = 1; i <= N; i++)
ret[i] = (long long) (c[i] + 0.5);

int testcase, L, R;
while (scanf("%d", &testcase) == 1) {
while (testcase--) {
scanf("%d %d", &L, &R);
assert(L <= R);
assert(L >= 1 && R <= 500000);
int mxA = L;
for (int i = L; i <= R; i++) {
if (ret[i] > ret[mxA])
mxA = i;
}
printf("%d %lld\n", mxA, ret[mxA]);
}
}
return 0;
}

NTT

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <assert.h>
#include <algorithm>
#include <complex>
using namespace std;

typedef unsigned int UINT32;
typedef long long INT64;
class TOOL_NTT {
public:
#define MAXN (1048576<<1)
const INT64 P = (479 << 21) + 1LL; // prime m = kn+1
const INT64 G = 3;
INT64 wn[22];
INT64 s[MAXN], d1[MAXN], d2[MAXN], y[MAXN];
TOOL_NTT() {
for (int i = 0; i < 22; i++)
wn[i] = mod_pow(G, (P-1) / (1<<i), P);
}
INT64 mod_mul(INT64 a, INT64 b, INT64 mod) {
return (a*b - (long long)(a/(long double)mod*b+1e-3)*mod+mod)%mod;
// INT64 ret = 0;
// for (a = a >= mod ? a%mod : a, b = b >= mod ? b%mod : b; b != 0; b>>=1, a <<= 1, a = a >= mod ? a - mod : a) {
// if (b&1) {
// ret += a;
// if (ret >= mod)
// ret -= mod;
// }
// }
// return ret;
}
INT64 mod_pow(INT64 n, INT64 e, INT64 m) {
INT64 x = 1;
for (n = n >= m ? n%m : n; e; e >>= 1) {
if (e&1)
x = mod_mul(x, n, m);
n = mod_mul(n, n, m);
}
return x;
}
int NumberOfBitsNeeded(int PowerOfTwo) {
for (int i = 0;; ++i) {
if (PowerOfTwo & (1 << i)) {
return i;
}
}
}
inline UINT32 FastReverseBits(UINT32 a, int NumBits) {
a = ( ( a & 0x55555555U ) << 1 ) | ( ( a & 0xAAAAAAAAU ) >> 1 ) ;
a = ( ( a & 0x33333333U ) << 2 ) | ( ( a & 0xCCCCCCCCU ) >> 2 ) ;
a = ( ( a & 0x0F0F0F0FU ) << 4 ) | ( ( a & 0xF0F0F0F0U ) >> 4 ) ;
a = ( ( a & 0x00FF00FFU ) << 8 ) | ( ( a & 0xFF00FF00U ) >> 8 ) ;
a = ( ( a & 0x0000FFFFU ) << 16 ) | ( ( a & 0xFFFF0000U ) >> 16 ) ;
return a >> (32 - NumBits);
}
void NTT(int on, INT64 *In, INT64 *Out, int n) {
int NumBits = NumberOfBitsNeeded(n);
for (int i = 0; i < n; ++i)
Out[FastReverseBits(i, NumBits)] = In[i];
for(int h = 2, id = 1; h <= n; h <<= 1, id++) {
for(int j = 0; j < n; j += h) {
INT64 w = 1, u, t;
int block = h/2, blockEnd = j + h/2;
for(int k = j; k < blockEnd; k++) {
u = Out[k], t = mod_mul(w, Out[k+block], P);
Out[k] = u + t;
Out[k + block] = u - t + P;
if (Out[k] >= P) Out[k] -= P;
if (Out[k+block] >= P) Out[k+block] -= P;
w = mod_mul(w, wn[id], P);
}
}
}
if (on == 1) {
for (int i = 1; i < n/2; i++)
swap(Out[i], Out[n-i]);
INT64 invn = mod_pow(n, P-2, P);
for (int i = 0; i < n; i++)
Out[i] = mod_mul(Out[i], invn, P);
}
}
void convolution(INT64 *a, INT64 *b, int n, INT64 *c) {
NTT(0, a, d1, n);
s[0] = b[0];
for (int i = 1; i < n; ++i)
s[i] = b[n-i];
NTT(0, s, d2, n);
for (int i = 0; i < n; i++)
s[i] = mod_mul(d1[i], d2[i], P);
NTT(1, s, c, n);
}
} tool;

INT64 a[MAXN], b[MAXN], c[MAXN];
long long ret[1048576];
long long pr[1048576] = {};
int main() {
int m;
const int N = 500000;
for (m = 1; m < (N<<1); m <<= 1);
/*
(a + c)(b + d) - 2 * 0.5 * (a b) - 2 * 0.5 * (c * d)
= ad + bc = area
a, b, c, d > 0,
ways of area = #(a, b, c, d) sat. ad + bc = area
*/

for (int i = 1; i <= N; i++) {
for (int j = i; j <= N; j += i)
pr[j]++;
}

memset(a, 0, sizeof(a[0]) * m);
memset(b, 0, sizeof(b[0]) * m);

for (int i = 1; i <= N; i++) {
a[i] = pr[i];
b[m-i] = pr[i];
}

tool.convolution(a, b, m, c);

for (int i = 1; i <= N; i++)
ret[i] = c[i];

int testcase, L, R;
while (scanf("%d", &testcase) == 1) {
while (testcase--) {
scanf("%d %d", &L, &R);
assert(L <= R);
assert(L >= 1 && R <= 500000);
int mxA = L;
for (int i = L; i <= R; i++) {
if (ret[i] > ret[mxA])
mxA = i;
}
printf("%d %lld\n", mxA, ret[mxA]);
}
}
return 0;
}