Euler Projesi 226. Soru

Pelte Kepçesi

Bir pelte eğrisi, s(x) = "x'ten en yakın tam sayıya olan uzaklık" olmak üzere $0\le x\le 1$ için $$y=\sum_{n=0}^{\infty} [\frac{s(2^nx)}{2^n}],$$ eşitliğini sağlayan (x,y) noktalarının kümesidir.

Pelte eğrisinin altındaki alan 1/2 dir (şekilde pembe alan).

Merkezi (1/4,1/2) ve yarıçapı 1/4 olan çember C olsun (şekilde siyah çember).

Pelte eğrisinin altındaki alanın ne kadarı C'nin içinde kalır? Cevabınızı 0,abcdefgh şeklinde 8 ondalık basamağa kadar veriniz.

Cevap: 0,11316017

C++:

#include <stdio.h>
#include <math.h>
 
int main()
{
 int i,p,n,x,y,res,P,depth=26,pow2=1<<depth;
        double dx,diff,sq,y1,y2,fx,area=0.0;
 
 for(p=0;2*p<pow2;p++)  {
  x=p;
  y=0;
  P=pow2;
  for(i=0;i<depth;i++)  {
   res=x&(P-1);
   if(res>P-res)  res=P-res;
                        y+=res;
                        P>>=1;
 
                }
                        fx=(double) y/pow2;
                        dx=(double) p/pow2-0.25;
   sq=sqrt(1.0/16.0-dx*dx);
                        y1=0.5+sq;
                        y2=0.5-sq;
                        if(fx>=y2)  {
                           if(fx<y1)  area+=(double) (fx-y2)/pow2;
                           else       area+=(double) (y1-y2)/pow2;
                        }
 
 }
 printf("area=%.8lf\n",area);
 return 0;
}

Java:

public static void main(String[] args) throws InterruptedException, ExecutionException {
    int k = 26;
    final double step = 1. / (1 << k);

    int threadCount = Runtime.getRuntime().availableProcessors();
    List<Callable<Double>> tasks = new ArrayList<Callable<Double>>();
    for (int i = 0; i < threadCount; i++) {
        final double start = 0.5 * i / threadCount;
        final double end = 0.5 * (i + 1) / threadCount;
        tasks.add(new Callable<Double>()
        {
            public Double call() throws Exception {
                return integrate(start, end, step);
            }
        });
    }
    
    ExecutorService exec = Executors.newFixedThreadPool(threadCount);
    List<Future<Double>> results = exec.invokeAll(tasks);
    exec.shutdown();
    exec.awaitTermination(1, TimeUnit.DAYS);
    
    double sum = 0;
    for (Future<Double> result : results) {
        sum += result.get();
    }

    System.out.println("sum = " + sum);
}

private static double integrate(double start, double end, double step) {
    double sum = 0;
    for (double x = start; x < end; x += step) {
        double y = curve(x, 26);
        double y1 = 0.5 - Math.sqrt(0.5 * x * (1 - 2 * x));
        if (y > y1) {
            sum += (y - y1);
        }
    }
    sum *= step;
    return sum;
}

private static double curve(double x, int k) {
    double sum = 0;
    for (int j = 0; j < k; j++) {
        sum += s((1 << j) * x) / (1 << j);
    }
    return sum;
}

private static double s(double v) {
    return Math.abs(v - Math.round(v));
}

Python:

from math import sqrt
def p226():
    def f(x, n):
 y = 0.0
 for n in xrange(0, n):
     y += abs(x * 2**n - round(x * 2**n))/2**n
 return y
    step = 0.000005
    n = 40
    x = 0.0
    area = 0.0
    while x <= 0.5:
        y_circle = 0.5 - sqrt(0.0625 - (x - 0.25)**2)
        y = f(x, n)
        if y_circle < y:
            area += (y - y_circle) * step
        x += step
    print round(area, 8)
p226()

Ruby:

xm=0; ym=0
x0=0; y0=0
x1=0.5; y1=0.5
sd=0.5
area=0.0; da=0.125
70.times{
  sd/=2
  da/=4
  xm=(x0+x1)/2; ym=(y0+y1)/2+sd
  y=0.5-Math.sqrt(0.0625-(xm-0.25)**2)
  if ym>y
    area+=da    # adding bump
    area+=(0.5-x1-xm)*(y1-ym)/2    # adding trapezoid
    x1=xm; y1=ym
  else
  x0=xm; y0=ym
  end
}
dx=0.25-xm; dy=0.5-ym
area-=dx*dy/2    # triangle from intersection to cicle center
area+=Math.atan2(dy,-dx)/32   # sector of circle
printf("%10.8f\n", area)

Matlab:

clear

tic
nx  = 1e5;
nn  = 100;

Hx  = 0.5;
dx  = Hx/(nx-1);
x   = 0:dx:Hx;
xc  = 0.25;
yc  = 0.5;
r   = 0.25;

y = 0;
for n=0:nn
    y = y + abs((2^n*x)-round((2^n*x)))/(2^n);
end

dycirc  = sqrt(r^2-(x-xc).^2);
ycirc   = yc - dycirc;
dycirc  = y - ycirc;
cut     = dycirc(dycirc>0);

intcirc = sum(cut)*dx;

disp(num2str(intcirc,8))
toc

Fortran:

int main(){
    ldouble x(0.0),integral(0.0),y;
    const ldouble nula(0.0);
    const ldouble dx=0.00000005;

    x=dx/2.0;

    while (x<0.5){
        y=max(nula,blanc(x)-ispodkruga(x));
        integral+=y*dx;
        x+=dx;
        i++;
    }

    cout<<"Rezultat je "<<setprecision(8)<<integral;

    system("PAUSE");
    return 0;
}