// approximsum -- approximate number as sum of integer multiples of two others
// David Eppstein, UC Irvine, 18 Jun 1996
//
// given a,b,c find x,y >= 0 s.t. ax+by is as close as possible to c

#include "approximsum.h"

static void approxbelow(long a, long b, long c, long & x, long & y, long & d);

// substitute possible better soln
static inline void testsoln(long & x, long & y, long & d,
		       				long xx, long yy, long dd)
{
	if (dd < d) {
		x = xx;
		y = yy;
		d = dd;
	}
}

// main entry point
void approximsum(long a, long b, long c, long & x, long & y)
{
	long d = c;					// recursion wants an extra argument
	x = 0;
	y = 0;
	approxbelow(a, b, c, x, y, d);
}

// main recursion
void approxbelow(long a, long b, long c, long & x, long & y, long & d)
{
	if (a < b) approxbelow(b, a, c, y, x, d);
	else if (c < 0) {
		// base case: we want to approximate a negative number
		// the best approximation is always the origin.
		//
		testsoln(x, y, d, 0, 0, -c);

	} else if (a <= 0) {
		// another bad case: stuck at the origin
		//
		testsoln(x, y, d, 0, 0, c);

	} else {
		// main case.  divide points into two regions qx+y <= r and qx+y >= r+1
		// where r is max value s.t. whole first region satisfies ax+by <= c
		long q = a/b;
		long r = (q*c)/a;

		// solution 1: qx + y <= r
		// find qx+y=r, y as small as possible
		long xx = c/a;
		long yy = r - q*xx;
		testsoln(x, y, d, xx, yy, c - a*xx - b*yy);

		// solution 2: qx + y >= r + 1
		//
		// replace with new vars xx, yy
		// xx = x, yy = y + qx - (r + 1)
		// to make the inequality qx+y >= r+1 turn into yy >= 0
		//
		// so y = yy - qx + (r + 1)
		// and (ax+by) = (a xx + b yy - bq xx + b(r+1)
		c -= b*(r+1);
		y += q*x - (r+1);
		approxbelow(b, a - b*q, c, y, x, d);
		y -= q*x - (r+1);
	}
}