Project Euler Problem 276, 268, 347

lastwinner

我认为Problem 370应从b找起
只有b是非素数平方或立方的合数，才有可能是满足题设的b（素数平方或立方一定不满足题设，易证）

〇〇

tree_new_bee 发表于 2012-2-8 23:53
采用上面的算法， 计算与OO相同的参数(每边不超过10000， 周长小于1e6), 时间只要2秒左右了。

with targ ...

p=q=>a=b=c恒成立，不用判断

tree_new_bee

〇〇 发表于 2012-2-9 07:47
p=q=>a=b=c恒成立，不用判断

恩。
我说的是你那个结果8341其实不完全。没把a=b=c的算进来。

tree_new_bee

改成pl/sql 后：

1. create or replace procedure Q370(arg number default 1e6) is
   
2. peri number;
   
3. cnt number :=0;
   
4. begin
   
5. 
   
6.   for q in 1.. sqrt(arg) loop
   
7.     for p in ceil(q/2) .. q loop
   
8.       if pkg_prim2.gcd(p,q) = 1 and p+q > q*q/p then
   
9.       peri := p*p + p*q + q*q;
   
10.       if peri <=arg then
    
11.             cnt:= cnt+ floor(arg/peri);
    
12.       end if;
    
13.       end if;
    
14.     end loop;
    
15.   end loop;
    
16.   
    
17.   dbms_output.put_line(cnt);      
    
18. end Q370;

复制代码

SQL> exec q370(1e6)

861805

PL/SQL procedure successfully completed

Executed in 0.827 seconds

SQL> exec q370(1e7)

9712598

PL/SQL procedure successfully completed

Executed in 9.157 seconds

这样性能又提升了一个数量级， 但比起题目要求还差的远。

对上面的代码做profile, 可以看出时间消耗主要在gcd上。

其实，回头看94#列出的p/q, 就可以发现q要使用3以上的所有整数， 所以关键是要高效的找到与q互质的p,
这个问题也许要需要一些以前的关于真分数的题目的结论。


也因为这个原因， 我觉得project euler的题目，还是从前往后挨个做更可取一点。

〇〇

tree_new_bee 发表于 2012-2-9 08:47
改成pl/sql 后：SQL> exec q370(1e6)

861805

去掉相等的，基本没有影响

create or replace procedure Q3702(arg number default 1e6) is
peri number;
cnt number :=0;
begin

  for q in 1.. sqrt(arg) loop
    for p in ceil(q/2) .. q-1 loop
      if pkg_prim2.gcd(p,q) = 1 and p+q > q*q/p then
      peri := p*p + p*q + q*q;
      if peri <=arg then
            cnt:= cnt+ floor(arg/peri);
      end if;
      end if;
    end loop;
  end loop;
  
  dbms_output.put_line(cnt+trunc(arg/3));      
end Q3702;
/
exec q3702(1E6)

〇〇

〇〇 发表于 2012-2-9 09:17
去掉相等的，基本没有影响

create or replace procedure Q3702(arg number default 1e6) is

改变q的上界，q最大不超过1/2周长
create or replace procedure Q3702(arg number default 1e6) is
peri number;
cnt number :=0;
begin

  for q in 1.. sqrt(arg/2) loop
    for p in ceil(q/2) .. q-1 loop
      if pkg_prim2.gcd(p,q) = 1 and p+q > q*q/p then
      peri := p*p + p*q + q*q;
      if peri <=arg then
            cnt:= cnt+ floor(arg/peri);
      end if;
      end if;
    end loop;
  end loop;
  
  dbms_output.put_line(cnt+trunc(arg/3));      
end Q3702;
/

〇〇

〇〇 发表于 2012-2-9 09:58
改变q的上界，q最大不超过1/2周长
create or replace procedure Q3702(arg number default 1e6) is
per ...

其实p也有上界，不能超过1/3周长，否则就不是最短边了,c程序验证通过
但改到pl/sql却出错了
#include <stdio.h>
#include <math.h>
#define M 10000000
int gcd(int a,int b)
{
int temp;if(a<b)/*交换两个数，使大数放在a上*/{temp=a;a=b;b=temp;}
while(b!=0)/*利用辗除法，直到b为0为止*/{temp=a%b;a=b;b=temp;}
return a;
}
int main()
{
int a;
int c;
int s=M/3;
for(int a0=2;a0<sqrt(M/3);a0++)
{a=a0*a0;
for(int c0=a0+1;c0<=sqrt(M/2);c0++)
{
if (gcd(a0,c0)==1)
{
c=c0*c0;
if( a+a0*c0>c && /*a<a0*c0 && a0*c0<c &&*/ a+a0*c0+c<=M)
s+=M/(a+a0*c0+c);
}
if(a%500000==1)
printf("%d,%d\n",a,s);
}
}
printf("%d",s);
return 1;
}

SQL> exec q3702(1E6)
853262

〇〇

〇〇 发表于 2012-2-9 09:58
改变q的上界，q最大不超过1/2周长
create or replace procedure Q3702(arg number default 1e6) is
per ...

条件p+q > q*q/p 可以解不等式，得到p/q>(sqrt(5)-1)/2=0.618
create or replace procedure Q3702(arg number default 1e6) is
peri number;
cnt number :=0;
begin

  for q in 1.. sqrt(arg/2) loop
    for p in 0.618*q .. q-1 loop
      if pkg_prim2.gcd(p,q) = 1 and p+q > q*q/p then
      peri := p*p + p*q + q*q;
      if peri <=arg then
            cnt:= cnt+ floor(arg/peri);
      end if;
      end if;
    end loop;
  end loop;
  
  dbms_output.put_line(cnt+trunc(arg/3));      
end Q3702;
/

SQL> exec q3702(1E6)
861805

PL/SQL 过程已成功完成。

已用时间:  00: 00: 00.43
SQL> exec q3702(1E7)
9712598

PL/SQL 过程已成功完成。

已用时间:  00: 00: 04.86

tree_new_bee

可以细调， 但调不出数量级的变化来。 这样是不可能做出2.5e13的。

sqrt(2.5e13/2) = 3535533.9
用plsql 仅仅简单循环3535533次，也要秒级的时间。

所以要想在合理的时间内作出， 要么有什么数学方法回避循环， 要么就控制好只做一层循环， 把复杂度控制在O(n)

〇〇

tree_new_bee 发表于 2012-2-9 11:22
可以细调， 但调不出数量级的变化来。 这样是不可能做出2.5e13的。

sqrt(2.5e13/2) = 3535533.9

如果q为质数，那么范围内的所有p都是可以的，也可以不计算
SQL> exec q3702(1000)
4:6:9
9:12:16
16:20:25
25:30:36
25:35:49
36:42:49
25:40:64
49:56:64
49:63:81
64:72:81
49:70:100
81:90:100

49:77:121
64:88:121
81:99:121
100:110:121

121:132:144
81:117:169
100:130:169
121:143:169
144:156:169
81:126:196
121:154:196
169:182:196
121:165:225
169:195:225
196:210:225
121:176:256
169:208:256
225:240:256

121:187:289
144:204:289
169:221:289
196:238:289
225:255:289
256:272:289

169:234:324
289:306:324

144:228:361
169:247:361
196:266:361
225:285:361
256:304:361
289:323:361

169:260:400
169:273:441
532