2013年6月18日 星期二

一道神祕題

給你兩個數字,A, B(A, B <= 10000)

你可以把他們各自拆成乘積:

33 -> 3, 11
24 -> 2, 3, 4
18 -> 2, 9

拆完的數字,由小排到大,計算差的平方和

EX:

33 40

-> 3, 11, 5, 8 -> 3, 5, 8, 11 -> (3-5)^2+(5-8)^2+(8-11)^2

= 22

而剛好這也就是所有拆法的「差的平方和」的最小值

給兩數字,求「差的平方和」的最小值

===========================

這種,近乎nHm的怪異題型,碰到就要先試試看DP

以某個數字X結尾的序列,A還剩A'沒拆,B還剩B'沒拆

最少的cost 為 DP[X][A'][B'],這樣。

可知A'會是A的因數,B'也是,X則是其一的因數

在1 ~ 10000裡,因數最多的數字也只有64個

所以DP的空間複雜度是64^3

至於要怎麼轉移,其實就看倒數第二是誰就行了

全部的複雜度是64^4。

但是,問題是他會問1000次:

1000 * 64^4 = 16777216000 = 10^10,就GG了




那要怎麼辦呢?

其實只要壓掉一個64就可以了

不是什麼性質,也沒有什麼怪東西

這裡不要想太多,簡單的DP優化而已 :)

http://ideone.com/hCStyp

2013年5月9日 星期四

APIO 2011 Table Coloring

題目乍看很像狀態壓縮DP,但其實有深遠的含意在裡頭

N, M, K <= 10^5, 看起來超級無敵瘋狂

Observation 1: 假設知道第一行、第一列,全部就都可以推出來了(2^(N+M-1)種)

a11 a12 a13
a21
a31

假設前三個是a, b, c,這一個 d 就要是a ^ b ^ c ^ 1,如此才會讓 a^b^c^d = 1
(也就是題目的定義)

Observation 2:

在5*4的表格上,會長的如下圖所示,恰好是 X ^ 該行 ^ 該列 ^ (bool)(偶行偶列)


因此,假設已經有填過了,就會變成 a 與 b 之間的條件限制(如果一個是1另一個就是0之類的),如同2-sat一般,而且邊是雙向的,就像以前討論過的題目,可以直接用disjoint set去實現。但比較麻煩的是有可能a1, a2這些單個直接被設限制,雖然依照目前為止的公式就可以求解,但編程複雜度會較高。

Observation 3:

接下來,我們希望把規則化簡。我想的作法是這樣,把第一行全部 xor x,如此後面一大堆的 xor x 都消失了,當然這不會改變任何東西,只是當b1是 1 時可能會變 0 罷了。再者,可以發現大家都是兩個東西在 xor ,為了統一化,把第一列全部再 xor 一個 y,變成如下圖所示。答案則是依此計算岀符合條件的變量種類後,除以二。


 如此一來規則就很簡單,直接給出code。
http://ideone.com/o9XFTY

2013年5月6日 星期一

Primality Testing // from CLRS

Prime Number Theorem: // page 965

pi(N) ~= N / ln(N)

According to Terence Tao, the proof is something like the below:

      Listening to the "music" of the primes. We start with a "sound wave" that is "noisy" at the prime numbers and silent at other numbers; this is the von Mangoldt functionThen we analyze its notes or frequencies by subjecting it to a process akin to Fourier transform; this is the Mellin transform. Then we prove, and this is the hard part, that certain "notes" cannot occur in this music. This exclusion of certain notes leads to the statement of the prime number theorem.

---------------------------------------------------------------------------------------------

Chinese Remainder Theorem: // page 950

n = PRODUCT( ni ) { ni is pairwise relatively prime }
a <---> (a1 ... ak) { ai = a mod ni } have bijection.

PF:
a --> (a1 ... ak), QED
(a1 ... ak) --> a, a = SUM( ai * mi * ( mi^(-1) mode ni ) ), mi = PRODUCT( nj except ni ), QED

Practice: X = 2 (mod 3), X = 3 (mod 5), X = 2 (mod 7), Find X


---------------------------------------------------------------------------------------------

Primality Testing 1: // page 956

if n is prime, only x = 1 or -1, x^2 = 1 (mod n)

PF:
if n = 2, 0^2 = 0, 1^2 = 1, 2^2 = 1, QED
if n > 2, p | (x-1)(x+1). p | (x-1) or p | (x+1), but not both.(Use Contradiction) If p | (x-1), x = 1(mod n). if p | (x+1), x = -1(mod n), QED

Corollary: if exist x != 1 & x != -1, x^2 = 1 (mod n), n is composite.(contrapositive)

---------------------------------------------------------------------------------------------

Primality Testing 2: // page 967

if n is prime, for all x != kn, x^(n-1) = 1 (mod n)

PF:

define Ai = x * i (mod n); P = PRODUCT( Ai ) { i = 1 ~ (n - 1) } ; Q = PRODUCT( i ){ i = 1 ~ (n-1) }

All Ai are distinct. ( Use Contradiction, if Ai = Aj (mod n), then x * (i - j) = 0 (mod n), then n | x * (i - j), then n | (i - j), but i - j < n, Contradict. )

Thus, P = Q (mod n)

x^(n-1) * Q = Q (mod n)

Q( x^(n-1) - 1 ) = 0 (mod n)

x^(n-1) = 1 (mod n)

Corollary: if exist x^(n-1) != 1 (mod n), n is composite.(contrapositive)
But there's something called Carmichael Number, that looks just like Prime under this test.

---------------------------------------------------------------------------------------------

Primality Testing 3: // page 968

For PT2, We use a^(n-1) = 1 (mod n). But here let a^(n-1) = a^(u*2^t), check a^u, (a^u)^2, ((a^u)^2)^2 ... 

So we won't be fooled by Carmichael Number.

And we can prove that if n is composite, then there are (n-1)/2 number a, that can discover it is composite.

http://ideone.com/FzgRvN

2013年4月29日 星期一

Πυθαγόρας

陳力!!!!!!!!!!!!!!!!!可以來看一下XD
之前那個神祕的東東的證明~

a^2 + b^2 = c^2 (a, b, c 是整數)

=> (a/c)^2 + (b/c)^2 = 1
---------------------------------
Let: x = a/c,y = b/c (x, y 是有理數)

=> x^2 + y^2 = 1

=> y * y = (1+x)(1-x)

=> y / (1+x) = (1-x) / y
---------------------------------
Let: t = y / (1+x) = (1-x) / y (t 是有理數)

x = 1 - yt,y = t + xt

=> x = (1-t^2) / (1+t^2), y = (2t) / (1+t^2)
---------------------------------
Let:t = u / v (u, v 是整數)

x = (v^2 - u^2) / (v^2 + u^2)
y = ( 2 * u * v ) / (v^2 + u^2)
---------------------------------
因為是比例,所以乘一個 r

a = (v^2 - u^2) * r
b = (2 * u * v) * r
c = (v^2 + u^2) * r
---------------------------------
放進任意正整數的 v, u, r 便可以得到所有的triple。(v > u)

但如果要找的是素勾股數(gcd(a, b, c) = 1):

不需要 r 這個數字,v 和 u 也要互質且必定一奇一偶
(這部份很容易證明,若符合任一條件則不是素勾股數)
---------------------------------
但要保證,剩下的全部都是素勾股數就比較困難

a = (v^2 - u^2) 奇
b = (2 * u * v) 偶
c = (v^2 + u^2) 奇

gcd(b, c)
= gcd(b, b + c) // b=gx, c=gy, x互質於y, 則x互質於x+y(cont.)
= gcd(2 * u * v, (v + u) ^ 2) // 用反證法證明上述論點
= gcd(u * v, (u + v) ^ 2) // 可以知道b + c是奇數

引理:A互質於C 且 B互質於C => A * B互質於C
Pf: 用反證,A=pa, B=qb, C=gc, p*q = g, 矛盾

u互質於v => u互質於u+v => u互質於(u+v)^2
同理可以得證  =========> v互質於(u+v)^2

因此 gcd(u * v, (u + v) ^ 2) = gcd(b, c) = 1
因此 gcd(a, gcd(b, c)) = gcd(a, b, c) = 1
---------------------------------
最後順便證「對於不同的(u, v)不可能會有相同的(a, b, c)」

先觀察一下
a = (v^2 - u^2)
b = (2 * u * v)
c = (v^2 + u^2)
可以知道c > a 且 c > b,但不保證b > a

因為我們希望的是a, b, c由小排到大,所以要證的如下:
若 (u, v) != (u', v'),則 (a, b, c) != (a', b', c') 且 (a, b, c) != (b', a', c')

這邊很簡單,亂搞一下就出來了~XD
話說,證岀來還滿爽的XDDDDDDDD

2013年4月13日 星期六

有趣的題目

自己亂創的題目:

       給你N堆石頭,a1, a2, a3, a4, a5, a6...,你可以選一群大小相同的石堆集合,一起拿掉1顆石頭,如: 1,3,5,5,5 -> 1,3,4,4,5 或 1,3,5,5,5 -> 1,3,4,5,5。兩個人在玩這個遊戲,問先手勝利還是後手。

2013年4月1日 星期一

最小圓覆蓋

好久沒有學習新東西了,今天因為剛好要用到 (poi的plot)

所以就順便學一下了,如果有敘述不清,歡迎糾正>_____<

=====================================


Lv 1. Naiwny Algorytm(樸素算法) O(n^4)


枚舉三個點,做出一個圓圈圈,看有沒有覆蓋所有的點點點點


Lv 2. Zhang Jiao(張角法) O(n^3)



       想法其實也很簡單,就是枚舉兩個點(如圖中的A, B)當圓周上的點,接下來,因為離這個圓越遠的角一定越小,所以找出在AB線段上面跟下面兩種的最小角點(也就是像圖中的C與D),如果這兩個角alpha + beta >= pi 的話,你就得到一個包含圓了(太小就變成一個怪怪的橢圓了),雖然要特別判一下,就是如果max(alpha, beta) >= pi,就可以把AB當成直徑,如果 < pi,就拿這三個點做圓(必定唯一)。

接下來就要進行大進化!

Lv 9999999. Venusaur(廟哇花) O(n)

       其實叫隨機增量法 A___A ,利用一個神定理和Random的特性,以最樸素作法少三個冪次的超神演算法。(雖然常數頗大)

神定理:假設{p1, p2 ... pi-1},形成一個最小包含圓,若pi不在這個圓內,那麼{p1, p2 ... pi}所產生的最小包含圓一定有pi在圓周上。

{p1 ~ pi-1} 產生一個最小包含圓:
加入 pi 之後(不是M_PI,只是p的第i項XDDD)
--------------------------------------------------------------------
1. pi 在圓內or圓上:爽!什麼毛都不用做
2. pi 在圓外:對於{p1 ~ pi}的最小包含圓,pi一定在圓上(依據神定理)

=====================================

至於他發生2.的這件事的機率是多少呢!因為1~i ,普通情況只會有三個點在圓上,所以機率是 3 / i !!!!!!!!!!! 接下來你輕鬆亂推廣神定理。

亂推神定理:假設{p1, p2 ... pi-1},且點 x 在圓周上,形成一個最小包含圓,若pi不在這個圓內,那麼{p1, p2 ... pi}且點 x 在圓周上,所產生的最小包含圓一定有pi在圓周上。

{p1 ~ pi-1} 且圓周上有點x,產生一個最小包含圓:
加入 pi 之後
--------------------------------------------------------------------
1. pi 在圓內or圓上:爽!什麼毛都不用做
2. pi 在圓外:對於{p1 ~ pi}且x在圓周上的最小包含圓,pi一定在圓上(依據亂推神定理)


=====================================


相同的,發生2.的機率還是3 / i 或更低。這時你已經確定有兩個點在圓周上了。
接下來,反正就再做一次亂推神定理,就有三個確定的點啦~~~~~~~~~

三個點形成一個圓~所以就結束了。

至於為什麼是O(n),請自行做點小小的計算吧。

2013年3月27日 星期三

2013年3月22日 星期五

2013年3月4日 星期一

USACO 2012 三月 金組

第一題:

關鍵:


這個數量約等於 nlgn,因為你在建表的時候,大致像下面那樣



for(int i = 1; i < N; i++)
        for(int j = i; j < N; j += i) P[j].PB(i); 

這樣子的數量級差不多是O(nlgn),因為 sum of 1/x 是 lgn。

所以:2^質因數個數 <= 因數的數量 <= nlgn

至於第一項是哪裡來的呢?是排容原理。

總之因數數量的總和其實是很小的,要注意這一個地方。

第二題:

相信大家都寫過,這題直鏈的版本,想通就解決了,較為簡單的一個題目。

第三題:

       我寫的是一個暴力搜的方法,並用目前答案去做剪枝,但這種作法在暴力方式是1*2*3*4*5*6*7...這樣子會比較有效。

       我一開始寫的暴力方式 n*n-1*n-2...,就會顯得沒有什麼意義,以後要多想由小到大的方法,會比較好Cut。

      官方題解不知道是在寫什麼毛,明天再去看看。

2013年2月28日 星期四

Codeforces 198 C. Delivering Carcinogen

Codeforces 的題目還蠻有品質保證的

很不錯的幾何題。不知道怎麼做就,二分搜

假設答案是 t 秒,那麼 > t 秒的都可以(先走到跟他一樣,而且V > Vp),利用這個性質來二分搜。

所以現在那個一直在繞圈圈的行星,被我們固定住

你只要找到所在位置與該行星的最短距,就解決了

直線永遠是最短距離,所以假如你不會穿越中間的圓圈,就直接拿他們的直線距離。

假設有穿越中間的圓圈,可以知道走切線再走圓弧再走切線,會是最短距離,

所以把他們一一算出來之後,就成功了。

=======================================

子問題-

1. 判斷有沒有穿越中間:直線與中心的距離 <= r 且OAB和OBA都是銳角

2. 直線與中心的距離:算外積後(平行四邊形面積),再除掉底 = 高 

3. 判斷兩個向量是否形成銳角:內積後,正的為銳角,負的為鈍角

4. 切線的長度:OA^2 - R^2

5. 圓弧的角度:設切點是C, D,設高是M,
若BOA是銳角: (BOM-BOD-AOC),若BOA是鈍角: (180 - BOM - BOD - AOC)

Codeforces 198E. Gripping Story

線段樹的變形真的很經典

這題要做的就是我們希望可以把 距離 <=  r 的區段

重量都小於等於力量的物質挑出來,放入queue裡

每一個物質都有兩個變量,距離以及重量

要對 <= 某距離內,找出 <= 某重量 的物質

這時,我們要用到一個很基本的概念:區段樹

就是一個線段樹,但是他的每一個節點,不是單單包含一個值

而是包含了這個區段的所有元素!

這麼看的話,感覺還挺崩潰的,但是可以知道每個元素只會被lgn個包含

所以元素數量是 NlgN ,非常的好用

有了這個方向之後,是否就很簡單了呢?答案是:Yes

只要讓那個區段的值都是由小排到大的,之後在query時

就從前面開始取出即可

(我是用Set去讓他由小到大的,整體複雜度nlg^2n,但只要質量小的放到質量大的,就可以用nlgn的時間,解出來了

想要看code,可以直接去codeforces,上面找。

2013年2月8日 星期五

Codeforces 266D BerDonalds

奇怪的題目名字,應該是學 Mcdonalds 的吧

如果他說只能在點上,就好辦了

直接做一次Floyd,就可以得到答案了

但!問題就在於你可以位於其中的無限多個點的其中一個

前面都很無腦,接下來就需要用到好腦袋惹~
(要了解很簡單,但要想到有點難==)

===================================

其實可以這樣想,先想對於一個邊你要到所有人的最短距離

就是,這邊的兩端為 a, b,min( x + dis[ a ][ i ], c - x + dis[ b ][ i ] )對所有 i 的max要最小

直接做枚舉可能答案的動作,是最簡單也是最方便的(比起二分搜或三分搜)

所以就先朝著枚舉答案的方向思考(也就是哪些中間點要試放看看這樣)

假設所有點的集合是這樣(a, b, c, d, e...),將他分成兩個集合A, B

A代表要讓這個邊的 a 去連會是min, B代表要讓這個邊的 b 去連會是min

也就是我們假定對A集合min(剛剛那團) 會是 a 的比較小,B也是同樣的道理

而若是這麼一回事的話,拿A的最大值,B的最大值,就可以組出想要的 x 了。

枚舉A的最大值,然後求B(補集)的最大值,弄出所有可能的X,再取最小的即可。

//By momo
#include
#include
#include
#define N 210
#define INF 999999999
using namespace std;

double ans = INF;
int n, m, dis[N][N];
int a[N*N], b[N*N], c[N*N];
struct pairr{ int a, b; }d[N];
bool comp(pairr x, pairr y){ return x.a > y.a; }

int main (){
    scanf("%d%d", &n, &m);
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++) dis[i][j] = INF;
        dis[i][i] = 0;
    }
    for(int i = 0; i < m; i++){
        scanf("%d%d%d", &a[i], &b[i], &c[i]);
        a[i]--, b[i]--;
        dis[a[i]][b[i]] = dis[b[i]][a[i]] = c[i];
    }
    for(int k = 0; k < n; k++)
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                dis[i][j] = min(dis[i][j],
                            dis[i][k] + dis[k][j]);
    for(int i = 0; i < m; i++){
        int C = c[i];
        for(int j = 0; j < n; j++){
            d[j].a = dis[a[i]][j];
            d[j].b = dis[b[i]][j];
        }
        sort(d, d+n, comp);
        int M = d[0].b; ans = min(ans, 1.0*d[0].a);
        for(int j = 1; j < n; j++){
            double val;
            if(d[j].a > M - C && M + C > d[j].a)
                val = (d[j].a+M+C) / 2.0;
            else val = max(d[j].a, M);
            ans = min(ans, val);
            M = max(M, d[j].b);
        }
    }
    printf("%lf\n", ans);
}

2013年2月5日 星期二

SGU 219 Synchrograph 6/10

想了一想,一個有>0值的環,或一個沒有入度的都可以不斷不斷的燃燒

環是可以循環燃燒,沒有入度的則可以無限燃燒

唯一會燒不起來的好像就是全部都是0的環所以就用0邊做出新的圖,用SCC判斷一下誰不能燒

而且如果有一個人的上游是無法亂燒的,那他就也不能亂燒(也就是所謂的alive)所以要再DFS一下,他的下游是誰

我想到一個比較好的解釋方式:
1. 對零的邊的SCC縮點,縮起來的那個點代表他無法燃燒
2. 因為他沒辦法燃燒,所以他下游的點也都不能燃燒
3. 但假設你把1,2,所說的點去除掉之後,剩下的點再做縮點,(因為環可以不斷的循環,所以可以視為是一個點)
4. 最後變成了DAG,一定存在一個無入度點,他可以無限燃燒,所以剩下的那些點,都是alive的。

//By momo
#include
#include
#include
#define N 1010
#define PB push_back
#define SZ(x) ((int)(x).size())
#define FOR(it, c) for(__typeof((c).begin())it=(c).begin();it!=(c).end();it++)
using namespace std;

int n, m, hit[N], cnt, vst[N];
vector<int> G[N], rG[N], tG[N], scc[N];
void dfs(int p){
    vst[p] = 1;
    FOR(it, G[p])
        if(!vst[*it]) dfs(*it);
    hit[cnt++] = p;
}
void rdfs(int p){
    vst[p] = 1;
    FOR(it, rG[p])
        if(!vst[*it]) rdfs(*it);
    scc[cnt].PB(p);
}
void kosaraju(){
    cnt = 0; fill(vst, vst + n, 0);
    for(int i = 0; i < n; i++)
        if(!vst[i]) dfs(i);
    cnt = 0; fill(vst, vst + n, 0);
    for(int i = n-1; i >= 0; i--)
        if(!vst[hit[i]]) rdfs(hit[i]), cnt++;
}

int dead[N];
void ddfs(int p){
    vst[p] = 1;
    FOR(it, tG[p])
        if(!vst[*it]) ddfs(*it);
    dead[p] = 1;
}

int main (){
    scanf("%d%d", &n, &m);
    for(int i = 0; i < m; i++){
        int a, b, c; scanf("%d%d%d", &a, &b, &c);
        a--, b--; if(c == 0) G[a].PB(b), rG[b].PB(a);
        if(a == b && c == 0) dead[a] = 1; tG[a].PB(b);
    }
    kosaraju();
    for(int i = 0; i < cnt; i++){
        if(SZ(scc[i]) == 1) continue;
        FOR(it, scc[i]) dead[*it] = 1;
    }
    fill(vst, vst+n, 0);
    for(int i = 0; i < n; i++)
        if(!vst[i] && dead[i]) ddfs(i);
    for(int i = 0; i < n; i++)
        printf("%d\n", 1 - dead[i]);
}

2013年2月4日 星期一

SGU 213 Strong Defence 6/10

還不錯的題目

先發現:
答案小於等於S-T的最短路長度L
(證明很簡單)

把從最短距X到X-1的設為一組邊集合

這樣總共有L個集合

因為X只能走到X或X-1,所以

這樣是一組合理解!

科科,證完了~

SGU 217 Two Cylinders 5/10


先做出積分式


接下來對其做數值模擬

因為精度的關係,這邊要使用一個稱為

『辛普森公式』的數值模擬方法

選三個點,用二次函數去fit他

接著就去算這曲線的面積,最後把這些東東全部加起來

寫的漂亮一點,可以弄成公式:

dx/3*(y0+4y1+2y2+…+2yn-2+4yn-1+yn)

//By momo
#include
#include
#include
#include
#define N 1000000
using namespace std;

double r1, r2, ans = 0;
double f(double x){
    return 8*sqrt((r1*r1-x*x)*(r2*r2-x*x));
}

int main (){
    scanf("%lf%lf", &r1, &r2);
    if(r1 > r2) swap(r1, r2);
    double dx = r1/N;
    for(int i = 0; i <= N; i++){
        double coe;
        if(i == 0 || i == N) coe = 1;
        else coe = (i&1)?4:2;
        ans += coe*dx*f(dx*i)/3;
    }
    printf("%.4lf\n", ans);
}

2013年1月29日 星期二

Codeforces 258B Little Elephant and Election 4/10

非常經典的一個題目

重點在他的第一個部分:

給一個  有N位 的數字 M,問1~M
有1個4 or 7的有幾個
有2個4 or 7的有幾個
有3個4 or 7的有幾個

這樣。

比較好的作法應該是

從最高位開始,從那位以前都一樣
但枚舉那一位的大小(小於原本的值)
再去枚舉剩下的位置有幾個4
最後再用nCm去計算~

這樣編程複雜度低,時間複雜度也很低

scanf("%d", &m), m++, cnt[0]--;
while(m) a[++n] = m % 10, m /= 10;
for (int i = 0; i <= 15; i++)
    for (int j = 0; j <= i; j++)
        C[i][j] = i && j? (C[i-1][j-1]+C[i-1][j]) % M: 1;
for (int i = n; i; i--){
    for (int j = 0; j < a[i]; j++){
        LL now = 1LL << i-1;
        for (int k = i-1; k >= 0; k--){
            cnt[prv+(j == 4 || j == 7)+k] += now * C[i-1][k];
            now *= 8/2;
        }
    }
    prv += (a[i] == 4 || a[i] == 7);
}