標籤:

Cody Note 001:除數數量的判斷

Cody Note 001:除數數量的判斷

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

馬良、祁彬彬


序言

準備做個系列:Cody寫代碼的筆記心得,一方面因為Cody精巧甚至異想天開的代碼方案層出不窮,對提高MATLAB編程水平、擴展思路和了解常用演算法的MATLAB實現,確實有很多好處;另一方面,在Cody苦逼解決問題,是個脫胎換骨淬鍊自身的過程,MATLAB代碼寫多了,常在遇到新問題時,常發現有些內容似曾相識,一時又不易搜到,所以想著幹嘛不把認為值得記錄的經典問題、或問題系列,以及多種代碼求解方案的分析寫下來,攢到一定程度,用TeX編成Cody筆記?他山之石可以攻玉,這些內容可能對初學MATLAB編程,或到達一定程度卻陷入瓶頸的人也提供一定的啟發。

關於Cody

1)什麼是「Cody」?

Cody (cn.mathworks.com/matlab) 是Mathworks公司主頁單獨列出的問題求解社區。存有大量MATLAB編程練習題目,世界各個國家的MATLAB高手基本都在此有自己的賬號,你可以解答各種各樣的編程問題,或創建自己的問題供大家求解,同時Mathworks公司專門組織Cody Team,專門維護這一社區的良性發展,其中所提出的一些有價值和啟發性的問題,甚至在明裡暗裡影響著MATLAB軟體的更新趨勢和走向,因為Mathworks代碼編程工程師基本都在此悄悄潛水。和其他論壇,或者Mathworks其他如:「Answer」、「File Exchange」等社區不同處在於:Cody不解答基礎編程疑問,也不是單純給出代碼以供下載,可以把Cody看做一個問題挑戰類的遊戲,需要做的就是在線提交相關問題的MATLAB程序,然後在問題下方給定的驗證代碼上,Pass或Fail。

2)Cody有沒有什麼門檻?

基本沒有,只需在MATLAB主頁註冊賬號,進入Cody,挑選能做的題目。

Done!

3)題目難不難?

難易不等,涉及數學、工程、生活等很多方面的問題,都有人把它們歸結為某種模型,需要按照背景知識提出見解,當然,要用MATLAB代碼體現和表達。每個人都能找到適合自己水平的題目。所謂「簡單」永遠是個相對和動態的概念,Cody最精彩和激動人心的地方就在於:能夠看到世界各地頂級MATLAB高手對一個問題的多種代碼實現方法,這對於開闊思路,啟發和提高MATLAB編程水平非常有益。因此,有時候簡單問題也有大開眼界、別開生面的求解方案,玩兒過Cody,有過比較的人都知道,一般無論編程天分有多高,學習MATLAB有多用心刻苦,但光在家拿本書閉門造車,甚至就是經常泡論壇的,也都一定會遇到某個無法突破的瓶頸階段,然後水平就停滯不前,而Cody就是突破瓶頸的最佳渠道。甚至現在我一看代碼,就知道這人是否常駐Cody、因為玩兒過的沒玩兒過的,代碼水平通常差一截,很多技巧並不出現在任何教材或者書中,都是Cody獨一份。

4)Cody有什麼規則嗎?

先談基本,一般只要按Problem順序往下解就是了,沒有什麼順序,看哪道題順眼,或看哪道題目不順眼,都可以求解。需要說明的是:解不出某道題目,則暫時不能看到其他人對這一問題的答案,很多時候得絞盡腦汁,費盡心思寫出正確代碼通過下方的testsuite;而且,即使寫出某道題目的答案,也不一定能馬上看到所有人的solutions,因為Cody Team設計了稱為「size」的計算方法,提交正確答案後,可以看到自己所寫代碼的size,比如你所寫代碼的size是24,那麼,你只能看到比24這個size值更高的其他人解決方案,如果想看到所有人的代碼,你要麼得寫出在size上更小的solution,要麼,就再做一道題目,這樣,就可以看到前一道題目的所有答案。

另外,Cody聚集著全世界當下火力最猛的MATLAB高手,地球上MATLAB編程能力最強的人基本都有賬號,並時不時偷偷上來算幾道題目,標識自己的存在感,如同猛獸圈畫自己的領地。Cody有所有Players的排名,排名靠前的兄弟們,私下都比較在意自己的Rank…

一個數有幾個正整數除數?

第1道題目:關於「獲取某個(有可能比較大的)正整數所有除數」的題目。

鏈接: cn.mathworks.com/matlab

【原題】:Inspired by Problem 1025 and Project Euler 12. Given n, return the number y of integers that divide N. For example:

with n = 10, the divisors are [1 2 5 10],

so y = 4.

Its easy with normal integer but how to proceed with big number?

【解釋】:受Project Euler(註:這是另一個受到歡迎和關注的網路解題項目,也有很多很出色的題目,經常被移植到Cody尋求MATLAB代碼方案) 第12題的啟發,要求通過編寫MATLAB代碼,獲得某個自然數所有正整數除數(divisor)的數量,但不要求把所有這些正整數都列出來。如果n=10,那麼,所有正整數除數是[1,2,5,10],一共4個。所以,我們關心的是「有幾個?」,而不是「到底哪幾個?」

求解代碼及分析

可能有人感覺整除類問題在計算機上比較簡單:既然要正整數除數,遍歷一下不就得了?例如:

f=@(n)sum(~mod(n./(1:n),1));>> f(10)ans = 4

再比如:

>> f(784)ans = 15

抬走!下一題!!!

呵呵…等等,先看看test suite:

Test

Code Input

1

%%x = 10;y_correct = 4;assert(isequal(divisors_Big(x),y_correct))

2

%%x = 28;y_correct = 6;assert(isequal(divisors_Big(x),y_correct))

3

%%x = 28;y_correct = 6;assert(isequal(divisors_Big(x),y_correct))

4

%%x = 784;y_correct = 15;assert(isequal(divisors_Big(x),y_correct))

5

%%x = 1452637;y_correct = 2;assert(isequal(divisors_Big(x),y_correct))

6

%%x = 5452637;y_correct = 4;assert(isequal(divisors_Big(x),y_correct))

7

%%x = 16452637;y_correct = 2;assert(isequal(divisors_Big(x),y_correct))

8

%%x = 116452637;y_correct = 8;assert(isequal(divisors_Big(x),y_correct))

9

%%x = 416452638;y_correct = 32;assert(isequal(divisors_Big(x),y_correct))

10

%%x = 12250000;y_correct = 105;assert(isequal(divisors_Big(x),y_correct))

11

%%x = 2031120;y_correct = 240;assert(isequal(divisors_Big(x),y_correct))

12

%%x = 76576500;y_correct = 576;assert(isequal(divisors_Big(x),y_correct))

13

%%x = 816452637;y_correct = 32;assert(isequal(divisors_Big(x),y_correct))

14

%%x = 103672800;y_correct = 648;assert(isequal(divisors_Big(x),y_correct))

15

%%x = 842161320;y_correct = 1024;assert(isequal(divisors_Big(x),y_correct))

總計15個Test,要用所寫代碼全部通過,每個Test最後一行assert命令讓求解代碼輸出和給定正確結果y_correct相符,或相差在某個閾值範圍內。從test 5,所有正整數都非常大,如果用向量除法的代碼遍歷,要麼內存不夠,要麼cpu高佔用率,目前Cody Sever規定:提交答案測試的運行總時間是1min,所以遍歷做法肯定不可行,但失敗方案也非一無是處:代碼中出現了判定何為整數的常見方案:

~mod(x,1)

數組比較小的時候,這個判斷方法很實用。當然並非唯一方案,關於整數判定的多種代碼實現方式,放在後一部分:延伸討論中詳細探討。言歸正傳,既然不能遍歷,那怎麼搞?

整除理論結合定數for或不定數while循環遍歷

容易想到的是循環,既然向量除法不方便一次放進內存處理,改成循環每次處理一個數判定,於是:

function y = divisors_Big(x) y = 0; for i = 1:round(x/2) y = y + ~mod(x,i); end y = y + 1end

循環次數不必對小於x的所有數遍歷,最小質因數是2,理論上循環一半即可。for循環的求解思路比較多,如開根號開方是否等於其本身:

function ans = divisors_Big(x)y=1:floor(sqrt(x));if y(end)^2==x k=1;else k=0;end2*sum(mod(x,y)==0)-k;end

開根號向下捨去小數部分,然後開方如果等於本身,則為除數。

下面是兩個數學上相對嚴謹的方案,第1個是利用整除理論求解:

function c = divisors_Big(x)s = num2str(x) - 0;len = length(s);c = 2;for i = 2:sqrt(x) carry = 0; for j = 1:len carry = mod( 10*carry + s(j) , i ); end if carry == 0 c = c+2; endendif (floor(sqrt(x)))^2 == x c = c - 1;endend

倒也不難理解,印度人rifat提交的演算法相當於對數字除法的模擬模擬,這是用內層循環,從最高位數字開始、求余借位迭代實現的:

carry = 0; for j = 1:len carry = mod( 10*carry + s(j) , i ); end

初值c=2代表整數(x>1)能被本身和1整除,因此結果至少是2;當對本次循環i逐位相除得到的變數迭代結果carry=0時,意味著x能夠對i整除,下方的第1個判斷中,c=c+2中的「+2」也隱藏了一個很妙的演算法,首先注意外層循環體是從2:sqrt(x),這意味著在這個範圍內找到某個數a能被x整除,則在[sqrt(x),x]範圍內必有另一個數b也能被x整除,且滿足x=a*b,所以計算1次,得到的是2個除數。這個演算法把數字除法替換成對各個位數的單獨除法以降低計算大數除法時單次運算的負擔。但整體循環次數方面沒有顯著改善,畢竟外層「for i=2:sqrt(x)」仍在對i<=floor(sqrt(x))的除數做遍歷。其他做法還有:

function ans = divisors_Big(x)i = 2;[];while x ~= 1 while rem(x, i) ~= 0 i = i + 1; end x = x / i; ans = [ans i];endbsxfun(@(X, Y) X .* Y, dec2bin(0:2^numel(ans) - 1) - 48, ans);ans(ans == 0) = 1;numel(unique(prod(ans, 2)));end

這是個先難後易、更全面的演算法:不僅算出能被x整除整數的個數,且所有符合條件的整數也通過dec2bin,10進位轉2進位函數給全部枚舉出來。演算法分成兩個部分,第1部分是嵌套2重while循環,這是尋找x分解質因數的標準演算法,當時不是很明白,為什麼沒直接用factor函數等效代替這部分的計算過程。後來發現作者又補了一個替換代碼:

function ans = divisors_Big(x)factor(x);bsxfun(@(X, Y) X .* Y, dec2bin(0:2^numel(ans) - 1) - 48, ans);ans(ans == 0) = 1;numel(unique(prod(ans, 2)));end

看來是一開始不知道有factor這個函數所致。第2部分的關鍵是用bsxfun構造了2進位數和所有質因數的單一維擴展點乘。其中利用質因數2進位轉換間接構造除數的做法非常新穎,很值得借鑒。

順便說句題外話,這種利用2進位,字元和數值類型無縫轉換構造某種條件的做法相當考驗模型構造能力,在許多骨灰級高手的代碼中經常見到,例如下面的小題目:

【鏈接】:ww2.mathworks.cn/matlab

【原題】:Draw a x-by-x matrix E using 1 and 0. (x is odd and bigger than 4)

Example:

x=5ans=[1 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 1 1 1 1 1]

x=7ans=[1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1]

這是個給定一個大於4的奇數,構造0-1矩陣,且其中的1元素形成類似「E」字母排布形式的簡單題目,例如能夠利用MATLAB彈性擴維的規則:

function ans = your_fcn_name(x) ans([1 ceil(x/2) x],1:x)=1;ans(:,1)=1;end

但J.S.Kowontan卻另闢蹊徑,按照dec2bin給出新奇的構造思路:

function ans = your_fcn_name(x) pow2(x-1); dec2bin(ans)+dec2bin(ans+sqrt(ans)+1)>96; end

要對不同進位數據轉換具有很深刻了解後,才能夠寫出這樣的代碼。

圍繞factor質因數尋找除數的快速解法

言歸正傳,回到原問題,前面提到整除問題中,和數字n的分解質因數密切相關,因此,快速的演算法都需要先把整數n所有質因數提取出來,MATLAB中獲得分解質因數的命令是「factor」,所以:

>> factor(10)ans = 2 5>> factor(784)ans = 2 2 2 2 7 7>> factor(1452637)ans = 1452637

對於n=10,分解出兩個唯一的質數2和5,加上它本身和1,結果為4;第2個n=784則分解出4個2和2個7的答案;最後一個分解後是其本身,說明是質數,此時結果是2。但第2個n=784就有點麻煩了,分解的質因數不唯一,所以又是個排列組合的問題,一共6個數,nchoosek函數依次取1,2,3,…,6個質因數,再以prod命令做乘積,再用unique命令取得唯一的除數:

function y = divisors_Big(n) f=factor(n); ps=[1;n;f] k=2; while k<length(f) ps=[ps;prod(nchoosek(f,k),2)]; k=k+1; end y=length(unique(ps));end

我覺得nchoosek這種精細的定位組合,搞法上有點煩,所以想了個簡單粗暴的辦法,數字範圍內隨機打靶:

function ans = divisors_Big(n) factor(n); numel(unique(arrayfun(@(i)prod(ans(randperm(numel(ans),randi([1 numel(ans)])))),1:100000)))+1;end

也就是說,分解質因數後,隨機選擇其中的k個質因數相乘,隨機打靶100000次,得到100000個乘積,再以unique排重,總能囊括所有整除的數字吧?!通常千萬以下的數字隨機打靶100000次左右,基本都能找到解。

當然,前述並不是分解質因數的理想方案,國棟提供了一個非常理想的矢量化代碼:

function y = divisors_Big(n) y = prod(sum(bsxfun(@eq,factor(n),unique(factor(n))))+1);end

利用bsxfun和element-wise的隱式擴展,包含對這一問題求解的最優思路。

首先分析輸入的自然數n內部包含的所有除數,究竟有什麼樣的規律,剛才提到這些除數是用factor命令計算得到的質因數組合遍歷,因此不同除數包含有不同組合質因數,因此問題變成:一個數n,究竟有多少個質因數組合?

以n=784為例,質因數為[2,2,2,2,7,7],4個2和2個7,共計6個質因數,對2,每個除數中最多可能含有:0,1,2,3,4個2,也就是說,從不包含質因數2到含有全部4個2,同樣道理,每個除數中最多含有0個、1個和2個質因數7,按照乘法定律,除數只可能有(4+1)×(2+1)=5×3=15個。

下面就談談國棟的代碼,首先把這行代碼分解,最裡面的bsxfun得到如下結果:

>> n=784;>> bsxfun(@eq,factor(n),unique(factor(n)))ans = 6×2 logical 數組 1 0 1 0 1 0 1 0 0 10 1

沒用過bsxfun估計覺得稍稍玄乎,其實適應了會感覺到該函數在實際代碼中無處不在的便利,這是用於單一維擴展,也就是隱式擴展的點對位操作函數。題外話:如果對bsxfun函數的應用感興趣,可參考我和彬彬寫的書《MATLAB向量化編程基礎精講》,目前貌似只有我倆的書結合實例詳細講過bsxfun函數的用法。

繼續討論:通過內部調用是否相等的句柄「@eq」,來逐個判斷第2個參數「factor(n)』」,也就是[2;2;2;2;7;7]和第3個參數「unique(factor(n))」的元素是否相等,如果擴展開,應該是如圖1所示:

圖1 單一維(隱式)擴展原理

左側框內是所有factor(n)的質因數,右側每行都是分解質因數經過排重之後的最簡,也就是各不相同的質因數,本例n=784時,顯然只有2和7兩個,中間紅色部分是判斷相同的句柄「@eq」對每行數據的比較結果,返回的是邏輯值,例如第1行左側框內的2和右側2和7進行比較,2和2相同,所以中間紅色部分第1個數是1,代表二者相等的邏輯「TRUE」,和右側第2個數7不等,返回0,也就是邏輯「FALSE」,其餘類推。此時發現質因數組合中,可能出現4次相同的2和2次相同的7,這裡面少了0次2和0次7,分別加上這倆1,中間紅色部分用sum求和就變成了最最需要的(4+1)和(2+1),最後再用5×3這麼一合…

DONE!

解法、思路,雙料牛逼!對MATLAB函數命令的深度了解和對演算法的熟悉,二者缺一不可。

理解了演算法,發現問題本質是對輸入n的所有質因數做頻數統計,這引出另一類採用頻數統計命令求解的方案。

MATLAB提供多種頻數統計的命令,例如histcounts,例如accumarray、例如histc等等…

為熟悉函數用法,舉個簡單例子:

>> str=aabbccccedabtttfffffffstr = aabbccccedabtttfffffff>> histcounts(+str,[sort(unique(+str)) 1000])ans = 3 3 4 1 1 7 3

或者使用accumarray:

>> accumarray((+str-96),1)ans = 3 3 4 1 1 7 0 0 0 0 0 0 0 0 0 0 0 0 0 3

或者histc命令:

>> histc(+str,sort(unique(+str)))ans = 3 3 4 1 1 7 3

結果都一致,在字元串str中,共有3個a,3個b,4個c,1個d,1個e以及7個f和3個t,不同之處是accumarray是自動把a-t分成了20個簇,或者20個垛(tokens),然後挨個統計,中間的0代表沒有出現的字母,「g-s」,但非零部分三種命令都是一樣的,必須說明,accumarray函數用法靈活,遠不止做個頻數統計這麼簡單,如果感興趣,也可以參考我和彬彬的書。之後的筆記中也將經常提到這個函數在各種問題中的應用。

言歸正傳,除數的問題就可以用這個命令搞定,例如:

prod(accumarray(factor(n)』,1)+1)

目前為止,我覺得該方案為最優代碼。當然histc得到的結果也類似:

function ans = divisors_Big(x)factor(x);prod(histc(ans,unique(ans)) + 1);end

除了用histc+factor或者factor+accumarray的,劉鵬使用新的分類函數countcats+categorical+factor提交了一個美妙的代碼:

function y = divisors_Big(x)y = prod(countcats(categorical(factor(x)))+1);end

countcats和categroical都是在R2013b版本中出現的新函數,categorical用於對數值進行categorical類型的轉換,外層接入countcats統計數量。

利用mod的矢量化操作特性

mod在整除判定時,可把第2參數設為1,如果輸出結果為0說明為整數,最美妙的是,命令支持矢量化操作,下面是劉鵬給出的圍繞求余函數的矢量化代碼方案:

function y = divisors_Big(n)y = 2*sum(~mod(n,1:n^.5)) - ~rem(n^.5,1);end

代碼分成2部分,sum(~mod(n,1:n^.5))很明顯是統計i<=sqrt(n)的數字有多少能被n整除,最開始的地方「2*」翻倍,其含義在前面印度人rifat的for循環代碼中已經提及:如果i<=sqrt(n)某個數a能被n整除,則在[sqrt(n),n]範圍內必有另一數b是n的除數,且a*b=n;第2部分「~rem(n^.5,1)」是對中間節點的是否會多加1次的處理,比如當n=15時,sqrt(15)=3.873,所以後半部分的結果為0不起作用,減了等於沒減,而前面部分的1:n^.5相當於1:floor(sqrt(n)),即到達小於3.873的最大整數3停止,因此mod的計算結果為1×3數組:

>> n=15;>> ~mod(n,1:n^.5)ans = 1×3 logical array 1 0 1

1代表能被整除的數字,所以在小於3.873的整數1:3中,1和3滿足條件,因此在3.87<=x<=15的範圍內,還有15和5滿足除數條件所以除數數量應該在sum([1 0 1])基礎上翻倍,即除數數量的輸出結果為4;但在n=16時情況發生了變化:

>> n=16;>> ~mod(n,1:n^.5)ans = 1×4 logical array 1 1 0 1

最後一位因為sqrt(n)為整數,因此結果為1,可是sqrt(16)=4的數,其另外一半的除數也是4,這樣就多算了一次,為避免這種情況,後半部分的「~rem(n^.5,1)」發揮作用,此時結果因為sqrt(n)為整數,輸出必然為1,中間節點位的除數中多餘的被扣除掉。這個方案顯然優於rifat用if條件構造的判斷:

if (floor(sqrt(x)))^2 == x c = c - 1;end

延伸討論

整數判定的幾種辦法

題目介紹完,還是有些意猶未盡,比如前文提到的整數判定,也就是說:MATLAB中如何判定某個數是整數,這裡說的整數不是整型變數,是浮點數里的數學概念,因此用isinteger來弄這事兒就錯了。比如:

>> isinteger(4)ans = logical 0

對數字「4」的判斷返回結果0,這隻能說明4不是整型變數,但不是我們通常需要的結果,判定整數可採用如下幾種思路:

開方平方法

其實前面說過了:

>> f=@(n) isequal(floor(sqrt(n))^2,n)f = 包含以下值的 function_handle: @(n)isequal(floor(sqrt(n))^2,n)>> f(4)ans = logical 1>> f(4.1)ans = logical 0

餘數為零法

這個也簡單,就是用mod或者rem命令,求一下數字n對1做除的餘數,有餘數就不是,沒有就是整數

>> g=@(n)~mod(n,1);>> g(4)ans = logical 1>> g(4.1)ans = logical 0

邏輯取反是因為如果除以1餘數為零,取反變為邏輯值1,說明是整數。

整型變數法

既然提到整型變數,不妨就利用uint8、uint16、uint64等命令的特徵實現整數的判定,如:

>> isinteger(4)ans = logical 0>> isinteger(uint8(4))ans = logical 1

雖然二者都是4,但類型不同:前者返回0是因為那是雙精度浮點數,後者才是真正的整型變數,不過,MATLAB中,對二者是否相等的判斷可從來不管是屬於哪種變數類型,只看數值在數學意義上的大小,利用這個規則,其中就包含一個實用的技巧:

>> t=@(n)uint8(n)==n;>> t(4.1)ans = logical 0>> t(4)ans = logical 1

就是比較二者是否相等,如果是uint(4.1),其返回結果是4,數值上相當於round(4.1),當然和其本身不相等。

有人可能問:我幹嘛要知道這麼多整數判定的技巧呢?有一個不就得了?來看看下面這個Cody題目:

【鏈接】: cn.mathworks.com/matlab

【原題】:Given a vector of elements, determine if each element is an integer and return true or false accordingly.

【注釋】:這是一道,紅果果地、不加掩飾地要求你判定給出的數組,其元素是否為整數的題目。

看似簡單,testsuite則有點兒棘手:

Test

Code Input

1

%%x = [1 2.2 3.3 4];y_correct = [true false false true];assert(isequal(checkIfinteger(x),y_correct))

2

%%x = [-1 2 3.3 4 -7.54];y_correct = [true true false true false];assert(isequal(checkIfinteger(x),y_correct))

3

%%x = [pi 2.2 single(3.78) -4.34 eps inf] ;y_correct = [false false false false false false];assert(isequal(checkIfinteger(x),y_correct))

4

%%x = [1.0 0 nan -6.0]; % infact nan is not a numbery_correct = [true true false true];assert(isequal(checkIfinteger(x),y_correct))

5

%%x = [true false logical(0) logical(1)]y_correct = [false false false false];assert(isequal(checkIfinteger(x),y_correct))

6

%%x = [1 2 a b !]y_correct = [false false false false false];assert(isequal(checkIfinteger(x),y_correct))

發現其中混進了相當奇怪的東西:單精度類型的single(3.78)、精度閾值eps、無窮大inf,非數NaN、邏輯值logical(0)、字元串…五花八門千奇百怪,此時求餘數、開根號之類的辦法就比較難用,首先得把前述這些元素排除,當然這也不是不能夠,MATLAB中給出許多針對性的「isa*」類型的命令,例如isfinite、isnan、ischar、isnumeric等等,但如此多的情況,不得「if」上好半天嗎?比如:

function y = checkIfinteger(x)y = false(1,length(x));for i = 1 : length(x) try if not(ischar(x(i))) && isfinite(x(i)) && floor(x(i)) == x(i) y(i) = true; end endend

再比如:

function x = checkIfinteger(x)find(isfinite(x).*(~ischar(x)).*(~islogical(x)));if ~isempty(ans) x(ans)=~mod(x(ans),1); x(setdiff(1:numel(x),ans))=0;else x=zeros(1,numel(x));endend

再比如:

function y = checkIfinteger(x)lx=length(x);y=[];for flag=1:lx j=x(flag); if strcmp(class(j),double) if j==round(j) y=[y 1]; else y=[y 0]; end else y=[y 0]; endend

各種比如:

function ans = checkIfinteger(x) x(isnan(x) | isinf(x)) = .1; if islogical(x) | ischar(x), x = x+.1; end ~mod(x,1);end

繼續各種比如:

function ans = checkIfinteger(x)try arrayfun(@(y) isfinite(y) & isnumeric(y) & floor(y) == y, x);catch str2num([0 0 0 0]);end

總之就是各種花樣百出的比如:

function y = checkIfinteger(x) y = arrayfun(@(x) isnumeric(x) && str2num(sprintf(%1.f, x)) - x == 0, x);end

而知道符號和無符號整型的命令後,並了解isnumeric的用法代碼就簡單多了。

function ans = checkIfinteger(x) int8(x)==x&isnumeric(x);end

再如:

function ans = checkIfinteger(x) isnumeric(x)&mod(x,1)==0;end

隱式擴展的進一步探討

前面提到bsxfun判定除數數量的精妙解法:

bsxfun(@eq,factor(n),unique(factor(n)))

自2016b版本後,Mathworks公司技術人員仿照Python,把單一維,也就是隱式擴展寫進普通矩陣運算,因此上述代碼可以在新版本中得到等效寫法:

~(factor(n)-unique(factor(n)))

這種方式在老版本會提示錯誤的,因為兩個減數和被減數的維數不同,一個是7×1,一個是1×2,而新版中則自動對後者按前者的行維度進行了擴展,擴展為7×2。

這件事情告訴我們,學習知識要與時代同行,不能刻舟求劍,死抱著老知識不放。

小結

本篇是CodyNote系列筆記中的第1篇,開篇的話並沒有,也不需要討論很難理解或者很難寫出解決代碼的問題,但就這樣,其中仍包含了對基本for、while循環、if判斷的基礎代碼能力的考查,還有bsxfun隱式擴展、元素數量統計和計算函數histc、accumarray、進位轉換dec2bin等命令在代碼中的綜合應用,這些代碼的應用,通過test suite中的幾個超大數的例子,已隱隱指向了對演算法效率的深入思考。用到的函數如下圖所示:

圖2 CodyNote001代碼方案中用到的主要函數

圖2看出同一個問題能夠選擇的MATLAB函數是豐富的。說明演算法是一條線,函數的理解運用是另一條線,這2條線索隱藏在每個精心構思的代碼中,而操縱這2條線索的,是屏幕後面那顆不服輸而耐心的頭腦。

Cody的精髓在於多種思路的碰撞,其中包含演算法,也有命令組合方面的創新,這屬於MATLAB編程深水區的保留節目。真正鑽研頂級高手的優秀代碼並深入分析,是讓MATLAB的代碼水平能獲得快速的提高的必由之路。


推薦閱讀:

MATLAB應用舉例
CodyNote007:超大整數的整除
Matlab編程實踐(一)
CodyNote002:富有彈性的數組長度(Part.1)
MATLAB筆記(3.2)循環

TAG:MATLAB |