CodyNote006：惡意滿滿的什一抽殺率問題

09-13

來自專欄 Cody習題7 人贊了文章

馬良，祁彬彬

序言

MATLAB代碼寫成垃圾的理由變化萬千，但想寫成經典，大體滿足如下條件：

邏輯縝密演算法合適，充分挖掘與利用問題的潛在條件；
函數選擇合理恰當，規避昏昏欲睡的常規套路；
高度精鍊濃縮，銜接絲絲入扣。

滿足上述條件時，for循環出現的頻率通常較低，Cody的優秀代碼看多了，往往下意識地，就留下：「好代碼和for/while循環互斥」的印象。

其實大謬不然。代碼中流程和函數選擇，歸根結底為解決問題服務，什麼對主旨有利就用什麼。不是for/while流程簡單和枯燥，我看根本是對這流程沒理解透。前天知乎看到關於for loop是否一定比矢量化代碼速度慢的問題，一個回答深得我心，關鍵就兩點：「矢量化代碼內存讀寫時間的消耗」、「新版循環速度並不慢」。對了，比較諷刺的是：回答問題的這位，個人感覺是地球上MATLAB矢量化代碼寫得最好的幾個人之一，我可以隨手擺出100條他在不同問題中留下的代碼，證明他絕非對矢量化代碼惱羞成怒，因愛生恨...

言歸正傳，這次專門舉例，討論以for/while循環為求解代碼主體結構的經典問題。

Pro.1092.Decimation

鏈接：https://ww2.mathworks.cn/matlabcentral/cody/problems/1092

【原題】：When dealing to the Roman Army, the term decimate meant that the entire unit would be broken up into groups of ten soldiers, and lots would be drawn. The person who was unlucky enough to draw the short straw would be executed by the other nine members of his group.

The bloodthirsty Roman Centurion Carnage Maximus decided to apply this to his prisoners, with a few gruesome differences. Rather than kill every tenth prisoner and allow the rest to live, he is going to leave only one prisoner alive and kill all of the others. Instead of killing every tenth prisoner, he chooses a number (kill_every). If kill_every=3, he kills every third prisoner. If kill_every=5, he kills every fifth prisoner. He always chooses a number between 2 and the number of prisoners he has, and this process will be repeated until there is only one prisoner left. For example, if there are 10 prisoners, and kill_every=3

First iteration: 1 2 3 4 5 6 7 8 9 10

1-2-3 4-5-6 7-8-9 10

Prisoners 3, 6 and 9 will be killed.

Second iteration:

1 2 4 5 7 8 10

Because Prisoner 10 was counted during the first iteration, the executions will proceed as such: 10-1-2 4-5-7 8-10, so prisoners 2 and 7 will be killed

Third iteration: 1 4 5 8 10 8-10-1 4-5-8 10, so prisoners 1 and 8 executed.

Fourth Iteration: 10-4-5 10 Prisoner 5 is executed.

Fifth iteration: 10-4 10 Prisoner 10 is executed

Since the sole survivor is prisoner 4, he is released.

You are an unlucky prisoner caught by Carnage Maximum. Prior to lining up the prisoners, he reveals the number of prisoners he has and his value of kill_every for the day. Your job is to figure out which prisoner you need to be in order to survive. Write a MATLAB script that takes the values of num_prisoners and kill_every. The output will be survivor, which is the position of the person who survives. If you write your script quickly enough, that person will be you.

Good luck!

【解釋】：內容是基於古羅馬惡名昭著的「什一抽殺率」的數學問題，為防止士兵逃跑，把士兵分成10人一組，抽籤抽出一個，讓其他9人殺死此人，此為「什一抽殺」。似乎古代羅馬人在節約糧食、降低管理成本時，「開源」和「節流」都做得非常徹底...

本題中根據「嗜血的羅馬百夫長Carnage Maximus決定對戰俘實施抽殺」的傳說。搜索Wiki未發現這位羅馬百夫長的相關條目，猜測題目敘述系野史杜撰。另一個與戰俘抽殺有關的數學問題也與羅馬人有關，即著名的「Josephus環」問題，這在Cody上也有題目，指的是羅馬皇帝尼祿於公元67年，指派韋斯巴薌將軍進攻猶太人聚集的高地之城約塔帕塔，44日後攻陷該城，歷史學家Josephus作為加利利軍團的指揮官，在其他39名猶太士兵掩護下退入山洞，並被羅馬軍隊包圍，士兵決定自殺，但Josephus不想死，因此想出主意，稱自殺有悖教典，提出讓士兵圍成一圈，每隔3人抽取1個，由別人動手，得到大家的贊同，而他把自己安排在一個特定的數字位置，直至最後剩下他和另一名士兵，他突然偷襲殺死最後這名士兵，並向韋斯巴薌投降，之後寫出《猶太戰爭史》。這個數學問題與現代計算機數據結構中的隊列、循環鏈表類型的定義有一些關聯。本問題對Josephus環問題做了推廣：戰俘一字排開，每人一個獨立編號，指揮官先抽一個隨機數(Kill_every)，用於決定每隔多少處決一個戰俘，然後循環，直至剩下最後一名囚犯被釋放，以示仁慈。例如10個戰俘，指揮官抽取到的隨機數為3:

第1次迭代中，第3、6、9名囚犯被處決：

1 2 (3) 4 5 (6) 7 8 (9) 10

第2次迭代中，因第10名囚犯第1次迭代沒輪到，第2輪從這名囚犯開始，則第2和第7名囚犯退場：

10 1 (2) 4 5 (7) 8 10

第3次迭代，由於8號前次迭代沒輪到，因此下一輪迭代放首位，此輪1和8號退場。

8 10 (1) 4 5 (8) 10

第4次迭代，5號退場

10 4 (5) 10

第5次迭代，10號退場：

(10) 4 (10)

幸運的4號囚犯獲釋，迭代結束。

Test

Code Input

%%assert(isequal(decimate(10,3),4))

%%assert(isequal(decimate(1024,3),676))

%%assert(isequal(decimate(2012,50),543))

%%assert(isequal(decimate(30,5),3))

%%assert(isequal(decimate(10,10),8))

%%assert(isequal(decimate(2048,2),1))

【分析】：描述有點歧義，對每次抽取完畢，尾數上剩下的人如何處理，解釋不詳細，不過不影響問題解決，這個歧義可以用後面的代碼來解釋。根據給出的test suite，如果指揮官決定在10個人中每隔3個處死1人循環，最後的幸運兒是4號戰俘；如果1024人每隔3人，676號生還；2012人每隔50人，543號生還，30人每隔5人，3號生還；10人每隔10人，8號生還；2048人每隔2人，1號生還。

因此，完全沒有什麼顯著的數字變化特徵，就不要想著投機取巧找規律了。

求解代碼

不定次數的while循環

問題關鍵在於對前次迭代結果的處理。首先，從流程結構上講，似乎一眼過去並不清楚具體的循環次數，因此不定次數while循環是第一選擇；其次，每次移除被整除數據，剩下元素的排列重組是關鍵。為理解這兩個問題的解決辦法，不妨對算例逐次迭代，觀察運算中間結果。

第1種while循環：逐步移位的索引變換

方便說明問題，運行過程用sprintf增加循環迭代次數的提示行。

function survivor=decimate(num_prisoners,kill_every) survivor=1:num_prisoners;x=1; while numel(survivor)>1 sprintf(第 %d 次迭代,x) mod(kill_every-1,numel(survivor))+1 survivor=survivor([ans+1:end 1:ans-1]) x=x+1; endend

Command Windows中運行結果如下：

>> num_prisoners=10,kill_every=3;survivor=decimate(num_prisoners,kill_every)num_prisoners = 10ans = 第 1 次迭代ans = 3survivor = 4 5 6 7 8 9 10 1 2ans = 第 2 次迭代ans = 3survivor = 7 8 9 10 1 2 4 5ans = 第 3 次迭代ans = 3survivor = 10 1 2 4 5 7 8ans = 第 4 次迭代ans = 3survivor = 4 5 7 8 10 1ans = 第 5 次迭代ans = 3survivor = 8 10 1 4 5ans = 第 6 次迭代ans = 3survivor = 4 5 8 10ans = 第 7 次迭代ans = 3survivor = 10 4 5ans = 第 8 次迭代ans = 3survivor = 10 4ans = 第 9 次迭代ans = 1survivor = 4survivor = 4

觀察逐級迭代的運行結果發現：

(1).每次去掉的是前次剩餘戰俘中的第3個，並且以第3個為軸點，將第1、2這兩個起始編號移至末尾。移至末尾的原因是代碼中最關鍵是計算倖存戰俘，使用mod命令，其參數為：

mod(kill_every-1,numel(survivor))+1

第1參數用間隔數減1，相當於在第kill_every位被抽中，則kill_every -1確保前面3-1=2位戰俘此次迭代安全，為不干擾下一波迭代，挪到序列最後。作為後續迭代的倖存戰俘編隊序列。換句話說移位方式「碎片化」了原來整體抽取間隔數過程，即每次去掉1個數，而不像題意中那樣，一次迭代把本次編隊中符合間隔數量的編號都去掉，所以稱為「逐步移位」；

(2). 當前迭代剩餘戰俘數量不足指定的間隔數量，本例在第9次(倒數第2次)迭代時出現這種情況，mod命令剛好把不夠抽取的尾數重新算到前面。好處在於無需每次迭代把尾數從前往後、從後往前反覆挪。

為加深印象，不妨用n=k=10的例子說明問題，即：編隊初始總人數為10，間隔數也設為10，如下為程序迭代結果：

>> num_prisoners=10,kill_every=10;survivor=decimate(num_prisoners,kill_every)num_prisoners = 10ans = 第 1 次迭代ans = 10survivor = 1 2 3 4 5 6 7 8 9ans = 第 2 次迭代ans = 1survivor = 2 3 4 5 6 7 8 9ans = 第 3 次迭代ans = 2survivor = 4 5 6 7 8 9 2ans = 第 4 次迭代ans = 3survivor = 7 8 9 2 4 5ans = 第 5 次迭代ans = 4survivor = 4 5 7 8 9ans = 第 6 次迭代ans = 5survivor = 4 5 7 8ans = 第 7 次迭代ans = 2survivor = 7 8 4ans = 第 8 次迭代ans = 1survivor = 8 4ans = 第 9 次迭代ans = 2survivor = 8

同理，第1次迭代，命令mod(9,10)+1得到結果為10，代表第一次被拿掉的編號是10，剩餘戰俘數量變為9，而第2次迭代時，mod(9,9)+1=1，代表拿掉編號為1，變成1的原因是總體為9，一次編隊不夠排序，重新挪回隊伍之首，所以1被抽中，余類推。

上述迭代過程關鍵在於Jan Orwat的逆向思維，用mod命令求余時，一般本能選擇讓總編隊放在mod第1參數位，間隔數kill_every放第2參數位，這樣第1次迭代可一次去掉3個數，但下次迭代構造新序列時比較麻煩。此外，把迭代編隊放在前面對間隔數求余，往往不用mod命令，改為circshift函數做整體輪轉移位。

第2種while循環：while+circshift整體移位

輪轉移位的circshift命令前一篇筆記已經介紹，直接看實現代碼。

function ans =decimate(n,k) 1:n; while numel(ans)>1 circshift(ans,-k+1); ans(2:end); endend

此while循環解法已相當簡練，結束while循環判斷條件仍是最終剩1人。以n=10,k=3為例，運行查看中間迭代運算結果：

Loop = 第 1 次迭代ans = 3 4 5 6 7 8 9 10 1 2ans = 4 5 6 7 8 9 10 1 2Loop = 第 2 次迭代ans = 6 7 8 9 10 1 2 4 5ans = 7 8 9 10 1 2 4 5Loop = 第 3 次迭代ans = 9 10 1 2 4 5 7 8ans = 10 1 2 4 5 7 8Loop = 第 4 次迭代ans = 2 4 5 7 8 10 1ans = 4 5 7 8 10 1Loop = 第 5 次迭代ans = 7 8 10 1 4 5ans = 8 10 1 4 5Loop = 第 6 次迭代ans = 1 4 5 8 10ans = 4 5 8 10Loop = 第 7 次迭代ans = 8 10 4 5ans = 10 4 5Loop = 第 8 次迭代ans = 5 10 4ans = 10 4Loop = 第 9 次迭代ans = 10 4ans = 4ans = 4

類似工廠流水線上的機械手，生產線由帶傳動或者鏈傳動單位間隔時間移動確定的步距，機械手拿掉確定位置上的某個物件。例如第1次迭代，把首位的兩個在第1次迭代肯定安全的數字挪到尾部，3號暴露在第1個位置，用提取ans(2:end)去掉首位3，第1次迭代結束，第2次迭代同樣挪走本次迭代安全的4和5，6至首位，以此類推，到了最後第8次迭代時，只剩元素10和4，circshift運算結果如下：

>> circshift([10 4],-1)ans = 4 10>> circshift([10 4],-2)ans = 10 4

迭代中固定的移動步距-2相當於原序列原封不動進入最後一次迭代，再抽掉第1位元素，最後剩下的還是4。

還有主體思路類似，細節安排則異常精緻的優化代碼，比如：

function ans = decimate(s,ke) 1:s; while diff(ans) circshift(ans, -ke); ans(end) = []; endend

個人認為，這段代碼大幅while循環，屬於不定數循環最優解法之一，注意circshift第2參數和之前有所不同，之前是1-k，上述代碼去掉前面的1，運行查看中間結果：

>> decimate(10,3)Loop = 第 1 次迭代ans = 4 5 6 7 8 9 10 1 2 3ans = 4 5 6 7 8 9 10 1 2Loop = 第 2 次迭代ans = 7 8 9 10 1 2 4 5 6ans = 7 8 9 10 1 2 4 5Loop = 第 3 次迭代ans = 10 1 2 4 5 7 8 9ans = 10 1 2 4 5 7 8Loop = 第 4 次迭代ans = 4 5 7 8 10 1 2ans = 4 5 7 8 10 1Loop = 第 5 次迭代ans = 8 10 1 4 5 7ans = 8 10 1 4 5Loop = 第 6 次迭代ans = 4 5 8 10 1ans = 4 5 8 10Loop = 第 7 次迭代ans = 10 4 5 8ans = 10 4 5Loop = 第 8 次迭代ans = 10 4 5ans = 10 4Loop = 第 9 次迭代ans = 4 10ans = 4

該解法和本節前1種先清除元素再輪轉移位不同，是反過來先把整個1-3環節移位到尾部，再通過去掉最後一個數來達到移除的效果；此外while循環結束判定語句不同，採用「diff(survivor)」，最後只剩1個數時，diff執行結果是空矩陣[]，可能有人奇怪空矩陣怎麼能作為循環結束的判定條件呢？

在MATLAB中，當非0(~=0)時邏輯判定值即為真值1(TRUE)，假值為0(FALSE)，空矩陣為何能當FALSE判定？見如下代碼執行結果：

>> numel([])ans = 0

這就很清楚了：空矩陣的元素數量為0，這也在MATLAB中視作結束循環的FALSE條件。從這裡能夠窺得MATLAB語言本身的靈活性。

第3種while循環：反推型的索引增長序列

直接都是逐步縮減倖存者人數，下面介紹1種非常新穎的逆推代碼。

function ans = decimate(s,ke) 1;x=1; while numel(ans)<s Loop=sprintf(第 %d ′次迭代,x) circshift([0 ans], ke-1) x=x+1; end find(ans)end

這段代碼初看有點兒費解，通過迭代的中間過程，看看變化趨勢：

>> decimate(10,3)Loop = 第 1 次迭代ans = 0 1Loop = 第 2 次迭代ans = 0 1 0Loop = 第 3 次迭代ans = 1 0 0 0Loop = 第 4 次迭代ans = 0 0 0 1 0Loop = 第 5 次迭代ans = 1 0 0 0 0 0Loop = 第 6 次迭代ans = 0 0 0 1 0 0 0Loop = 第 7 次迭代ans = 0 0 0 0 0 0 1 0Loop = 第 8 次迭代ans = 1 0 0 0 0 0 0 0 0Loop = 第 9 次迭代ans = 0 0 0 1 0 0 0 0 0 0ans = 4

迭代整個過程都是處理0和1，不停增加0並改變1的索引位置，乍一看不大容易理解，不過比較前面解法對同一參數的執行結果，發現歷次迭代1的位置就是編號4倖存者在歷次迭代中所處的位置！前述解法while循環最後1次迭代從[10,4]中剔除10，由於本解法是逆序倒推，所以變成第1次迭代，1的位置正好在索引位2，再看最後的移位距離，是3-1=2，這正好是之前正向縮減序列移位的相反數！於是倒推就通過初值設定而成立。這個初始條件看似輕鬆，實則需要對移位循環問題的數字變化規律有極其深刻的理解才可以。這個反推代碼及之前的「先移位、後移除」做法，還有下面馬上要介紹的遞歸方法，都由Cody上ID為yurenchu的兄弟提交，猜測是中國人，之前名不見經傳，今天看到他的代碼，真是驚到了！

遞歸

看到上述流程，也可能會想到遞歸，如下是yurenchu提交的遞歸代碼：

function s = decimate(s,ke) if s>1, s = mod(decimate(s-1,ke)-1+ke,s)+1 endend

一樣極致簡潔。遞歸流程實際是調用自身，對for/while循環逆序回溯，n=10,k=3執行結果：

>> s = decimate(10,3)s = 2s = 2s = 1s = 4s = 1s = 4s = 7s = 1s = 4

中間結果這段數字貌似讓人一頭霧水，其實就是前面反推型索引增長序列歷次迭代的執行結果。有這個結果反推就比較容易理解：遞歸流程使用編隊遞減的mod代碼，相當於反推了遞增型代碼的結果。再次對這位yurenchu兄弟的精彩構思擊節讚歎！

給定次數的for循環

下面主角定數循環for流程登場！

粗看貌似for循環不適合出現在循環節不明顯或者不易控制的問題，而這恰恰是精彩的部分：

從答案的提交時間上看，彬彬是第1個使用for循環完成這道題目的，大神Alfonso還「like」了這個解法，由於彬彬解法可能令人感到費解，所以我先列出彬彬的原始做法，最後用Alfonso的優化代碼解釋原理。

彬彬的解法：

function ans = decimate(num_prisoners,kill_every)0;for i = 2 : num_prisoners ans = mod(ans + kill_every,i);endans + 1;end

Alfonso的優化代碼

如下是Alfonso在彬彬思路基礎上優化過的代碼：

function ans = decimate(ans,k)for i=1:ans 1+mod(ans+k-1,i)endend

n=10,k=3的執行結果：

>> decimate(10,3)ans = 1ans = 2ans = 2ans = 1ans = 4ans = 1ans = 4ans = 7ans = 1ans = 4

從Alfonso的代碼發現他被Cody所有人敬服，不是沒有原因的，許多代碼堪稱藝術作品，經常提出一些經典解法的標準套路。從這寥寥3行的代碼中得到充分體現。首先，ans既是輸入又是函數的輸出，在輸入時，它作為編隊的總數，順次進入循環節和冒號操作符無縫對接，變成定數循環的循環總次數；其次，彬彬代碼中需要在結尾做「+1」處理，在Alfonso的代碼中被濃縮進了每個循環節內部，當然，從執行效率上講，彬彬的更優，因為「+1」運算在彬彬代碼中只需1次，且循環數量小1，但Alfonso代碼本身更簡潔，size更小，畢竟算例給出的數字都不大，效率並非主要矛盾。

最後，也是最重要的：循環節內的迭代。

(1). 第1次迭代，ans=10，故：ans=1+mod(ans+k-1,i)相當於ans=1+mod(10+3-1,1)=1+mod(12,1)，整數對1求余，結果為0，第1次迭代結果等於1；運算不複雜，可是1的含義有講究！它代表最後倖存者的編號索引，從時間上講，這是逆序的，第1次迭代相當於移除、或者是處決其他編號囚犯的最終結果，也就是打開牢門，把這貨放出去那一刻，從下一次迭代開始，就開始倒推了。

(2). 第2次迭代，前次迭代ans=1，故ans=1+mod(ans+k-1,i)=1+mod(1+3-1,2)=1+1=2，意味著幸運先生在最終處決誰的生死競賽決賽中，處於最後一個編號位置2，換句話說，另一個倒在最後一次的兄弟排在第1位，所以這次迭代倆人實際位置是[0 1]；

(3). 第3次迭代，再回溯一步：剩3個人時，幸運先生在這3個人中的索引位，前次迭代結果為2，所以本次ans=1+mod(ans+k-1,i)=1+mod(2+3-1,3)=2，本次移除或處決時，老哥站在第2位，此時3人的站位是[0 1 0]；

...

同理類推，1+mod(ans+k-1,i)中，i代表本次迭代還剩多少人參加抽殺;k-1代表安全的移位距離3-1=2，誰幸運地移動距離2呢？當然是ans，也就是「幸免於難」先生。最後，為和在mod命令前用1來加？實際不是+1，剛才說過，for循環是抽殺的逆過程，所以從正常時間順序，應該是減1才對！因為這正代表從序列中移除1個編號、或者說處決了某個戰俘。

這個for循環就是我今天主要想說的，這是目前位置，我見過最簡潔，也是最讓我覺得燒腦的for循環。整整想了1天才完全明白其功用。

下面問題來了：把答案攤開在面前，都看了一頭汗，這兩位怎麼想到如此精彩的演算法和構思的？

更進一步，我覺得不能因為看了些矢量化解法，就錯誤以為for/while循環都是菜鳥用的。「摘花飛葉」、「草木皆是利器」，說的就是到一定境界，並不光比其他人知道更多的函數命令與組合，而是能把簡單函數用出不簡單的味道。解決問題心中沒有各種套路約束，俯拾仰看，皆是生花妙筆。

這段代碼就是出色的演算法和簡潔的代碼，二者完美結合的範例。

延伸討論

說到for循環，不得不提它的姊妹函數arrayfun，二者一定程度上有互換性，也有一定區別，下面通過一個問題，談談arrayfun函數結合匿名函數發揮出的功效。

Pro2528. Expand Intervals

【鏈接】：https://ww2.mathworks.cn/matlabcentral/cody/problems/2528

【原題】：Youre given a row vector of an even number of monotonically increasing integers. Each pair of consecutive integers is the lower and upper bound (included) of an integer interval. Return a row vector which consists of all the integers, ordered, within these intervals.

【解釋】：給定整數端點序列數組，其元素數量為偶數，數組內每對數構成單調增加的區間，由上、下邊界構成，要求返回區間內按次序的所有間隔整數。

Test

Code Input

%%bounds = [1 5 7 9 24 32];elements = [1 2 3 4 5 7 8 9 24 25 26 27 28 29 30 31 32];assert(isequal(ExpandIntervals(bounds),elements))

%%bounds = [9 9 11 11];elements = [9 11];assert(isequal(ExpandIntervals(bounds),elements))

%%bounds = [100 200 300 400];elements = [100:200 300:400];assert(isequal(ExpandIntervals(bounds),elements))

%%temp = [-10:10; -10:10];bounds = temp(:);elements = -10:10;assert(isequal(ExpandIntervals(bounds),elements))

%%bounds = [-10 10];elements = -10:10;assert(isequal(ExpandIntervals(bounds),elements))

解釋拗口，直接看算例。

Test1給定3個區間：[1,5],[7,9]和[24,32]，要求返回數組是區間包含的所有整數(含端點)。區間[1,5]擴展為1:5、區間[7,9]擴展為7:9、區間[24,32]擴展為24:32。結果為：

[1 2 3 4 5 7 8 9 24 25 26 27 28 29 30 31 32]

輸入數組bounds中，區間統一放在行向量中，以首、末、首、末……的規律循環，這些端點要單獨提取出來，兩個常規方案：

(1). 索引號提取，例如：採用bounds(1:2:end)和bounds(2:2:end)

(2). 用reshape命令。例如：

>> a=[1 5 7 9 24 32]a = 1 5 7 9 24 32>> reshape(a,2,[])ans = 1 7 24 5 9 32

所以for循環以1:size(ans,2)設定循環次數，寫成：

function j = ExpandIntervals(bounds)a=bounds(1:2:end)b=bounds(2:2:end)j=[]for flag=1:numel(a) j=[j a(flag):b(flag)]endend

或：

function elements = ExpandIntervals(bounds) bounds = reshape(bounds,2,[]) elements = []; for i=1:size(bounds,1) elements=[elements bounds(i,1):bounds(i,2)]; endend

乾脆不在開頭分首末端點，降低循環次數：

function ans = ExpandIntervals(x) [];for i=1:2:length(x) [ans x(i):x(i+1)];endend

定數循環for改成while：

function elements = ExpandIntervals(bounds) elements=[]; i=1; while i<length(bounds) elements=[elements bounds(i):bounds(i+1)]; i=i+2; endend

如下我認為是本問題leading solution，用的仍然是for循環：

function ans = ExpandIntervals(ans) for k=reshape(ans,2,) union(ans,k:k(2)); endend

看似不起眼，實則體現作者Jan Orwat高超的代碼編寫能力，及對MATLAB流程和函數命令理解方面深邃的洞察力！

喂！吹流弊不必吹這麼凶吧？

呃...分析下就知道了：

ans既是輸入又是輸出，節省一個變數名，這點不說了，Alfonso前面問題中也用過；
省去第1行初始空矩陣做底擴展，改用union動態合併數組向量；
循環節充分利用MATLAB中，列排布優先的特徵，省掉類似索引「x(i)」這樣的中間引用步驟，直接用低維索引擴展序列併合並。聽起來似乎有點啰嗦和費解，不妨分解運行查看中間過程：

ans = 1 2 3 4 5 7 9 24 32ans = 1 2 3 4 5 7 8 9 24 32ans = 1 2…… 5 7 8 9 24 25 26…… 32

第1次循環，ans還是原輸入數組，也就是[1,5,7,9,24,32]，執行到循環內部時，變成讓這個數組和k:k(2)合併，什麼是循環內的「k:k(2)」呢？原來在MATLAB中，默認循環是執行列的，第1次循環中的k，代表著reshape的數組，即：

reshape(a,2,[])ans = 1 7 24 5 9 32

其中[1;5]當然按低維索引，k(1)=1,k(2)=5，k:k(2)就變成了1:5，讓它和ans合併，二者中間有相同元素1和5，union自動剔重，變成[1 2 3 4 5 7 9 24 32]』。

貌似簡單、人人以為自己會用的for循環，就這麼讓Jan Orwat寫出不簡單的味道。

與for類似的arrayfun也是常用命令，和for循環不同，arrayfun關鍵在於處理對象以「行數組」形式出現，處理方法必須以匿名函數形式，好處是書寫形式簡潔，符合對象編程主旨，缺點是：既然匿名函數，就不能在arrayfun內部完成更複雜的變數賦值。

不過，對Pro.2528而言，arrayfun是適用的。先說說函數具體用法，再度吐槽arrayfun函數的幫助文件，寫得簡直神頭鬼腦，直接限制了閱讀熱情。

通常這個函數有兩種用法，一是每次調用匿名函數操作array數組得到結果是1×1維度的數值(scalar)，此時其後尾隨的「統一輸出」開關默認是打開的，所有返回結果將排成一個行序列，統一輸出，例如：

>> arrayfun(@(i)i^2,1:5)ans = 1 4 9 16 25>> x=randi(10,1,10)x = 2 10 10 5 9 2 5 10 8 10>> arrayfun(@(i)nnz(x(1:i)>=3),1:numel(x))ans = 0 1 2 3 4 4 5 6 7 8

它也支持多輸出，只要匿名函數內調用方法有這個選項：

>> x=randi(10,5)x = 3 5 8 10 9 7 10 3 6 3 7 4 6 2 9 2 6 7 2 3 2 3 9 3 10>> [iX,maxX]=arrayfun(@(i)max(x(i,:)),1:size(x,1))iX = 10 10 9 7 10maxX = 4 2 5 3 5

上述命令執行結果代表5×5方陣x每行最大值iX以及其對應的索引位maxX。注意結果逐行執行，每次只處理一行數據，得到一個最大值和該最大值索引位。還有種情況是當每次輸出結果不是一個數，而是矩陣、字元串等，就需要關閉默認的統一輸出開關返回結果，有人也許問：既然返回不是數字，排列的不就不是行序列了嗎？

未必。舉個例子，接前述5×5的隨機整數矩陣x，當關閉統一輸出開關後，出現結果：

>> [ix,jx]=arrayfun(@(i)find(x(i,:)>=mean(x(:))),1:size(x,1),un,0)ix = 1×5 cell 數組 {1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double}jx = 1×5 cell 數組{1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double}

命令用於查找x每一行中，大於等於x均值數字的行列索引。顯然，每行滿足要求的索引位數量不定，且可能是1也可能不是1，而把』un』設為0，關閉統一輸出開關後，得到結果是用cell形式存儲，整體看，它們仍然是行序列，只是改成cell類型。

說到這裡，大致上函數arrayfun算說明白了，看看它在Pro.2528中的應用，首先：

function ans = ExpandIntervals(x) cell2mat(arrayfun(@(i)x(i):x(i+1),1:2:length(x), un,0));end

中規中矩，得到結果是cell，用cell2mat對結果橫向排布，一行搞定。

再看逗號操作符實現的新穎方案：

function ans = ExpandIntervals(x) cell2mat(arrayfun(@colon,x(1:2:end),x(2:2:end),uni,0));end

都是cell2mat+arrayfun，兩者有何不同？

後者內部調用的是基本操作符「@colon」，多數人用慣了「1:n」這種方式，可能都不記得還有這種等價寫法了

>> 1:3ans = 1 2 3>> colon(1,3)ans = 1 2 3

最後看調用對象，並未在第2參數出現熟悉的索引格式，例如：「1:k:numel(x)」這種方式，而直接用x(2:2:end)的實際數組代替，這仍然充分利用了循環節默認按列的順序形式，非常巧妙，同時實錘了，arrayfun和for循環之間的密切關聯：畢竟前面剛剛提及Jan Orwat那個精妙的for循環，裡面也使用了類似技巧。

小結

別瞅不上循環，for/while循環體構思和設計同樣是個技術活兒，對流程有充分認識，還得知道哪些函數能夠搭配。例如Jan大神用for+union+reshape組合得到序列擴展，彬彬和Alfonso用for循環對編號移位做回溯倒序，得到抽殺率問題的極簡解法，這些都是再具體問題中，把for循環優化到極致的經典案例。

能這樣寫for循環的人，又怎麼會寫不好矢量化代碼？！