HMPP演算法源碼

1. 什麼是algor ithm和程序

演算法是用來解決一個問題的程序或者公式。「演算法」（algorithm）這個詞源自生活在公元9世紀的波斯數學家Al-Khowarizmi。一個電腦程序可以被看成是一個精心設計的演算法。在數學和電腦學科，一個演算法通常意味著一個可以解決循環問題的小程序。
程序序（program）是為實現特定目標或解決特定問題而用計算機語言編寫的命令序列的集合。
程序（港台稱之為程式） program(me)或procere
chénɡ xù
概念1.：為進行某活動或過程所規定的途徑。
概念2.程序是由序列組成的，告訴計算機如何完成一個具體的任務。由於現在的計算機還不能理解人類的自然語言，所以還不能用自然語言編寫計算機程序。

2. 泰坦的軟體支持

參見：NVIDIA CUDA
不同於以往以中央處理器作主要數據處理單元或以圖形處理器輔助處理數據的超級電腦，泰坦超級電腦以圖形處理器為主要的數據處理單元。在初期，不少專案的程序源碼一般重新編寫或編譯以方便泰坦的通用圖形處理器運行，而且還要求這些源碼能兼容於以中央處理器為主的系統，以便專案可以不僅可以由泰坦進行處理，還能供其它仍使用傳統中央處理器架構的超級電腦運行。Oak Ridge Leadership Computing Facility（OLCF）還成立了Center for Accelerated Application Readiness（加速申請中心，CAAR）來幫助科學研究者們重寫或重新編譯他們所寫的代碼以供泰坦運作之，並且在英偉達總部舉辦開發者研習會，讓研究人員和程序員學習並使用泰坦和其它類似的超級電腦的架構、編譯器和應用程序。
以使研究者們可以和以前一樣，以Fortran、C或C++來編寫他們的代碼而無須學習新的編程語言，並且編譯器可以將這些代碼編譯為通用圖形處理器可以高效運作的程序。
「加速編譯器指令」包括OpenACC、PGI Accelerator、CAPS HMPP Workbench；「低級圖形處理器編程語言」有NVIDIA C for CUDA、PGI CUDA Fortran、OpenCL。這些工具中有部分甚至允許開發者在不了解關於NVIDIA CUDA的知識的情況下使用。在泰坦上適用的編譯器有PGI（Portland Group Compiler Suite，默認使用這個）、GCC、CCE（Cray Compiling Environment）以及英特爾提供的編譯器。
即所謂的多核多線程優化，只是更深入、更適合通用圖形處理器的架構特性（通用圖形處理器擁有比中央處理器多得多的線程和更好的運算並發性）。不過這樣的程序源碼優化也適用於傳統以中央處理器為基礎的機器上，這些機器也可以從中獲得性能提升。天體物理學家兼程序員Bronson Messer，表示「一個應用程序要將泰坦的性能發揮至極限，必須尋求方法來使通用圖形處理器滿負載運作，需要明確的是通用圖形處理器盡管處理速度很快，但靈活度顯然不如中央處理器。」作業管理者也發現，一些計劃專案的代碼在重寫後，在那些不是基於圖形處理器的機器上的效率也變得大有提升，「Denovo」專案的程序在基於中央處理器的機器上甚至還比原先獲得雙倍的性能提升表現。
要求調用圖形處理器的進行運算作業的代碼變更量因不同專案而變。根據負責NRDF專案的梅瑟博士（Dr. Messer）的意見，他們的代碼中調用圖形處理器的演算法源碼只佔整個代碼中很小的百分比，因為這些程序演算法的運算都相對簡單，但處理進程重復而且具有平行性。NRDF的程序用CUDA Fortran寫成，CUDA Fortran是Fortran語言鏈接CUDA指令擴充庫（NVIDIA CUDA庫）給予圖形處理器使用的Fortran改進版本。

3. 求教：蟻群演算法選擇最短路徑問題

這個例子其實是當初數模比賽時用來完成碎片拼接的，但其所用到原理還是求解最短路徑的原理。但這里的最短路徑和數據結構中最短路徑有一定的區別。在數據結構中，對於最短路徑的求解常用的一般有Dijkstra演算法與Floyd演算法，但對於要求出一條經過所有的點的並且要求路徑最短，這些演算法還是有一定的局限性的。而蟻群演算法則很好地滿足了這些條件。話說回來，很想吐槽一下網路流傳的一些蟻群演算法的例子，當初學習這個時候，身邊也沒有相關的書籍，只好到網上找例子。網上關於這個演算法源代碼的常見的有2個版本，都是出自博客，但是在例子都代碼是不完整的，缺失了一部分，但就是這樣的例子，居然流傳甚廣，我很好奇那些轉載這些源碼的人是否真的有去學習過這些，去調試過。當然，我下面的例子也是無法直接編譯通過的，因為涉及到圖像讀取處理等方面的東西，所以就只貼演算法代碼部分。但是對於這個問題蟻群演算法有一個比較大的缺點，就是收斂很慢，不過對於數量小的路徑，效果還是很好的。function bestqueue =aco1(nt,nc_max,m ,st, sd ,Alpha ,Beta ,Rho ,Q,gethead,getend)%參數解釋：%nt 路徑所經過的點的個數；%nc_max 迭代的次數；%m 螞蟻的個數；%st 起點序號；%sd 終點序號；%Alpha 信息素系數；�ta 啟發因子系數；%Rho 蒸發系數；% Q 信息量；%gethead getend 是用來求距離矩陣的，可根據實際情況修改
% nt = 209;%碎片個數full = zeros(nt,nt);tic;%初始化距離矩陣for i =1:nt for t = 1:nt if i ~= t full(i,t) = sum(abs(getend(:,i) - gethead(:,t))); else full(i,t) = inf; end endend% a =full(156,187)eta = 1./full;%啟發因子，取距離的倒數% eta% e = eta(4,2)tau = ones(nt,nt);%信息素矩陣% tabu = zeros(nt,nt);%禁忌矩陣，取螞蟻數量和碎片數量一致，以減少迭代次數nc =1;%初始化迭代次數;rbest=zeros(nc_max,nt);%各代最佳路線rbest(:,1) = (linspace(st,st,nc_max))';rbest(:,nt) =(linspace(sd,sd,nc_max))'; lbest=zeros(nc_max,1);%各代最佳路線的長度pathlen = 0;%臨時記錄每代最佳路線長度stime = 1;%記錄代數進度for i = 1:nc_max % 代數循環 delta_tau=zeros(nt,nt);%初始化改變數 stime for t = 1:m % 對螞蟻群體的循環， tabu=zeros(1,nt);%禁忌向量，標記已訪問的碎片，初試值設為0，訪問之後則變為1; viseted = zeros(1,nt);%記錄已訪問的元素的位置 tabu(st) = 1;%st為起點，在此表示為碎片矩陣的編號，因為已經將蟻群放在起點，故也應將禁忌向量和位置向量的狀態進行修改 tabu(sd) =1;%同上 visited(nt) = sd ;%同上； visited(1) = st;%同上； ht = 0; for r = 2:nt-1 %記錄了還沒訪問的圖片編號 vp = 1;%visited指示量 pp = [];%置空的概率向量 jc = 0; %獲取尚未訪問的位置的向量。 wv = zeros( nt -2 - ht ); for k =1 : nt if tabu(k) == 0 jc = jc +1; wv(jc) = k; end end% a =(tau(visited(end),ju(3))^Alpha)*(eta(visited(end),ju(3))^Beta)% visited(end) %計算選擇的概率 for k=1:length(wv) pp(k)=(tau(visited(vp),wv(k))^Alpha)*(eta(visited(vp),wv(k))^Beta);%下一張碎片的選擇概率計算，p =（信息素^信息素系數）*（啟發因子^啟發因子系數） end pp=pp./(sum(pp));%歸一化 pcum =cumsum(pp); psl = find(pcum >= rand);%輪盤賭法 to_visit= wv(psl(1)) ;%完成選點 tabu(to_visit) =1; visited(r) = to_visit; ht =ht +1;%已訪問碎片個數變化 vp =vp+1; end %路徑變化信息 %對單個螞蟻的路徑進行統計 sum1 =0; for pr = 1:nt -1 x = visited(pr); y = visited(pr+1) ; sum1 =sum1 + full(x,y); end% vcell{t} =visited;%元胞記錄每個螞蟻的路徑，即碎片順序;% msum(t) = sum1; %信息素變化； for ww=1:(nt-1) delta_tau(visited(ww),visited(ww+1))=delta_tau(visited(ww),visited(ww+1)) + Q/sum1; end% delta_tau(visited(end),visited(1))=delta_tau(visited(end),visited(1))+Q/(sum1/100);% if t == m & i == nc_max % bestqueue = visited% end if t == m bestqueue = visited end end tau=(1-Rho).*tau+delta_tau; %完成信息素的更新，找出現有的最新的最佳路徑,即信息素最多的路徑； stime =stime +1;end toc;

4. 求五子棋AI演算法

在本次「五子棋「程序的編寫中，只編寫了人機對弈部分，運用了博弈樹進行搜索，在選取最優的走步時使用極大極小分析法，考慮到搜索的時間復雜度和空間復雜度，在程序中只進行了2步搜索，即計算機在考慮下一步的走法時，只對玩家進行一步的推測。（程序中的棋盤規格為15*15）
下面對具體做法進行描述：
1． 數據結構定義：
棋盤定義：int board[15][15]；
在15*15的棋盤上，獲勝的情況總共有572種，
如：
* * * * * ……
…… …… …… …… …… ……
…… …… …… …… …… ……
…… …… …… …… …… ……
…… …… …… …… …… ……
…… …… …… …… …… ……
中的第一行「*「所代表的格子就是一種獲勝組合。
計算機和玩家的獲勝組合情況bool ctable[15][15][572],
bool ptable[15][15][572]，來表示棋盤上的各個位置都在那種獲勝組合中。
計算機和玩家在各個獲勝組合中所填入的棋子數int win[2][572],如有一方在某一獲勝組合的棋子數達到5個，該方即獲勝。
Bool player:是否輪到玩家下棋
Bool computer:是否輪到計算機下棋
Bool start:游戲是否開始
Bool pwin:玩家是否獲勝
Bool cwin:計算機是否獲勝
CPoint m_pplastpos;//玩家走的前一步棋
CPoint m_pclastpos;//計算機走的前一步棋
為便於說明程序的主要演算法，這里先說本程序中估價函數的選取方法：
e=p1+p2;
p1為下完當前這步棋時計算機的得分；p2為下完當前這步棋時玩家的得分（p2其實為負），這樣做即考慮了進攻的因數，由考慮了防守的因數，兩個方面都進行了考慮，防止計算機只考慮進攻而忽略防守，同時也防止計算機只考慮防守而忽略進攻，從而達到比較好的情況。
2．主要流程描述
其程序流程圖如下：

（1） 初始化棋盤：判斷哪方先開始，（2） 初始化計算機和玩家的獲勝組合情況
bool ctable[15][15][572],bool ptable[15][15][572]
void CMyChessDlg::InitializeBoard()
{
//初始時雙方都還沒下子
int i,j,count=0,k;
m_pclastpos.x=-1;
m_pclastpos.y=-1;
m_pplastpos.x=-1;
m_pplastpos.y=-1;
start=true;
//判斷哪方先開始
if(m_bwfirst)
{
player=false;
computer=true;
}
else
{
player=true;
computer=false;
}
pwin=cwin=false;
//初始化計算機和玩家的獲勝組合情況
for(i=0;i<15;i++)
for(j=0;j<15;j++)
for(k=0;k<572;k++)
{
ptable[i][j][k]=false;
ctable[i][j][k]=false;
}
for(i=0;i<2;i++)
for(j=0;j<572;j++)
win[i][j]=0;
for(i=0;i<15;i++)
for(j=0;j<15;j++)
board[i][j]=2;
for(i=0;i<15;i++)
for(j=0;j<11;j++)
{
for(k=0;k<5;k++)
{
ptable[j+k][i][count]=true;
ctable[j+k][i][count]=true;
}
count++;
}
for(i=0;i<15;i++)
for(j=0;j<11;j++)
{
for(k=0;k<5;k++)
{
ptable[i][j+k][count]=true;
ctable[i][j+k][count]=true;
}
count++;
}
for(i=0;i<11;i++)
for(j=0;j<11;j++)
{
for(k=0;k<5;k++)
{
ptable[j+k][i+k][count]=true;
ctable[j+k][i+k][count]=true;
}
count++;
}
for(i=0;i<11;i++)
for(j=14;j>=4;j--)
{
for(k=0;k<5;k++)
{
ptable[j-k][i+k][count]=true;
ctable[j-k][i+k][count]=true;
}
count++;
}
}
（3） 給出下了一個子後的分數：
int CMyChessDlg::GiveScore(int type, int x, int y)
{
int i,score=0;
for(i=0;i<572;i++)
{
//計算機下
if(type==1)
{
if(ctable[x][y][i])
{
switch(win[1][i])
{
case 1:
score+=5;
break;
case 2:
score+=50;
break;
case 3:
score+=100;
break;
case 4:
score+=10000;
break;
default:
break;
}
}
}
//人下
else
{
if(ptable[x][y][i])
{
switch(win[0][i])
{
case 1:
score-=5;
break;
case 2:
score-=50;
break;
case 3:
score-=500;
break;
case 4:
score-=5000;
break;
default:
break;
}
}
}
}
return score;
}
（4） 核心程序 ，（5） 即計算機如何運用極小極大法分析，（6） 選取最優走法，（7）
在程序中極小極大法即體現在兩個while循環前後及之間的內容，其估價函數的體現在ctemp+pscore中
void CMyChessDlg::ComTurn()
{
//bestx,best為當前最佳位置，i，j是人能下的各種位置；pi，pj是計算機能下的各種位置
Int bestx,besty,i,j,pi,pj,ptemp,ctemp,pscore=10,cscore=-10000,
ctempboard[15][15],ptempboard[15][15];
int m,n,temp1[20],temp2[20];//暫存第一步搜索的信息
if(start)
{
if(board[7][7]==2)
{
bestx=7;
besty=7;
}
else
{
bestx=8;
besty=8;
}
start=false;
}
else
{//尋找最佳位置
GetBoard(ctempboard,board);
while(SearchBlank(i,j,ctempboard))
{//進行第一步查找
n=0;
pscore=10;
GetBoard(ptempboard,board);//獲取當前棋盤狀態
ctempboard[i][j]=3;//標記已被查找
ctemp=GiveScore(1,i,j);
for(m=0;m<572;m++)
{//暫時更改玩家信息
if(ptable[i][j][m])
{
temp1[n]=m;
ptable[i][j][m]=false;
temp2[n]=win[0][m];
win[0][m]=7;
n++;
}
}
ptempboard[i][j]=1;
/***********************************************/
/* 是否勝利,是返回1 */
int victory(int player)
{
int column=E_column,row=E_row;
int left_sum=0,right_sum=0;
int top_sum=0,down_sum=0;
int left_top_sum=0,right_down_sum=0;
int right_top_sum=0,left_down_sum=0;
int i,j;

for (i=column;i>=column-4 && i>0;i--)
if (chess_map[i-1][row-1]==player)
left_sum++;
else
break;
for (i=column;i<=column+4 && i<17;i++)
if (chess_map[i-1][row-1]==player)
right_sum++;
else
break;
if ((left_sum+right_sum)==6)
return 1;

for (j=row;j>=row-4 && j>0;j--)
if (chess_map[column-1][j-1]==player)
top_sum++;
else
break;
for (j=row;j<=row+4 && j<17;j++)
if (chess_map[column-1][j-1]==player)
down_sum++;
else
break;
if ((top_sum+down_sum)==6)
return 1;

for (i=column,j=row;i>=column-4 && j>=row-4 && i>0 && j>0;i--,j--)
if (chess_map[i-1][j-1]==player)
left_top_sum++;
else
break;
for (i=column,j=row;i<=column+4 && j<=row+4 && i<17 && j<17;i++,j++)
if (chess_map[i-1][j-1]==player)
right_down_sum++;
else
break;
if ((left_top_sum+right_down_sum)==6)
return 1;

for (i=column,j=row;i<=column+4 && j>=row-4 && i<17 && j>0;i++,j--)
if (chess_map[i-1][j-1]==player)
right_top_sum++;
else
break;
for (i=column,j=row;i>=column-4 && j<=row+4 && j<17 && i>0;j++,i--)
if (chess_map[i-1][j-1]==player)
left_down_sum++;
else
break;
if (right_top_sum+left_down_sum==6)
return 1;
return 0;
}

/**********************AI**********************/
int AI(int *m,int *n)
{
int bestx=8,besty=8;
long f_score=0,score=0;
int i,j;
int a[8],b[8];
for (i=0;i<8;i++)
{a[i]=0,b[i]=0;}

for (i=0;i<16;i++)
for (j=0;j<16;j++)
{
score=0;
if (chess_map[i][j]==0)
{
/********************進攻分數**********************/
chess_map[i][j]=computer;
for (i=column;i>=column-4 && i>0;i--)
if (chess_map[i-1][row-1]==computer)
a[0]++;
else
break;
for (i=column;i<=column+4 && i<17;i++)
if (chess_map[i-1][row-1]==computer)
a[1]++;
else
break;
for (j=row;j>=row-4 && j>0;j--)
if (chess_map[column-1][j-1]==computer)
a[2]++;
else
break;
for (j=row;j<=row+4 && j<17;j++)
if (chess_map[column-1][j-1]==computer)
a[3]++;
else
break;
for (i=column,j=row;i>=column-4 && j>=row-4 && i>0 && j>0;i--,j--)
if (chess_map[i-1][j-1]==computer)
a[4]++;
else
break;
for (i=column,j=row;i<=column+4 && j<=row+4 && i<17 && j<17;i++,j++)
if (chess_map[i-1][j-1]==computer)
a[5]++;
else
break;
for (i=column,j=row;i<=column+4 && j>=row-4 && i<17 && j>0;i++,j--)
if (chess_map[i-1][j-1]==computer)
a[6]++;
else
break;
for (i=column,j=row;i>=column-4 && j<=row+4 && j<17 && i>0;j++,i--)
if (chess_map[i-1][j-1]==computer)
a[7]++;
else
break;
for (i=0;i<8;i++)
{
switch(a[i]-1)
{
case 1:
score+=5;
break;
case 2:
score+=50;
break;
case 3:
score+=100;
break;
case 4:
score+=10000;
break;
default:
break;
}
}
/********************進攻分數**********************/
/********************防守分數**********************/
chess_map[i][j]=elva6401;
for (i=column;i>=column-4 && i>0;i--)
if (chess_map[i-1][row-1]==elva6401)
b[0]++;
else
break;
for (i=column;i<=column+4 && i<17;i++)
if (chess_map[i-1][row-1]==elva6401)
b[1]++;
else
break;
for (j=row;j>=row-4 && j>0;j--)
if (chess_map[column-1][j-1]==elva6401)
b[2]++;
else
break;
for (j=row;j<=row+4 && j<17;j++)
if (chess_map[column-1][j-1]==elva6401)
b[3]++;
else
break;
for (i=column,j=row;i>=column-4 && j>=row-4 && i>0 && j>0;i--,j--)
if (chess_map[i-1][j-1]==elva6401)
b[4]++;
else
break;
for (i=column,j=row;i<=column+4 && j<=row+4 && i<17 && j<17;i++,j++)
if (chess_map[i-1][j-1]==elva6401)
b[5]++;
else
break;
for (i=column,j=row;i<=column+4 && j>=row-4 && i<17 && j>0;i++,j--)
if (chess_map[i-1][j-1]==elva6401)
b[6]++;
else
break;
for (i=column,j=row;i>=column-4 && j<=row+4 && j<17 && i>0;j++,i--)
if (chess_map[i-1][j-1]==elva6401)
b[7]++;
else
break;
for(i=0;i<8;i++)
{
switch(b[i]-1)
{
case 1:
score-=5;
break;
case 2:
score-=50;
break;
case 3:
score-=500;
break;
case 4:
score-=5000;
break;
default:
break;
}
}
/********************防守分數**********************/
chess_map[i][j]=0;
if (score>f_score)
{
bestx=i+1;
besty=j+1;
f_score=score;
}
}
}
E_column=bestx,* m=bestx;
E_row=besty,* n=besty;
chess_map[bestx-1][besty-1]=computer;
}
/**********************AI***********************/

/***********************************************/
/* main */
int main()
{
int gdriver=DETECT,gmode;
initgraph(&gdriver,&gmode," ");
setbkcolor(BLUE);
init();
column=9,row=9;
chess_map[column-1][row-1]=computer;
draw_chess(computer);
while(1)
{
init_col_row();
flag=1;
while(flag==1)
{
move_chess(elva6401);
getch();
}
draw_chess(elva6401);
if(victory(elva6401))
{
outtextxy(550,400,"You win the game\n");
getch();
return 0;
}
AI(&column,&row);
draw_chess(computer);
if(victory(computer))
{
outtextxy(550,400,"You lost the game\n");
getch();
return 0;
}
}
}
// ChangeStatus(ptempboard);
pi=i;
pj=j;
while(SearchBlank(i,j,ptempboard))
{//進行第二不查找
ptempboard[i][j]=3;//標記已被查找
ptemp=GiveScore(0,i,j);
if(pscore>ptemp)//此時為玩家下子，運用極小極大法時應選取最小值
pscore=ptemp;
}
for(m=0;m<n;m++)
{//恢復玩家信息
ptable[pi][pj][temp1[m]]=true;
win[0][temp1[m]]=temp2[m];
}
//ctemp+pscore為實際估價函數
if(ctemp+pscore>cscore) //此時為計算機下子，運用極小極大法時應選取最最大值
{
cscore=ctemp+pscore;
bestx=pi;
besty=pj;
}
}
}
board[bestx][besty]=1;
if(m_pclastpos.x!=-1&&m_pclastpos.y!=-1)
{//畫前一棋子
if(!m_bwfirst)
DrawBlackChess(m_pclastpos.x,m_pclastpos.y);
else
DrawWhiteChess(m_pclastpos.x,m_pclastpos.y);
}

if(!m_bwfirst)
DrawNowBlack(bestx,besty);
else
DrawNowWhite(bestx,besty);
m_pclastpos.x=bestx;
m_pclastpos.y=besty;
for(i=0;i<572;i++)
{//修改計算機下子後，棋盤的變化狀況
if(ctable[bestx][besty][i]&&win[1][i]!=7)
win[1][i]++;
if(ptable[bestx][besty][i])
{
ptable[bestx][besty][i]=false;
win[0][i]=7;
}
}
computer=false;
player=true;
}
（8） 玩家下棋：（其中坐標（9） 加減目的是調整圖片的位置和棋盤對齊）
void CMyChessDlg::OnLButtonDown(UINT nFlags, CPoint point)
{
int x,y,tx,ty;
// TODO: Add your message handler code here and/or call default
if(player&&point.x<=535&&point.y<=535)//判斷是否在有效區域
{
tx=x=point.x-24;
ty=y=point.y-25;
while(tx>=36)
tx-=36;
while(ty>=36)
ty-=36;
tx+=x/36;
ty+=y/36;
if(tx>18)
x=x/36+1;
else
x=x/36;
if(ty>18)
y=y/36+1;
else
y=y/36;
if(board[x][y]==2)
{
board[x][y]=0;//設為玩家的棋子
if(m_pplastpos.x!=-1&&m_pplastpos.y!=-1)
{
if(!m_bwfirst)
DrawWhiteChess(m_pplastpos.x,m_pplastpos.y);
else
DrawBlackChess(m_pplastpos.x,m_pplastpos.y);
}

if(!m_bwfirst)
DrawNowWhite(x,y);
else
DrawNowBlack(x,y);
m_pplastpos.x=x;
m_pplastpos.y=y;
for(int i=0;i<572;i++)//修改玩家下子後棋盤狀態的變化
{
if(i==80)
i=80;
if(ptable[x][y][i]&&win[0][i]!=7)
win[0][i]++;
if(ctable[x][y][i])
{
ctable[x][y][i]=false;
win[1][i]=7;
}
}
player=false;
computer=true;
}
}
CDialog::OnLButtonDown(nFlags, point);
}

5. 張澤民教授課題組再發Cell，用到的單細胞技術了解一下

4月16日，百奧智匯創始人、科學顧問張澤民教授在北京大學的課題組與合作者在國際頂級學術期刊《Cell》上發表了題為＂Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer＂的文章，利用單細胞轉錄組測序技術對結直腸癌患者的腫瘤微環境，特別是浸潤髓系細胞類群首次進行了系統性的刻畫，同時利用小鼠模型，對anti-CSF1R抑制劑和anti-CD40激動劑兩種靶向髓系細胞的免疫治療策略潛在的作用機理給出了解釋。

該研究建立了結合腫瘤患者及小鼠模型的單細胞轉錄組來研究腫瘤免疫治療的範例，為人們研究其他疾病中免疫細胞以及開發新的治療方案提供了思路。

對於該研究，上海市免疫學研究所蘇冰所長、上海交通大學醫學院葉幼瓊研究員等相關專家點評道：該工作全面地解析結腸癌的腫瘤微環境細胞圖譜，闡明細胞與細胞之間的相互作用，對靶向腫瘤免疫微環境為基礎的腫瘤治療提供了更詳細的理論基礎，為今後靶定髓系細胞的精準治療提供助力，具有重要的臨床轉化意義。

在該研究中，研究者共使用了Smart-seq2、10x 3′ Gene Expression、10x V(D)J + 5′ Gene Expression等 3種單細胞測序技術 ，以及tSNE、UMAP、PCA、RNA velocity、URD、PAGA、STARTRAC、GSVA富集分析、熱圖、小提琴圖、氣泡熱圖、軌跡圖、Circos圖、火山圖、生存分析等 十多種生物信息學分析及展示方法。

值得一提的是， 這些實驗技術和分析方法，百奧智匯皆可實現。

下面，百奧我就為大家詳細解析，帶大家了解這些技術和方法是如何在研究中應用的。

1.10x Genomics和Smart-seq2技術比較

該研究首先評估並比較了10x Genomics和Smart-seq2兩種單細胞測序技術的結果。結果顯示，與10x scRNA-seq平台相比，Smart-seq2捕獲了更多的基因，包括細胞因子，CD分子，配體/受體和轉錄因子，並且顯示出較弱的批次效應 ，從而可以對調節途徑進行更深入的分析（圖2），而10x平台則獲得了更多的分群。因此，作者將兩種技術同時使用，從而能夠最大化地確定細胞類型或稀有種群的數量，改善細胞聚類的結果，得出更准確的結論。

2.tSNE降維——展示分群、基因表達模式、組織分布等多層信息

該研究使用無監督聚類、PCA、CGA等方法對來自18位CRC患者腫瘤、鄰近組織和血液樣本的10× 3′ Gene Expression （ 43,817個細胞）和Smart-seq2（10,468）單細胞測序結果進行整合聚類分析，然後用tSNE進行降維展示，分別得到38個和36個群，包括6個內皮細胞群，2個成纖維細胞群，13個髓系細胞群，4個ILC群，18個T細胞群和5個B細胞群等（圖3）。

對於淋巴細胞，研究通過特異的的免疫球蛋白重鏈特徵基因加以區分（圖4）。

同時，研究還在tSNE結果中展示了各細胞的組織分布情況（圖5）。

對於Smart-seq2非免疫細胞的測序結果，該研究利用tSNE將其分為12個亞群，包括4個惡性細胞亞群和8個非惡性細胞亞群（圖6）。其中，由推斷的拷貝數變異（CNV）定義的惡性細胞表現出高度的基因表達異質性，形成了患者特異性的分群（圖7）。

3.熱圖、氣泡熱圖——展示細胞間基因表達模式的差異

單個tSNE圖或組圖可展示單個或數個基因在不同細胞群中的表達情況，但在展示多個樣本大量基因的全局表達情況時就相對吃力。此時就需要用到熱圖。而氣泡熱圖則是在熱圖的基礎上加入了基因表達細胞在特定細胞群中佔比的信息，從而對細胞群進行更詳細的表徵。

在該研究中，作者用熱圖展示了10x測序分析得到的38個白細胞群中的差異基因表達差異（圖8），以及整合10x和Smart-seq2測序分析得到的48個群的基因表達差異（圖9）。

研究還分析了13個髓系細胞中的特徵基因表達情況（其中9個以氣泡熱圖展示），並據此將它們區分為肥大細胞（hM01），樹突狀細胞(hM02-04)，單核細胞(hM05-07，hM11)，組織駐留性巨噬細胞（hM08-10）和 腫瘤相關巨噬細胞（TAMs，hM12-13）。

4.擴散圖、RNA velocity、URD、PAGA——推斷和展示細胞發育軌跡

為進一步探索 TAMs 的來源，該研究利用擴散圖中嵌入的RNA velocity對隨機選擇的單核細胞和巨噬細胞的發育軌跡進行推斷和展示，鑒定出了從表達CD14的單核細胞向FCN1+ 單核樣細胞和不同的巨噬細胞群體的強烈定向流動（圖11），體現二者在發育上的前後順序。

進一步利用URD和PAGA這兩種正交演算法分析巨噬細胞的轉錄軌跡，發現巨噬細胞發育成TAMs（圖12，圖13）。

綜合上述軌跡推斷結果，研究發現TAM主要來自獨特的腫瘤浸潤性單核樣細胞前體。其中， C1QC+ 和 SPP1+ TAM都從浸潤腫瘤的單核細胞樣前體形成，而 SPP1+ TAM也可能源自 NLRP3+ RTM（圖14）。

5.火山圖、富集分析、熱圖、URD圖——多角度展示兩類TAMs差異

對於兩類在發育軌跡上存在顯著差異的TAMs，作者進一步使用熱圖、火山圖、富集分析等方法進行了分析，從多角度揭示了它們的差異。

由於細胞的發育軌跡受轉錄調控網路控制，作者先分析了兩類TAMs轉錄因子的表達差異，結果以熱圖顯示。其中，C1QC+ TAMs顯示出MAF/MAFB和FOS/JUNS高表達，而PP1+ TAMs則高表達CEBPB和ZEB2。

然後，作者又用火山圖展示兩類TAMs中顯著差異表達的基因。該圖包含兩個維度，其中縱軸P value體現顯著性，橫軸fold change展示差異性。

結果顯示，C1QC+ TAMs顯示出補體C1Q、TREM2、MERTK和CD80等基因的高表達， 而SPP1+ TAMs顯示出SPP1、MARCO和VEGFA的特異性表達。

隨後，作者又使用基因集合變異分析（GSVE）分析了兩類TAMs在通路上的差異。GSVE是一種以非監督方式對一個群體評估通路活性差異的基因集富集（GSE）分析方法。

該研究的結果顯示，SPP1+ TAMs顯示出腫瘤血管生成，ECM受體相互作用、腫瘤脈管系統、大腸腺瘤和轉移性肝癌等通路的特異性富集，而C1QC+ TAMs則顯示出補體激活以及抗原加工和呈遞途徑的富集（圖17）。

此外，SPP1 + TAMs還顯示了大腸腺瘤和轉移性肝癌通路的特異性富集和相關基因的表達（圖17和圖18），表明在它們CRC中有促癌/促轉移作用。

6.相互作用網路圖和Circos圖——展示細胞間相互作用、受體配體相互作用

更進一步地，作者在CRC中建立了細胞間相互作用網路圖。將該研究中的單細胞測序數據集與GTEx、TCGA的組織整體RNA測序數據集進行整合分析，發現TAM和cDC作為預測網路的核心，與其他細胞類型的聯系最多（圖19，圖20）。其中，C1QC+ TAM和兩組cDC主要與其他免疫細胞（尤其是T細胞亞群）相互作用（圖20），提示其在抗腫瘤T細胞應答中起調節作用。

再進一步的細胞群間的受體——配體相互作用分析顯示，在與髓系細胞和T細胞有關的配體-受體對中，CXCL10-CXCR3在C1QC + TAM中富集，暗示了C1QC+ TAM具有募集或激活T細胞的潛在作用。而SPP1+ TAMs中富集的SPP1-ITGAV、SDC2-MMP2、FN1-ITGA5等配體-受體對，則暗示了其在可能與某些整合蛋白相互作用，以促進CRC的腫瘤發生。

7.相似性分析、熱圖——展示跨物種的細胞群相似性

為了將上述對人髓系細胞異質性的研究發現與臨床應用相結合，作者接下來將相同的實驗和分析方法用於兩種對腫瘤免疫治療的小鼠模型中，其中Renca對CSF1R阻斷抗體敏感，而MC38對CD40激動劑抗體敏感（圖21）。

作者使用10x Genomics平台對免疫治療後小鼠的腫瘤中分離出的免疫細胞進行了單細胞轉錄組測序，並與人髓系細胞群進行了相似性分析，確定了多個跨物種相對應的髓系種群，包括兩種cDC群和兩種TAM群（圖22）。

此外，對小鼠TAM群進行與人類TAM群相同的途徑分析發現，小鼠TAM群體也基於它們的血管生成，低氧和T細胞相互作用基因特徵而分離（圖23）。這些數據表明人類CRC患者和小鼠腫瘤模型之間存在功能相似的TAM群體。

8.細胞豐度、生存分析——體現不同細胞類群的抗葯性

進一步的耐葯性研究顯示，抗CSF1R治療後F4/80高表達的巨噬細胞優先減少（圖24），說明不同的巨噬細胞群體對抗CSF1R治療的敏感性不同。同時，治療後小鼠細胞群中mC12和mM14簇幾乎完全丟失，TAM簇mM11，mM13和mM15的減少最小（圖24），說明TAMs對CSF1R阻斷的治療具有抗性。同時，抗CSF1R的TAM亞群優先表達參與血管生成和免疫抑制的基因，如Vegfa，Cd274和Arg1。

為了將在小鼠中的發現與人類CRC相關聯，作者使用生存分析比較了具有不同水平C1QC+ TAM和SPP1+ TAM基因特徵患者的生存率，發現低C1QC+ TAM、高SPP1 +TAM組合與CRC患者的預後更差相關（圖26）。

這些發現表明，抗CSF1R治療可能不足以耗盡所有具有促進腫瘤生長潛力的巨噬細胞，這一特性可能是其單葯療效差的原因。

類似的分析應用於抗CD40治療則發現，抗CD40治療能夠激活cDC1細胞，CCL22等激活的基因特徵與CRC患者的總體生存期呈正相關（圖27），這可能部分解釋了抗CD40激動劑治療CRC的機制。

9.STARTRAC——分析T細胞的遷移、克隆擴增和發育轉變

眾所周知，T細胞在腫瘤免疫治療中發揮重要作用。為進一步探索抗CD40激動劑治療CRC的機制，作者基於TCRα和β鏈序列，應用STARTRAC演算法分析抗CD40治療後T細胞的遷移、克隆擴增和發育轉變，以了解抗CD40激動劑治療對腫瘤浸潤性T細胞功能的影響。

結果顯示，抗CD40激動劑治療後，Ccl5+ Tem CD8+ T細胞具有比其他CD8+ T細胞亞群更高的遷移指數（圖28）。進一步剖析Ccl5+ Tem亞群內的TCR克隆型表明，克隆擴增較高的細胞在腫瘤和腫瘤引流淋巴結之間表現出更多的TCR共享，表明這些細胞在抗CD40處理後具有更強的遷移能力（圖29）。

同時，Ccl5+ Tem和Cxcr6+ Trm CD8+ T細胞之間的過渡指數在基線時顯著高於CD8+ T細胞亞群中其他對的過渡指數，表明某些Cxcr6+ 腫瘤中的Trm細胞可能在某一過程中由浸潤的Ccl5+ Tem細胞發育轉變而來，這一過程在抗CD40治療後得到了進一步增強（圖30）。

這些數據表明，抗CD40治療對腫瘤浸潤性T細胞的擴增，遷移和發育轉變具有獨特影響，能夠加強Ccl5+ Tem的遷移和克隆擴增能力，及其向Cxcr6+ Trm CD8+ T細胞轉變的能力。

10.相關性分析——揭示細胞間存在相互作用

在研究中作者發現，Th1類細胞與成熟和未成熟的cDC1細胞均顯示正相關，表明這兩類細胞之間存在相互作用。在此基礎上，結合該研究中發現的Bhlhe40 + Th1類細胞在抗CD40治療後的佔比以及特異性擴增、Bhlhe40 + CD4 + T細胞能夠產生更多IFNγ，以及之前研究發現的表達IFNG的BHLHE40 + Th1樣細胞富集於MSI CRC患者人群中且顯示出對ICB治療的良好反應等結果，作者認為抗CD40治療導致增加的Bhlhe40 + Th1樣細胞可能為在該模型中CD40激動劑治療能夠成功與抗PD1協同的機制提供了解釋。

總的來說，該研究的思路是：

先通過單細胞轉錄組、免疫組測序等技術，對CRC患者的腫瘤微環境進行細胞分群、基因差異表達、細胞發育軌跡、細胞間相互作用等多水平、多維度的詳細表徵，發現了腫瘤浸潤髓系細胞類群中兩類特殊的細胞類群TAM和cDC可能在CRC的抗腫瘤T細胞應答、腫瘤發生等過程起調節作用。然後利用小鼠模型，將上述發現應用於對anti-CSF1R抑制劑和anti-CD40激動劑兩種靶向髓系細胞的免疫治療策略的潛在作用機理的解釋中。

參考文獻：

Lei Z. et al., Single-Cell Analyses Inform Mechanisms of Myeloid-Targeted Therapies in Colon Cancer.  Cell . 2020.

6. 遊程編碼源代碼

這個...........樓上的諸位說的都是什麼啊。今天剛好看到這個問題，把你的E-mail給我把，我有純c的源碼（RLC）。

算了，直接貼關鍵部分吧。這個有一點C++成分，很容易改的。
bool CompressRunLength(BYTE *pSrc, int nSrcLen, BYTE *&pDes, int &nDesLen, int nBitsPerSample, void* nRuns, int nRunCount, int nRunSize);
bool DecompressRunLength(BYTE *pSrc, int nSrcLen, BYTE *&pDes, int &nDesLen);

#define GetDWORD(buf,bit,mask) ((*(DWORD*)(((BYTE*)buf)+((bit)>>3)))>>((bit)&7)&(mask))
#define GetWORD(buf,bit,mask) ((*(WORD*)(((BYTE*)buf)+((bit)>>3)))>>((bit)&7)&(mask))

int GetBitCount(int n)
{
int nBitCount = 0;
while(n)
n >>= 1, nBitCount++;
return nBitCount;
}

int BinarySearch(void* pValue, int nVlaueSize, void* pArray, int nCount)
{
int nIndex, nResult, nStart = 0, nEnd = nCount-1;
while(nStart <= nEnd)
{
nIndex = (nEnd+nStart)/2;
if((nResult = memcmp((BYTE*)pArray+nIndex*nVlaueSize, pValue, nVlaueSize)) == 0)
return nIndex;
if(nResult > 0)
nEnd = nIndex-1;
else
nStart = nIndex+1;
}
return -1;
}

bool CompressRunLength(BYTE *pSrc, int nSrcLen, BYTE *&pDes, int &nDesLen, int nBitsPerSample, void* nRuns, int nRunCount, int nRunSize)
{
pDes = (BYTE*)malloc(nSrcLen*2);
memset(pDes, 0, nSrcLen*2);

nDesLen = sizeof(DWORD);
*(DWORD*)pDes = nSrcLen; // save source length
*(pDes+nDesLen++) = nBitsPerSample; // save bits per sample
*(pDes+nDesLen++) = nRunCount; // save runs count
*(pDes+nDesLen++) = nRunSize; // save run bytes
memcpy(pDes+nDesLen, nRuns, nRunCount*nRunSize); // save runs
nDesLen += nRunCount*nRunSize;
nDesLen <<= 3; // bytes to bits
if(nRunCount == 0)
nRunCount = 256, nRunSize = 1, nRuns = NULL;

int nBitsPerTypeIndex = GetBitCount(nRunCount-1);
int nMaxRunLength = (1 << nBitsPerSample)-1, nRunLength, nRunIndex, nByte = 0;
// loop in the source buffer
while(nByte < nSrcLen)
if((nRuns && (nRunIndex = BinarySearch(pSrc+nByte, nRunSize, nRuns, nRunCount)) != -1 &&
memcmp(pSrc+nByte+nRunSize, (BYTE*)nRuns+nRunIndex*nRunSize, nRunSize) == 0) ||
(!nRuns && (nRunIndex = *(pSrc+nByte)) == *(pSrc+nByte+1)))
{ // set bit to 1 to indicate type found
*(pDes+(nDesLen>>3)) |= 1 << (nDesLen&7);
*(DWORD*)(pDes+(++nDesLen>>3)) |= nRunIndex << (nDesLen&7);
nDesLen += nBitsPerTypeIndex;
// skip the two repeated runs
nByte += nRunSize*2;
// get run length - 2 (without the two repeated runs)
nRunLength = 0;
while(nRunLength < nMaxRunLength && nByte < nSrcLen &&
((nRuns && memcmp(pSrc+nByte, (BYTE*)nRuns+nRunIndex*nRunSize, nRunSize) == 0) || (!nRuns && (BYTE)nRunIndex == *(pSrc+nByte))))
nRunLength++, nByte += nRunSize;
// save run length and increment destination length by bits per sample
*(DWORD*)(pDes+(nDesLen>>3)) |= nRunLength << (nDesLen&7);
nDesLen += nBitsPerSample;
}
else // one byte
*(WORD*)(pDes+(++nDesLen>>3)) |= *(pSrc+nByte++) << (nDesLen&7), nDesLen += 8;
nDesLen = (nDesLen+7)/8; // bits to bytes
pDes = (BYTE*)realloc(pDes, nDesLen);

return true;
}

bool DecompressRunLength(BYTE *pSrc, int nSrcLen, BYTE *&pDes, int &nDesLen)
{
if(nSrcLen == 0)
return true;

// allocate destination buffer
nDesLen = *(DWORD*)pSrc;
pDes = (BYTE*)malloc(nDesLen);
memset(pDes, 0, nDesLen);

// compression information
int nSrcIndex = sizeof(DWORD);
int nBitsPerSample = *(pSrc+nSrcIndex++);
int nRunCount = *(pSrc+nSrcIndex++);
int nRunSize = *(pSrc+nSrcIndex++);
void* nRuns = pSrc+nSrcIndex;
nSrcIndex += nRunSize*nRunCount;
nSrcIndex <<= 3; // bytes to bits
if(nRunCount == 0)
nRunCount = 256, nRunSize = 1, nRuns = NULL;

int nBitsPerTypeIndex = GetBitCount(nRunCount-1);
int nMaxTypeIndex = (1 << nBitsPerTypeIndex)-1;
int nMaxRunLength = (1 << nBitsPerSample)-1;
int nDesIndex = 0, nRunLength, nRunIndex, nRun, nByte;

nSrcLen <<= 3; // bytes to bits
while(nSrcIndex < nSrcLen-8)
if((*(pSrc+(nSrcIndex>>3)) >> (nSrcIndex++&7)) & 1)
{
nRunIndex = GetDWORD(pSrc, nSrcIndex, nMaxTypeIndex), nSrcIndex += nBitsPerTypeIndex;
nRunLength = GetDWORD(pSrc, nSrcIndex, nMaxRunLength)+2, nSrcIndex += nBitsPerSample;
for(nRun = 0; nRun < nRunLength; nRun++)
for(nByte = 0; nByte < nRunSize; nByte++, nDesIndex += 8)
*(WORD*)(pDes+(nDesIndex>>3)) |= nRuns ? GetWORD(nRuns+nRunSize*nRunIndex, nByte<<3, 0xff) : (BYTE)nRunIndex;
}
else // one byte
*(WORD*)(pDes+(nDesIndex>>3)) |= GetWORD(pSrc, nSrcIndex, 0xff), nDesIndex += 8, nSrcIndex += 8;

return true;
}

7. 求32位MD5加密c語言源碼

#include<stdio.h>
#include<string.h>
#define F(x, y, z) (((x) & (y)) | ((~x) & (z)))
#define G(x, y, z) (((x) & (z)) | ((y) & (~z)))
#define H(x, y, z) ((x) ^ (y) ^ (z))
#define I(x, y, z) ((y) ^ ((x) | (~z)))

#define RL(x, y) (((x) << (y)) | ((x) >> (32 - (y)))) //x向左循環移y位

#define PP(x) (x<<24)|((x<<8)&0xff0000)|((x>>8)&0xff00)|(x>>24) //將x高低位互換,例如PP(aabbccdd)=ddccbbaa

#define FF(a, b, c, d, x, s, ac) a = b + (RL((a + F(b,c,d) + x + ac),s))
#define GG(a, b, c, d, x, s, ac) a = b + (RL((a + G(b,c,d) + x + ac),s))
#define HH(a, b, c, d, x, s, ac) a = b + (RL((a + H(b,c,d) + x + ac),s))
#define II(a, b, c, d, x, s, ac) a = b + (RL((a + I(b,c,d) + x + ac),s))

unsigned A,B,C,D,a,b,c,d,i,len,flen[2],x[16]; //i臨時變數,len文件長,flen[2]為64位二進製表示的文件初始長度
char filename[200]; //文件名
FILE *fp;

void md5(){ //MD5核心演算法,供64輪

a=A,b=B,c=C,d=D;
/**//* Round 1 */
FF (a, b, c, d, x[ 0], 7, 0xd76aa478); /**//* 1 */
FF (d, a, b, c, x[ 1], 12, 0xe8c7b756); /**//* 2 */
FF (c, d, a, b, x[ 2], 17, 0x242070db); /**//* 3 */
FF (b, c, d, a, x[ 3], 22, 0xc1bdceee); /**//* 4 */
FF (a, b, c, d, x[ 4], 7, 0xf57c0faf); /**//* 5 */
FF (d, a, b, c, x[ 5], 12, 0x4787c62a); /**//* 6 */
FF (c, d, a, b, x[ 6], 17, 0xa8304613); /**//* 7 */
FF (b, c, d, a, x[ 7], 22, 0xfd469501); /**//* 8 */
FF (a, b, c, d, x[ 8], 7, 0x698098d8); /**//* 9 */
FF (d, a, b, c, x[ 9], 12, 0x8b44f7af); /**//* 10 */
FF (c, d, a, b, x[10], 17, 0xffff5bb1); /**//* 11 */
FF (b, c, d, a, x[11], 22, 0x895cd7be); /**//* 12 */
FF (a, b, c, d, x[12], 7, 0x6b901122); /**//* 13 */
FF (d, a, b, c, x[13], 12, 0xfd987193); /**//* 14 */
FF (c, d, a, b, x[14], 17, 0xa679438e); /**//* 15 */
FF (b, c, d, a, x[15], 22, 0x49b40821); /**//* 16 */

/**//* Round 2 */
GG (a, b, c, d, x[ 1], 5, 0xf61e2562); /**//* 17 */
GG (d, a, b, c, x[ 6], 9, 0xc040b340); /**//* 18 */
GG (c, d, a, b, x[11], 14, 0x265e5a51); /**//* 19 */
GG (b, c, d, a, x[ 0], 20, 0xe9b6c7aa); /**//* 20 */
GG (a, b, c, d, x[ 5], 5, 0xd62f105d); /**//* 21 */
GG (d, a, b, c, x[10], 9, 0x02441453); /**//* 22 */
GG (c, d, a, b, x[15], 14, 0xd8a1e681); /**//* 23 */
GG (b, c, d, a, x[ 4], 20, 0xe7d3fbc8); /**//* 24 */
GG (a, b, c, d, x[ 9], 5, 0x21e1cde6); /**//* 25 */
GG (d, a, b, c, x[14], 9, 0xc33707d6); /**//* 26 */
GG (c, d, a, b, x[ 3], 14, 0xf4d50d87); /**//* 27 */
GG (b, c, d, a, x[ 8], 20, 0x455a14ed); /**//* 28 */
GG (a, b, c, d, x[13], 5, 0xa9e3e905); /**//* 29 */
GG (d, a, b, c, x[ 2], 9, 0xfcefa3f8); /**//* 30 */
GG (c, d, a, b, x[ 7], 14, 0x676f02d9); /**//* 31 */
GG (b, c, d, a, x[12], 20, 0x8d2a4c8a); /**//* 32 */

/**//* Round 3 */
HH (a, b, c, d, x[ 5], 4, 0xfffa3942); /**//* 33 */
HH (d, a, b, c, x[ 8], 11, 0x8771f681); /**//* 34 */
HH (c, d, a, b, x[11], 16, 0x6d9d6122); /**//* 35 */
HH (b, c, d, a, x[14], 23, 0xfde5380c); /**//* 36 */
HH (a, b, c, d, x[ 1], 4, 0xa4beea44); /**//* 37 */
HH (d, a, b, c, x[ 4], 11, 0x4bdecfa9); /**//* 38 */
HH (c, d, a, b, x[ 7], 16, 0xf6bb4b60); /**//* 39 */
HH (b, c, d, a, x[10], 23, 0xbebfbc70); /**//* 40 */
HH (a, b, c, d, x[13], 4, 0x289b7ec6); /**//* 41 */
HH (d, a, b, c, x[ 0], 11, 0xeaa127fa); /**//* 42 */
HH (c, d, a, b, x[ 3], 16, 0xd4ef3085); /**//* 43 */
HH (b, c, d, a, x[ 6], 23, 0x04881d05); /**//* 44 */
HH (a, b, c, d, x[ 9], 4, 0xd9d4d039); /**//* 45 */
HH (d, a, b, c, x[12], 11, 0xe6db99e5); /**//* 46 */
HH (c, d, a, b, x[15], 16, 0x1fa27cf8); /**//* 47 */
HH (b, c, d, a, x[ 2], 23, 0xc4ac5665); /**//* 48 */

/**//* Round 4 */
II (a, b, c, d, x[ 0], 6, 0xf4292244); /**//* 49 */
II (d, a, b, c, x[ 7], 10, 0x432aff97); /**//* 50 */
II (c, d, a, b, x[14], 15, 0xab9423a7); /**//* 51 */
II (b, c, d, a, x[ 5], 21, 0xfc93a039); /**//* 52 */
II (a, b, c, d, x[12], 6, 0x655b59c3); /**//* 53 */
II (d, a, b, c, x[ 3], 10, 0x8f0ccc92); /**//* 54 */
II (c, d, a, b, x[10], 15, 0xffeff47d); /**//* 55 */
II (b, c, d, a, x[ 1], 21, 0x85845dd1); /**//* 56 */
II (a, b, c, d, x[ 8], 6, 0x6fa87e4f); /**//* 57 */
II (d, a, b, c, x[15], 10, 0xfe2ce6e0); /**//* 58 */
II (c, d, a, b, x[ 6], 15, 0xa3014314); /**//* 59 */
II (b, c, d, a, x[13], 21, 0x4e0811a1); /**//* 60 */
II (a, b, c, d, x[ 4], 6, 0xf7537e82); /**//* 61 */
II (d, a, b, c, x[11], 10, 0xbd3af235); /**//* 62 */
II (c, d, a, b, x[ 2], 15, 0x2ad7d2bb); /**//* 63 */
II (b, c, d, a, x[ 9], 21, 0xeb86d391); /**//* 64 */

A += a;
B += b;
C += c;
D += d;

}

main(){
while(1){
printf("Input file:");
gets(filename); //用get函數,避免scanf以空格分割數據,
if (filename[0]==34) filename[strlen(filename)-1]=0,strcpy(filename,filename+1); //支持文件拖曳,但會多出雙引號,這里是處理多餘的雙引號
if (!strcmp(filename,"exit")) exit(0); //輸入exit退出
if (!(fp=fopen(filename,"rb"))) {printf("Can not open this file!\n");continue;} //以二進制打開文件
fseek(fp, 0, SEEK_END); //文件指針轉到文件末尾
if((len=ftell(fp))==-1) {printf("Sorry! Can not calculate files which larger than 2 GB!\n");fclose(fp);continue;} //ftell函數返回long,最大為2GB,超出返回-1
rewind(fp); //文件指針復位到文件頭
A=0x67452301,B=0xefcdab89,C=0x98badcfe,D=0x10325476; //初始化鏈接變數
flen[1]=len/0x20000000; //flen單位是bit
flen[0]=(len%0x20000000)*8;
memset(x,0,64); //初始化x數組為0
fread(&x,4,16,fp); //以4位元組為一組,讀取16組數據
for(i=0;i<len/64;i++){ //循環運算直至文件結束
md5();
memset(x,0,64);
fread(&x,4,16,fp);
}
((char*)x)[len%64]=128; //文件結束補1,補0操作,128二進制即10000000
if(len%64>55) md5(),memset(x,0,64);
memcpy(x+14,flen,8); //文件末尾加入原文件的bit長度
md5();
fclose(fp);
printf("MD5 Code:%08x%08x%08x%08x\n",PP(A),PP(B),PP(C),PP(D)); //高低位逆反輸出
}
}