CodyNote007：超大整數的整除

08-23

CodyNote007：超大整數的整除

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

馬良，祁彬彬

序言

按realmax和intmax兩個命令幫助的解釋，有：

>> format long g>> [t1,t2]=deal(realmax(double),realmax(single))t1 = 1.79769313486232e+308t2 = single 3.402823e+38

符號64位整型int64和無符號64位整型的最大數限定：

>> [i1,i2]=deal(intmax(uint64),intmax(int64))i1 = uint64 18446744073709551615i2 = int64 9223372036854775807

那麼有問題：假如整數大到超過realmax限定，不妨稱其為超大整數。這樣的數字如何判斷是否能被某個整數(很小)整除？

本篇打算通過Cody題目，探討如何判定超大整數是否能被其他數字整除的問題。

Pro.42455. Divisible by n, prime divisors-11,13,17,&19

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

【原題】：Divisibility checks against prime numbers can all be accomplished with the same routine, applied recursively, consisting of add or subtract x times the last digit to or from the remaining number. For example, for 13, add four times the last digit to the rest:

2392: 239 + 4*2 = 247: 24 + 4*7 = 52: 5 + 4*2 = 13 -> 2392 is divisible by 13.

For 17, subtract five times the last digit from the rest:

3281: 328 - 5*1 = 323: 32 - 5*3 = 17 -> 3281 is divisible by 17.

For 19, add two times the last digit to the rest:

16863: 1686 + 2*3 = 1692: 169 + 2*2 = 173: 17 + 2*3 = 23: 2 + 2*3 = 8 -> 16863 is not divisible by 19.

And, for 11, subtract the last digit from the rest:

269830: 26983 - 0 = 26983: 2698 - 3 = 2695: 269 - 5 = 264: 26 - 4 = 22: 2 - 2 = 0 -> 269830 is divisible by 11.

Write a function to return a true-false vector for the prime numbers in the 11:20 range ([11 13 17 19]) based on a number supplied as a string.

Restrictions on Java, mod, ceil, round, and floor are still in effect.

【解釋】：給定一個整數n，判斷n是否能被質數11、13、17和19整除。

Test

Code Input

%%filetext = fileread(prime_divisors_11_to_20.m);assert(isempty(strfind(filetext, mod)),mod() forbidden)assert(isempty(strfind(filetext, round)),round() forbidden)assert(isempty(strfind(filetext, rem)),rem() forbidden)assert(isempty(strfind(filetext, ceil)),ceil() forbidden)assert(isempty(strfind(filetext, floor)),floor() forbidden)assert(isempty(strfind(filetext, java)),java forbidden)

%%n = 143;tf = [1 1 0 0]; %[11 13 17 19]assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 187;tf = [1 0 1 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 221;tf = [0 1 1 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 247;tf = [0 1 0 1];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 46189;tf = [1 1 1 1];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 2133423721;tf = [1 1 1 1];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 233296158667;tf = [1 1 1 1];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 1011001000101010101010110101001010101001010101001001011010101000101010101010101010010101010010101010100101010101001100101010010101;tf = [0 0 0 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 1011001000101010101010110101001010101001010101001001011010101000101010101010101010010101010010101010100101010101001100101010010103;tf = [0 1 0 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 1011001000101010101010110101001010101001010101001001011010101000101010101010101010010101010010101010100101010101001100101010010107;tf = [0 0 0 1];assert(isequal(prime_divisors_11_to_20(n),tf))

n = 14300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000;tf = [1 1 0 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 14300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001;tf = [0 0 0 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%%n = 22100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000;tf = [0 1 1 0];assert(isequal(prime_divisors_11_to_20(n),tf))

%% anti-cheating testind = randi(4);switch ind case 1 n = 221; tf = [0 1 1 0]; case 2 n = 233296158667; tf = [1 1 1 1]; case 3 n = 46189; tf = [1 1 1 1]; case 4 n = 247; tf = [0 1 0 1];endassert(isequal(prime_divisors_11_to_20(n),tf))

%% anti-cheating testind = randi(4);switch ind case 1 n = 187; tf = [1 0 1 0]; case 2 n = 143; tf = [1 1 0 0]; case 3 n = 221; tf = [0 1 1 0]; case 4 n = 233296158667; tf = [1 1 1 1];endassert(isequal(prime_divisors_11_to_20(n),tf))

%% anti-cheating testind = randi(4);switch ind case 1 n = 2133423721; tf = [1 1 1 1]; case 2 n = 46189; tf = [1 1 1 1]; case 3 n = 187; tf = [1 0 1 0]; case 4 n = 247; tf = [0 1 0 1];endassert(isequal(prime_divisors_11_to_20(n),tf))

【分析】：前8個測試還算正常，第9個開始瘋，Test11和Test12的整數位甚至達到503位，數位多到被迫用word字數統計功能來統計的地步，遠超realmax定義。

對了，Cody規定不能使用符號數學工具箱，這2個限制條件堵死了拿計算器、甚至計算機一般四則運算判定整除的傳統路數。也迫使代碼轉從其他角度，比如數字整除理論或其他一些特殊分析手段完成這項工作。但整除判定是個變化多端的bitch，有的很簡單，例如2的整除查看末尾數是否為偶數；3的整除看數字所有的位數相加是否能被3整除；4的整除看末兩位數是否是0:4:99中的任何一個，這種情況用ismember+str2num通常可以搞定；5的整除則看末尾是否0或者5……

需要判斷整除的數字太多個，似乎還各有各的判斷規則，例如被7、11、13整除的數，存在共同特徵：把原數從千分位處分開，二者相減，如所得數能被三數整除，則原數也可以，當然，這個原則是可以連續迭代的，比如判斷71858332:

71858-332=71526

526-71=455

數字455就很容易判斷：，能被7和13，而不能被11整除。

因此，goc3利用這點，一口氣出了大約20來個整除類型的題目(有些似乎功能重複未算在列)，這其中有prime number的，有composite numer的，本例要求的是能被11、13、17、19這4個質數整除的統一解。

求解代碼

先說明為什麼不能直接用四則運算操作，先看如下除法的執行結果：

>> 2416577611356487953168413234558596937377133784871649154573851275246273/13ans = 1.8589e+68

在mod、round等函數在Test1中被明確禁止使用的情況下，整除判定有一定難度，通過format再增加幾個位數呢？

>> format long g>> 2416577611356487953168413234558596937377133784871649154573851275246273/13ans = 1.85890585488961e+68

呃...還是不行。

Cleve Moler專門撰寫過關於MATLAB在數值精度方面的文章，把這種浮點形式的整數稱之為「flint」，意思是「float integer」，比如數字5，在MATLAB中判定其距離真正的「5」到底有多遠，這是個問題，因為在數學意義上的「5」和MATLAB中的「5」之間，有著細微的差別，多細微呢：

>> eps(5)ans = 8.88178419700125e-16

換句話說，在：「5」和「5+eps(5)」間，已經沒有其他雙精度數了。這樣的話，整除判定的意義就顯得有些複雜。更噁心的是，測試代碼中給出的輸入變數n，全都是字元串形式，因為雙精度格式數字無法顯示前述長度動輒以「行」為單位的超級大數，用str2num、str2double這類字元串轉換數字的命令，都無法給出超過16位的整數，例如：

>> a=sprintf(%d,randi(10,1,80))a =48478514532959484984381047598291069622589846122752521961076991019710659351069863717789>> str2num(a)ans = 4.84785145329595e+85>> str2double(a)ans = 4.84785145329595e+85

Cody中不允許使用符號運算工具箱，貌似最後一條輕鬆的路也被堵死。

怎麼辦？

分情況討論

不少人(包括我自己在內)，被題干中那一套煞有介事的整除理論給帶進溝里，半天沒爬出來。吭哧吭哧分4種不同情況討論，被11整除，blahblah…被13整除，blahblah…，流程沒問題，貌似唯一的可行方案，問題是：面對4種不同規則，程序實在太長了…。

function tf = prime_divisors_11_to_20(n)[n11 n13 n17 n19]=deal(n-0);while numel(n11)>4 n11=[n11(1:end-2) n11(end-1)-n11(end)]; while any(n11<0) smaller=find(n11<0); n11(smaller)=n11(smaller)+10; n11(smaller-1)=n11(smaller-1)-1; endendp11=polyval(n11,10);i11=0:11:10000;while numel(n13)>4 n13=[n13(1:end-2) n13(end-1)+4*n13(end)]; while any(n13>9) bigger=find(n13>9); n13(bigger)=n13(bigger)-10; n13(bigger-1)=n13(bigger-1)+1; endendp13=polyval(n13,10);i13=0:13:10000;while numel(n17)>4 n17=[n17(1:end-2) n17(end-1)-5*n17(end)]; while any(n17<0) smaller=find(n17<0); n17(smaller)=n17(smaller)+10; n17(smaller-1)=n17(smaller-1)-1; endendp17=polyval(n17,10);i17=0:17:10000;while numel(n19)>4 n19=[n19(1:end-2) n19(end-1)+2*n19(end)]; while any(n19>9) bigger=find(n19>9); n19(bigger)=n19(bigger)-10; n19(bigger-1)=n19(bigger-1)+1; endendp19=polyval(n19,10);i19=0:19:10000;tf=~[isempty(find(i11==p11,1)) isempty(find(i13==p13,1)) ... isempty(find(i17==p17,1)) isempty(find(i19==p19,1))];end

James寫的這個程序其實算夠簡練了，就這樣A4紙也放不下，畢竟架不住容納4個不同的整除規則。幸運的是：總有人不樂意屈從於重複勞動的低級趣味，於是諸多通用和半通用解法應聲覺醒。

與最大公約數有關的解法

J.S.Kowontan提出利用最大公約數命令gcd來構造通用解法，以質數11的整除為例：

function ans = divisible_by_11(n_str) n_str-0; gcd(11,sum(ans(1:2:end))-sum(ans(2:2:end)))==11;end

第1句把字元串形式數字換成數值型，但注意，不是變成一個整體的數，而是每個位數上只有一個單一的[0-9]整數，例如：

>> 123-0ans = 1 2 3

下面一句中，可能是利用了某個整除數理論的內容，讓奇數數字之和減去偶數數字之和，看它和11的最大公約數是不是11，也就是能否整除的判斷。再看15的：

function ans = divisible_by_15(n_str) gcd(3,sum(n_str))==3&ismember(n_str(end),05);end

15是合數，分成兩個公約數5和3分別處理，3的部分用了通用解法，5的部分比較簡單，還有末尾為0的整除已經被前面的通用部分包括了。

再看7的：

function ans = divisible_by_7(n_str) str2num(fliplr(char(regexp(flip(n_str),[0-9]{1,3},match)))); gcd(7,sum(ans(1:2:end))-sum(ans(2:2:end)))==7;end

最後一句核心代碼沒變，但前面的處理方式有點複雜：實際上就是利用正則表達式，把源數字按照3個一組，千分位的形式，分成許多「垛」，按照7整除的理論，用gcd來處理整除判定。13的也是類似：

function ans = divisible_by_13(n_str) str2num(fliplr(char(regexp(flip(n_str),[0-9]{1,3},match)))); gcd(13,sum(ans(1:2:end))-sum(ans(2:2:end)))==13;end

應當指出，涉及gcd的解法並非通用，原作者J.S.Kowontan在某個解法下專門以注釋形式給出了自己的解法在某些情況下不適用的說明，個人認為使用gcd是mod等函數被禁止使用的從權之舉，核心仍然是之前諸多不同形式的整除理論。

那麼到底有沒有徹底的通用解法呢？

右移消減

Jan Orwat按除法原則，提出針對不同數字判定整除特性的通用解法。

function ans = divisible_by_11(n_str) 0; for digit = n_str-0 10.*ans + digit; while ans >= 11 ans - 11; end end ~ans;end

寥寥幾行，事情安排得明明白白，注意：儘管還沒有測試，但感覺這個演算法應該對所有兩位以內數字整除都有效。

自左至右，從輸入數字最大位數開始，用它做10位，和相鄰位構成單獨數字，不定數while內部次第減去需要判斷是否整除的數字n，以數11的整除判定為例，用去掉多個11後的剩餘數m(m<11)乘以10做十位數和右側相鄰位數相加繼續上述循環，首行的0是起始補位，即初值0*10+d1，n_str為最左側第1位數，循環到最末尾，消減掉所有能消減的11後，剩餘為0說明能整除，邏輯取反得到最終結果；

可能有人不明白為什麼每次參與去掉11的都是2位或者個位數，這麼大的數字每次只用2-3位逐步消去被整除的數字，比如本例中的11，這樣做靠譜嗎？

很靠譜！相當靠譜！

為什麼？首先說明，每回while中去掉的11，有時候是11，有時則是110，1100,11000，… 比如找個小點兒的數字，2797來判定是否能整除11，上述程序的運行結果如下：

(1). 第1次循環，左端補0，10*ans+digit的結果是10*0+2=2，小於11，因此跳過下面的while，直接進入第2次循環；

(2). 第2次循環，左端為10*2+7=27，注意這裡的27其實不是27，從其在原輸入中的位數來看，應該是2700！然後再while循環中，由於大於11，先去掉11，這個11在整體數字中是1100，肯定可以被11整除，所以相當於從原數字中減去1100，剩下16，代表1600，再次while，剩下的是下次循環中需要的ans=5(5<11)，即代表原來輸入中的500，它和下面十位上的數字「9」，組成新數字590，同理，通過while逐次減11，此時代表從原數字逐次減去110，循環5次去掉5個110，變為ans=4，代表原輸入還剩40，再回到for循環和最後的個位數7結合，進入while，去掉4個11，剩下最後的數字3；

(3). 此時數字已經全數用完，因此脫離for循環進入最後一行，最後一行就是個邏輯取反，數字3代表的是2797對11的餘數，取反結果為0，因此原數2797不能被11整除；反之如果for循環中最後出來的數是0，取反就是1，意味著原輸入數字2797可以整除11.

是不是很妙的做法？把一個長長的大數字，三下五除二，連削帶砍，碎屑橫飛，變成兩位或者三位數對除數的減法迭代運算，這真是相當牛叉的思路！

字元串消去法

大神alfonso提出另一方案，思路和Jan Orwat類似，但每次循環用字元串重組方式消減輸入數字：

function ans = divisible_by_7(str)1;while ans 0; for n=find(str) if str2num(str(1:n))>=7, str=[num2str(str2num(str(1:n))-7) str(n+1:end)]; 1; break; end endendall(str==0);

以對7的整除為例，內層的for循環用到if條件，當大於7時，相當於前面說法中的減11，一直減不能減為止。例如2797是否能被7整除，就是先不斷地從原來數字中扣除700，扣到不夠扣為止，也就是變成2797-3*700=697，開始從這個數基礎上不斷地扣70，一直扣到不能扣為止，也就是697-9*70=67，然後再繼續不停地扣7，一直扣到67-9*7=4，等到從while循環中出去的時候，結果是』4』，最後一行「all(str==0)」的執行結果當然為0，也就是不能整除。這段代碼中，while循環結束判定是個幌子，重點是裡面循環位數和if語句中的跳出循環break，它是當用完所有數字後，直接從for循環跳到外面while循環，如果不進if，ans的結果是0，所以for循環執行完畢時，ans不變化，這是再次跳出while循環的flag，因此str順順噹噹被最後一句all(…)給接收了。

個人感覺這是段充滿金屬質感的古典味兒代碼。

擴展

Test 1禁止使用java，菜鳥會感到費解，好像禁不禁止就那麼回事兒。其實MATLAB和java兼容性良好，甚至可以在MATLAB的命令窗口中直接使用java代碼，Java具有一些十分實用的庫，比如java的GUI窗口控制項就是其中之一，由於Cody沒有引入符號計算工具箱，Java數學庫里的超大整數(BigDecimal 類型)計算就算是有了用武之地：

function ans=JavaInteger(n_str,n)import java.math.*t1=BigDecimal(n_str);t2=BigDecimal(n);char(t1.divide(t2,3,BigDecimal.ROUND_HALF_UP))end

第1句引入java的math庫，裡面有超大整數計算程序包，然後定義數字類型和計算。看幾組對一個74位大數做除法的運算結果，最後保留了3位小數，當然想保留更多位數也沒問題：

n_str = 10100101010110101010010110101001011010100101001011010100110101001010010110;>> JavaInteger(n_str,13)ans = 776930846931546231539239238538539308469238538539308469239238538539231546.923>> JavaInteger(n_str,19)ans = 531584263690005316316321584263211105794742157947947900005794789526842637.368>> JavaInteger(n_str,7)ans = 1442871572872871572858587157285858715728585857287287157158585857287144301.429>> JavaInteger(n_str,11)ans = 918191000919100091819100918272819182736372818273728190919100091000910010.000

結果發現大整數對11用3位小數判斷是不夠的，7、9和13肯定不行。

小結

整除特性實際上又是利用for/while循環，逐步迭代的過程，終值為0則判定為整除，這個思路徹底越過各種振振有辭的整除理論，充分利用計算機完成大劑量重複工作的特有能力，連消帶打地，獲得了最終結果。

問題求解難點在於構造合理的迭代格式，把原來很大的數字，逐步分塊消除，露出最後的真相。之所以找到了Jan Orwat的解法，是當時我挖空心思琢磨各種數字除法，苦逼之餘，無意中翻起goc3整除系列的其他題目，發現JO這廝所有題目，竟一套代碼通殺！接著就好奇，然後就分析，再接著擊節讚歎，最後長嘆一聲：還是丫夠雞賊啊...