Shapley-Shubik Power Index Calculator

shapley-shubik.cc

#include <iostream>
#include <sstream>
#include <vector>
#include <algorithm>
using namespace std;
using Num = long long;
using Array = vector<Num>;

Num math_combination(Num total, Num select)
{
	if (select > total / 2) select = total - select;

	Num ret = 1;
	for (Num i = 0; i < select; ++i) ret *= (total - i);
	for (Num i = 2; i <= select; ++i) ret /= i;
	return ret;
}

Num math_power(Num n)
{
	if (n == 0) return 1;
	Num ret = n;
	for (Num i = 2; i < n; ++i) ret *= i;
	return ret;
}

/*
 * Solve by generating all permutation and check the key element.
 * Time Complexity: O(n!)
 * 5.4h / ops
 */
Array solve_perm(const Array &weight, int threshold)
{
	vector<int> perm;
	for (Array::size_type p = 0; p < weight.size(); ++p) {
		perm.push_back(perm.size());
	}

	Array result(weight.size());
	do {
		Num sum = 0;
		for (const auto i : perm) {
			auto v = weight[i];
			sum += v;
			if (sum >= threshold) {
				++result[i];
				break;
			}
		}
		static int i;
	} while (next_permutation(perm.begin(), perm.end()));

	return result;
}

/*
 * Solve by generating all combination and infer the key time for each element.
 * Time Complexity: O(n * 2 ^ n)
 * 400 ops/s
 */
Array solve_comb(const Array &weight, int threshold)
{
	using Bit = unsigned char;
	vector<Bit> choice(weight.size());
	vector<Num> sums(weight.size());

	for (;;) {
		Bit c = 1; // carry
		for (auto &i : choice) {
			Bit pre_i = i;
			Bit pre_c = c;
			i = pre_i ^ pre_c;
			c = pre_i & pre_c;
		}
		if (c) break; // overflow

		int sum = 0;
		int num = 0;
		for (Array::size_type i = 0; i < weight.size(); ++i) {
			if (choice[i]) {
				sum += weight[i];
				++num;
			}
		}

		if (sum < threshold) continue;

		Num perm_num = math_power(num - 1)
			* math_power(weight.size() - num);
		for (Array::size_type i = 0; i < weight.size(); ++i) {
			if (!choice[i]) continue;
			if (sum - weight[i] >= threshold) continue;
			sums[i] += perm_num;
		}
	}
	return sums;
}

/*
 * Solve by generating all combination and infer the key time for each element.
 * Optimize by combining the same weights.
 * Time Complexity: O(sum(k) ^ 2 * prod(k + 1)!)
 * 32k ops/s
 */
Array solve_comb_optim(const Array &weight, int threshold)
{
	struct Weight {
		Num val;
		int num;
		Weight(Num val, int num): val(val), num(num) {};
	};

	// Combine the same weights
	vector<Weight> weights;
	using WeightsSize = vector<Weight>::size_type;
	Num max_weight = *max_element(weight.begin(), weight.end());
	for (Num i = 1; i <= max_weight; ++i) {
		Num num = count(weight.begin(), weight.end(), i);
		if (num) weights.emplace_back(i, num);
	}

	Array comb(weights.size()); // number of each weight are selected
	Array sums(weights.size()); // number of key role for each weight
	for (;;) {
		// Calculate next combination
		bool carry = true;
		for (WeightsSize i = 0; i < weights.size(); ++i) {
			comb[i] += carry;
			if (comb[i] > weights[i].num) {
				carry = true;
				comb[i] = 0;
			} else {
				carry = false;
			}
		}
		if (carry) break;

		// Calculate number of weights selected and sum of the weights
		int sum = 0;
		int num = 0;
		for (WeightsSize i = 0; i < weights.size(); ++i) {
			num += comb[i];
			sum += comb[i] * weights[i].val;
		}
		if (sum < threshold) continue;

		// For each choice of equivalent weights
		// there are perm_num possible ways
		Num perm_num = math_power(num - 1) * math_power(weight.size() - num);

		for (WeightsSize i = 0; i < weights.size(); ++i) {
			if (!comb[i]) continue;
			if (sum - weights[i].val >= threshold) continue;

			// Calculate the number of equivalent weight schemes
			Num tmp_perm_num = perm_num;
			for (WeightsSize j = 0; j < weights.size(); ++j) {
				// When the item being counted now is chosen,
				// the available weights and the needed weights
				// both decreases by 1.
				tmp_perm_num *= math_combination(
					weights[j].num - (i == j),
					comb[j] - (i == j));
			}
			sums[i] += tmp_perm_num;
		}
	}

	// Convert the representation method back
	Array result(weight.size());
	for (WeightsSize i = 0; i < weights.size(); ++i) {
		for (Array::size_type j = 0; j < weight.size(); ++j) {
			if (weight[j] == weights[i].val) result[j] = sums[i];
		}
	}
	return result;
}

int generate_votes(Array &votes, int mode)
{
	switch (mode) {
		case 1:
			// Euro
			for (int i = 0; i < 4; ++i) votes.push_back(10);
			for (int i = 0; i < 1; ++i) votes.push_back(8);
			for (int i = 0; i < 4; ++i) votes.push_back(5);
			for (int i = 0; i < 2; ++i) votes.push_back(4);
			for (int i = 0; i < 3; ++i) votes.push_back(3);
			for (int i = 0; i < 1; ++i) votes.push_back(2);
			return 62;
		case 2:
			// Test 1
			votes.push_back(8);
			votes.push_back(8);
			votes.push_back(4);
			votes.push_back(2);
			return 20;
		case 3:
			// United Nations
			for (int i = 0; i < 5; ++i) votes.push_back(7);
			for (int i = 0; i < 10; ++i) votes.push_back(1);
			return 37;
		case 4:
			// Australia
			for (int i = 0; i < 6; ++i) votes.push_back(1);
			for (int i = 0; i < 1; ++i) votes.push_back(3);
			return 5;
		default:
			return -1;
	}
}

int main(int argc, char **argv)
{
	Array votes;
	int threshold = generate_votes(votes, 1);

	int n = 1;
	if (argc == 2) {
		istringstream s(argv[1]);
		s >> n;
	}

	Array ret;
	for (int i = 0; i < n; ++i) ret = solve_comb_optim(votes, threshold);
	Num total = math_power(votes.size());

	for (Array::size_type i = 0; i < ret.size(); ++i) {
		cout << "[" << i + 1 << "] ";
		cout << "#=" << votes[i] << " ";
		cout << "%=" << ret[i] * 100.0 / total << " ";
		cout << "count=" << ret[i] << " ";
		/*
		cout << ret[i] << " ";
		cout << corr[i] << " ";
		cout << corr[i] - ret[i];
		*/
		cout <<endl;
	}

	return 0;
}