Efficient algorithm to generate all solutions of a linear diophantine equation with ai=1

Question

I am trying to generate all the solutions for the following equations for a given H.

With 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

For my problem, there are always 4 equations to solve (independently from the others). There are a total of 2^(H-1) solutions. For the previous one, here are the solutions :

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

Here is an R algorithm which solve the problem.

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,])
}

However, this algorithm makes more calculations than needed. I am sure it is possible to go faster. By that, I mean not generating the permutations for which the sums is > H

Any idea of a better algorithm for a given H ?

Answer 1

Here's an implementation in 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;
}

commandline stuff:

$ 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
$

Here's an implementation in 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

Here's an implementation in 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)

Answer 2

As with many problems, the solution becomes much easier to find/research when some terminology is known.

The solutions to these problems are known as Integer Compositions, which are in term a generalisation of Integer Partitions (where the order doesn't matter, i.e. only answers that are unique under permutation are considered).

For example, the integer partitions of 4 are: 1+1+1+1, 1+1+2, 1+3, 2+2, 4, whereas the integer compositions of 4 are: 1+1+1+1, 1+1+2, 1+2+1, 2+1+1, 1+3, 3+1, 2+2, 4.

There are a few implementations readily available (references to language-agnostic algorithms follows):

• Since you are working in R, the partitions package can generate partitions for you. You would need to find unique permutations of each partition to get the compositions (see this SO question).
• If you are able to use another language (either by interfacing with R, or by precomputing the answer) then Mathematica has a function to compute Compositions: Compositions.
• Sage is free (unlike Mathematica), and also has a function to generate Compositions built in: Compositions. It is worth noting that this one is implemented using generators, which may use memory more efficiently.
• Python: see this Stack Overflow question (which generates partitions, which you can then permute). I did something similar here (although it uses itertools to find permutations, which I then need to filter for unique permutations, so this could be done more efficiently by using a permutation algorithm specifically for multisets).

In order to understand the algorithms better (or implement them yourself), you could check out this incomplete but useful ebook: Combinatorial Generation by Frank Ruskey, which shows how to generate partitions in constant amortized time (CAT). Since you want compositions, you can also use a CAT algorithm for generating permutations (also in the book) to generate permutations of each integer partition.

Ruskey also explains how to rank and unrank them, which can be handy for storing/hashing the results.

I believe these are also nicely covered in Knuth's The Art of Computer Programming Volume 4A, if you happen to have it handy.

ElKamina's suggestion to solve it recursively is a good one, but I wouldn't use this approach for large H; since R (as well as Python) doesn't optimise tail-calls, you may end up with a stack overflow.

Answer 3

I assume you are not trying to simultaneously solve the equations.

You can either use recursion or dynamic programming to solve this.

If you are using recursion, just assign a valid value to the first variable and solve the rest of it recursively.

Here n is the number of variables and sum is the sum. cursol is the partial solution (initially set to [] )

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])

If you want to use dynamic programming, you have to remember the set of solutions for each combination of n and sum.