User login

You are here

Home » Blogs » st's blog

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

Fri, 24/10/2014 - 01:46 - 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 (у меня древняя версия 7). Для Си программу надо слегка переделать.

#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 порядков.

Математическое программирование
»

Comments

Submitted by Vlad_AB on Fri, 24/10/2014 - 11:12.

Приведённая Вами задача не является задачей линейного программирования, во-первых, потому, что отсутствует целевая функция, а во-вторых, потому, что система ограничений не является линейной системой уравнений (присутствуют произведения переменных).

Вашу программу можно немного оптимизировать, если учесть имеющуюся в условии информацию о ценах до и после обеда, а именно, что x>y.

Методов решения классической задачи линейного программирования несколько, но все они требуют наличия целевой функции.

P.S. Запустил программку на HP Prime, считает уже полчаса. Результаты выложу позже.

»
st's picture

Submitted by st on Fri, 24/10/2014 - 11:39.

Была опечатка в третьем операторе сравнения, я исправил, добавил оптимизацию X>Y. Решение получается одно. ЦФ здесь - полная распродажа товара :)

»
Vitasam's picture

Submitted by Vitasam on Fri, 24/10/2014 - 11:56.

А вообще налицо сговор. Антимонопольный комитет дремлет :)

»

Submitted by Vlad_AB on Fri, 24/04/2015 - 05:17.

.

»

Submitted by Vlad_AB on Fri, 24/10/2014 - 17:32.

Для HP Prime пришлось ещё раз оптимизировать программку, т.к. в указанном варианте время счёта ОЧЕНЬ велико (задачка не для калькулятора). Теперь программка отбрасывает часть вариантов, среди которых гарантированно нет решения. Предполагаемое время счёта от 15 до 30 часов.

Текст программы:

EXPORT KURKI()
   BEGIN
      PRINT();
      FOR A FROM 0 TO 10 DO
         FOR B FROM 0 TO 16 DO
            FOR C FROM 0 TO 26 DO
               FOR X FROM 1 TO 3499 DO
                  Y:=1;
                  IF (A*X+(10-A)*Y)>3500 OR (B*X+(16-B)*Y)>3500 OR (C*X+(26-C)*Y)>3500 THEN
                     BREAK;
                  END;
                  FOR Y FROM 1 TO X-1 DO
                     IF (A*X+(10-A)*Y)>3500 OR (B*X+(16-B)*Y)>3500 OR (C*X+(26-C)*Y)>3500 THEN
                        BREAK;
                     END;
                     IF (A*X+(10-A)*Y)==3500 AND (B*X+(16-B)*Y)==3500 AND (C*X+(26-C)*Y)==3500 THEN
                        PRINT("X="+STRING(X)+"; Y="+STRING(Y));
                        PRINT("a1="+STRING(A)+"; b1="+STRING(10-A));
                        PRINT("a2="+STRING(B)+"; b2="+STRING(16-B));
                        PRINT("a3="+STRING(C)+"; b3="+STRING(26-C));
                        PRINT("------------------------");
                     END;
                  END;
               END;
            END;
         END;
      END;
   END;
END;
»

Submitted by Ardo_79 on Thu, 12/03/2015 - 22:09.

Для оптимизации можно использовать условие B

»

Submitted by Vlad_AB on Fri, 24/10/2014 - 18:11.

HP 50g по оптимизированной программе посчитал примерно за 15-16 минут.

Текст программы KURKI:

#include <hpgcc49.h>
 
int main(void)
 {
	int a1, a2, a3;
	int x, y;
	clear_screen();
	sys_slowOff();
	for (a1 = 0; a1 <= 10; a1++)
		for (a2 = 0; a2 <= 16; a2++)
			for (a3 = 0; a3 <= 26; a3++)
				for (x = 1; x < 3500; x++)
				{
				y=1;
				if((a1 * x + (10 - a1) * y > 3500) && (a2 * x + (16 - a2) * y > 3500) && (a3 * x + (26 - a3) * y > 3500))
					break;
				for (y = 1; y < x; y++)
					{
					if((a1 * x + (10 - a1) * y > 3500) && (a2 * x + (16 - a2) * y > 3500) && (a3 * x + (26 - a3) * y > 3500))
						break;
					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");
	sys_slowOn();
	beep();
	WAIT_CANCEL;
	return(0);
 }
»
st's picture

Submitted by st on Fri, 24/10/2014 - 18:15.

Добавил в текст оптимизированный вариант, на HP 50g должен считать за секунды.

»

Submitted by Vlad_AB on Fri, 24/10/2014 - 18:35.

Проверил Ваш вариант оптимизированной программки. HP 50g посчитал за 18 секунд: решение выдаёт через секунду, а остаток вариантов проверяет 17 секунд.

»
st's picture

Submitted by st on Fri, 24/10/2014 - 20:27.

Добавил еще одну оптимизацию :)

»

Submitted by Vlad_AB on Sat, 25/10/2014 - 15:55.

У HP Prime время расчёта по программе, которая на HP 50g выполняется 18 секунд, составило примерно 75 минут.

»
st's picture

Submitted by st on Sat, 25/10/2014 - 18:45.

Наверняка какая-то ошибка в программе или программы разные.

»

Submitted by Vlad_AB on Sat, 25/10/2014 - 20:14.

Ошибку не нашёл, а программы такие:

EXPORT KURKI3()
BEGIN
 
FOR X FROM 1 TO 3499 DO
	FOR Y FROM 1 TO X-1 DO
		FOR A FROM 0 TO 10 DO
			IF A*X+(10-A)*Y==3500 THEN
				FOR B FROM 0 TO 16 DO
					IF B*X+(16-B)*Y==3500 THEN
						FOR C FROM 0 TO 26 DO
							IF C*X+(26-C)*Y==3500 THEN
								PRINT("X="+STRING(X)+"; Y="+STRING(Y));
								PRINT("a1="+STRING(A)+"; b1="+STRING(10-A));
								PRINT("a2="+STRING(B)+"; b2="+STRING(16-B));
								PRINT("a3="+STRING(C)+"; b3="+STRING(26-C));
								PRINT("----------------");
							END;
						END;
					END;
				END;
			END;
		END;
	END;
END;
 
PRINT("END.");
 
END;
#include <hpgcc49.h>
 
int main()
 {
	int a1, a2, a3;
	int x, y;
	clear_screen();
	sys_slowOff();
	for (x = 1; x < 3500; x++)
		for (y = 1; y < x; y++)
			for (a1 = 0; a1 <= 10; a1++)
				if (a1 * x + (10 - a1) * y == 3500)
					for (a2 = 0; a2 <= 16; a2++)
						if (a2 * x + (16 - a2) * y == 3500)
							for (a3 = 0; a3 <= 26; a3++)
								if (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");
	sys_slowOn();
	beep();
	WAIT_CANCEL;
	return(0);
 }

Запустил программку на HP Prime ещё раз, результат прежний (75 минут плюс/минус пара секунд). Решение появляется через 52 секунды. На эмуляторе результат появляется примерно через 5 секунд, а общее время счёта около 6 с половиной минут.

»
st's picture

Submitted by st on Sat, 25/10/2014 - 22:22.

А, так это программа на Си для HP 50g, я-то думал, что на UserRPL. Тогда результаты объяснимые. Со второй оптимизацией время должно снизиться еще раз в 100.

»

Submitted by Vlad_AB on Fri, 24/04/2015 - 05:18.

.

»

Submitted by Vlad_AB on Mon, 27/10/2014 - 13:16.

Сегодня после работы попробую написать и запустить на HP 50g (опять на C) и на HP Prime программку со второй оптимизацией.

P.S. Если бы это была задача линейного программирования, то можно было попробовать некоторые алгоритмы. Упоминание линейного программирования в названии этой записи некорректно (прошу прощения за нахальство). Правильнее говорить либо о целочисленном программировании, либо о математическом или нелинейном программировании.

P.P.S. По следующей ссылке, вероятно, один из возможных первоисточников задачи (есть интересные комментарии).

»
st's picture

Submitted by st on Mon, 27/10/2014 - 15:14.

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

»

Submitted by pmk_emu_user on Wed, 18/02/2015 - 22:52.

Заинтересовался. Набросал программку на питоне. Алгоритм таков: для первого продавца, у которого 10 кур, методом перебора находится совпадение суммы выручки до обеда и выручки после обеда с числом 3500. Затем эти цены подставляются двум другим продавцам. Если и у них сумма 3500, то выдается распечатка, кто, сколько и почем продавал. Программа не оптимизирована, находит первый верный вариант, а затем ещё 4 варианта, когда до обеда продавцы продавали своих кур, а после обеда плюнули на всё и раздали оставшихся бесплатно, т.е. цена у = 0 :)

#! /usr/bin/env python
# -*- coding: utf8 -*-

def other(x, y, n):
for a in range(n):
if a * x + (n - a) * y == 3500: return (True, a, n - a)
return (False, 0, 0)

for x in range(3499):
for y in range (x):
for a1 in range(10):
if a1 * x + (10 - a1) * y == 3500 :
b = other(x, y, 16)
c = other(x, y, 26)
if b[0] and c[0]:
print 'a1 = ', a1, 'x = ', x, 'b1 = ', 10 - a1, 'y = ', y
print 'a2 = ', b[1], 'x = ', x, 'b2 = ', b[2], 'y = ', y
print 'a3 = ', c[1], 'x = ', x, 'b3 = ', c[2], 'y = ', y
print '\n\n'

На моём HP550, linuxmint 17 первый вариант найден примерно за 1 секунду. Весь перебор занимает 27.641 сек.

»

Submitted by Ajaja on Fri, 13/03/2015 - 01:30.

Задача заинтересовала. Написал решение для HP-42s. Использовал Free42 на Андроиде. На стареньком телефоне Samsung Galaxy S5660 находит решение и перебирает все остальные варианты почти мгновенно. Думаю, на реальном калькуляторе работала бы тоже не очень долго.
00 { 214-Byte Prgm }
01>LBL "CHICK"
02 135
03 STO "Y"
04>LBL "Y" @ начало цикла Y: 135..1 (3500/26 ~ 135 - максимальный Y)
05 10
06 STO "A1"
07>LBL "A1" @ начало цикла A1: 10..1 (если A1 = 0, то Y = 350, что > 135)
08 3500 @ X = (3500 - 10*Y + A1*Y)/A1
09 10
10 RCL "Y"
11 ×
12 -
13 RCL "A1"
14 RCL "Y"
15 ×
16 +
17 RCL ST X
18 RCL "A1"
19 MOD
20 X!=0?
21 GTO 99
22 Rv
23 LASTX
24 ÷
25 X<=0?
26 GTO 99
27 RCL "Y"
28 X>=Y?
29 GTO 99
30 Rv
31 STO "X"
32 3500 @ A2 = (3500-16*Y)/(X-Y)
33 16
34 RCL "Y"
35 ×
36 -
37 RCL ST X
38 RCL "X"
39 RCL "Y"
40 -
41 MOD
42 X!=0?
43 GTO 99
44 Rv
45 LASTX
46 ÷
47 X<=0?
48 GTO 99
49 16
50 X<Y?
51 GTO 99
52 Rv
53 STO "A2"
54 3500 @ A3 = (3500-26*Y)/(X-Y)
55 26
56 RCL "Y"
57 ×
58 -
59 RCL ST X
60 RCL "X"
61 RCL "Y"
62 -
63 MOD
64 X!=0?
65 GTO 99
66 Rv
67 LASTX
68 ÷
69 X<=0?
70 GTO 99
71 26
72 X<Y?
73 GTO 99
74 Rv
75 STO "A3"
76 RCL "A2"
77 RCL "A1"
78 STOP @ остановка при нахождении решения, продолжение вычислений по R/S
79>LBL 99
80 DSE "A1" @ конец цикла A1
81 GTO "A1"
82 DSE "Y" @ конец цикла Y
83 GTO "Y"
84 99999999 @ конец вычислений
85 STOP
86 END

»
Vitasam's picture

Submitted by Vitasam on Fri, 13/03/2015 - 15:54.

Как можно загнать проще всего загнать эту программу во Free42?

»

Submitted by Ajaja on Fri, 13/03/2015 - 16:59.

Проще всего через импорт. Залил программу сюда:chick.raw

»

Submitted by Ardo_79 on Sat, 14/03/2015 - 14:06.

A*X + (10 - A)*Y = 3500
B*X + (16 - B)*Y = 3500
C*X + (26 - C)*Y = 3500
Если учесть что X>Y и всегда будет соблюдаться условие (A>B>C при A<=10) - Это условие не может не соблюдаться,как в свою очередь X не может быть меньше 350 (X>=350),то можно написать программу на CASIO 9860
For 350→X To 3499
For 1→Y To X-1
X-Y→Z
For 1→A To 10
If A*Z+10*Y=3500
Then
For 1→B To 9
If B*Z+16*Y=3500
Then
For 1→C To 8
If C*Z+26*Y=3500
Then
X¸
Y¸
A¸
B¸
C¸
IfEnd
Next
IfEnd
Next
IfEnd
Next
Next
Next
Время до появления результатов 7:50 , при разогнанном процесоре в режиме (Quadruple)
время 4:10
Есть еще старая задача 3-го класса церковно-приходской школы ,когда надо было на 100 копеек
купить 100 животных. Корова стоила 5 копеек. овца 3 копейки и десяток гусей стоил 1 копейку.

»

Submitted by Ajaja on Sat, 14/03/2015 - 20:51.

[quote=Ardo_79]Есть еще старая задача 3-го класса церковно-приходской школы ,когда надо было на 100 копеек
купить 100 животных. Корова стоила 5 копеек. овца 3 копейки и десяток гусей стоил 1 копейку.[/quote]Ответ: 16 коров, 4 овцы и 80 гусей :)
Решил в Free42, находит мгновенно:

00 { 205-Byte Prgm } 26 |-" "
01>LBL "POP"         27 RCL "SHEEP"
02 CLST              28 AIP
03 STO "GOOSE"       29 |-" "
04>LBL "GOOSE"       30 RCL "GOOSE"
05 CLST              31 10
06 STO "COW"         32 ×
07>LBL "COW"         33 AIP
08 100               34 AVIEW
09 10                35 STOP
10 RCL× "GOOSE"      36>LBL "NEXT"
11 RCL+ "COW"        37 ISG "COW"
12 -                 38 DEG
13 STO "SHEEP"       39 RCL "COW"
14 3                 40 20
15 ×                 41 X>=Y?
16 RCL+ "GOOSE"      42 GTO "COW"
17 5                 43 ISG "GOOSE"
18 RCL× "COW"        44 DEG
19 +                 45 RCL "GOOSE"
20 100               46 10
21 X!=Y?             47 X>=Y?
22 GTO "NEXT"        48 GTO "GOOSE"
23 "SOLUTION "       49 "THE END"
24 RCL "COW"         50 AVIEW
25 AIP               51 .END.

Программа тут

»
Vitasam's picture

Submitted by Vitasam on Tue, 17/03/2015 - 22:06.

Многобуквенные Labels и переменные радуют после HP35s

»