Euler Projesi 283. Soru

Alan/çevre oranının integral olduğu tamsayı kenarlı üçgenler

Kenarları 6, 8 ve 10 olan üçgeni düşünün. Hem çevrenin hem de alanın 24'e eşit olduğu görülebilir. Yani alan/çevre oranı 1'e eşittir.

Ayrıca 13, 14 ve 15 kenarlı üçgeni de düşünün. Alan 84'e eşitken çevre 42'ye eşittir. Yani bu üçgen için alan/çevre oranı 2'ye eşittir.

Alan/çevre oranlarının 1000'i geçmeyen pozitif tam sayılara eşit olduğu tüm tamsayı kenarlı üçgenlerin çevrelerinin toplamını bulun.

Cevap: 28038042525570324

C/C++: #include <iostream> #include <vector> #include <cmath> #include <fstream> using namespace std; #define root3 1.7320508075688772935274463415059 vector<__int64> prosti; void init(){ ifstream fajl("prosti.txt"); __int64 pom; while (1){ fajl>>pom; if (pom>8000000) break; prosti.push_back(pom); } fajl.close(); } void factors(__int64 broj, vector<int> &faktori, int prethodni){ bool ima_delilaca=0; if (broj>1){ for (int i=prethodni;prosti.at(i)<=sqrt(broj);i++) if (broj/prosti.at(i)*prosti.at(i)==broj) { faktori.push_back(prosti.at(i)); factors(broj/prosti.at(i), faktori, i); ima_delilaca=1; break; } if (!ima_delilaca) faktori.push_back(broj); } } vector<__int64> delioci(__int64 broj){ vector<int> faktori; factors(broj, faktori, 0); vector<__int64> delioci, aux, pom; delioci.push_back(1); aux.push_back(1); for (int i=0;i<faktori.size();){ if (aux.size()==1||aux.at(1)==faktori.at(i)) { aux.push_back(aux.at(aux.size()-1)*faktori.at(i)); i++; } else { for (int j=0;j<delioci.size();j++) for (int k=0;k<aux.size();k++) pom.push_back(delioci.at(j)*aux.at(k)); delioci=pom; pom.clear(); aux.clear(); aux.push_back(1); } } for (int j=0;j<delioci.size();j++) for (int k=0;k<aux.size();k++) pom.push_back(delioci.at(j)*aux.at(k)); delioci=pom; sort(delioci.begin(), delioci.end()); return delioci; } unsigned int gcd(unsigned int u, unsigned int v){ int shift; if (u == 0 || v == 0) return u | v; for (shift = 0; ((u | v) & 1) == 0; ++shift) { u >>= 1; v >>= 1; } while ((u & 1) == 0) u >>= 1; do { while ((v & 1) == 0) v >>= 1; if (u <= v) { v -= u; } else { int diff = u - v; u = v; v = diff; } v >>= 1; } while (v != 0); return u << shift; } int main() { init(); __int64 result = 0; for (__int64 m=1;m<=1000;m++) { cout<<m<<endl; for (__int64 u=1;u<=2*m;u++) { if (2*m/u*u==2*m) { for (__int64 v=1;v<=root3*u;v++) { if (gcd(u,v)==1) { __int64 pom = (__int64)4*m*m*(u*u+v*v); vector<__int64> del = delioci(pom); for (int i=0;del.at(i)<=sqrt(pom);i++) { __int64 d1 = del.at(i); __int64 d2 = pom/d1, q1 = d1+2*m*u, q2 = d2+2*m*u; if (q1/v*v==q1&&q2/v*v==q2) { __int64 a = q1/v + 2*m*v/u, b = q2/v + 2*m*v/u, c = (d1 + d2 + 4*m*u)/v; if (a<=b&&b<=c) result+=a+b+c; //cout<<a<<" "<<b<<" "<<c<<endl; } } } } } } } cout<<result<<endl; return 0; } Python: #!/usr/bin/env python # -*- coding: utf-8 -*- from math import sqrt import zm import zdecs import psyco psyco.full() def process(m): s, process.count, sqrt3 = 0, 0, sqrt(3) for u in xrange(1, 2*m+1): if (2*m)%u: continue u2, mu2 = u**2, 2*m*u for v in xrange(1, int(sqrt3 * u)+1): if zm.gcd(u, v) == 1: v2, mv2 = v**2, 2*m*v t2 = 4 * m**2 * (u2+v2) t = int(sqrt(t2)) mn = (2*m)*(v2-u2)//u # get separate factorization factorsm = dict([(p, 2*e) for (p, e) in zm.factorize(2 * m, 4)]) factorsuv = dict(zm.factorize(u2+v2)) for p in factorsuv: factorsm[p] = factorsm.setdefault(p, 0) + factorsuv[p] for d1 in zm.divisors(t2, factorsm.items()): if not (mn <= d1 <= t): continue d2 = t2 // d1 if (d1+mu2)%v == 0 and (d2 + mu2)%v == 0: p = 2*(d1 + d2 + 2*mu2)//v + (2*mv2)//u process.count += 1 s += p print('found {0} results for m={1}'.format(process.count, m)) return s process.count = 0 @zdecs.timing def p283(): s, n = 0, 1000 for m in xrange(1, n+1): s += process(m) print 'Solution: {0}'.format(s) p283() ''' Solution: 28038042525570324 p283 took 8 m 59.136420 s ''' Mathematica: Timing[sm = 0; For[n = 1, n <= 1000, ++n, For[c = 1, c <= 4 Sqrt[3] n, c++, Q = 16 n^2 (c^2 + 16 n^2); ar = Divisors[Q]; For[i = 1, i <= Length[ar], ++i, W1 = ar[[i]] + 16 n^2; W2 = Q/ar[[i]] + 16 n^2; If[Mod[W1, c] != 0, Continue[]]; If[Mod[W2, c] != 0, Continue[]]; a = W1/c; b = W2/c; If[a < b, Continue[]]; If[b < c, Continue[]]; If[Mod[a + b, 2] != 0, Continue[]]; If[Mod[a + c, 2] != 0, Continue[]]; sm += a + b + c]]]; sm] Delphi: procedure TForm1.Button3Click(Sender: TObject); const target=1000; var x,k,f:integer;y,sm,k2:integer;z,bm:int64;big,big2,h:int64;maxfactor,minfactor:int64;count,maxf:integer; divisors:array[0..10000] of integer; top,curtop,multi,multi2,multi3,i,j,d:integer; fp:int64; perimeter:int64; area:extended; sq12:extended; smpf:array of integer; sievemax,siv:integer; u,v,uf,vf,aux:integer; s:string; procedure add(f,multi:integer); var i,j:integer;h:int64; begin curtop:=top; for j:=0 to top do begin h:=divisors[j]; for i:=1 to multi do begin h:=h*f; if h>maxfactor then break; inc(curtop); divisors[curtop]:=h; end; end; top:=curtop; end; procedure process(d,m:integer); begin if (m<0) or (d>maxfactor)then exit; if d>=minfactor then begin sm:=d+k2; y:=sm div x; if sm=x*y then begin bm:=trunc(big/d)+k2; z:=trunc(bm/x); if bm=x*z then begin inc(count);inc(perimeter,x+y+z);area:=area+k*(x+y+z) end; end; end; process(d shl 1,m-1); end; begin starttimer; maxfactor:=0; perimeter:=0; area:=0; count:=0; divisors[0]:=1; sievemax:=8*target*target; siv:=4*target; setlength(smpf,sievemax+1); for i:=1 to sievemax do smpf[i]:=2*i+1; for i:=1 to siv do begin j:=i shl 1 ; if smpf[i]=j+1 then begin f:=j+1; j:=j*(i+1); while j<=sievemax do begin if smpf[j]>f then smpf[j]:=f; inc(j,f); end end; end; for k:=1 to target do begin k2:=(k*k) shl 2; for x:=1 to trunc(sqrt(12)*k) do begin multi2:=2; u:=k;v:=k2+x*x; big:=int64(k2)*v; maxfactor:=trunc(sqrt(big+0.0)); minfactor:=x*x-k2; top:=0; while u and 1=0 do begin inc(multi2,2);u:=u shr 1;end; while v and 1=0 do begin inc(multi2);v:=v shr 1;end; uf:=smpf[u shr 1];vf:=smpf[v shr 1]; while (u>1) and (v>1) do begin if uf=vf then begin multi:=3; u:=u div uf; aux:=u div uf; while u-aux*uf=0 do begin inc(multi,2);u:=aux;aux:=u div uf; end; v:=v div uf; aux:=v div uf; while v-aux*uf=0 do begin inc(multi);v:=aux;aux:=v div uf;end; add(uf,multi); uf:=smpf[u shr 1];vf:=smpf[v shr 1]; end else begin if (uf<vf) then begin multi:=2; u:=u div uf; aux:=u div uf; while u-aux*uf=0 do begin inc(multi,2);u:=aux;aux:=u div uf; end; add(uf,multi); uf:=smpf[u shr 1]; end else begin multi:=1; v:=v div vf; aux:=v div vf; while v-aux*vf=0 do begin inc(multi);v:=aux;aux:=v div vf;end; add(vf,multi); vf:=smpf[v shr 1]; end; end; end; while u>1 do begin multi:=2; u:=u div uf; aux:=u div uf; while u-aux*uf=0 do begin inc(multi,2);u:=aux;aux:=u div uf; end; add(uf,multi); uf:=smpf[u shr 1]; end; while v>1 do begin multi:=1; v:=v div vf; aux:=v div vf; while v-aux*vf =0 do begin inc(multi);v:=aux;aux:=v div vf; end; if vf<=maxfactor then add(vf,multi) else break; vf:=smpf[v shr 1]; end; for i:=0 to top do begin d:=divisors[i]; if odd(x) then process(d,multi2) else process(d shl 1,multi2-1); end; end; end; stoptimer; if area<intpower(2,62) then s:=inttostr(trunc(2*area)) else s:=floattostrf(2*area,ffgeneral,18,18); memo1.lines.add(inttostr(target)+' '+inttostr(count)+' '+inttostr(perimeter*2)+' '+s+' '+elapsedtime); end;