2016-04-30

UVa 1718 - Tile Cutting

1. Problem
2. Sample Input
3. Sample Output
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

#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

#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;
}

Morris' Blog

UVa 1718 - Tile Cutting

contents

Problem

Sample Input

Sample Output

Solution

FFT

NTT