Эффективный алгоритм генерации всех решений линейного диофантова уравнения с ai=1
Я пытаюсь сгенерировать все решения для следующих уравнений для данного H.
С Н =4:
1) ALL solutions for x_1 + x_2 + x_3 + x_4 =4
2) ALL solutions for x_1 + x_2 + x_3 = 4
3) ALL solutions for x_1 + x_2 = 4
4) ALL solutions for x_1 =4
Для моей задачи всегда нужно решить 4 уравнения (независимо от других). Всего существует 2^(H-1) решения. Для предыдущего, вот решения:
1) 1 1 1 1
2) 1 1 2 and 1 2 1 and 2 1 1
3) 1 3 and 3 1 and 2 2
4) 4
Вот алгоритм R, который решает проблему.
library(gtools)
H<-4
solutions<-NULL
for(i in seq(H))
{
res<-permutations(H-i+1,i,repeats.allowed=T)
resum<-apply(res,1,sum)
id<-which(resum==H)
print(paste("solutions with ",i," variables",sep=""))
print(res[id,])
}
Однако этот алгоритм выполняет больше вычислений, чем необходимо. Я уверен, что можно идти быстрее. Под этим я подразумеваю не генерировать перестановки, для которых суммы>H
Любая идея лучшего алгоритма для данного H?
3 ответа
Вот реализация в C++
blah.cpp:
#include <stdlib.h>
#include <iostream>
#include <vector>
using namespace std;
vector<int> ilist;
void diophantine(int n)
{
size_t i;
if (n==0)
{
for (i=0; i < ilist.size(); i++) cout << " " << ilist[i];
cout << endl;
}
else
{
for (i=n; i > 0; i--)
{
ilist.push_back(i);
diophantine(n-i);
ilist.pop_back();
}
}
}
int main(int argc, char** argv)
{
int n;
if (argc == 2 && (n=strtol(argv[1], NULL, 10)))
{
diophantine(n);
}
else cout << "usage: " << argv[0] << " <Z+>" << endl;
return 0;
}
материал командной строки:
$ g++ -oblah blah.cpp
$ ./blah 4
4
3 1
2 2
2 1 1
1 3
1 2 1
1 1 2
1 1 1 1
$
Вот реализация в bash:
blah.sh:
#!/bin/bash
diophantine()
{
local i
local n=$1
[[ ${n} -eq 0 ]] && echo "${ilist[@]}" ||
{
for ((i = n; i > 0; i--))
do
ilist[${#ilist[@]}]=${i}
diophantine $((n-i))
unset ilist[${#ilist[@]}-1]
done
}
}
RE_POS_INTEGER="^[1-9]+$"
[[ $# -ne 1 || ! $1 =~ $RE_POS_INTEGER ]] && echo "usage: $(basename $0) <Z+>" ||
{
declare -a ilist=
diophantine $1
}
exit 0
Вот реализация в Python
blah.py:
#!/usr/bin/python
import time
import sys
def output(l):
if isinstance(l,tuple): map(output,l)
else: print l,
#more boring faster way -----------------------
def diophantine_f(ilist,n):
if n == 0:
output(ilist)
print
else:
for i in xrange(n,0,-1):
diophantine_f((ilist,i), n-i)
#crazy fully recursive way --------------------
def diophantine(ilist,n,i):
if n == 0:
output(ilist)
print
elif i > 0:
diophantine(ilist, n, diophantine((ilist,i), n-i, n-i))
return 0 if len(ilist) == 0 else ilist[-1]-1
##########################
#main
##########################
try:
if len(sys.argv) == 1: x=int(raw_input())
elif len(sys.argv) == 2: x=int(sys.argv[1])
else: raise ValueError
if x < 1: raise ValueError
print "\n"
#diophantine((),x,x)
diophantine_f((),x)
print "\nelapsed: ", time.clock()
except ValueError:
print "usage: ", sys.argv[0], " <Z+>"
exit(1)
Как и во многих проблемах, решение становится намного легче найти / исследовать, когда известна некоторая терминология.
Решения этих проблем известны как целочисленные композиции, которые в терминах являются обобщением целочисленных разделов (где порядок не имеет значения, т.е. рассматриваются только ответы, которые являются уникальными при перестановке).
Например, целочисленные разделы 4: 1+1+1+1, 1+1+2, 1+3, 2+2, 4, тогда как целочисленные составы 4: 1+1+1+1, 1+1+2, 1+2+1, 2+1+1, 1+3, 3+1, 2+2, 4.
Есть несколько доступных реализаций (ссылки на не зависящие от языка алгоритмы приведены ниже):
- Поскольку вы работаете в R,
partitionsПакет может генерировать разделы для вас. Вам нужно будет найти уникальные перестановки каждого раздела, чтобы получить композиции (см. Этот вопрос SO). - Если вы можете использовать другой язык (либо взаимодействуя с R, либо предварительно вычисляя ответ), то в Mathematica есть функция для вычисления Compositions:
Compositions, - Sage бесплатен (в отличие от Mathematica), а также имеет функцию для генерации встроенных композиций:
Compositions, Стоит отметить, что этот реализован с использованием генераторов, которые могут использовать память более эффективно. - Python: посмотрите этот вопрос переполнения стека (который генерирует разделы, которые вы затем можете переставлять). Я сделал что-то подобное здесь (хотя он использует
itertoolsчтобы найти перестановки, которые мне затем нужно отфильтровать для уникальных перестановок, чтобы это можно было сделать более эффективно, используя алгоритм перестановок специально для мультимножеств).
Чтобы лучше понять алгоритмы (или реализовать их самостоятельно), вы можете обратиться к этой неполной, но полезной книге: " Комбинаторная генерация " Фрэнка Руски, в которой показано, как генерировать разделы в постоянном амортизированном времени (CAT). Поскольку вам нужны композиции, вы также можете использовать алгоритм CAT для генерации перестановок (также в книге) для генерации перестановок каждого целочисленного раздела.
Руски также объясняет, как их ранжировать и отменять, что может быть полезно для сохранения / хеширования результатов.
Я полагаю, что они также хорошо освещены в книге Кнута " Искусство компьютерного программирования", том 4А, если она вам пригодится.
Предложение Эль-Камины рекурсивно решить ее хорошо, но я бы не стал использовать этот подход для больших Н; поскольку R (как и Python) не оптимизирует хвостовые вызовы, вы можете получить переполнение стека.
Я предполагаю, что вы не пытаетесь одновременно решить уравнения.
Вы можете использовать рекурсию или динамическое программирование, чтобы решить эту проблему.
Если вы используете рекурсию, просто присвойте правильное значение первой переменной и решите остальную часть рекурсивно.
Здесь n - количество переменных, а sum - сумма. cursol - это частичное решение (изначально установлено в [])
def recSolve(n,sum, cursol):
if n==1:
print cursol + [sum]
return
if n == sum:
print cursol + [1 for i in range(n)]
return
else:
for i in range(1, sum-n+2):
recsolve(n-1, sum-i, cursol+[i])
Если вы хотите использовать динамическое программирование, вы должны запомнить набор решений для каждой комбинации n и суммы.