Посчитаем цыплят по осени или немного матпрограмирования

| рубрика «Программы» | автор st
Метки:

Друзья подкинули задачку из кружка математики для школьников. Инженерный калькулятор для неё не годится. А программируемому может не хватить процессорных силёнок.

Три фермера продавали кур на местном рынке. У одного было 10 кур, у второго — 16, у третьего — 26. Чтобы не конкурировать между собой, они договорились продавать кур по одной цене. К обеду они решили, что продажи идут не так уж хорошо, поэтому они все одинаково понизили цену. К концу дня они продали всех кур. Оказалось, что каждый из фермеров за этот день выручил 35 долларов. Какова была цена за курицу до обеда и после обеда?

Обозначим цены как X и Y, а количества проданных кур, соответственно, до и после обеда i-м фермером как Ai и Bi. Получается система линейных уравнений:

A1*X + (10 - A1)*Y = 3500
A2*X + (16 - A2)*Y = 3500
A3*X + (26 - A3)*Y = 3500

Дополнительные условия:

  • точность до центов;
  • 0 > X > 3500;
  • 0 > Y > 3500;
  • X > Y;
  • куры продаются только целиком.

Это позволяет немного упростить задачу и решить её в целых числах.

Небольшой объем позволяет решить задачу полностью путем перебора. Вот соответствующая программа на FreePascal.

program Chickens;

var
  a1, a2, a3: integer;
  x, y: integer;
begin
  for x := 1 to 3499 do
    for y := 1 to x - 1 do
      for a1 := 0 to 10 do
        for a2 := 0 to 16 do
          for a3 := 0 to 26 do
            if (a1 * x + (10 - a1) * y = 3500) and
               (a2 * x + (16 - a2) * y = 3500) and
               (a3 * x + (26 - a3) * y = 3500)
            then
              writeln('x = ', x, ', y = ', y,
                ', a1 = ', a1,
                ', b1 = ', 10 - a1,
                ', a2 = ', a2,
                ', b2 = ', 16 - a2,
                ', a3 = ', a3,
                ', b3 = ', 26 - a3
                );
  writeln('Finished');
  readln;
end.

Программа считает достаточно долго, несколько минут на двуядерном Intel 2,2 ГГц (что хорошо иллюстрирует полную неразрешимость задач из обычной логистики даже средствами современных суперкомпьютеров). Ответ:

x = 375, y = 125, a1 = 9, b1 = 1, a2 = 6, b2 = 10, a3 = 1, b3 = 25

Вопросы любителям вычислений:

  • Сколько считает программа на вашем калькуляторе?
  • Какие методы линейного программирования позволят сократить счет, чтобы получить хотя бы одно решение за минимальное время?

Дополнение: тесты производительности

Тот же текст программы можно без изменений запустить в Delphi. Для Си программу надо слегка переделать.

#include <stdio.h>
#include <stdlib.h>

int main()
{
  int a1, a2, a3;
  int x, y;
  for (x = 1; x < 3500; x++)
    for (y = 1; y < x; y++)
      for (a1 = 0; a1 <= 10; a1++)
        for (a2 = 0; a2 <= 16; a2++)
          for (a3 = 0; a3 <= 26; a3++)
            if ((a1 * x + (10 - a1) * y == 3500) &&
                (a2 * x + (16 - a2) * y == 3500) &&
                (a3 * x + (26 - a3) * y == 3500))
              printf("x = %d, y = %d, a1 = %d, b1 = %d, a2 = %d, b2 = %d, a3 = %d, b3 = %d\n",
                     x, y, a1, 10 - a1, a2, 16 - a2, a3, 26 - a3);
  printf("Finished\n");
  return 0;
}

Результаты на однопроцессорной виртуальной машине под Windows XP:

  • FreePascal 2.6 : 75 сек
  • Delphi 7 : 55 сек
  • GNU C 4.6 : 38 сек

Для C# программу также надо немного переделать

using System;

namespace chicken_csharp
{
  class MainClass
  {
    public static void Main (string[] args)
    {
      DateTime d1 = DateTime.Now;
      int a1, a2, a3;
      int x, y;
      for (x = 1; x < 3500; x++)
        for (y = 1; y < x; y++)
          for (a1 = 0; a1 <= 10; a1++)
            for (a2 = 0; a2 <= 16; a2++)
              for (a3 = 0; a3 <= 26; a3++)
                if ((a1 * x + (10 - a1) * y == 3500) &&
                   (a2 * x + (16 - a2) * y == 3500) &&
                   (a3 * x + (26 - a3) * y == 3500))
                  Console.WriteLine ("x = {0}, y = {1}, a1 = {2}, b1 = {3}, a2 = {4}, b2 = {5}, a3 = {6}, b3 = {7}",
                                   x, y, a1, 10 - a1, a2, 16 - a2, a3, 26 - a3);
      Console.WriteLine ("Finished");
      Console.WriteLine ("Elapsed sec {0}", DateTime.Now.Subtract (d1).TotalSeconds);
    }
  }
}

Результаты на 4-ядерном ПК под Ubuntu 12.04 x64 (используется только 1 ядро, судя по монитору):

  • FreePascal 2.6 : 58 сек
  • C# Mono/.NET 4 : 53 сек
  • GNU C 4.6 : 29 сек

Во всех случаях включена оптимизация, отключена отладочная информация.

Си даже в его GNU-той реализации показывают впечатляющее качество генерируемого кода.

Есть желающие сравнить с другими интерпретаторами и байткод-компиляторами типа Явы?

Оптимизация 1

В первоначальной программе можно сделать простую оптимизацию: в следующий уровень цикла входим только если условие "все куры i-го продавца проданы за 3500 центов".

program Chickens;

var
  a1, a2, a3: integer;
  x, y: integer;
begin
  for x := 1 to 3499 do
    for y := 1 to x - 1 do
      for a1 := 0 to 10 do
        if (a1 * x + (10 - a1) * y = 3500) then
          for a2 := 0 to 16 do
            if (a2 * x + (16 - a2) * y = 3500) then
              for a3 := 0 to 26 do
                if (a3 * x + (26 - a3) * y = 3500) then
                  writeln('x = ', x, ', y = ', y,
                    ', a1 = ', a1,
                    ', b1 = ', 10 - a1,
                    ', a2 = ', a2,
                    ', b2 = ', 16 - a2,
                    ', a3 = ', a3,
                    ', b3 = ', 26 - a3);
  writeln('Finished');
  readln;
end.

Время счета снижается с 58 до примерно 0,5 секунды, то есть на 2 порядка.

Оптимизация 2

Цену Y можно вычислять, как функцию от переменных (a1, x). Для наглядности я закомментировал "лишний" цикл, а значение цену Y=-1 используется как любое (весь товар был продано i-м торговцем ещё до обеда). Необходимо в этих случаях добавить проверки деления на ноль.

program Chickens;

var
  a1, a2, a3: integer;
  x, y: integer;
begin
  for x := 1 to 3499 do
    //    for y := 1 to x - 1 do
    for a1 := 0 to 10 do
      if ((a1 * x) <= 3500) and ((a1 = 10) or ((3500 - (a1 * x)) mod (10 - a1) = 0)) then
      begin
        if a1 = 10 then
          y := -1 // цена может быть любой
        else
          y := (3500 - (a1 * x)) div (10 - a1);
        if (y = 0) or (x <= y) then
          continue;
        for a2 := 0 to 16 do
        begin
          if y = -1 then
            if a2 = 16 then
              y := -1
            else if (((a2 * x) <= 3500) and ((3500 - (a2 * x)) mod (16 - a2) = 0)) then
              y := (3500 - (a2 * x)) div (16 - a2)
            else
              continue;
          if (y = 0) or (x <= y) then
            break;
          if ((a2 * x) <= 3500) and ((y = -1) or (a2 * x + (16 - a2) * y = 3500)) then
            for a3 := 0 to 26 do
            begin
              if y = -1 then
                if a3 = 26 then
                  y := -1
                else if (((a3 * x) <= 3500) and ((3500 - (a3 * x)) mod (26 - a3) = 0)) then
                  y := (3500 - (a3 * x)) div (26 - a3)
                else
                  continue;
              if (y = 0) or (x <= y) then
                break;
              if (a3 * x + (26 - a3) * y = 3500) then
                writeln('x = ', x, ', y = ', y,
                  ', a1 = ', a1,
                  ', b1 = ', 10 - a1,
                  ', check S1 = ', a1 * x + (10 - a1) * y,
                  ', a2 = ', a2,
                  ', b2 = ', 16 - a2,
                  ', check S2 = ', a2 * x + (16 - a2) * y,
                  ', a3 = ', a3,
                  ', b3 = ', 26 - a3,
                  ', check S3 = ', a3 * x + (26 - a3) * y
                  );
            end;
        end;
      end;
  writeln('Finished');
end.

Время счета снижается с 58 до примерно 0,001 секунды, то есть на 4-5 порядков.