Euler Projesi 222. Soru

Küre Paketleme

Yarıçapları 30 mm, 31 mm, ..., 50 mm olan 21 topu tam olarak içine alabilecek iç yarıçapı 50 mm olan en kısa borunun uzunluğu nedir?

Cevabınızı en yakın tamsayıya yuvarlanmış halde mikrometre ($10^{-1}$ m) cinsinden veriniz?

Cevap: 159093

Java:

static final double R = 50.;
static final int n = 21;
static final double sizes[] = new double[n];

public static void main(String[] args) {
    for (int i = 0; i < sizes.length; i++) {
        sizes[i] = (30. + i) / R;
    }

    Permutation p = new Permutation(n);
    p.shuffle();
    double cur_best = solution_cost(p);
    double best_cost = cur_best;

    int iters = 0, restarts = 0;
    while (true) {
        double min = cur_best;
        int m_i = -1, m_j = -1;
        for (int i = 0; i < n - 1; i++) {
            for (int j = i + 1; j < n; j++) {
                p.swap(i, j);
                double s = solution_cost(p);
                if (s < min) {
                    min = s;
                    m_i = i;
                    m_j = j;
                }
                p.swap(i, j);
            }
        }

        if (min < cur_best) {
            p.swap(m_i, m_j);
            cur_best = min;
        } else {
            best_cost = Math.min(cur_best, best_cost);
            restarts++;
            p.shuffle();
            cur_best = solution_cost(p);
        }

        iters++;
        if (iters % 100000 == 0) {
            System.out.println("best_cost = " + best_cost * R * 1000 + ", iters = " + iters + ", restarts = " + restarts);
        }
    }
}

private static double solution_cost(Permutation p) {
    double s = sizes[p.get(0)] + sizes[p.get(n - 1)];
    for (int k = 1; k < n; k++) {
        s += 2 * Math.sqrt(-1 + sizes[p.get(k - 1)] + sizes[p.get(k)]);
    }
    return s;
}

C++:

#include <cstdio>
#include <cmath>
#include <vector>


double preCost = 1680*1000;

double R = 50*1000;

double cost(double r1, double r2)
{
    double r = r1 + r2;
    //return R*std::sqrt(2/R*r - 1) - r;
    return std::sqrt(r*r - (2*R - r)*(2*R - r)) - r;
}

double vecCost(std::vector<unsigned> const &v)
{
    double ret = preCost;
    for(unsigned i = 0; i < v.size()-1; ++i)
        ret += cost((v[i]+30)*1000, (v[i+1]+30)*1000);

    return ret;  //unsigned(ret + 0.5);
}

int main()
{
    std::vector<unsigned> v;
    //unsigned order[21] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
    //unsigned order[21] = {0, 20, 1, 19, 2, 18, 3, 17, 4, 16, 5, 15, 6, 14, 7, 13, 8, 12, 9, 11, 10};
    unsigned order[21] = {20, 18, 16, 14, 12, 10, 8, 6, 4, 2, 0, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19};
    for(unsigned r = 0; r < 21; ++r)
        v.push_back(order[r]);
    std::printf("%lf\n", vecCost(v));
    //std::printf("%lf\n", cost(30000, 30000));
    /*unsigned sum = 0;
    for(unsigned r = 30; r <= 50; ++r)
        sum += r*2;
    std::printf("%u\n", sum);*/
}

Delphi:

const
  MinR=30;
  MaxR=50;

procedure TForm1.Button1Click(Sender: TObject);
var
  LR, R, S: integer;
  L: extended;
begin
  TheTimer.Start;

  L:= MaxR;
  R:= MaxR;
  S:= -2;
  repeat
    LR:= R;
    if R=MinR then
      begin
        S:= 2;
        R:= MinR+1;
      end
    else
      R:= R+S;
    L:= L+Sqrt((LR+R)*(LR+R)-(100-LR-R)*(100-LR-R))
  until R=MaxR-1;
  L:= L+MaxR-1;

  TheTimer.Stop;
  Caption:= 'Time= '+FloatToStrF(TheTimer.TimeElapsed,ffGeneral,5,0)+' sec';
  Edit1.Text:= IntToStr(Round(L*1000));
end;

Python:

private void Run() {
 double[] r = new double[] {30};
 double min=0;
 for (double i = 31; i <= 50; i++) {
  r = Calc2(r, i, out min);
 }
 Console.WriteLine("Solution: {0}", min);
}

private double Calc(double[] r) {
 double h=0;
 double rp = 50.0;
 double rads=0, b=0;
 for (int i = 1; i < r.Length; i++) {
  rads = r[i]+r[i-1];
  b = 2*rp-rads;
  h += rads*Math.Sin(Math.Acos(b/rads));
 }
 h+=r[0]+r[r.Length-1];   
 return (double)Math.Round(h*1000);
}

private double[] Calc2(double[] r, double ball, out double best) {
 List<double> original = new List<double>(r);
 double min=double.MaxValue, h=0;
 double[] ret=r;
 for (int i = 0; i <= original.Count; i++) {
  List<double> tmp = new List<double>(original);
  tmp.Insert(i, ball);

  double[] res = tmp.ToArray<double>();
  h = Calc(res);
  if (h<min) {
   min = h;
   ret = res;
  }
 }

 best = min;
 return ret;   
}

Perl:

use strict;
use warnings 'all';
  
my @r = (50000,48000,46000,44000,42000,40000,38000,36000,34000,32000,30000,31000,33000,35000,37000,39000,41000,43000,45000,47000,49000);
  
my $s = $r[0]+$r[$#r];
my @s = ();
my $r1 = shift @r;
my $r2 = shift @r;
my $differenz = 100000-$r1-$r2;
while () {
   push @s, (($r1+$r2)**2-$differenz**2);
   $r1 = $r2;
   last if $#r == -1;
   $r2 = shift @r;
   $differenz = 100000-$r1-$r2;
}

foreach (@s) {
   $s += sqrt($_);
}

print "\n$s";

Matlab:

y = 30;
y = input('what is the radius of the first ball\n');
r2 = y;
control = y*2;


for i = 31:1:50
    r = input('what is the diameter of the  next ball?\n');
    tx = (50-r)+(50-(r2));
    th = (r2)+r;
    dy = sqrt(th^2-tx^2);
    y = y+dy;
    control = control + 2*r;
    r2 = r;
end


y = (y+r)*10^3
round(y)