Используйте алгоритм редукции базиса решетки Ленстры–Ленстры–Ловаша, чтобы найти нормальную форму Эрмита вашей системы.

Это просто в Matlab http://fr.mathworks.com/help/symbolic/mupad_ref/linalg-hermiteform.html

Проверьте NTL на c++ http://www.shoup.net/ntl/index.html

Следуйте этому примеру: предположим, что уравнения:

7 = x + 3y + 4z + 9w
12 = 4x + 2y + 3z + 7w

Есть 2 уравнения и 4 неизвестных. Вы можете установить 2 неизвестных в качестве параметров, поэтому в системе будет столько же уравнений, сколько и неизвестных:

7 - (4z + 9w) = x + 3y
12 - (3z + 7w) = 4x + 2y

Мы будем использовать x и y в качестве неизвестных. Можно решить его и оставить w и z в качестве параметров в линейном решении:

x = (22 - 3w - z)/10
y = (16 - 29w - 13z)/10

Теперь вам нужно сделать так, чтобы числители делились на 10, а знаменатель. Если есть решение, можно протестировать все параметры от 1 до 10.

В общем, вы делаете это:

• Подберите параметры так, чтобы неизвестных было столько же, сколько уравнений. Попробуйте оставить неизвестные, которые генерируют наименьшее абсолютное значение для определителя (в примере это было 10, но выбор y и z был бы лучше, так как это было бы |det|=3)
• Решите линейную систему и сгенерируйте ответ в зависимости от параметров
• Протестируйте значения параметров от 1 до абсолютного значения det (если есть решение, вы найдете его в этой точке), пока не будет целочисленное решение для всех неизвестных (поэтому вы выбираете наименьший определитель абсолютное значение) и проверьте, является ли оно положительным (это не было показано в примере).

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

РЕДАКТИРОВАТЬ1

Существует способ избежать перебора последней части. Опять же, в примере вы должны сделать:

22 = 3w + z (congruent, mod 10)
16 = 29w + 13z (congruent, mod 10)

Применение модуля:

2 = 3w + z (congruent, mod 10)
6 = 9w + 3z (congruent, mod 10)

Вы можете сделать систему сравнений треугольной (исключение Гаусса), умножая сравнение на обратное по модулю 10 и суммируя с другими. Обратное число 9 по модулю 10 равно -1, поэтому мы умножаем последнее сравнение:

2 = 3w + z (congruent, mod 10)
-6 = -9w + -3z (congruent, mod 10)

Эквивалентно:

2 = 3w + z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Затем умножьте на -3 и прибавьте к первому:

2 - 3*4 = 3w + z -3w - 21z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Эквивалентно:

-10 = -20z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Затем вы решаете сверху вниз (этот пример кажется верным для любого значения z). Может быть некоторая степень свободы, если у вас больше параметров, чем сравнений, но они эквивалентны, и вы можете установить лишние параметры на любое значение, предпочтительно наименьшее, равное 1.

Дайте мне знать, если есть что-то еще не ясно!

Я бы попробовал следующий подход (обратите внимание, что мне никогда не приходилось сталкиваться с этой проблемой, поэтому воспринимайте этот ответ как попытку решить проблему вместе с вами, а не настоящий аналитический ответ).

Вы просто находите решения, как если бы это была обычная система, поэтому вы можете найти бесконечно много решений:

Пример:

y = x + 3

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

Конечно, у вас есть ограничения, в нашем случае решение имеет только 1 свободную переменную, поэтому мы можем написать все решения следующим образом:

(x, x+3)

Сюрприз:

Если где-то есть иррациональное число, целочисленных решений нет:

(x, x+pi)  => neither 1st or 2nd number here can be whole at same time

Таким образом, вы, вероятно, можете найти целочисленные решения тогда и только тогда, когда ваши «бесконечно много решений» ограничены целыми или рациональными числами.

Предположим, у вас есть следующий вектор:

 ( 3x, (2/5)y, y, x, x+y)

Тогда целое решение может быть:

 x=3
 y=10/2

Если вы чувствуете, что ответ вас не удовлетворяет, просто скажите: я удалю его, чтобы не получать очки за вознаграждение.

Возможно, вы захотите взглянуть на решатели ограничений, такие как Z3. Также есть JavaExample.

Полезное руководство, объясняющее некоторые возможности Z3, можно найти здесь.

РЕДАКТИРОВАТЬ: Также взгляните на алгоритм решения систем линейных неравенств

Это взято из соответствующего раздела Википедии.

1. Перепишите уравнения с матричной записью AX=C.
2. Вычислить нормальную форму Смита числа A. Грубо говоря, это нахождение B=UAV, где B[i][j] == 0, если i != j.
3. B(inverse(V))X=UC
4. Так как B[i][j] == 0 если i != j, найти значение (inverse(V))X и, следовательно, X несложно.

Если я правильно понимаю проблему, у вас есть несколько посылок, каждая из которых имеет разную стоимость доставки. Вы хотите оплатить эти почтовые расходы из имеющегося у вас пула марок. Можно решить проблему с помощью целочисленного программирования. Приведенное ниже решение подходит для всех пакетов одновременно. У вас будет количество переменных, равное numPackages*numStampDenominations (вероятно, неудобно для большого количества пакетов).

Ограничение равенства выглядит так: Aeqx = beq. Матрица Aeq для двух пакетов с 4 номиналами марок (numVarsnumPackages) имеет вид:

.21 .68 .47 .01 .00 .00 .00 .00           .89
                                   * x = 
.00 .00 .00 .00 .21 .68 .47 .01           .50 

Первая строка — это коэффициенты ограничений (значения штампа) для пакета 1. Ненулевые коэффициенты — это значения штампа. Переменная с нулевыми коэффициентами, связанная с другими пакетами, не имеет значения.

Второй набор ограничений (неравенство) связан с имеющимся у меня набором марок. Ограничение неравенства выглядит как A*x = b. Матрица A для 4 номиналов марок и 8 переменных (numPackages * numStampDenominations) выглядит так:

1 0 0 0 1 0 0 0         10

0 1 0 0 0 1 0 0         10
                 * x <=  
0 0 1 0 0 0 1 0         10

0 0 0 1 0 0 0 1         20

Функция стоимости — это f'*x, которая представляет собой просто общее количество марок. Его коэффициенты - всего лишь единицы в виде вектор-столбца.

Ниже приведен сценарий, запускаемый в Matlab, который решает проблему описанным способом. Вероятно, это будет сформулировано примерно так же в октаве, предложении GNU, похожем на Matlab. Matlab или octave могут быть неподходящим решением для вас, но они, по крайней мере, говорят вам, что работает, и дают вам песочницу для выработки решения.

% The value of each stamp available as a 4x1 matrix
sv = [.21; .68; .47; .01];
% The number of each stamp available as a 4x1 matrix
sq = [10;10;10;40];
% Number of demominations
m = size(sv, 1);
% The postage required for each package as a 4x1 matrix
% -- probably read from a file
postage = [.89;.50;1.01;.47;.47];
% Number of packages -- just a count of the postage rows
packageCount = size(postage, 1);
% Number of variables is m*packageCount
numVar = m*packageCount;
% lower bound -- zero stamps of a given denomination
lb = zeros(numVar,1);
% upper bound -- use constraints for upper bound
ub = [];
% The cost function -- minimize the number of stamps used
% min(f'*x) 
f = ones(numVar,1);
% integer constraints
intcon = 1:numVar;
% Construct postage constraint matrix
Aeq = zeros(numVar, packageCount);

for i = 1:packageCount
    first = i*m - 3;
    last = first + 3;
    Aeq(first:last,i) = sv(:);
end

% Construct stamp count constraint matrix
A = zeros(numVar, m);

for r = 1:m
    for j = 1:m
        c = (j-1)*m + 1;
        A(c,r) = 1;
    end
end

x = intlinprog(f, intcon, A', b, Aeq', beq, lb, ub)