欧拉计划747:三角披萨

〇〇 · 发表于 2021-2-17 22:42

你点页面底下的电脑版就能看见对齐的格式了

〇〇 · 发表于 2021-2-18 08:55

用了公式

In[1]:= Solve[(1+n+c+d)^2-4*n(c+1)*(d+1)==0,d]
2 2 2 2
Out[1]= {{d -> -1 - c + n + 2 c n - 2 Sqrt[-(c n) - c n + c n + c n ]},
2 2 2 2
> {d -> -1 - c + n + 2 c n + 2 Sqrt[-(c n) - c n + c n + c n ]}}
>>> def f1(n):
... s=0
... for c in range(1,n+1):
... v=0
... for d in range(1,n+1):
... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
... if x<0: break #print(c,x,v);break
... else:
... if x==0:s+=3;v+=1
... else:s+=6;v+=1
... return s
...
>>> def f2(n):
... s=0
... for c in range(1,n-1):
... #for d in range(1,n+1):
... w=(-(c* n) - c*c* n + c* n*n + c*c*n* n )
... if w>=0:
... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
... #print(c,x,d)
... if x<0: break
... else:
... if x==0:s+=6*d-3 #;print('****')
... else:s+=6*d
... return s
import time
def test(m):
t=time.time()
s=0
for n in range(3,1001):s+=f1(n)
print('f1',s,time.time()-t)
t=time.time()
s=0
for n in range(3,1001):s+=f2(n)
print('f2',s,time.time()-t)
>>> test(1001)
f1 3008595 1.7994389533996582
f2 3008595 2.2176740169525146

复制代码

〇〇 · 发表于 2021-2-18 09:10

本帖最后由〇〇于 2021-2-18 09:23 编辑

根据c==1时求出的最大d退出c循环，快了不少。这次公式f2b终于比列举f1b快了

>>> def f1b(n):
... s=0
... up=0
... for c in range(1,n+1):
... if c==1:w=(-(c* n) - c*c* n + c* n*n + c*c*n* n );up=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
... if c>up:break
... v=0
... for d in range(1,n+1):
... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
... if x<0: break #print(c,x,v);break
... else:
... if x==0:s+=3;v+=1
... else:s+=6;v+=1
... return s
...
>>> def f2b(n):
... s=0
... up=0
... for c in range(1,n-1):
... #for d in range(1,n+1):
... w=(-(c* n) - c*c* n + c* n*n + c*c*n* n )
... if w>=0:
... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
... if c==1:up=d
... if c>up:break
... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
... #print(c,x,d)
... if x<0: break
... else:
... if x==0:s+=6*d-3 #;print('****')
... else:s+=6*d
... return s
...
>>> import time
>>> def test(m):
... t=time.time()
... s=0
... for n in range(3,m):s+=f1(n)
... print('f1',s,time.time()-t)
... t=time.time()
... s=0
... for n in range(3,m):s+=f2(n)
... print('f2',s,time.time()-t)
... t=time.time()
... s=0
... for n in range(3,m):s+=f2b(n)
... print('f2b',s,time.time()-t)
... t=time.time()
... s=0
... for n in range(3,m):s+=f1b(n)
... print('f1b',s,time.time()-t)
...
>>> test(1001)
f1 3008595 1.7506561279296875
f2 3008595 2.199603796005249
f2b 3008595 0.3860177993774414
f1b 3008595 0.9080519676208496
>>> test(2001)
f1 14100711 7.668433427810669
f2 14100711 9.103666543960571
f2b 14100711 1.5756027698516846
f1b 14100711 4.127683877944946

复制代码

〇〇 · 发表于 2021-2-18 09:32

在pypy中，反而是f1b快

>>>> test(1001)
('f2b', 3008595, 0.17499995231628418)
('f1b', 3008595, 0.032000064849853516)
>>>> test(2001)
('f2b', 14100711, 0.24900007247924805)
('f1b', 14100711, 0.04699993133544922)

复制代码

〇〇 · 发表于 2021-2-18 13:17

加速版

>>>> def f4c(n):
.... s=0
.... up=0
.... a=[0 for i in range(100)]
.... for c in range(1,n-1):
.... #for d in range(1,n+1):
.... w=(-(c* n) - c*c* n + c* n*n + c*c*n* n )
.... if w>=0:
.... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
.... if c==1:up=d
.... if c>up:break
.... if c>d:
.... a[c]=d
.... #print('&&&&',a[slice(10)])
.... for i in range(1,c):
.... if a[i+1]==0 or a[i]==i: break
.... s+=6*i*(a[i]-a[i+1]) #;print(i,s)
.... break
.... #if d0>d:s+=6*(d0-d);break
.... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
.... #print(c,x,d)
.... if x<0: break
.... else:
.... a[c]=d
.... if x==0:s+=6*d-3 #;print('****',6*d-3)
.... else:s+=6*d
.... #if d0>d:s+=6*(st)*(d0-d) #;print(6*d)
.... return s
....
>>>> import time
>>>> def test(m):
.... t=time.time()
.... s=0
.... for n in range(3,m):s+=f2b(n)
.... print('f2b',s,time.time()-t)
.... t=time.time()
.... s=0
.... for n in range(3,m):s+=f1b(n)
.... print('f1b',s,time.time()-t)
.... t=time.time()
.... s=0
.... for n in range(3,m):s+=f4c(n)
.... print('f4c',s,time.time()-t)
....

复制代码

〇〇 · 发表于 2021-2-18 15:16

改正了一个bug

>>>> def f4c(n):
.... s=0
.... up=0
.... x0=0
.... a=[0]# for i in range(100)]
.... for c in range(1,n-1):
.... #for d in range(1,n+1):
.... w=(-(c* n) - c*c* n + c* n*n + c*c*n* n )
.... if w>=0:
.... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
.... if c==1:up=d
.... if c>up:break
.... if c>d:
.... if x0>0:s-=3*x0
.... a.append(d) #a[c]=d
.... #print('&&&&',a[slice(10)])
.... for i in range(1,c):
.... if a[i]==i: break
.... s+=6*i*(a[i]-a[i+1]) #;print(i,s)
.... break
.... #if d0>d:s+=6*(d0-d);break
.... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
.... #print(c,x,d)
.... if x<0: break
.... else:
.... a.append(d) #a[c]=d
.... if x==0:
.... s+=6*d-3
.... if c!=d:x0+=1 #;print('****',6*d-3)
.... else:s+=6*d
.... #if d0>d:s+=6*(st)*(d0-d) #;print(6*d)
.... return s
....

复制代码

〇〇 · 发表于 2021-2-18 21:10

越算越慢
testget2(2000001,3000001)
1882.57299995
804155548
testget2(3000001,4000001)
2632.04400015
774864131
testget2(4000001,5000001)
2993.84100008
463780742

〇〇 · 发表于 2021-2-18 21:28

本帖最后由〇〇于 2021-2-19 21:20 编辑

不能一个n算一次，要把它们合在一起。比如+6的几个，+12的几个
(c,d)==(1,1)的个数*12
+(1,2)的个数*24
+(1,3)的个数*36
-(x==0)的个数*3
-(c==d)的个数*3
-n最后一个d是1的个数*(12*1-6)
-n最后一个d是2的个数*(12*2-6)
-n最后一个c平方*6*个数

>>>> def f3(n):
.... s=0
.... for c in range(1,n-1):
.... w=(c* n)*(-1 - c + n + c* n )
.... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
.... if c>d:s-=d*d*6;break
.... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
.... print(n,c,x,d)
.... s+=12*d
.... if x==0:
.... s-=3
.... if c!=d:s-=3
.... if c==d:s-=12*d-6
.... return s
....
>>>> s=0
>>>> for i in range(9,26):s+=f3(i)
(9, 1, 0, 1)
(10, 1, 9, 1)
(11, 1, 20, 1)
(12, 1, 33, 1)
(13, 1, 48, 1)
(14, 1, 65, 1)
(15, 1, 1, 2)
(16, 1, 16, 2)
(17, 1, 33, 2)
(18, 1, 52, 2)
(19, 1, 73, 2)
(20, 1, 96, 2)
(21, 1, 4, 3)
(22, 1, 25, 3)
(23, 1, 48, 3)
(24, 1, 73, 3)
(25, 1, 100, 3)
(25, 2, 0, 2)
>>>> print(s)
294
>>>> m=12*1*6 +12*2*7 +12*3*5 -3 -3 -6*(12*1-6) -(12*2-6) -0*0*6*5 -1*1*6*(6+5)
>>>> print(m)
294

复制代码

〇〇 · 发表于 2021-2-19 12:23

f4c简化为不需要数组

>>>> def f3b(n):
.... s=0
.... up=0
.... for c in range(1,n-1):
.... w=(c* n)*(-1 - c + n + c* n )
.... d=int( -1 - c + n + 2* c* n - 2* pow(w,0.5))
.... if c>d:s-=d*d*6;break
.... x=(-1-c-d-n)**2-4*n*(1+d)*(1+c)
.... #print(c,x,d)
.... s+=12*d
.... if x==0:
.... s-=3
.... if c!=d:s-=3
.... if c==d:s-=12*d-6
.... return s

复制代码

〇〇 · 发表于 2021-2-19 13:08

newkid 发表于 2021-2-17 22:24
方程我看得懂，就是答案看不懂，因为在手机上分数显示不出来。你把那个不等式展开看看有没有希望简化得到公 ...

In[17]:= Reduce[{
c==1,
2> -1 - c + n + 2* c* n - 2* Sqrt[(c* n)*(-1 - c + n + c* n )]>=1
},n]
Out[17]= n == 1 || 9 <= n < 4 (2 + Sqrt[3])
In[18]:= Reduce[{
c==1,
3> -1 - c + n + 2* c* n - 2* Sqrt[(c* n)*(-1 - c + n + c* n )]>=2
},n]
Out[18]= 4 (2 + Sqrt[3]) <= n < 11 + 4 Sqrt[6]
In[19]:=
In[19]:= Reduce[{
c==1,
4> -1 - c + n + 2* c* n - 2* Sqrt[(c* n)*(-1 - c + n + c* n )]>=3
},n]
Out[19]= 11 + 4 Sqrt[6] <= n < 2 (7 + 2 Sqrt[10])

复制代码

欧拉计划747:三角披萨

浏览过的版块