[Python_4_President]

Ok, now it seems to work properly. I assume that if you go through an entire series of some particular tens (ie, 70-80) without finding any non-solutions, you're probably done. At least, my tests would indicate 10 is a good number. This is my first mathy combinatorial problem =)

Spoiler

def calc(packages, looping_vals):
    packlen = len(packages)
    return sum([packages[i]*looping_vals[i] for i in range(packlen)])

def gen_for(packages, maxVal):
    iNeed = len(packages)
    c = """def calc(packages, looping_vals):
    packlen = len(packages)
    return sum([packages[i]*looping_vals[i] for i in range(packlen)])\n\n\n"""
    d = ""
    c += "nope = []\n"
    c += "ld = []\n"
    c += "maxVal = {0}\n".format(maxVal)
    c += "threshold = 10\n"
    c += "lastVal = 0\n"
    c += "single_increments = 0\n"
    c += "try:\n"
    c += "    for n in xrange(1, maxVal):\n"
    for i in range(0, iNeed):
        c += (" "*4)*(i+2)+"for x{0} in xrange(n):\n".format(i)
        d += "x{0},".format(i)
    d = "["+d.rstrip(",")+"]"
    c += (" "*4)*(i+3)+"comb = calc({0}, {1})\n".format(packages, d)
    c += (" "*4)*(i+3)+"if n == comb:\n"
    c += (" "*4)*(i+4)+"if n not in ld:\n"
    c += (" "*4)*(i+5)+"print 'You can buy {0}. n-lastVal: {1}, single_increments: {2}'.format(n, n-lastVal, single_increments)\n"
    c += (" "*4)*(i+5)+"if n-lastVal == 1:\n"
    c += (" "*4)*(i+6)+"single_increments += 1\n"
    c += (" "*4)*(i+5)+"else:\n"
    c += (" "*4)*(i+6)+"single_increments = 0\n"
    c += (" "*4)*(i+5)+"lastVal = n\n"
    c += (" "*4)*(i+5)+"ld.append(n)\n"
    c += (" "*4)*(i+5)+"if single_increments == threshold: raise StopIteration()\n"
    c += (" "*4)*(i+3)+"else:\n"
    c += (" "*4)*(i+4)+"if n not in nope and n not in ld: nope.append(n)\n"
    c += "except StopIteration:\n"
    c += "    pass\n"
    c += "nope = [x for x in nope[:-1] if x not in ld]\n"
    c += "print 'The biggest number (up to {0}) that you cannot buy is: {1}'.format(maxVal, sorted(nope, reverse=True)[0])"
    return c
 
def fixMyPackingProblem(packages, maxVal):
    exec(compile(gen_for(packages, maxVal), "<String>", "exec"))
    
    
#print gen_for((6, 9, 20), 150)
fixMyPackingProblem((6, 9, 20), 150)

This post has been edited by Python_4_President: 01 November 2012 - 09:09 AM

[blackcompe]

I did mine in Scala. Purely functional baby!! Can we get some test cases?? Am I correct in saying the equation must consist of co-prime coefficients only?

Updated code

Spoiler

/**
 * Computes the Frobenius number for a linear equation with co-prime coefficients.
 * A Frobenius number is the biggest integer for which the equation doesn't have an
 * exact solution for.
 */
object Frobenius extends App {

	/**
	 * Compute the greatest common denominator of two integers.
	 *
	 * @param a		The first integer.
	 * @param b		The second integer.
	 * @return 		The gcd value.
	 */
	def gcd(a: Int, b: Int): Int =
		if (b == 0) a
		else gcd(b, a % B)/>

	/**
	 * Check whether a list of integers is co-prime; that is, their greatest
	 * common denominator is 1.
	 *
	 * @param nums	The list of numbers.
	 * @return 		<code>True</code> if the integers are co-prime, otherwise returns
	 * 				<code>false</code>.
	 */
	def isCoprime(nums: List[Int]): Boolean = {
		val allPairs = nums.map(num => nums.map(otherNum => (num, otherNum)))
		val allPairsGCD = allPairs.map(list => list.map(tuple => gcd(tuple._1, tuple._2)))
		allPairsGCD.map(list => list.exists(_ == 1)).exists(_ == true)
	}

	/**
	 * Evaluate a linear equation specified by a list of coefficients and values
	 * for the unknown variables. The coefficients should match up with the values.
	 * E.g. To evaluate the expression <code>1*x + 2*y</code>, where <code>x = 3</code>
	 * and <code>y = 4</code> the call is as follows:
	 *
	 * 				<br/>
	 * 				<code>
	 * 				val coeffs = List(1, 2)
	 * 				val varVals = List(3, 4)
	 * 				val result = eval(coeffs, varVals)
	 * 				</code>
	 * 				<br/>
	 *
	 * @param coeffs	A list of linear coefficients.
	 * @param varVals	The values of the equation variables.
	 */
	def eval(coeffs: List[Int], varVals: List[Int]): Int =
		coeffs.zip(varVals).map(p => p._1 * p._2).sum

	/**
	 * Check whether a solution for a given equation represented as a list of
	 * linear coefficients exists. Note that the coefficients must be co-prime,
	 * i.e. <code>gcd(coeffs) = 1</code>, or this method will throw an
	 * {@link IllegalArgumentException}. To check whether a solution exists
	 * for the expression <code>1*x + 2*y = 2</code>, the call is as follows:
	 *
	 * 					<br/>
	 * 					<code>
	 * 					val coeffs = List(1, 2)
	 * 					val N = 2
	 * 					val exists = solutionExists(coeffs, N)
	 * 					</code>
	 * 					<br/>
	 *
	 * @param coeffs	A list of linear co-prime coefficients.
	 * @param N			The value of the known side of the equation.
	 * @throws IllegalArgumentException
	 * 					If the coefficients aren't co-prime.
	 * @return 			<code>True</code> if a solution was found, otherwise returns <code>false</code>.
	 */
	def solutionExists(coeffs: List[Int], N: Int): Boolean = {

		def _solutionExists(vals: List[Int]): Boolean = {
			val res = eval(coeffs, vals)
			if (res > N) false
			else if (res == N) true
			else {
				val retVals = for (i <- 0 until coeffs.size) yield _solutionExists(vals.updated(i, vals(i) + 1))
				if (retVals.exists(_ == true)) true
				else false
			}
		}
		if (!isCoprime(coeffs)) throw new IllegalArgumentException("Coefficients must be co-prime!")
		else _solutionExists((for (i <- 1 to coeffs.size) yield 0).toList)
	}

	/**
	 * Find the Frobenius number for a given equation represented as a list of
	 * linear coefficients. To find the Frobenius number for the expression
	 * <code>1*x + 2*y</code>, the call is as follows:
	 * goog
	 * 	 				<br/>
	 * 					<code>
	 * 					val coeffs = List(1, 2)
	 * 					val frobNum = frobenius(coeffs)
	 * 					</code>
	 * 					<br/>
	 *
	 * @param coeffs	A list of linear co-prime coefficients.
	 * @return 			The Frobenius number or -1 if the coefficients aren't co-prime.
	 */
	def frobenius(coeffs: List[Int]): (Int, Int) = {

		def solve(N: Int, solns: Int, candidate: Int, iters: Int, solvable: List[Int]): (Int, Int) = {
			if (solvable.contains(N)) solve(N + 1, solns + 1, candidate, iters + 1, (N + coeffs(0)) :: solvable)
			else if (!solutionExists(coeffs, N)) solve(N + 1, 0, N, iters + 1, solvable)
			else if (solns == (coeffs(0) - 1) && candidate != -1) (candidate, iters)
			else solve(N + 1, solns + 1, candidate, iters + 1, (N + coeffs(0)) :: solvable)
		}
		if (!isCoprime(coeffs)) (-1, 0)
		else solve(0, 0, -1, 0, List())
	}

	//Valid
	println(frobenius(List(6, 9, 20)))
	println(frobenius(List(9, 20)))
	println(frobenius(List(12, 8, 23)))

	println()

	//Invalid
	println(frobenius(List(2, 2)))
	println(frobenius(List(2, 4)))
	println(frobenius(List(24, 48, 100)))
}

This post has been edited by blackcompe: 02 November 2012 - 06:51 PM

[Simown]

They are not necessarily co-prime, but if they aren't there is no solution. So return 0, -1, or "none" or some indicator of that.

Test cases are a good idea. I will have to update the first post later, so check back!

[blackcompe]

Has anyone been able to get the correct output for 9x + 20y = N? It's 151, but using the algorithm given in the OCW handout gives 115. Below is the debugging output from my code. I get 115 as having no solutions and then immediately after I get 6 consecutive numbers with solutions. So 115, according to MIT is correct.

115:false
116:true
117:true
118:true
119:true
120:true
121:true
(115,120)

It appears that this algorithm doesn't generalize at all. The only problem that's correctly solved is this the particular case of 6, 9, and 20. Heh, maybe that's what was intended. I also found a calculator for the Frobenius problem (McNuggets generalized).

This post has been edited by blackcompe: 02 November 2012 - 03:19 PM

[jon.kiparsky]

Yeah, 151 is what I get.

If you're talking about the six consecutive numbers, the generalization for that is, taking x to be the least term in the set, if you have x consecutive numbers then any number above that point must be a solution. Easy to prove: if p is a solution, then p+x is a solution. If xa + yb +c = p, then x(a+1) +yb+zc = p+x.

So in this case you need to look for 9 consecutive solutions.

[blackcompe]

Quote

If you're talking about the six consecutive numbers, the generalization for that is, taking x to be the least term in the set, if you have x consecutive numbers then any number above that point must be a solution. Easy to prove: if p is a solution, then p+x is a solution. If xa + yb +c = p, then x(a+1) +yb+zc = p+x.

Ah thanks Jon!

Edit:

I got it working.

This post has been edited by blackcompe: 02 November 2012 - 03:39 PM

[jon.kiparsky]

So I think we've had plenty of time for people to come up with solutions. Here's the one that I came up with, using a variation on the idea of the sieve: we don't need to test numbers which are congruent modulo x to a known solution. This limits the search space somewhat, resulting in improved running times.

There's an obvious improvement to make, of course - testing for GCD at the outset would tidy things up a little. However, the question of what number constitutes a reasonable upper bound is an interesting one, so I wanted to leave it in place while I think it over.

Spoiler

#!/usr/bin/python


# Find maximal non-solution integer for diophantine equation xa + yb + zc +...  = n using 
# a sieve. 
# Coeffiecients are taken as a tuple, though the first thing we do with them is turn them into a non-repeating sorted list. 


def diophantine_sieve (coeffs):
	coeffs = list(set(coeffs))
	coeffs.sort()
	hugenum = coeffs[-1] * 10000  	## to prevent overruns, for example coeffs=(3,9,27)
									## To do: what is the correct value for hugenum? 
	x= coeffs[0]				
	coeffs = coeffs[1:]				## don't need to check x, since it's built in to the solution
	try_list = [0] 					## start by assuming zero is a solution
	solution_grid = list([hugenum]*x) 			
	solution_grid[0] = x			## solution grid[n] contains the minimal non-zero solution p 
									## where p%x==n
	outer_count = 0	
	inner_count = 0	
	while len(try_list)>0:			## try_list is the list of potential solutions which have yet to be tested
		outer_count = outer_count +1 
		(c, try_list) = try_list[0], try_list[1:]	## remove an item to test 
		for i in coeffs:
			inner_count = inner_count +1 			
			row = (c +i ) % x				
			if solution_grid[row] > c + i:	## if this is the lowest solution for this row of the grid
				solution_grid[row] = c +i	## add it to the grid
				try_list.append(c + i)		## and put it in the list of values to try	
							
										## implicit else: if not, we don't need to try any further with this one
										## because we've been here already: if p1 is congruent mod x to p, then 
										## p1 + y is congruent mod x to p + y. 

	max_s = max(solution_grid)			## max_s is the "largest least" solution: solution_grid[n] contains the least solution p
										## such that p%x == n. Of these p1,p2,...pn, we want the largest. 
	if max_s == hugenum:				## if max_s == hugenum, then an overrun occurred: coefficients were co-prime
		return 
	else:
		return  max_s - x, inner_count 		## otherwise, the number we are looking for is the largest solution 
											## we've found, call it max_s, less x. 
											## Proof: If max_s - x were a solution, it would be in the solution grid, 
											## therefore it is not a solution. 
											## If any number t, larger than max_s - x, were a non-solution, 
											## max_s would not be the largest number in the
											## grid since there would necessarily be a larger solution, 
											## or else our coefficients were co-prime
											
											## Therefore max_s -x is the largest integer which cannot be a solution for the 
											## equation described. 	
											## "Neat, huh?"


print "Some test values (solution and number of times through the inner loop)"
print "6, 9, 20:" +				str(diophantine_sieve((6,9, 20)))
print "691, 935, 2137: "+ 		str(diophantine_sieve((691,935,2137)))
print "6913, 9353, 21374" + 	str(diophantine_sieve((6913,9353,21374)))
print "69135, 94712, 213542" +	str(diophantine_sieve((69135, 94712, 213542)))  
print "123, 83, 233" +			str(diophantine_sieve((123, 83, 233)))
print "9, 20:" +				str(diophantine_sieve((9, 20)))

[blackcompe]

Cutting down the search space would definitely help my algorithm's efficiency. As of now, it's taking quite a while for even small coefficients.

Quote

we don't need to test numbers which are congruent modulo x to a known solution.

Can you expand on this? Perhaps give an example? When you say known solution, are you referring to a number that we find can (or can't) find a solution for? For instance, in the McNugget problem, the numbers that don't have solutions are:

1,2,3,4,5,7,8,10,11,13,14,16,17,19,22,23,25,28,31,34,37,43

So are the known solutions the complement of that set?

Secondly, what's x referring to here when you say "modulo x". Is it the a coefficient in the equation (e.g. 6 for the McNugget problem) or is it the value of the unknown x in the equation for that particular solution.

E.g. In finding a solution for N = 42 in the McNugget problem, we have:

6 * x + 9 * y + 20 * z = 42

A solution is [x = 7, y = 0, z = 0]

So the x you're referring to is then 7. So we'd have something like this:

42 = test_number (mod 7)

Obviously, there's no test_number that satisfies that congruence, so I know my logic is incorrect, but is this close to what you mean?

[jon.kiparsky]

Quote

When you say known solution, are you referring to a number that we find can (or can't) find a solution for

By "a known solution" I mean some number which is known to be a possible solution to that equation. So 42, as you say, is known to be a possible solution (and therefore 49 is as well)
The set of known solutions would be, as you say, the complement of the set of numbers which cannot be solutions.

Quote

Secondly, what's x referring to here when you say "modulo x". Is it the a coefficient in the equation (e.g. 6 for the McNugget problem) or is it the value of the unknown x in the equation for that particular solution.

Yes, I mean the "a" coefficient. Sorry, that was unclear. Particularly, I mean the least of the coefficients. So perhaps I should have said "a". That would have made more sense.

So from your example, if we know that 42 is a legitimate "move" in this game, then we can infer that 48, as well as 51, and 62, are legal moves. This is a simple rule of inference which in a way ignores the math of it, but it turns out to be correct, as we saw above.
In fact, we could reduce this to a simple nim-type "game", though not a very entertaining one. The rules are, you start from zero, and you can derive any number which you can reach by adding some member of the coefficient set. If you imagine a grid of "x" (a, really! ) rows, numbered 0..x-1 and filled in with the numbers in the obvious pattern (0..x-1, x..2x-1, etc), extending arbitrarily far to the right, you can imagine the goal being to "claim" rows by deriving the lowest possible number in that row, according to the simple rule of derivation.
My algorithm, basically, plays this game as a solitaire.

That game, by the way, is an example, or perhaps an illustration, of my assertion that "we don't need to test numbers which are congruent modulo x to a known solution." it should be obvious that once we reach any square on that grid, all of the squares to the right of it are known to be solutions. Therefore, finding the least value in each row, and then taking the largest of that set, is sufficient to solve the problem.

This post has been edited by jon.kiparsky: 06 November 2012 - 12:01 PM

[blackcompe]

Quote

Yes, I mean the "a" coefficient. Sorry, that was unclear. Particularly, I mean the least of the coefficients. So perhaps I should have said "a". That would have made more sense.

So from your example, if we know that 42 is a legitimate "move" in this game, then we can infer that 48, as well as 51, and 62, are legal moves.

I'm still unclear. Can you show the Math behind finding that 48, 51, and 62 are solutions. From what I understand, your saying (for the McNugget problem, where a = 6):

0 = 48 mod 6
3 = 51 mod 6
2 = 62 mod 6

0, 3, 2 are previously found known solutions, but for

1 = 43 mod 6

1 is a known solution, and therefore 43 is a known solution, which is obviously not true.

[jon.kiparsky]

1 is not a known solution - there's no way to get 1 to be a solution of 6x + 9y +20z. The rule of inference proceeds from zero. Zero can be a solution to any equation in this form, simply by setting x, y,z, ... to zero. And term in the list of coefficients can be a solution, simply by setting the appropriate x, y, z, ... to 1, and leaving the rest at zero. (That is, by adding that coeffeicient to zero.)


Okay, I admit that I've had a few beers at this point (happy re-election to all!) but there's three steps as I see it.

First of all, zero is a solution to any equation of this form. (see above)

Second, if the coefficients of the equation are listed as (a, b, c...) then any previously derived solution, increased by some member of that set, is also a solution.

Third, if we have a solution in row of this notional grid, described above, then we know that there is a point above which all integers are solutions to this equation.

I think each of these is proveable, and I think that

Fourth, taken together they establish that the algorithm is sound.

I would be happy to try to walk through any of these tomorrow evening, when I'm not at work and when I'm not just back from the bar. Let me know which you'd like.

[macosxnerd101]

I've got a basic solution cranked out. I optimized to ignore coefficients that are greater than the sum, and to check right off if the remaining capacity (mod coefficient[i]) == 0. I also started the lowest non-solution as the lowest coefficient-1. Beyond that, I basically took the Sieve approach.

I did play around with upper bounds. I tried to guestimate that since the solution is coprime to all the coefficients, and the Sieve of Eratosthenes is O(n^1.5), then the upper bounds for the largest non-solution is (sum of coefficients)^1.5. Unfortunately, that didn't scale. The lower bounds I found was sqrt(3 * product(coefficients)) - (Sum of coefficients). My lower bounds was more efficient. Given the equation 77x + 42y = c, sqrt(3*77*42) - (77+42) < 0, and my lower bounds is 41 there.

Spoiler

package test;
import java.util.*;

/**
 * @author Michael Levet
 * @date November 9, 2012
 ***/
class Diophantine{

    /**
     * 
     * @param coefficients The coefficients for the equation
     * @param capacity The target number
     * @return boolean true if the coefficients can be made to add to capacity, false otherwise
     ***/
    public static boolean isDiophantine(TreeSet<Integer> coefficients, int capacity){
        
        //if we are done, return true
        if(capacity == 0){
            return true;
        }
        
        Iterator<Integer> iterator = coefficients.iterator();
        
        while(iterator.hasNext()){
            int temp = iterator.next();
            
            //if a coefficient is a factor of capacity, then this equation is a solution
            if(capacity % temp == 0){
                return true;
            } 
            
            //if the coefficient is larger than capacity, prune it
            if(temp > capacity){
                break;
            }
        }
        
        //and recurse for the remaining coefficients
        TreeSet<Integer> remainingCoefficients = new TreeSet<Integer>(coefficients);
        
        for(int i:coefficients){
            if(i > capacity){
                break;
            }
            
            if(isDiophantine(remainingCoefficients, capacity - i)){
                return true;
            }
        }
        
        return false;
    }
    
    /**
     * @param coefficients The coefficients for the equation
     * @return int The largest number that doesn't satisfy the equation
     ***/
    public static int largestDiophantine(TreeSet<Integer> coefficients){
        int sum = 1;
        
        //add up the factors
        for(int i:coefficients){
            sum += i;
        }

        int lowest = coefficients.iterator().next();
        
        int k = 0;
        int lastNonSolution = lowest-1;
        
        ArrayList<Integer> solutions = new ArrayList<Integer>();
        
        while(k < lowest){
            
            if(zeroModSolutions(sum, solutions)){
                k++;
                solutions.add(sum);
            }
            
            else if(isDiophantine(coefficients, sum)){
                k++;
                solutions.add(sum);
            }
            
            else{
                lastNonSolution = sum;
                k = 0;
            }
            
            sum++;
        
        }
        
        return lastNonSolution;
    }
    
    public static boolean zeroModSolutions(int sum, ArrayList<Integer> solutions){
        
         for(int i:solutions){
             if(sum % i == 0){
                 return true;
             }
         }
         
         return false;
    }
    
    public static void main(String[] args){
        TreeSet<Integer> coefficients = new TreeSet<Integer>();
        coefficients.add(6);
        coefficients.add(9);
        coefficients.add(20);
       
        System.out.println("Largest number: " + largestDiophantine(coefficients));
    }
}

This post has been edited by Simown: 10 November 2012 - 07:14 PM
Reason for edit:: Code tags?!