|x| + |y| ≤ r denklemini sağlayan (x,y) tamsayı koordinatlara sahip noktalar kümesi S(r) ile gösterilsin. O(0,0) ve C(r/4,r/4) noktaları verilsin. S(r) kümesinde OBC üçgeni geniş açılı olacak şekilde alınabilen B noktalarının sayısı N(r) olsun (yani en geniş açı α, 90°<α<180° eşitsizliğini sağlar). Bu durumda, örneğin N(4)=24 ve N(8)=100 olur.
Şimdi, N(1,000,000,000) değeri kaçtır?
Cevap: 1598174770174689458
Java:
public class Euler210 {
Euler210(){
final long r = 1000000000;
long sum=0;
// whole square
sum += 2*r*(r+1)+1;
// remove band between 0 <= x+y <=r/2
sum -= (2*r+1)*r/4 + r+1;
// remove points on line x=y, not in previous band
sum -= r/2 + r/4;
// add all points strictly interior to circle around r/8,r/8, radius r/8 sqrt2
sum += countCircle( r*r/32 );
// remove points on line x=y, interior to circle
sum -= r/4 - 1;
System.out.println(sum);
}
// count all points strictly interior to circle with squared radius r2.
long countCircle( long r2 ){
long x = 0;
long x2 = 1;
while( x2<r2 ){
x2 += x+x+1;
x++;
}
long total = 0;
long y = 0;
long y2next = 1;
while( true )
{
x--;
x2 -= x+x+1;
// find top point of this column
while( x2 + y2next < r2){
y++;
y2next += y+y+1;
}
if( x<y) break;
total += 8*y+4;
}
total += (2*x+1)*(2*x+1);
return total;
}
} C++:#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#include "string.h"
long min(a, b)
{
return a > b ? b : a;
}
long max(a, b)
{
return a > b ? a : b;
}
long r = 1000000000;
long okpoint(double x, double y)
{
x -= r/8.0;
y -= r/8.0;
double l = sqrt(x*x+y*y);
double m = sqrt(2.0*(r/8.0)*(r/8.0));
return l < m;
}
int main(void)
{
long total = 0;
long o, a, b;
long x, y;
total += r*(2*r+1)/4;
total -= r/4;
total += r*(2*r+1)/2;
total -= r/2;
total += (r/4+1)*(r/4+1)-(r/4+1)-2;
a = 0;
b = r/4;
y = r/4+1;
long sub = 0;
while(a < b)
{
while(!okpoint(a, y))
{
a++;
if(a > b) break;
}
while(!okpoint(b, y))
{
b--;
if(a > b) break;
}
if(a > b) break;
sub += b-a+1;
y++;
}
printf("%li\n", sub);
total += sub*4;
printf("%li\n", total);
return 0;
} Mathematica:(* Compiled into C to make it run at a reasonable speed: *)
inzone2 =
Compile[{{n, _Integer}, {a, _Integer}, {b, _Integer}},
Module[{eye = 0.0, enn = 0.0, p = 0}, enn = n;
Table[eye = i; p = Ceiling[Sqrt[2 enn^2 - eye^2]];
2 p - 1, {i, a + 1, b}]], CompilationTarget -> "C"];
(* Counts the number of points in one quarter of the area produced by \
removing an inscribed square from the circle: *)
inzone4 = Function[{n}, Module[{q, b, x, y}, Sum[
Print[a];
b = Min[a + 10^6, Floor[n Sqrt[2]]];
q = inzone2[n, a, b];
(* Be careful with points on the circumference of the circle.
This deals with these cases separately,
since I don't trust the floating-point computations in inzone2: *)
Map[(x = #[[1]]; y = #[[2]];
If[a < x && x <= b, q[[x - a]] = 2 y - 1];
If[a < y && y <= b, q[[y - a]] = 2 x - 1]) &,
PowersRepresentations[2 n^2, 2, 2]];
Total[q]
(* Split it into chunks of 10^6 to avoid using ridiculous \
swathes of memory: *)
, {a, n, Floor[n Sqrt[2]], 10^6}]]];
(* Gives the number of points in the interior of a circle of radius n \
Sqrt[2]: *)
interior = Function[{n}, 4 inzone4[n] + (2 n + 1)^2 - 4];
(* Gives the number of obtuse-angled triangles: *)
obtuse = Function[{r}, interior[r/8] + 3 r^2/2 - (r/4 - 1)];
(* Actually do the computation: *)
AbsoluteTiming[obtuse[10^9]] Vb:Public Shared Function Problem210() As String
Const r = 1000L * 1000 * 1000
Const d = r \ 4
Dim total = 0L, innertotal = 0L
total += (r + 1) * (r + 1) + r * r '' add all points
total -= (d + 1) * (r + 1) + d * r '' remove points internal to (1,-1)-(-1,1) and (d,0)-(0,-d)
total -= (r - d) '' remove outer colinear non-triangle points
Dim radiussquare = r * r \ 32
Dim radius = Math.Sqrt(radiussquare)
Dim maxwidth = CLng(radius)
For height = 1 To CLng(Math.Floor(radius))
Do While maxwidth * maxwidth + height * height >= radiussquare AndAlso maxwidth > 0
maxwidth -= 1
Loop
innertotal += CInt(Math.Floor(maxwidth)) + 1
Next
innertotal *= 4 '' we calculated only a quarter of the circle
innertotal -= (r \ 4) - 2 '' remove inner colinear points, except those ignored in center, at O and at C
total += innertotal
Return total.ToString()
End Function c#:class Task210:Task
{
const long r = 1000000000;
const long r2 = r / 2;
public override void doTask()
{
//BruteForce();
long sum = 3 * r * r2+N(r)-r/4+((r%8==0)?1:0);
Print(sum.ToString());
}
private long N(long R)
{
double Ox = R / 8.0;
long Oy = (long)Math.Ceiling(R / 8.0);
long count = 0;
long y = Oy;
for (long x = (long)((1.0+Math.Sqrt(2))/8.0*R); x > Ox; --x)
{
while ((8*x-R) * (8*x-R) + (8*y-R) * (8*y-R) < 2* R * R) ++y;
count += (y-Oy);
}
return 4 * count + 1;
}
private void BruteForce()
{
int x0 = 0;
int y0 = 0;
int xC = (int)r / 4;
int yC = (int)r / 4;
ulong count = 0;
for (long xB = -r; xB <= r; ++xB)
{
long yMin = (xB > 0) ? xB - r : -xB - r;
long yMax = (xB > 0) ? -xB + r : xB + r;
for (long yB = yMin; yB <= yMax; ++yB)
{
if (isObtuse(x0 - xC, y0 - yC, xB - xC, yB - yC) ||
isObtuse(x0 - xB, y0 - yB, xC - xB, yC - yB) ||
isObtuse(xB - x0, yB - y0, xC - x0, yC - y0))
{
++count;
Print(xB.ToString() + "," + yB.ToString() + ",true");
}
}
}
Print(DateTime.Now.ToString() + ": " + count.ToString());
}
private bool isObtuse(long x1, long y1, long x2, long y2)
{
long dotProduct = x1 * x2 + y1 * y2;
long lenProduct = (x1*x1+y1*y1)*(x2*x2+y2*y2);
return (dotProduct < 0 && dotProduct*dotProduct < lenProduct);
}
}