2014年12月27日 星期六

最大團Maximum Clique

好久沒有po文了,來講一下最大團的一派作法,
更改自Bron-Kerbosch演算法。

主要想法是,首先先照度數從大排到小,並將他狀態壓縮成int128,
也就是128bit的一個整數,要#include <limits.h>才會有的東西。

最naive的算法就是從第一個可以取的人開始枚舉放進clique裡頭,
但放過的人之後就不要再放了,所以將他從candi裡頭移除,可以看Without Pivot那部分。

第一個cut是算一下最多可以再放多少人+已經放多少人(potential),連ans都不到就cut掉。
第二個cut來自Bron-Kerbosch演算法,選出一個pivot來加速,詳細可以看pivot前面的那段註解。
第三個cut則在枚舉的裡頭,如果有個點他連出去的點被pivot連出去的點覆蓋了,我們就不需要去枚舉那個點。

基本上就是這樣子。


 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
int128 one = 1;
int128 linkto[MXN];

int popcount(int128 val){
 return __builtin_popcountll(val >> 50)
   + __builtin_popcountll(val & ((one << 50) - 1));
}

int lowbit(int128 val){
 if((val & ((one << 50) - 1)) == 0)
    return __builtin_ctzll(val >> 50) + 50;
 return __builtin_ctzll(val & ((one << 50) - 1));
}

/* candi should be sorted from big deg. to small deg. */

int ans;

void maxclique(int elem_num, int128 candi){
 if(elem_num > ans) ans = elem_num;

 int potential = elem_num + popcount(candi);
 if(potential <= ans) return;
 
 /*  Since Maximum Clique must contain >= 1 node not in Neighbor(pivot),
  *  Otherwise pivot can be added into the Clique.
  */

 int pivot = lowbit(candi);
 int128 smaller_candi = candi & (~linkto[pivot]);

 /* Without pivot:
  *
  *  while(candi && potential > ans){
  *  int next = lowbit(candi);
  *   
  *   candi ^= one << next;
  *   potential --;
  * 
  *   maxclique(elem_num + 1, candi & linkto[next]);
  *  }
  */

 while(smaller_candi && potential > ans){
  int next = lowbit(smaller_candi);
  
  candi ^= one << next;
  smaller_candi ^= one << next;
  potential --;

  /* The "later if" is 0,
   * iff "next" is connected to a subset of Neighbor(pivot).
   * Then it is suffice to pick pivot only.
   */

  if(next == pivot || (smaller_candi & linkto[next]))
   maxclique(elem_num + 1, candi & linkto[next]);
 }
}

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。兩個人在玩這個遊戲,問先手勝利還是後手。