用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

这是一个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 2

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

这些问题的解决方案被称为整数组合，它是整数分区的推广(其中顺序并不重要，即只考虑在排列下唯一的答案)。

例如，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).
• If )。如果你能够使用另一种语言(通过与R接口，或者通过预先计算答案)，那么Mathematica就有一个计算组合的函数：Compositions.
• Sage是免费的(不像)，并且还具有一个生成内置组合的函数：Compositions。值得注意的是，这是使用实现的，它可以使用更多的内存efficiently.
• Python：查看this堆栈溢出问题(它生成分区，然后您可以对其进行置换)。我做了类似的here (尽管它使用itertools来查找排列，然后我需要过滤出唯一的排列，所以通过使用专门针对多集的排列算法可以更有效地完成这项工作)。

为了更好地理解算法(或自己实现它们)，您可以查看这本不完整但很有用的电子书： by Frank Ruskey，它展示了如何在常量摊销时间(CAT)内生成分区。因为您想要组合，所以您还可以使用CAT算法生成排列(也在本书中)，以生成每个整数分区的排列。

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

我相信这些在Knuth的计算机编程艺术卷4A中也有很好的介绍，如果你手头有的话。

ElKamina's suggestion to solve it recursively是一个很好的方法，但是我不会对大的H使用这种方法；因为使用R (as well as Python) doesn't optimise tail-calls，您可能会导致堆栈溢出。

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的每一个组合的解的集合。