求 a^2+ab+b^2=N^2的正整数解

waynebuaa

lin2009

用C比较快（与10#的M$，及7#的Matlab比较）。
这是我编的第一个稍微像样的C程序。速度应该还可以快一点。仅供参考。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<string.h>

int gcd(int x, int y);

void main()
{
  double N;
  double mm;
  int k;
  double NN;
  double a = 0;
  double b = 0;
  FILE *fid1;
  
  
  N  = 50*10*10;//*10*10*10*10*10*10*10; // 10^3 不等于 1000 ? C没有^运算符？
  mm = floor(sqrt(N*N/3));   // N^2 不等于 N*N ?
  
  fid1 = fopen("mydat5.txt","w+");
  
  k = 1;
  for (a = 1; a<=mm; a++)
    {
      double tmp1 = 4*N*N-3*a*a;
      double tmp2 = (-a+sqrt(tmp1))/2; //
      for (b = a ; b<=tmp2; b++)       //
        {
          if (gcd(a,b)==1 )
            {
              NN = sqrt(a*a+a*b+b*b);
              if (NN == floor(NN))
                {
                  if (fmod(k,5) == 0)
                    {
                      fprintf(fid1,"(%5g,%5g,%5g)\t\n",a,b,NN);
                    }
                  else
                    {
                      fprintf(fid1,"(%5g,%5g,%5g)\t",a,b,NN);
                    }
                  k++;
                }
            }
        }
    }
  printf("%d\n",k-1);
  fclose(fid1);
}


int gcd(int x, int y)
{
  int f;
  
  if (x>y)
    {
      int tmp = x;
      x = y;
      y = tmp;
    }
   
  if (x==1)
    {
      return(1);
    }
   
  f = y - x;
  if (f==1 || f==0)
    {
      return (f);
    }
  else
    {
      y = f;
      gcd(x,y);
    }      
}

lin2009

N = 10^8 = a^2+ab+b^2 >= 3*ab
ab<= 10^8/3 当 a=b时，a=b=5773.5027。
所以在a，b均小于5773时，所有互质的数对(a<b)都符合条件。

suancaibao

2# lin2009 你是来搞笑的吧

lin2009

本楼有误，下面已更正

楼上给的答案不对按照（3）1不是质数也不是合数，它和任何一个自然数在一起都是互质数。如1和9908。 http://zhidao.baidu.com/question/5591912.html 在N=10^6 的情况下，有183793组。详细的数据文件太大近3M。仅将最后的几组粘上来：
997077, 572, 581
986707, 573, 574
988429, 573, 575
991879, 573, 577
993607, 573, 578
997069, 573, 580
998803, 573, 581
990151, 574, 575
993603, 574, 577
997063, 574, 579
993601, 575, 576
995329, 575, 577
997059, 575, 578
998791, 575, 579
997057, 576, 577

lin2009

Matlab的参考程序：
clear all;
clc;
N  = 10^3;
mm = floor(sqrt(N^2/3));
tic
fid1 = fopen('mydat5.m','w+');
k = 1;
for  a  = 1:1:mm
    for b = a:1:floor(1/2*(-a+sqrt(4*N-3*a^2)))
        if gcd(a,b) == 1
            NN = sqrt(a^2+a*b+b^2);
            %             if NN == floor(NN) && gcd(a,NN)==1
            if NN == floor(NN)
                if mod(k,5) == 0
                    fprintf(fid1,'(%-5d,\t%-5d,\t%-5d)\n',a,b,NN);
                else
                    fprintf(fid1,'(%-5d,\t%-5d,\t%-5d)',a,b,NN);
                end
                k = k + 1;
            end
        end
    end
end
fclose(fid1);
toc

Elapsed time is 172.087603 seconds.

lin2009

8# waynebuaa
是误解了楼主的意思啦。
可以看出 a,b,N三个数互质：在N=10^6时只有137个。
(3   , 5   , 7  )    (5   , 16  , 19 )    ( 7   , 8   , 13 )    (7   , 33  , 37 )   (9   , 56  , 61 )
(11  , 24  , 31 )    (11  , 85  , 91 )    ( 13  , 35  , 43 )    (13  , 120 , 127)   (15  , 161 , 169)
(16  , 39  , 49 )    (17  , 63  , 73 )    ( 17  , 208 , 217)    (19  , 80  , 91 )   (19  , 261 , 271)
(21  , 320 , 331)    (23  , 120 , 133)    ( 23  , 385 , 397)    (24  , 95  , 109)   (25  , 143 , 157)
(25  , 456 , 469)    (27  , 533 , 547)    ( 29  , 195 , 211)    (29  , 616 , 631)   (31  , 224 , 241)
(31  , 705 , 721)    (32  , 45  , 67 )    ( 32  , 175 , 193)    (33  , 800 , 817)   (35  , 288 , 307)
(35  , 901 , 919)    (37  , 323 , 343)    ( 40  , 51  , 79 )    (40  , 77  , 103)   (40  , 279 , 301)
(41  , 399 , 421)    (43  , 440 , 463)    ( 47  , 528 , 553)    (48  , 407 , 433)   (49  , 575 , 601)
(53  , 675 , 703)    (55  , 57  , 97 )    ( 55  , 728 , 757)    (56  , 115 , 151)   (56  , 165 , 199)
(56  , 559 , 589)    (59  , 840 , 871)    ( 61  , 899 , 931)    (64  , 221 , 259)   (64  , 735 , 769)
(65  , 88  , 133)    (69  , 91  , 139)    ( 72  , 203 , 247)    (72  , 935 , 973)   (75  , 112 , 163)
(80  , 357 , 403)    (85  , 168 , 223)    ( 87  , 160 , 217)    (88  , 315 , 367)   (88  , 437 , 487)
(93  , 187 , 247)    (95  , 217 , 277)    ( 104 , 105 , 181)    (104 , 451 , 511)   (104 , 621 , 679)
(105 , 247 , 313)    (105 , 272 , 337)    ( 111 , 280 , 349)    (112 , 725 , 787)   (115 , 333 , 403)
(119 , 145 , 229)    (120 , 611 , 679)    ( 123 , 352 , 427)    (129 , 391 , 469)   (133 , 192 , 283)
(135 , 473 , 553)    (136 , 209 , 301)    ( 136 , 795 , 871)    (141 , 475 , 559)   (144 , 155 , 259)
(145 , 552 , 637)    (147 , 520 , 607)    ( 152 , 273 , 373)    (155 , 637 , 727)   (159 , 616 , 709)
(161 , 304 , 409)    (165 , 667 , 763)    ( 165 , 728 , 823)    (175 , 369 , 481)   (176 , 259 , 379)
(177 , 775 , 877)    (183 , 832 , 937)    ( 184 , 425 , 541)    (185 , 231 , 361)   (189 , 440 , 559)
(200 , 513 , 637)    (203 , 517 , 643)    ( 205 , 299 , 439)    (208 , 387 , 523)   (215 , 336 , 481)
(217 , 600 , 733)    (231 , 689 , 829)    ( 232 , 713 , 853)    (235 , 416 , 571)   (240 , 253 , 427)
(240 , 287 , 457)    (240 , 539 , 691)    ( 245 , 459 , 619)    (248 , 825 , 973)   (264 , 325 , 511)
(265 , 551 , 721)    (272 , 715 , 883)    ( 275 , 301 , 499)    (280 , 423 , 613)   (295 , 704 , 889)
(297 , 368 , 577)    (305 , 759 , 949)    ( 312 , 493 , 703)    (319 , 441 , 661)   (320 , 583 , 793)
(329 , 351 , 589)    (336 , 589 , 811)    ( 341 , 520 , 751)    (360 , 767 , 997)   (371 , 480 , 739)
(377 , 400 , 673)    (385 , 527 , 793)    ( 385 , 696 , 949)    (403 , 477 , 763)   (413 , 627 , 907)
(427 , 680 , 967)    (429 , 560 , 859)    ( 448 , 495 , 817)    (455 , 649 , 961)   (464 , 561 , 889)
(531 , 544 , 931)    (549 , 595 , 991)
Elapsed time is 137.031043 seconds.

ljwxhlzp

花50秒做到5000，离10^6还差很远。

1. m = Table[{n, n^2}, {n, 1, 10000}];
   
2. ClearAll[ss];
   
3. Evaluate[ss /@ m[[All, 2]]] = m[[All, 1]];
   
4. Timing@Length@Union@(Reap@Do[
   
5.        If[
   
6.         IntegerQ[
   
7.           ss[4 k^2 - 3 a^2]] && (b = (ss[4 k^2 - 3 a^2] - a)/2) > a &&
   
8.           IntegerQ[b], Sow[{a, b, k}/GCD[a, b, k]]],
   
9.        {k, 2500, 5000}, {a, 1, Floor[k/Sqrt[3]]}])[[2, 1]]
   

复制代码

lin2009

11# waynebuaa
算到{100000000, 13783174} 我的算法／程序可能不合适。
请楼主在合适的时候公布一下算法。

lin2009

算法不同，结构相似。也是用双重循环，花了多长的时间？M$在跑循环的效率怎样？

lin2009

算法才是关键，但算法也有快慢、繁简的问题。即效率的问题。
就这道题而言，算法很简单，但若要非要算出N<=10^8   13780042组那就不得不考虑算法／语言的效率和快慢。

lin2009

int gcd(int x, int y)中用
f = fmod(y,x)
代替
f = y - x;
步骤更少，速度更快。

FreddyMusic

很多时候我们是孤独的，我们要学会品味这种孤独和寂寞。

不要放弃，继续向你的方向前进、前进、再前进。

lin2009

21# waynebuaa

可以，算法如下，与1#问题解的算法差不多。
比较感兴趣的是，这个通式怎么得出来的，能用M$,或Maple等数学软件推导出来吗？


根据 u^2+uv+v^2 = N，v >= 1（ N = u^2+uv+v^2 > u^2+u + 1 ）确定最大的u值umax，确定u的范围： 1<= u <= umax; （第一层循环）；
再根据 u>v， N = u^2+uv+v^2  确定v 的最大值 vmax = (sqrt(-3*u^2+4*N)-u) / 2，进而确定v的范围。 （第二层循环）；
然后根据Gcd(u,v)=1, u≠v(mod3) 确定某一u值的v值的集合（循环主体）；

lin2009

根据通式编程：

#include<stdio.h>
#include<math.h>
#include<string.h>

int gcd(int x, int y);

void main()
{
  int N = 1000;
  
  int k = 1;
  int a = 0;  
  int b = 0;
  int amax = 0;
  int bmax = 0;
  
  amax = floor( (-1+sqrt(-3+4*N))/2 );
  // bmax = floor( (sqrt(-3*u^2+4*N)-u)/2 )
  
  FILE *fid1;
  fid1 = fopen("mydat6.txt","w+");
  
  for (a = 1; a <= amax; a++)
    {
      bmax = floor( (sqrt(-3*a*a+4*N)-a)/2 );
      for (b = a; b <= bmax; b++)
        {
          if ( gcd(a,b) == 1  &&  (a%3) != (b%3) )
            {
              if (k % 5 == 0)               // 每行5列（5组数据）
                {
                  fprintf(fid1,"(%5d,%5d,%5d)\t\n",a,b,a*a+a*b+b*b);
                }
              else
                {
                  fprintf(fid1,"(%5d,%5d,%5d)\t",a,b,a*a+a*b+b*b);
                }
              k++;
            }
        }
    }
  printf("\n\n\t\t总个数为：%d\n\n\n\n",k-1);
  fclose(fid1);
}

int gcd(int m, int n)
{
  int tmp;
  while (m != 0)
    {
      tmp = m;
      m   = n % m;
      n   = tmp;
    }
  return n;
}

fsm

哈哈,这个我也能看懂它的语法哦,