2-SAT問題-拓撲列可行解

由对称性解2-SAT问题 - 博客园

演算法筆記

一、建立有向圖。
二、尋找所有Strongly Connected Component。
三、收縮每個Strongly Connected Component，成為有向無環圖DAG。
四、尋找縮圖的Topological Ordering。
五、在縮圖上，以Topological Ordering的反序，列出解。

* N个集团，每个集团2个人(0<=ni<2n p="p">

* 且每个集团只能选出一个人。如果两人有矛盾，他们不能同时被选中

* 问最多能选出多少人並列出

sample input:
3 4 (N M) (N組人M組條件)
1 2 (第1組)
3 4
5 6 (第3組)
1 3 (第1條件: 1不與3同時選上)
1 5
2 3
2 4

//reference: 演算法筆記 http://www.csie.ntnu.edu.tw/~u91029/SimulationModeling.html#a6
#include<stdio.h>
#include<string.h>
#include<iostream>
#define DBG 1

using namespace std;
int visit[20],map[20][20];
int sat[20];        
int finish[20], scc[20], n, N;
int f[20],x[20],y[20];

void DFS1(int i)
{
    visit[i] = true;
    for (int j=0; j<N+N; ++j)
        if (map[i][j] && !visit[j])
            DFS1(j);
    finish[n++] = i;
}
void DFS2(int i)
{
    scc[i] = n;
    if(DBG)printf("scc[%d]=%d\n",i+1,n);
    visit[i] = true;
    for (int j=0; j<N+N; ++j)
        if (map[j][i] && !visit[j]){
            DFS2(j);
        }
}
void two_satisfiability()
{
    n = 0;  // 時刻
    memset(visit, false, sizeof(visit));
    for (int i=0; i<N; ++i)
        if (!visit[i])
            DFS1(i);
            
    memset(visit, false, sizeof(visit));
    for (int k=N-1; k>=0; --k)
    {
        
        int i = finish[k];
        if (!visit[i])
            DFS2(i), n++;
    }
    // 檢查是否有解，無解則立即結束。
    for (int i=0; i<N; ++i)
        if (scc[i] == scc[f[i]])
            return;

    // 找出一組解
    memset(sat, 0, sizeof(sat));
    for (int k=0; k<N; ++k)
    {
        int i = finish[k];
        if (!sat[scc[i]])
        {
            sat[scc[i]] = 1;
            sat[scc[f[i]]] = 2;
        }
    }
    // 印出一組解
    for (int i=0; i<N; ++i)
        if (sat[scc[i]] == 1)
            cout << i+1 <<endl;
}

int main(){
    int m;
    while (scanf("%d%d",&N,&m)!=EOF && N+m){
        int p,q;
        for (int i=0;i<N;i++){
            scanf("%d%d",&p,&q),p--,q--;
            f[p]=q,f[q]=p;
        }
        for (int i=0;i<m;i++) scanf("%d%d",&p,&q),p--,q--,x[i]=p,y[i]=q;
        for (int i=0;i<m;i++){
            map[x[i]][f[y[i]]]=map[y[i]][f[x[i]]]=1;
        }
        N*=2;
        two_satisfiability();
    }
    return 0;
}

SCC (Strong Connected Component) Tarjan解

Strongly Connected Component 


在無向圖當中，只要邊邊相連，即形成連通，不必在意方向。在有向圖當中，由於有了方向限制，因此細分為「強連通」和「弱連通」：所有兩點之間，雙向都有路可通，則是「強連通」；所有兩點之間，至少單向有路可通，則是「弱連通」。

一張有向圖的「強連通分量」，是所有兩點之間，雙向皆有路可通的連通分量。一張有向圖的強連通分量們，不可能互相重疊。

Tarjan解 - DFS紀錄當前節點所連接到最高的祖先，有環則輸出整個環，無環則個別為一個Component

//reference: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<vector>
#define V 100
#define DBG 0

using namespace std;

vector<int> vec[V];
int id[V], pre[V], low[V], s[V], stop, cnt, scnt;
void tarjan(int v, int n) // vertex: 0 ~ n-1
{
    int t, minc = low[v] = pre[v] = cnt++;
    vector<int>::iterator pv;
    s[stop++] = v;
    for (pv = vec[v].begin(); pv != vec[v].end(); ++pv) {
        if(-1 == pre[*pv]) tarjan(*pv, n);
        if(low[*pv] < minc) minc=low[*pv];
    }
    if(minc < low[v]) {
        low[v] = minc; return;
    }
    do {
        id[t = s[--stop]] = scnt; low[t] = n;
    } while(t != v);
    ++scnt;
}
int main(){
    int a,b,n,m,i;
    memset(pre,-1,sizeof(pre));
    scanf("%d%d",&n,&m);
    for(int i=0;i<m;i++)scanf("%d%d",&a,&b),a--,b--,vec[a].push_back(b);
    for(i=0; i<n; ++i) if(-1==pre[i]) tarjan(i, n);
    printf("scc: %d\n",scnt);
}

最大團問題- Clique

wikepedia
團（clique）是圖論中的用語。對於給定圖G=(V,E)。其中，V={1,…,n}是圖G的頂點集，E是圖G的邊集。圖G的團就是一個兩兩之間有邊的頂點集合。如果一個團不被其他任一團所包含，即它不是其他任一團的真子集，則稱該團為圖G的極大團（maximal clique）。頂點最多的極大團，稱之為圖G的最大團（maximum clique）。最大團問題的目標就是要找到給定圖的最大團。
dynamic programming解:

//reference: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<iostream>
#include<vector>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define V 100               //maximum point

using namespace std;

int g[V][V], dp[V], stk[V][V], mx;
int dfs(int n, int ns, int dep){
    if (0 == ns) {
        if (dep > mx) mx = dep;
        return 1;
    }
    int i, j, k, p, cnt;
    for (i = 0; i < ns; i++) {
        k = stk[dep][i]; cnt = 0;
        //if (dep + n - k <= mx) return 0;
        if (dep + dp[k] <= mx) return 0;
        for (j = i + 1; j < ns; j++) {
            p = stk[dep][j];
            if (g[k][p]) stk[dep + 1][cnt++] = p;
        }
        dfs(n, cnt, dep + 1);
    }
    return 1;
}
int clique(int n){
    int i, j, ns;
        for (mx = 0, i = n - 1; i >= 0; i--) {
        // vertex: 0 ~ n-1
            for (ns = 0, j = i + 1; j < n; j++)
                if (g[i][j]) stk[1][ ns++ ] = j;
            dfs(n, ns, 1);
            dp[i] = mx;
        }
    return mx;
}
int main(){
    int N,M,a,b,res;        //N是點數，M是邊數
    cin>>N>>M;              
    for(int i=0;i<b;i++)cin>>a>>b,g[a][b]=1,g[b][a]=1;
    res=clique(N);
    printf("%d\n",res);
}

比較高效的大數運算

10000進位大數

//source: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<string.h>

const int base = 10000; // (base^2) fit into int
const int width = 4; // width = log base
const int N = 1000; // n * width: 
struct bint{
    int ln, v[N];
    bint (int r = 0) {  //to use string use read
        for (ln = 0; r > 0; r /= base) v[ln++] = r % base;
    }
    bint& operator = (const bint& r) {
        memcpy(this, &r, (r.ln + 1) * sizeof(int));// !
        return *this;
    }
};
bool operator < (const bint& a, const bint& b){
    int i;
    if (a.ln != b.ln) return a.ln < b.ln;
    for (i = a.ln - 1; i >= 0 && a.v[i] == b.v[i]; i--);
    return i < 0 ? 0 : a.v[i] < b.v[i];
}
bool operator <= (const bint& a, const bint& b){
    return !(b < a);
}
bint operator + (const bint& a, const bint& b){
    bint res; int i, cy = 0;
    for (i = 0; i < a.ln || i < b.ln || cy > 0; i++) {
        if (i < a.ln) cy += a.v[i];
        if (i < b.ln) cy += b.v[i];
        res.v[i] = cy % base; cy /= base;
    }
    res.ln = i;
    return res;
}
bint operator - (const bint& a, const bint& b){
    bint res; int i, cy = 0;
    for (res.ln = a.ln, i = 0; i < res.ln; i++) {
        res.v[i] = a.v[i] - cy;
        if (i < b.ln) res.v[i] -= b.v[i];
        if (res.v[i] < 0) cy = 1, res.v[i] += base;
        else cy = 0;
    }
    while (res.ln > 0 && res.v[res.ln - 1] == 0) res.ln--;
    return res;
}
bint operator * (const bint& a, const bint& b){
    bint res; res.ln = 0;
    if (0 == b.ln) { res.v[0] = 0; return res; }
    int i, j, cy;
    for (i = 0; i < a.ln; i++) {
        for (j=cy=0; j < b.ln || cy > 0; j++, cy/= base) {
            if (j < b.ln) cy += a.v[i] * b.v[j];
            if (i + j < res.ln) cy += res.v[i + j];
            if (i + j >= res.ln) res.v[res.ln++] = cy % base;
            else res.v[i + j] = cy % base;
        }
    }
    return res;
}
bint operator / (const bint& a, const bint& b)
{
    bint tmp, mod, res;
    int i, lf, rg, mid;
    mod.v[0] = mod.ln = 0;
    for (i = a.ln - 1; i >= 0; i--) {
        mod = mod * base + a.v[i];
        for (lf = 0, rg = base -1; lf < rg; ) {
            mid = (lf + rg + 1) / 2;
            if (b * mid <= mod) lf = mid;
            else rg = mid - 1;
        }
        res.v[i] = lf;
        mod = mod - b * lf;
    }
    res.ln = a.ln;
    while (res.ln > 0 && res.v[res.ln - 1] == 0) res.ln--;
    return res; // return mod 就是%運算
}
int digits(bint& a) // 返回位數
{
    if (a.ln == 0) return 0;
    int l = ( a.ln - 1 ) * 4;
    for (int t = a.v[a.ln - 1]; t; ++l, t/=10) ;
    return l;
}
bool read(char s[],bint& b, char buf[]) //
{
    if (1 != sscanf(s,"%s", buf)) return 0;
    printf("buf: %s\n",buf);
    int w, u, ln = strlen(buf);
    memset(&b, 0, sizeof(bint));
    if ('0' == buf[0] && 0 == buf[1]) return 1;
    for (w = 1, u = 0; ln; ) {
        u += (buf[--ln] - '0') * w;
        if (w * 10 == base) {w=1, b.v[b.ln++] = u,u=0;}
        else w *= 10;
    }
    if (w != 1) b.v[b.ln++] = u;
    return 1;
}
void write(const bint& v){
    int i;
    printf("%d", v.ln == 0 ? 0 : v.v[v.ln - 1]);
    for (i = v.ln - 2; i >= 0; i--)
    printf("%04d", v.v[i]); // ! 4 == width
    printf("\n");
}
int main(){
    char s1[]="12345678901234567890",
         s2[]="1234567890123456789",buf[100];
    bint bn1,bn2;
    read(s1,bn1,buf);
    read(s2,bn2,buf);
    write(bn1);
    write(bn2);
    bint bn3 = bn1*bn2;
    printf("bn3 = ");
    write(bn3);
}

RMQ - Region Minimum Query

演算法筆記- Region Minimum Query -
complexity: nlgn

#include<stdio.h>
#include<stdlib.h>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define N 100

//RMQ離線算法 O(N*logN)+O(1)
// INIT: val[]為待查詢數組,  initrmq(n);
int st[20][N], val[N];
int min(int a,int b){return a<b?a:b;}
void initrmq(int n){
    int i, j, k, sk;
    for (i = 0; i < n; i++) st[0][i] = val[i];
    for (i = 1, k = 2; k < n; i++, k <<= 1) {
        for (j = 0, sk = (k >> 1); j+k-1 < n; ++j, ++sk)
            st[i][j] = min(st[i-1][j],st[i-1][sk]);
    }
}
int query(int x, int y) // min of { val[x] ... val[y] }
{
    int bl = 31-__builtin_clz(y-x+1);
    return min(st[bl][x], st[bl][y-(1<<bl)+1]);
}
int main(){
    rep(i,N)val[i]=i;
    initrmq(N);
    rep(i,20){
        int c=rand()%60,d=rand()%100;
        while(d<c)d=rand()%100;
        printf("(%d %d): %d\n",c,d,query(c,d));
    }
}

suffix array - 使用倍增算法

演算法筆記-String Matching

suffix array ( 後綴陣列 )

使用倍增算法(Prefix Doubling)构造后缀数组

如果采用对每个后缀排序的方法来生成后缀数组，即使采用快速排序，由于每次比较的对象是字符串（设输入字符串的长度为n），因此每一个比较操作的复杂度不再是常数，而是O(n)，快速排序本身的平均情况下时间复杂度为O(nlgn)，因此总的时间复杂度是O(n^2*lgn)，如果考虑到采用快速排序最坏情况下复杂度为O(n^2)，那么最坏时总的复杂度为O(n^3)，显然在输入串长度很大时这种做法不可取。

Prefix Doubling算法（倍增法）是构造后缀数组一个比较实用的算法。其基本思想是先计算出每个后缀的k-前缀的rank值，然后在此基础上计算每个后缀的2k-前缀rank值，k从1开始。直到每个后缀都排出先后顺序为止（任何两个后缀都不会相等，换句话说，每个后缀终究都能排出先后顺序）。在处理2k-前缀时，只需要使用基数排序（radix sort）算法，先后对两位数字排序（可以采用计数排序算法(counting sort)对每一位数字排序）。在最坏情况下，需要做lgn次基数排序，每一次基数排序的操作次数为2*O(n)，因此它的时间复杂度是O(nlgn)。倍增法虽然没有达到像DC3算法的线性复杂度，但是它的优点是实现比较简单，因此在实践中常被采用。

在以下实现中，当k=1时，由于只需要对输入串的每个字符排序，因此在采用基数排序时，只有一位数字需要排序。当k>1时，需要对两位数字排序（为考虑通用性，代码中Tuple.digits数组长度可以取>=1的任何整数，而不限于两位数字的基数排序）。如果rank数组中某个后缀具有最大rank值，而且该rank值等于输入串的长度，这时说明每一个后缀都排出了次序，因而可以提前终止程序。另外，假设字母表中最大的字符为MAX_CHAR。

注：rank数组中的rank值对应SA数组的下标+1。例如T=mississippi#,那么SA=(11,10,7,4,1,0,9,8,6,3,5,2),这里第7个后缀（即ippi$）在SA数组中的下标是2，因此Rank[7]=3。

//reference: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<string>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define DBG 0
#define N 512
using namespace std;
//////////////////////////////////////////////////////
// 后綴數組 O(N * log N), prefix doubling algorithm
// INIT: n = strlen(s) + 1;    (利用\0作為分隔符號)
// CALL: makesa(); lcp();
//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
//sa: suffix array
//height: LCP value
//rank: sa之反函數，i,e., sa[rank[i]]= i
//tmp: 暫存最新之sa
//top: top[i]在迴圈中表示目前計算中排序第i名之字串前方(第1~i-1名)共有幾個字串
//////////////////////////////////////////////////////
int n,sa[N], height[N], rank[N], tmp[N], top[N];
char s[N];      // N > 256
void makesa(){  // O(N * log N)
    int i, j, len, na;
    na = (n < 256 ? 256 : n);
    memset(top, 0, na * sizeof(int));
    for (i = 0; i < n ; i++) top[ rank[i] = s[i]]++;
    for (i = 1; i < na; i++) top[i] += top[i - 1];          //counting sort
    for (i = 0; i < n ; i++) sa[ --top[ rank[i] ] ] = i;
    for (len = 1; len < n; len <<= 1) {             //len每次乘2
        if(DBG){rep(k,n)printf("%d ",sa[k]);printf("\n");}
        for (i = 0; i < n; i++) {                   //計算最新之sa, radix sort
            j = sa[i] - len;
            if(DBG)printf("j= %d - %d = %d(%d)\n",sa[i],len,j,j<0?j+n:j);
            if (j < 0) j += n;
            tmp[ top[ rank[j] ]++ ] = j;
            if(DBG)printf("tmp[%d]=%d\n",top[ rank[j] ]-1,j);
        }
        //此處sa被拿來暫存新的rank
        sa[ tmp[ top[0] = 0 ] ] = j = 0;
        for (i = 1; i < n; i++) {                   //計算top和rank
            if (rank[ tmp[i] ] != rank[ tmp[i-1] ] || rank[ tmp[i]+len ]!=rank[ tmp[i-1]+len ])
                top[++j] = i;
            sa[ tmp[i] ] = j;
        }
        if(DBG)printf("j= %d\n",j);
        memcpy(rank, sa , n * sizeof(int));
        memcpy(sa , tmp, n * sizeof(int));
        if (j >= n - 1) break;
    }
}
//lcp: 計算從每個字串之LCP (from s[0] to s[n-1])
//note: rank[0]>0, 因sa[0]=n (\0)
void lcp(){ // O(4 * N)
    int i, j, k;
    for (j = rank[height[i=k=0]=0]; i < n - 1; i++, k++){
        while (k >= 0 && s[i] != s[ sa[j-1] + k ]){
            height[j] = (k--), j = rank[ sa[j] + 1 ];
        }
    }
}
int main(){
    bool ll=0;
    char arr[5][100]={"ttttttt", "mississippi#", "abcdefghijklmmnopqrstuvwxyz#",
        "yabbadabbado#", "DFDLKJLJldfasdlfjasdfkldjasfldafjdajfdsfjalkdsfaewefsdafdsfa#"
        };
    rep(i,5){
        strcpy(s,arr[i]);
        n=strlen(s)+1;              //+1!
        if(ll)printf("********************************\n");ll=1;
        printf("Text: %s\n",s);
        makesa();
        lcp();
        printf("suffix array: \n");
        rep(i,n)printf(" %d",sa[i]);printf("\n");
        printf("rank array: \n");
        rep(i,n)printf(" %d",rank[i]+1);printf("\n");
        //printf("lcp: ");
        //rep(i,n)printf("%d ",height[i]);printf("\n");
    }
}

左偏樹

图解数据结构（9）——左偏树

十三、左偏树（Leftist Tree）
树这个数据结构内容真的很多，上一节所讲的二叉堆，其实就是一颗二叉树，这次讲的左偏树（又叫“左翼堆”），也是树。
二叉堆是个很不错的数据结构，因为它非常便于理解，而且仅仅用了一个数组，不会造成额外空间的浪费，但它有个缺点，那就是很难合并两个二叉堆，对于“合并”，“拆分”这种操作，我觉得最方面的还是依靠指针，改变一下指针的值就可以实现，要是涉及到元素的移动，那就复杂一些了。
左偏树跟二叉堆比起来，就是一棵真正意义上的树了，具有左右指针，所以空间开销上稍微大一点，但却带来了便于合并的便利。BTW：写了很多很多的程序之后，我发觉“空间换时间”始终是个应该考虑的编程方法。:)
左偏左偏，给人感觉就是左子树的比重比较大了，事实上也差不多，可以这么理解：左边分量重，那一直往右，就一定能最快地找到可以插入元素的节点了。所以可以这样下个定义：左偏树就是对其任意子树而言，往右到插入点的距离（下面简称为“距离”）始终小于等于往左到插入点的距离，当然了，和二叉堆一样，父节点的值要小于左右子节点的值。

//reference: ACM/ICPC代碼庫  http://ir.hit.edu.cn/~lfwang/docs/MyLib(For%20ACM).pdf
#include<stdio.h>
#include<algorithm>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define DBG 0

#define N 100

using namespace std;

#define typec int // type of key val
const int na = -1;
struct node { typec key; int l, r, f, dist; } tr[N];
int iroot(int i){ // find i's root
    if (i == na) return i;
    while (tr[i].f != na) i = tr[i].f;
    return i;
}
int merge(int rx, int ry){ // two root: rx, ry
    if (rx == na) return ry;
    if (ry == na) return rx;
    if(DBG)printf("merge %d %d\n",rx,ry);
    if (tr[rx].key > tr[ry].key)swap(rx, ry);
    int r = merge(tr[rx].r, ry);
    tr[rx].r = r; tr[r].f = rx;
    if (tr[r].dist > tr[tr[rx].l].dist)
    swap(tr[rx].l, tr[rx].r);
    if (tr[rx].r == na) tr[rx].dist = 0;
    else tr[rx].dist = tr[tr[rx].r].dist + 1;
    if(DBG)printf("root %d(%d)\n",rx,tr[rx].key);
    return rx; // return new root
}
int del(int i) { // delete node i
    if (i == na) return i;
    int x, y, l, r;
    l = tr[i].l; r = tr[i].r; y = tr[i].f;
    tr[i].l = tr[i].r = tr[i].f = na;
    tr[x = merge(l, r)].f = y;
    if (y != na && tr[y].l == i) tr[y].l = x;
    if (y != na && tr[y].r == i) tr[y].r = x;
    for ( ; y != na; x = y, y = tr[y].f) {
        if (tr[tr[y].l].dist < tr[tr[y].r].dist)
        swap(tr[y].l, tr[y].r);
        if (tr[tr[y].r].dist + 1 == tr[y].dist) break;
        tr[y].dist = tr[tr[y].r].dist + 1;
        }
        if (x != na) return iroot(x); // return new root
        else return iroot(y);
    }
node top(int root){
    return tr[root];
}
node pop(int &root){
    node out = tr[root];
    int l = tr[root].l, r = tr[root].r;
    tr[root].l = tr[root].r = tr[root].f = na;
    tr[l].f = tr[r].f = na;
    root = merge(l, r);
    return out;
}
int ins(int i, typec key, int &root){ // add a new node(i, key)
    tr[i].key = key;
    tr[i].l = tr[i].r = tr[i].f = na;
    tr[i].dist = 0;
    if(DBG)printf("ins %d %d\n",i,key);
    return root = merge(root, i); // return new root
}
int add(int i, typec val) // tr[i].key += val
{
    if (i == na) return i;
    if (tr[i].l == na && tr[i].r == na && tr[i].f == na) {
        tr[i].key += val;
        return i;
    }
    typec key = tr[i].key + val;
    int rt = del(i);
    return ins(i, key, rt);
}
void init(int n){
    for (int i = 1; i <= n; i++) {
        scanf("%d", &tr[i].key); //%d: type of key
        tr[i].l = tr[i].r = tr[i].f = na;
        tr[i].dist = 0;
    }
    if(DBG)printf("end of init\n");
}
int main(){
    int root=1;
    node n;
    init(1);
    REP(i,2,11)ins(i,i,root);
    rep(i,10)printf("pop %d\n",pop(root).key);
}

中國餘數定理

wiki

在中國古代著名數學著作《孫子算經》卷下第28题，叫做“物不知數”，原文如下：

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

即，一個整數除以三餘二，除以五餘三，除以七餘二，求這個整數。

克罗内克符号

为了便于表述，对任意的正整数用一个常用函数表示，称之为克罗内克符号(Kronecker).定义：

使用该符号，即可给出上述一般同餘方程组的求解过程，分两步完成

• 对每个，先求出正整数满足，即所求的满足条件：，但被每个整除。其求法如下：记，根据条件两两互素，可知和也互素，故存在整数和使得.令，则对每个，相对应的显然整除，并且。由此表明即为所求。
• 对于前述所求的，令，则，这说明为上述同餘方程组的一个解，从而所有的解可表示为，其中的n可以取遍所有的整数。

#include<stdio.h>
#define DBG 1

int w[3],b[3];
int extgcd(int a, int b, int & x, int & y){
    if (b == 0) { x=1; y=0; return a; }
    int d = extgcd(b, a % b, x, y);
    int t = x; x = y; y = t - a / b * y;
    return d;
}
int china(int k){
    int i, d, x, y, m, a = 0, n = 1;
    for (i = 0; i < k; i++) n *= w[i];
    for (i = 0; i < k; i++) {
        m = n / w[i];
        d = extgcd(w[i], m, x, y);
        a += (y * m * b[i]) % n;
        if(DBG)printf("a y m b %d %d %d %d\n",a,y,m,b[i]);
    }
    if (a > 0) return a;
    else return (a + n);
}
int main(){
    w[0]=3,w[1]=5,w[2]=7,b[0]=0,b[1]=3,b[2]=4;
    int ans = china(2);
    printf("ans: %d\n",ans);
}

質數測試算法(基于Miller-Rabin的MC算法)

素數測試算法(基于Miller-Rabin的MC算法)
wiki 

強偽素數
設n是一個大于4的奇整數，s和t是使得(n-1)=2^s*t的正整數，其中t為奇數，設B(n)是如下定義的整數集合：
a屬于集合B(n)當且僅當2≤a≤n-2且
1: a^tmodn=1
2: 存在整數i，0<=i<s滿足a^((2^i)*t) mod n=n-1
當n為素數時， 任意a在2和n-1中，均有a屬于集合B(n)
當n為合數時，若a屬于集合B(n)，則稱n為一個以a為底(基)的強偽素數，稱a為n素性的強偽證據。
n為素數，說明它對所有底均為強偽素數

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define REP(i, b, n) for (int i = b; i < n; i++)
#define rep(i, n) REP(i, 0, n)
#define DBG 0

using namespace std;

bool witness(int a,int n){
    int t=n-1,s=0,x=1;
    while(t%2==0)s++,t/=2;
    rep(i,t)x=(x*a)%n;
    if(x==n-1||x==1)return 0;
    rep(i,s-1){
        x=(x*x)%n;
        if(x==1)return 1;
        if(x==n-1)return 0;
    }
    return 1;
}
bool miller(int n,int s){
    if(n==2)return 1;
    if(n%2==0)return 0;
    rep(i,s){
        int a=rand()*(n-2)/RAND_MAX+1;
        if(witness(a,n))return 0;
    }
    return 1;
}
int main(){
    float t1,t2;
    REP(i,2,10001){
        if(miller(i,10))printf("%d\n",i);
    }
    t1 = clock();
    printf("time: (%f seconds)\n",(t1)/CLOCKS_PER_SEC);
}

Network Flow

source: http://pisces.ck.tp.edu.tw/~peng/index.php?action=showfile&file=f3cec71910d4a0106624e839f2891b17198ef58be

網路流 Network Flow

(Network Flow)是近年來在圖論中相當熱門的問題，在
1955 年 ，T.E. Harris 在研究鐵路最大通量時，首先提出在一個給定
的網路上尋求兩點間最大運輸量的問題。1956 年，L.R. Ford 和 D.R.
Fulkerson 給出解決這類問題的演算法，從而建立了網路流理論。
最大流問題的研究密切了圖論和運籌學，特別是與線性規劃的聯繫，開闢了圖論應用的新途徑。
關於網路流的定義 Definitions of Network Flow
1. 網路(Network)：圖 G = ( V, A )為一有向圖，稱為網路
2. 源點與匯點(Source and Sink)：令一點 S 為源點、一點 T 為匯點，其餘點則為中間點
3. 容量(Capacity)：每條弧上定義一個非負數 C(u, v)
4. 流量(Flow)：每條弧上定義一個 F(u, v) ，所有 。
5. 剩餘容量(Residual Capacity)：每條弧上定義一個 Cf(u, v) = C(u, v) – F(u, v) 為該弧的剩餘
容量，而剩餘容量的集合則稱為剩餘網路(Residual Network)
6. 網路的流量(Flow of Network)：由源點發出，匯點匯集的總流量 ，，， ，若其為該網路能產生的
量，則稱其為 量量量 (Maximum Flow)。 ...