Herhangi Matematiksel Şeyler: Euler Projesi 264. Soru

Üçgen Merkezleri

Aşağıda verilen özelliklere sahip tüm üçgenleri düşünün:

Tüm köşeleri örgü noktaları üzerinde
Çevrel çemberin merkezi O orijin noktasında
Yüksekliklerin kesim noktası H(5,0)

Çevre uzunluğu ≤ 50 olan böyle 9 adet üçgen vardır.
Yarıçapları artan sırada olmak üzere bunlar:

A(-4, 3), B(5, 0), C(4, -3)
A(4, 3), B(5, 0), C(-4, -3)
A(-3, 4), B(5, 0), C(3, -4)

A(3, 4), B(5, 0), C(-3, -4)
A(0, 5), B(5, 0), C(0, -5)
A(1, 8), B(8, -1), C(-4, -7)

A(8, 1), B(1, -8), C(-4, 7)
A(2, 9), B(9, -2), C(-6, -7)
A(9, 2), B(2, -9), C(-6, 7)

Bunların çevre uzunlukları toplamı 4 ondalık basamağa kadar 291,0089'dur.

Çevre uzunluğu ≤ 10⁵ olan böyle üçgenleri bulun. Bunların çevre uzunlukları toplamını 4 ondalık basamağa kadar cevabınız olarak girin.
Cevap: 2816417.1055

Mathematica:
#include<iostream>
#include<cmath>
#include<cassert>
#include<set>
#include<complex>

using namespace std;

typedef long long Int;
typedef long double ld;
typedef complex<int> P;

const Int MAX=100000;
const ld EPS=1e-10;

namespace std {
  bool operator < (const P& a, const P& b) {
    return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
  }
  bool operator > (const P& a, const P& b) {
    return real(a) != real(b) ? real(a) > real(b) : imag(a) > imag(b);
  }

}

Int cross(const P& a, const P& b) {
  return imag(conj(a)*b);
}
Int dot(const P& a, const P& b) {
  return real(conj(a)*b);
}

ld absc(const P &a){
    return sqrt((ld)dot(a, a));
}

struct TRI{
    P a, b, c;
    void sort(){
        P temp;
        if(arg(a)>arg(b)){
            temp=a;
            a=b;
            b=temp;
        }
        if(arg(a)>arg(c)){
            temp=c;
            c=a;
            a=temp;
        }
        if(arg(b)>arg(c)){
            temp=b;
            b=c;
            c=temp;
        }
        return;
    }
    TRI(P a, P b, P c): a(a), b(b), c(c){sort();}
    Int area(){
        return cross(b-a, c-a);
    }
    void print(){
        cout<<a<<b<<c<<endl;
    }
    ld perm(){
        return absc(b-a)+absc(c-b)+absc(a-c);
    }
};

bool operator<(const TRI &l, const TRI &r){
    if(l.a>r.a)return false;
    if(l.a<r.a)return true;
    if(l.b>r.b)return false;
    if(l.b<r.b)return true;
    if(l.c>r.c)return false;
    if(l.c<r.c)return true;
    return false;
}

Int gcd(Int a, Int b){
    if(b==0)return a;
    else return gcd(b, a%b);
}

Int sqrt2(Int n){
    ld val=sqrt((ld)n);
    Int ret=(Int)(val+EPS);
    assert(ret*ret<=n and n<(ret+1)*(ret+1));
    return ret;
}

ld dist(Int x1, Int y1, Int x2, Int y2){
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

Int absol(Int n){
    if(n<0)return -n;
    else return n;
}

set<tri> tri;
ld sum=0;

void hantei(TRI &t){
    if(t.area()<=0)return;
    if(tri.find(t)==tri.end()){
        sum+=t.perm();
        t.print();
        tri.insert(t);
    }
    return;
}

int main(){
    for(Int x3=-MAX/4;x3<=MAX/4;x3++){
        for(Int y3=0;y3<=MAX/4;y3++){
            if(x3==0 and y3==0)continue;
            Int gs=5-x3, gt=-y3;
            if(gs==0 and gt==0)continue;
            Int g=gcd(absol(gs), absol(gt));
            Int s=gs/g, t=gt/g;
            Int d=s*s+t*t;
            Int val=((Int)25-gs*10)*4;
            if(val%d!=0)continue;
            Int k2=g*g*3+val/d;
            if(k2<0)continue;
            Int k=sqrt2(k2);
            if(k*k!=k2)continue;
            Int x1=(gs+k*t)/2, x2=(gs-k*t)/2;
            Int y1=(gt-k*s)/2, y2=(gt+k*s)/2;
            ld per=dist(x1, y1, x2, y2)+dist(x2, y2, x3, y3)+dist(x3, y3, x1, y1);
            if(per-EPS>MAX)continue;
            P p1=P(x1, y1), p2=P(x2, y2), p3=P(x3, y3);
            TRI t1=TRI(p1, p2, p3);
            hantei(t1);
            TRI t3=TRI(conj(p1), conj(p2), conj(p3));
            hantei(t3);

}
    }
    printf("%.5Lf\n", sum);
    return 0;
}

C/C++:
#include <stdio.h>
#include <math.h>

#define N 100000
#define LIMITR (N*5000)

#define PN 500000
int prime_search[PN+1];
int prime_t[PN/5];
int primen=0;

int prime5[PN/5];
int primen5=0;
int prime52[PN/5][2];

double answer = 0;

typedef struct _Point {
  int x;
  int y;
} Point;

Point pp[100000];
int pn;

void make_primetable();
void calc_radius0 ();
void calc_radius(int nth, int min, int r);
void calc_point(int nth, int r);
int comp_pp(Point *a, Point *b);
void search_tri();
void merge_ppoint(int nth);

void decompose_prime5();

int main(int argc, char **argv) {
  int i;

make_primetable();

for (i=1; i<primen; i++) {
    if (LIMITR < prime_t[i])
      break;
    if (prime_t[i] % 4 == 1)
      prime5[primen5++] = prime_t[i];
  }

decompose_prime5();

calc_radius0();

printf ("answer: %.4f  (%.7f)\n", answer, answer);
}

void decompose_prime5() {
  int i, j;
  int d, dsq;

for (i=0; i<primen5; i++) {
    for (j=1; j<prime5[i]; j++) {
      d = prime5[i]-j*j;
      dsq = sqrt(d);
      if (d == dsq*dsq) {
 prime52[i][0] = j;
 prime52[i][1] = dsq;
 break;
      }
    }
  }
}

int ii[20];
int iif[20];

// nth src/trg n element
int pset[20][2][10][2];
int ppoint[20][10000][2];
int ppointn[20];

void calc_radius0 () {
  int i;
  long long int nr;
  int trg, src;
  int nth = 0;
  int tmp;

ii[0] = 0; // 5
  {
    nr = prime5[ii[0]];

pset[0][1][0][0] = prime52[ii[0]][0];
    pset[0][1][0][1] = prime52[ii[0]][1];
    
    for (iif[0]=1;; iif[0]++) {
      ppointn[0] = iif[0]/2+1;
      for (i=0; i<ppointn[0]; i++) {
 ppoint[0][i][0] = pset[0][iif[0]%2][i][0];
 ppoint[0][i][1] = pset[0][iif[0]%2][i][1];
 // printf ("(%d %d)", ppoint[0][i][0],ppoint[0][i][1]);
      }
      //      printf ("\n");

calc_radius(1, 1, nr);

if (iif[0]!=1)
 calc_point(0, nr);

nr *= prime5[ii[0]];
      if (nr > LIMITR)
 break;

src = iif[0]%2;
      trg = 1-src;
      if (iif[0]==1) {
 pset[0][0][0][0] = 0;
 pset[0][0][0][1] = prime5[ii[0]];
 i = 1;
      } else {
 for (i=0; i<(iif[0]+1)/2; i++) {
   pset[0][trg][i][0] *= prime5[ii[0]];
   pset[0][trg][i][1] *= prime5[ii[0]];
 }
      }
 tmp = pset[0][src][iif[0]/2][0]*prime52[ii[0]][0]+pset[0][src][iif[0]/2][1]*prime52[ii[0]][1];
 if (tmp % prime5[ii[0]]) {
   pset[0][trg][i][0] = tmp;
   pset[0][trg][i][1] = abs(pset[0][src][iif[0]/2][0]*prime52[ii[0]][1]-pset[0][src][iif[0]/2][1]*prime52[ii[0]][0]);
 } else {
   pset[0][trg][i][0] = abs(pset[0][src][iif[0]/2][0]*prime52[ii[0]][0]-pset[0][src][iif[0]/2][1]*prime52[ii[0]][1]);
   pset[0][trg][i][1] = pset[0][src][iif[0]/2][0]*prime52[ii[0]][1]+pset[0][src][iif[0]/2][1]*prime52[ii[0]][0];
 }
    }
  }

}

void calc_radius(int nth, int min, int r) {
  int i;
  long long int nr;
  int trg, src;
  int tmp;

long long int limitrn;

limitrn = LIMITR/r;

for (ii[nth]=min; ii[nth]<primen5; ii[nth]++) {

if (prime5[ii[nth]] > limitrn)
      return;

nr = prime5[ii[nth]]*(long long int)r;

pset[nth][1][0][0] = prime52[ii[nth]][0];
    pset[nth][1][0][1] = prime52[ii[nth]][1];
    
    for (iif[nth]=1;; iif[nth]++) {
      merge_ppoint(nth);
      calc_point(nth, nr);
      calc_radius(nth+1, ii[nth]+1, nr);

nr *= prime5[ii[nth]];
      if (nr > LIMITR)
 break;

src = iif[nth]%2;
      trg = 1-src;
      if (iif[nth]==1) {
 pset[nth][0][0][0] = 0;
 pset[nth][0][0][1] = prime5[ii[nth]];
 i=1;
      } else {
 for (i=0; i<(iif[nth]+1)/2; i++) {
   pset[nth][trg][i][0] *= prime5[ii[nth]];
   pset[nth][trg][i][1] *= prime5[ii[nth]];
 }
      }
      tmp = pset[nth][src][iif[nth]/2][0]*prime52[ii[nth]][0]+pset[nth][src][iif[nth]/2][1]*prime52[ii[nth]][1];
      if (tmp % prime5[ii[nth]]) {
 pset[nth][trg][i][0] = tmp;
 pset[nth][trg][i][1] = abs(pset[nth][src][iif[nth]/2][0]*prime52[ii[nth]][1]-pset[nth][src][iif[nth]/2][1]*prime52[ii[nth]][0]);
      } else {
 pset[nth][trg][i][0] = abs(pset[nth][src][iif[nth]/2][0]*prime52[ii[nth]][0]-pset[nth][src][iif[nth]/2][1]*prime52[ii[nth]][1]);
 pset[nth][trg][i][1] = pset[nth][src][iif[nth]/2][0]*prime52[ii[nth]][1]+pset[nth][src][iif[nth]/2][1]*prime52[ii[nth]][0];
      }
    }
  }

}

void calc_point(int nth, int r) {
  int i;
  long long int r2;
  long long int d;
  int x, y;
  int xmax;

pn = 0;

// duplicate mirror images
  for (i=0; i<ppointn[nth]; i++) {
    pp[pn].x = ppoint[nth][i][0];
    pp[pn++].y = ppoint[nth][i][1];
    pp[pn].x = ppoint[nth][i][0];
    pp[pn++].y = -ppoint[nth][i][1];
    pp[pn].x = ppoint[nth][i][1];
    pp[pn++].y = ppoint[nth][i][0];
    pp[pn].x = -ppoint[nth][i][1];
    pp[pn++].y = ppoint[nth][i][0];
    if (ppoint[nth][i][0]) {
      pp[pn].x = -ppoint[nth][i][0];
      pp[pn++].y = ppoint[nth][i][1];
      pp[pn].x = -ppoint[nth][i][0];
      pp[pn++].y = -ppoint[nth][i][1];
      pp[pn].x = ppoint[nth][i][1];
      pp[pn++].y = -ppoint[nth][i][0];
      pp[pn].x = -ppoint[nth][i][1];
      pp[pn++].y = -ppoint[nth][i][0];
    }
  }

qsort (pp, pn, sizeof(Point), comp_pp);
  search_tri();
}

int comp_pp(Point *a, Point *b) {
  if (a->x == b->x)
    return 0;
  else if (a->x < b->x)
    return -1;
  else
    return 1;
}

//search triangle
void search_tri() {
  int i, j, k, ktmp;
  long long int tmp;
  int x3, y3, xk;
  double p2, p;
  for (i=0; i<pn-2; i++) {
    if (pp[i].x >=2)
      break;
    j = i+1;
    k = pn-1;
    while (j<k) {
      x3 = pp[i].x + pp[j].x + pp[k].x;
      if (x3 > 5)
 k--;
      else if (x3 < 5)
 j++;
      else {// (x3 == 5)
 xk = pp[k].x;
 ktmp = k;
 while(xk == pp[ktmp].x) {
   if (ktmp <= j)
     break;
   y3 = pp[i].y + pp[j].y + pp[ktmp].y;
   if (y3 == 0) {
     // get answer
     tmp = (pp[i].x - pp[j].x);
     tmp *= tmp;
     p2 = tmp;
     tmp = (pp[i].y - pp[j].y);
     tmp *= tmp;
     p2 += tmp;
     p = sqrt(p2);
     tmp = (pp[i].x - pp[ktmp].x);
     tmp *= tmp;
     p2 = tmp;
     tmp = (pp[i].y - pp[ktmp].y);
     tmp *= tmp;
     p2 += tmp;
     p += sqrt(p2);
     tmp = (pp[j].x - pp[ktmp].x);
     tmp *= tmp;
     p2 = tmp;
     tmp = (pp[j].y - pp[ktmp].y);
     tmp *= tmp;
     p2 += tmp;
     p += sqrt(p2);
     if (p <= N) {
       printf (">> (%d %d) (%d %d) (%d %d) ", pp[i].x, pp[i].y, pp[j].x, pp[j].y, pp[ktmp].x, pp[ktmp].y);
       answer += p;
       printf ("p %f %.7f\n", p, answer);
     }
   }
   ktmp--;
 }
 j++;
      }
    }
    
  }

}

void make_primetable() {
  int i, j;
  int SQRN = sqrt(PN);

primen = 0;
  
  prime_t[0] = 2;
  primen++;
  prime_t[1] = 3;
  primen++;

for (i=5; i<=SQRN; i+=2) { // sqrt(SQRN)
    if (i%3==0)
      continue;
    if (!prime_search[i/3]) {
      // prime
      prime_t[primen] = i;
      for (j=i; j<PN; j += (i+i)) {
 if (j%3==0)
   continue;
        prime_search[j/3] = primen;
      }
      primen++;
    }
  }

for (; i<=PN; i+=2) {
    if (i%3==0)
      continue;
    if (!prime_search[i/3]) {
      prime_search[i/3] = primen;
      prime_t[primen++] = i;
    }
  }

printf("prime %d %d  PN: %lld\n", primen, prime_t[primen-1], PN);
}

//pset[nth][iif[nth]%2][iif[nth]/2+1][] * ppoint[nth-1][- ppointn[nth-1]][]
void merge_ppoint(int nth) {
  int i;
  int j;
  int tmp1, tmp2, tmp3, tmp4;

ppointn[nth] = 0;
  for (i=0; i<ppointn[nth-1]; i++) {
    for (j=0; j<iif[nth]/2+1; j++) {
      if (pset[nth][iif[nth]%2][j][0]==0) {
 ppoint[nth][ppointn[nth]][0] = ppoint[nth-1][i][0] * pset[nth][iif[nth]%2][j][1];
 ppoint[nth][ppointn[nth]][1] = ppoint[nth-1][i][1] * pset[nth][iif[nth]%2][j][1];
 ppointn[nth]++;
      } else if (ppoint[nth-1][i][0]==0) {
 ppoint[nth][ppointn[nth]][0] = pset[nth][iif[nth]%2][j][0] * ppoint[nth-1][i][1];
 ppoint[nth][ppointn[nth]][1] = pset[nth][iif[nth]%2][j][1] * ppoint[nth-1][i][1];
 ppointn[nth]++;
      } else {
 tmp1 = ppoint[nth-1][i][0] * pset[nth][iif[nth]%2][j][0];
 tmp2 = ppoint[nth-1][i][1] * pset[nth][iif[nth]%2][j][1];
 tmp3 = ppoint[nth-1][i][0] * pset[nth][iif[nth]%2][j][1];
 tmp4 = ppoint[nth-1][i][1] * pset[nth][iif[nth]%2][j][0];
 ppoint[nth][ppointn[nth]][0] = tmp1 + tmp2;
 ppoint[nth][ppointn[nth]][1] = abs(tmp3 - tmp4);
 ppointn[nth]++;
 ppoint[nth][ppointn[nth]][0] = abs(tmp1 - tmp2);
 ppoint[nth][ppointn[nth]][1] = tmp3 + tmp4;
 ppointn[nth]++;
      }
    }
  }
}

Haskell:
import Data.List
import System.Environment
import Control.Monad

isSq n = (n==).(^2).iSqrt $ n
iSqrt = round.sqrt.fromIntegral

f 5 0 = []
f x y
    | not $ isSq e = []
    | mod n1 (2*a) == 0 && mod n2 (2*a) == 0 = [div n1 (2*a), div n2 (2*a)]
    | mod n1 (2*a) == 0 =  [div n1 (2*a)]
    | mod n2 (2*a) == 0 =  [div n2 (2*a)]
    | otherwise = []
    where a = ((x-5)^2+y^2)
          e = a*(3*x^2+10*x+3*y^2-25)
          n1 = -(x-5)*a - y * iSqrt e
          n2 = -(x-5)*a  + y * iSqrt e

triangles p = filter valid.concat $ do
  x <- [- p3 - 5..div p 3 + 5]
  y <- [0..p3- abs x]
  z <- f x y
  guard.isSq $ x^2+y^2-z^2
  let u = iSqrt $ x^2+y^2-z^2
  return $ nub [[x,y,z,u,5-x-z,-y-u],[x,y,z,-u,5-x-z,-y+u],[x,-y,z,u,5-x-z,y-u],[x,-y,z,-u,5-x-z,y+u]]
    where valid x = onCircle x && descend x
          p3 = div p 3

perimeter [x,y,z,u,v,w] = sum.map dist $ [(x,y,z,u),(z,u,v,w),(v,w,x,y)]
    where dist (a,b,c,d) = sqrt.fromIntegral $ (a-c)^2 + (b-d)^2

onCircle [x,y,z,u,v,w] = x^2+y^2 == z^2 + u^2 && x^2+y^2 == v^2 + w^2
descend [x,y,z,u,v,w] = (y,x) < (u,z) && (u,z) < (w,v)

p264 p = sum.filter ((<=p).ceiling).map perimeter.triangles $ p

main = do n <- liftM (read.head) getArgs
          print.p264 $ n

Python:
import math

k = 5
pmax = 10 ** 5

def sqrt(x):
    xx = int(math.sqrt(x))
    while xx ** 2 > x:
        xx -= 1
    while (xx + 1) ** 2 <= x:
        xx += 1
    return xx

def solve(x, y):
    a = 4*((x-k)**2+y**2)
    if a == 0:
        return
    b = y*a
    c = y**4-((x+k)*(3*x-k)+2*y**2)*(x-k)**2
    d = b**2-4*a*c
    if d >= 0:
        dd = sqrt(d)
        if dd ** 2 == d:
            r = x**2+y**2
            def test(x1,y1):
                x2,y2 = k-x-x1,-y-y1
                if x2**2+y2**2 == r and (x1 != x or y1 != y) and (x2 != x or y2 != y) \
                   and (x2 != x1 or y2 != y1):
                    return x,y,x1,y1,x2,y2
                return None
            if (-b+dd) % (2*a) == 0:
                y1 = (-b+dd) // (2*a)
                xx1 = x**2+y**2-y1**2
                if xx1 >= 0:
                    x1 = sqrt(xx1)
                    if x1 ** 2 == xx1:
                        t = test(x1,y1)
                        if t:
                            yield t
                        if dd != 0:
                            t = test(-x1,y1)
                            if t:
                                yield t

def d(a,b):
    return math.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2)
def p(t):
    return d(t[0:2],t[2:4]) + d(t[2:4],t[4:6]) + d(t[4:6],t[0:2])

tr = set()
s = 0

def test(tt):
    global tr
    global s
    ttt = frozenset(tt)
    if ttt not in tr:
        tr.add(ttt)
        s += pp

for xy in range(pmax // 3):
    for x in range(xy + 1):
        for t in solve(x, xy - x):
            pp = p(t)
            if pp <= pmax:
                test(((t[0], t[1]), (t[2], t[3]), (t[4], t[5])))
                test(((t[0], -t[1]), (t[2], -t[3]), (t[4], -t[5])))

print('{0:.4f}'.format(s))

ml:
module S = Set.Make(struct type t = (int * int * int) let compare = compare end)

let next s =
  let (a,b,c) = S.min_elt s in
  let rec f (a',b',c') l s =
    if a'>a then (a,l,s) else
    let s = S.add (a'+2*b'+1,b'+1,c') (S.remove (a',b',c') s) in
    let s = if b'=c'+1 then S.add (a'+2*c'+1,b',c'+1) s else s in 
    f (S.min_elt s) ((a',b',c')::l) s in f (a,b,c) [] s

let full_circle l = let l = List.map (fun (a,b,c) -> (b,c)) l in
  let tl = List.rev (match l with (a,b)::tl -> if b=0 then tl else l |[] -> []) in
  let tl = match tl with (a,b)::l -> if a=b then l else tl |[] -> [] in
  let l3 = (List.map (fun (a,b) -> (-a,-b)) l)@(List.map (fun (a,b) -> (-b,-a)) tl)@
  (List.map (fun (a,b) -> (b,-a)) l)@(List.map (fun (a,b) -> (a,-b)) tl) in
  let l2 = (List.map (fun (a,b) -> (-b,a)) l)@(List.map (fun (a,b) -> (-a,b)) tl)
    @(List.map (fun (a,b) -> (-a,-b)) l)@(List.map (fun (a,b) -> (-b,-a)) tl) in
  let l1 =   l@(List.map (fun (a,b) -> (b,a)) tl)@(List.map (fun (a,b) -> (-b,a)) l)
    @(List.map (fun (a,b) -> (-a,b)) tl)
  in (l1,l2,l3)

let test = let f (x1,y1) (x2,y2) (x3,y3) = x1*x2+y1*y2+5*x3 in fun a b c -> (f a b c, f b c a, f c a b)

let rawrrawr k = let rec f (l1,l2,l3) z = match l1,l2,l3 with
  | (a::k1, b::k2, c::k3) ->
    let (a',b',c') = test a b c in
    let d = min a' (min b' c') in
    f ((if a'=d then k1 else l1),(if b'=d then k2 else l2),(if c'=d then k3 else l3)) 
      (if a'=b' && a'=c' then k a b c z else z)
 |_ -> z in f

let p = let d (a,b) (c,d) = sqrt (float_of_int (c-a) ** 2. +. float_of_int (d-b) ** 2.) in
  fun a b c -> d a b +. d b c +. d c a

let answer lim =
  let rlim = lim*lim/27 and lim = float_of_int lim in
  let rec f s z =
    let (r,l,s) = next s in
    if r>rlim then z else
    if r mod 5 > 0 then f s z else
    f s (rawrrawr (fun a b c d -> let p = p a b c in if p > lim then d else (p+.d)) 
    (full_circle l) z) in
  f (S.singleton (1,1,0)) 
    (2.*.(p (-4,3) (5,0) (4,-3) +. p (-3,4) (5,0) (3,-4)) +. p (0,5) (5,0) (0,-5))

let _ = answer 100000

Herhangi Matematiksel Şeyler

Sayfalar

27 Nisan 2019

Euler Projesi 264. Soru