2013-05-20 81 views
2

我試圖編寫代碼來測試n^2 + (n+1)^2是否完美。由於我在編程方面沒有太多的經驗,所以我只能使用Matlab。 到目前爲止,這是我曾嘗試快速測試n^2 +(n + 1)^ 2是否完美的方法

function [ Liste ] = testSquare(N) 

     if exist('NumberTheory') 
      load NumberTheory.mat 
     else 
      MaxT = 0; 
     end 

     if MaxT > N 
      return 
     elseif MaxT > 0 
      L = 1 + MaxT; 
     else 
      L = 1; 
     end 


    n = (L:N)';   % Makes a list of numbers from L to N 
    m = n.^2 + (n+1).^2; % Makes a list of numbers on the form A^2+(A+1)^2 
    P = dec2hex(m);  % Converts this list to hexadecimal 

    Length = length(dec2hex(P(N,:))); %F inds the maximum number of digits in the hexidecimal number 
    Modulo = ['0','1','4','9']';  % Only numbers ending on 0,1,4 or 9 can be perfect squares in hex 

    [d1,~] = ismember(P(:,Length),Modulo); % Finds all numbers that end on 0,1,4 or 9 

    m = m(d1);        % Removes all numbers not ending on 0,1,4 or 9 
    n = n(d1);        % -------------------||----------------------- 
    mm = sqrt(m);       % Takes the square root of all the possible squares 

    A = (floor(mm + 0.5).^2 == m);   % Tests wheter these are actually squares 
    lA = length(A(A>0));     % Finds the number of such numbers 

    MaxT = N; 
    save NumberTheory.mat MaxT; 

if lA>0 

    m = m(A);        % makes a list of all the square numbers 
    n = n(A);        % finds the corresponding n values 
    mm = mm(A);        % Finds the squareroot values of m 

    fid = fopen('Tallteori.txt','wt');  % Writes everything to a simple text.file 
     for ii = 1:lA 
      fprintf(fid,'%20d %20d %20d\t',n(ii),m(ii),mm(ii)); 
      fprintf(fid,'\n'); 
     end 
    fclose(fid); 

end 

end 

這將寫有與對應的n值到文件的平方。現在我看到使用十六進制是一種在C++中找到完美正方形的快速方法,並試圖在matlab中使用它。不過,我不確定這是否是最好的方法。

由於十六進制轉換,上述代碼在m > 2^52時發生故障。

是否有另一種方法/更快地將n^2 + (n+1)^2上的所有完美正方形寫入從1到N的文本文件?

回答

9

有一種更快的方式,甚至不需要測試。你需要一點點初等數論中找到這種方式,但這裏有雲:

如果n² + (n+1)²是一個完美的正方形,這意味着有一個m,使得該類型的

 m² = n² + (n+1)² = 2n² + 2n + 1 
<=> 2m² = 4n² + 4n + 1 + 1 
<=> 2m² = (2n+1)² + 1 
<=> (2n+1)² - 2m² = -1 

方程迎刃而解中,從 「最小的」(正)

1² - 2*1² = -1 

的溶液起始

x² - 2y² = -1 

對應於數1 + √2,則乘以與的

a² - 2b² = 1 

原始溶液是(1 + √2)² = 3 + 2*√2的功率獲得的所有進一步的解決方案。

寫的是矩陣形式,你獲得的所有x² - 2y² = -1解決方案,

|x_k| |3 4|^k |1| 
|y_k| = |2 3| * |1| 

和所有x_k必然是奇數,因此可以寫爲2*n + 1

最初的幾個解決方案(x,y)

(1,1), (7,5), (41,29), (239,169) 

通過

(n_(k+1), m_(k+1)) = (3*n_k + 2*m_k + 1, 4*n_k + 3*m_k + 2) 

對應(n,m)

(0,1), (3,5), (20,29), (119,169) 

你可以得到下一個(n,m)解決方案對從(n_0, m_0) = (0,1)開始。

快速Haskell代碼,因爲我不說話MatLab的:

Prelude> let next (n,m) = (3*n + 2*m + 1, 4*n + 3*m + 2) in take 20 $ iterate next (0,1) 
[(0,1),(3,5),(20,29),(119,169),(696,985),(4059,5741),(23660,33461),(137903,195025) 
,(803760,1136689),(4684659,6625109),(27304196,38613965),(159140519,225058681) 
,(927538920,1311738121),(5406093003,7645370045),(31509019100,44560482149) 
,(183648021599,259717522849),(1070379110496,1513744654945),(6238626641379,8822750406821) 
,(36361380737780,51422757785981),(211929657785303,299713796309065)] 
Prelude> map (\(n,m) -> (n^2 + (n+1)^2 - m^2)) it 
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] 

編輯由EitanT

這裏的MATLAB代碼來計算第一N數字:

res = zeros(1, N); 
nm = [0, 1]; 
for k = 1:N 
    nm = nm * [3 4; 2 3] + [1, 2]; 
    res(k) = nm(1); 
end 

The結果數組res應該保持滿足完美平方條件的值n

+3

+1:我冒昧地將MATLAB代碼添加到您的答案中。 –

+0

謝謝。我可以理解它,但我無法寫出它;) –

+0

+1:非常好的解釋。很好的答案。 – Schorsch