生成 ai=1 的线性丢番图方程的所有解的高效算法

Question

我正在尝试为给定 H 生成以下方程的所有解。

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 有更好的算法吗？

Answer 1

与许多问题一样，当某些术语已知时，解决方案会变得更容易找到/研究。

这些问题的解决方案被称为 Integer Compositions，在术语上是 Integer Partitions 的概括（其中顺序无关紧要，即只回答在排列下是唯一的）。

例如，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 中工作，the partitions package 可以为您生成分区。您需要找到每个分区的唯一排列才能获得组合（请参阅this SO question）。
• 如果您能够使用另一种语言（通过与 R 交互，或通过预先计算答案），那么 Mathematica 具有计算合成的功能：Compositions。
• Sage 是免费的（与 Mathematica 不同），并且还具有生成内置合成的功能：Compositions。值得注意的是，这个是使用generators实现的，可能会更有效地使用内存。
• Python：请参阅this Stack Overflow 问题（生成分区，然后您可以对其进行置换）。我做了类似here 的事情（虽然它使用itertools 来查找排列，然后我需要过滤独特的排列，所以这可以通过使用专门用于multisets的排列算法更有效地完成)。

为了更好地理解算法（或自己实现它们），您可以查看这本不完整但有用的电子书：Combinatorial GenerationFrank Ruskey，它展示了如何在 恒定摊销时间 (CAT)。由于您想要组合，您还可以使用 CAT 算法生成排列（也在书中）来生成每个整数分区的排列。

Ruskey 还解释了如何对它们进行排名和unrank，这对于存储/散列结果非常方便。

我相信 Knuth 的 计算机编程艺术第 4A 卷 中也很好地介绍了这些内容，如果你碰巧得心应手的话。

ElKamina's suggestion to solve it recursively 是一个不错的方法，但我不会将这种方法用于大 H；因为R (as well as Python) doesn't optimise tail-calls，你最终可能会出现堆栈溢出。

Answer 2

这是 C++ 中的implementation

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 中的implementation：

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 中的 implementation

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)

Answer 3

我假设您没有尝试同时解方程。

您可以使用递归或动态编程来解决这个问题。

如果您使用递归，只需为第一个变量分配一个有效值，然后递归求解其余部分。

这里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 和 sum 的每个组合的解集。