(* :Title: Binomial Coefficients *)

(* :Context: BinomialCoefficients` *)

(* :Author: Eric Rowland *)

(* :Date: {2010, 3, 8} *)

(* :Package Version: 1.10 *)

(* :Mathematica Version: 7.0 *)

(* :Discussion:
	BinomialCoefficients is a package containing implementations
	of several theorems on the residues and divisibility properties of binomial coefficients.
	Documentation can be found at  http://math.tulane.edu/~erowland/packages.html .
*)


BeginPackage["BinomialCoefficients`", {"Combinatorica`", "SequenceIdentification`"}]

{
	BinomialComplementMod,
	BinomialExponent,
	BinomialMod,
	BinomialResidueCount,
	Carries,
	CoprimeFactorial
}


Unprotect["BinomialCoefficients`*"]

Begin["`Private`"]

var[variables___String] :=
	StringJoin[Riffle[("\*StyleBox[\"" <> # <> "\", \"TI\"]" &) /@ {variables}, ", "]]

sub[variable_String, indices : {___}] :=
StringJoin[Riffle[
	Replace[indices, {
		"..." :> var["\[Ellipsis]"],
		Null :> var[variable],
		index_String :> "\!\(\*SubscriptBox[StyleBox[\"" <> variable <> "\", \"TI\"], StyleBox[\"" <> index <> "\", \"TI\"]]\)",
		index_ :> "\!\(\*SubscriptBox[StyleBox[\"" <> variable <> "\", \"TI\"], StyleBox[\"" <> ToString[index] <> "\", \"TR\"]]\)"
	},
	{1}
], ", "]]
sub["...", _] := var["..."]
sub[variable_String, index_] := sub[variable, {index}]
sub[variables_List, index : Except[_List]] := sub @@ ({#, index} &) /@ variables
sub[variables : {_String, _} ...] := StringJoin[Riffle[sub @@@ {variables}, ", "]]


BinomialComplementMod::usage =
"BinomialComplementMod[" <> var["n", "m", "p"] <> "] gives the residue of Binomial[" <> var["n", "m"] <> "] modulo " <> var["p"] <> " after dividing by the largest power of " <> var["p"] <> " that divides Binomial[" <> var["n", "m"] <> "].
BinomialComplementMod[" <> var["n", "m", "p"] <> "^" <> var["\[Alpha]"] <> "] gives the residue of Binomial[" <> var["n", "m"] <> "] modulo " <> var["p"] <> "^" <> var["\[Alpha]"] <> " after dividing by the largest power of " <> var["p"] <> " that divides Binomial[" <> var["n", "m"] <> "]."

BinomialExponent::usage =
"BinomialExponent[" <> var["n", "m", "k"] <> "] gives the exponent of the highest power of " <> var["k"] <> " that divides Binomial[" <> var["n", "m"] <> "]."

BinomialMod::usage =
"BinomialMod[" <> var["n", "m", "k"] <> "] gives the residue of Binomial[" <> var["n", "m"] <> "] modulo " <> var["k"] <> "."

BinomialResidueCount::usage =
"BinomialResidueCount[" <> var["n", "k", "r"] <> "] gives the number of binomial coefficients Binomial[" <> var["n", "m"] <> "] that are congruent to a nonzero residue " <> var["r"] <> " modulo " <> var["k"] <> ".
BinomialResidueCount[" <> var["n", "k"] <> "] gives the number of binomial coefficients Binomial[" <> var["n", "m"] <> "] that are nonzero modulo " <> var["k"] <> "."

Carries::usage =
"Carries[" <> var["m", "r", "b"] <> "] gives the number of carries involved in adding " <> var["m"] <> " and " <> var["r"] <> " in base " <> var["b"] <> ".
Carries[" <> var["m", "r", "b", "j"] <> "] gives the number of carries onto or beyond the " <> var["b"] <> "^" <> var["j"] <> " place."

CoprimeFactorial::usage =
"CoprimeFactorial[" <> var["n", "p"] <> "] computes the product of all natural numbers less than or equal to " <> var["n"] <> " that are not divisible by a prime " <> var["p"] <> "."


PositivePrimeQ[n_] := PrimeQ[n] && Positive[n]

PositivePrimePowerQ[n_] := PrimePowerQ[n] && Positive[n]

SetAttributes[Case, HoldRest]
Case[list : _[___], pattern_, alt_ : Null, n : (_Integer?Positive) : 1] :=
If[Length[#] >= n,
	Last[#],
	alt
] & @ Cases[list, pattern, {1}, n]

SetAttributes[Occurrence, HoldRest]
Occurrence[list : _[___], p_, alt_ : Null, n : (_Integer?Positive) : 1] :=
If[Length[#] >= n,
	#[[-1,1]],
	alt
] & @ Position[list, p, {1}, n, Heads -> False]

coefficient[p_, w_List] :=
	(p - First[w] - 1) Times @@ (p - Take[w, {2, -2}]) Last[w] / Times @@ (w + 1)

DisjointIntervalListQ[intervallist : {{_, _} ...}] :=
MatchQ[
	Apply[IntervalIntersection, Subsets[Interval /@ intervallist, {2}], {1}],
	{Interval[] ...}
]

NonOverlappingSubwordSum[f_, wordlengths : {___Integer?Positive}, word_List] :=
Total[(f @@ (word[[Span @@ #]] &) /@ # &) /@
	Select[
		Join @@@ Tuples[
			Replace[
				Tally[wordlengths],
				{length_, count_} :> Select[
					Subsets[Table[{i, i + length - 1}, {i, Length[word] - length + 1}], {count}],
					DisjointIntervalListQ
				],
				{1}
			]
		],
		DisjointIntervalListQ
	]
]

MyConnectedComponents[pairs_] :=
	(pairs[[#]] &) /@
		ConnectedComponents[FromAdjacencyMatrix[Outer[Boole[Intersection[##] != {}] &, pairs, pairs, 1]]]

FindPath[edges : {{_, _} ..}] :=
Module[{edge, orderededges = {}, remainingedges = edges},
	edge = Case[Tally[Join @@ edges], {v_, 1} :> Case[edges, {v, _} | {_, v}], First[edges]];
	AppendTo[orderededges, edge];
	remainingedges = DeleteCases[remainingedges, edge, {1}, 1];
	While[remainingedges != {},
		edge = Case[remainingedges, _?(Intersection[#, Join @@ orderededges] != {} &), $Failed];
		AppendTo[orderededges, edge];
		remainingedges = DeleteCases[remainingedges, edge, {1}, 1]
	];
	orderededges
]

RestrictedClusters[subwordlengths_, requiredoverlaps_] :=
With[{wordindices = DeleteDuplicates[Join @@ FindPath[requiredoverlaps]]},
	(# + 1 - Min[#] &) /@ Fold[
		Function[{clumplist, newwordindex},
			Join @@ Function[clump,
				Function[i, Append[clump, Range[i, i + subwordlengths[[newwordindex]] - 1]]] /@ Replace[
					IntervalIntersection @@ Cases[
						requiredoverlaps,
						{newwordindex, overlappingwordindex_} | {overlappingwordindex_, newwordindex} /;
								Occurrence[wordindices, overlappingwordindex] < Occurrence[wordindices, newwordindex] :>
							Interval[
								{1 - subwordlengths[[newwordindex]], subwordlengths[[overlappingwordindex]] - 1} +
									First[clump[[Occurrence[wordindices, overlappingwordindex]]]]
							]
					],
					{Interval[interval_] :> Range @@ interval, Interval[] -> {}}
				]
			] /@ clumplist
		],
		{{Range[subwordlengths[[First[wordindices]]]]}},
		Rest[wordindices]
	]
]


SetAttributes[BinomialComplementMod, Listable]
BinomialComplementMod[n_Integer?NonNegative, m_Integer?NonNegative, p_?PositivePrimeQ] /; m <= n :=
With[{k = BinomialExponent[n, m, p], r = n - m},
	Mod[
		(-1)^k Times @@ 
			(#1! PowerMod[#2! #3!, -1, p] &) @@@
				Transpose[PadLeft[IntegerDigits[{n, m, r}, p]]],
		p
	]
]
BinomialComplementMod[n_Integer?NonNegative, m_Integer?NonNegative, k_?PositivePrimePowerQ] /; m <= n :=
Module[{p, q, r = n - m},
	{p, q} = NumberTheory`PrimePower[k];
	Mod[
		If[p == 2 && q >= 3, 1, (-1)^Carries[m, r, p, q-1]]
		*
		Times @@
			(CoprimeFactorial[#1, p] PowerMod[CoprimeFactorial[#2, p] CoprimeFactorial[#3, p], -1, k] &) @@@
				Transpose[PadLeft[
					Mod[Floor[#/p^Range[IntegerLength[#, p], 0, -1]], k] & /@
						{n, m, r}
				]],
		k
	]
]
BinomialComplementMod[n_Integer?NonNegative, m_Integer?NonNegative, _?PositivePrimePowerQ] /; m > n := 0
BinomialComplementMod[_Integer?NonNegative, _Integer?Negative, _?PositivePrimePowerQ] := 0
SyntaxInformation[BinomialComplementMod] = {"ArgumentsPattern" -> {_, _, _}}

SetAttributes[BinomialExponent, Listable]
BinomialExponent[n_Integer?NonNegative, m_Integer?NonNegative, p_?PositivePrimeQ] /; m <= n :=
	Carries[m, n - m, p]
BinomialExponent[n_Integer?NonNegative, m_Integer?NonNegative, k_Integer /; k >= 2] /; m <= n :=
With[{primepowers = FactorInteger[k]},
	Floor[Min[BinomialExponent[n, m, First /@ primepowers] / (Last /@ primepowers)]]
]
BinomialExponent[_Integer?NonNegative, _Integer, k_Integer /; k >= 2] := Infinity
SyntaxInformation[BinomialExponent] = {"ArgumentsPattern" -> {_, _, _}}

SetAttributes[BinomialMod, Listable]
BinomialMod[n_Integer?NonNegative, m_Integer?NonNegative, p_?PositivePrimeQ] :=
	Mod[Times @@ Binomial @@@ Transpose[PadLeft[IntegerDigits[{n, m}, p]]], p]
BinomialMod[n_Integer?NonNegative, m_Integer?NonNegative, k_?PositivePrimePowerQ] /; m <= n :=
With[{p = First[NumberTheory`PrimePower[k]]},
	Mod[
		p^BinomialExponent[n, m, p] BinomialComplementMod[n, m, k],
		k
	]
]
BinomialMod[n_Integer?NonNegative, m_Integer?NonNegative, k_Integer /; k >= 2] :=
With[{primepowers = Power @@@ FactorInteger[k]},
	ChineseRemainder[BinomialMod[n, m, primepowers], primepowers]
]
BinomialMod[_Integer?NonNegative, _Integer?Negative, k_Integer /; k >= 2] := 0
SyntaxInformation[BinomialMod] = {"ArgumentsPattern" -> {_, _, _}}

SetAttributes[BinomialResidueCount, Listable]
BinomialResidueCount[n_Integer?NonNegative, p_?PositivePrimeQ, r_Integer] /; !Divisible[r, p] :=
Module[{permutation, x},
	permutation = Ordering[PowerMod[PrimitiveRoot[p], Range[0, p - 2], p]] - 1;
	Coefficient[PolynomialMod[
		Product[
			Total[#2 x^permutation[[#1]] & @@@
				Tally[DeleteCases[Mod[Binomial[j, Range[0, j]], p], 0]]
			]^DigitCount[n, p, j],
			{j, p - 1}
		],
		x^(p - 1) - 1
	], x, permutation[[Mod[r, p]]]]
]
BinomialResidueCount[n_Integer?NonNegative, k_Integer /; k >= 2, r_Integer] /; 1 <= r <= k - 1 :=
	2 Count[BinomialMod[n, Range[0, (n - 1)/2], k], r] + Boole[EvenQ[n] && BinomialMod[n, n/2, k] == r]
BinomialResidueCount[n_Integer?NonNegative, 1] := 0
BinomialResidueCount[n_Integer?NonNegative, k_?PositivePrimePowerQ] :=
Module[{p, alpha, digits},
	{p, alpha} = NumberTheory`PrimePower[k];
	digits = DigitList[n, p];
	Times @@ (digits + 1)
	*
	Sum[
		Total[Function[partition,
			NonOverlappingSubwordSum[Times @@ (coefficient[p, #] &) /@ {##} &, partition, digits]
		] /@ IntegerPartitions[gamma, {Max[0, gamma - (alpha - 1)], Infinity}, Range[2, gamma]]],
		{gamma, 0, 2 (alpha - 1)}
	]
]
BinomialResidueCount[n_Integer?NonNegative, k_Integer /; k >= 2] :=
	2 Count[BinomialMod[n, Range[0, (n - 1)/2], k], Except[0]] + Boole[EvenQ[n] && BinomialMod[n, n/2, k] != 0]
BinomialResidueCount[n_, k_?PositivePrimePowerQ] :=
Module[{p, alpha},
	{p, alpha} = NumberTheory`PrimePower[k];
	Product[(w + 1)^DigitsCount[n, p, w], {w, 0, p - 1}]
	*
	Expand[Sum[
		Total[Function[partition,
			Total[Function[pairs,
				With[{connectedcomponents = MyConnectedComponents[pairs]},
					(-1)^Length[pairs]
					*
					Times @@ (Total[Function[cluster,
							Total[Apply[
								Evaluate[
									DigitsCount[n, p, {##}]
									*
									Times @@ (coefficient[p, #] &) /@ (cluster /. i_Integer :> Slot[i])
								] &,
								Tuples[Range[0, p - 1], Max[cluster]],
								{1}
							]]
						] /@ #] &) /@
							Join[
								RestrictedClusters[partition, #] & /@ connectedcomponents,
								{{Range[#]}} & /@
									partition[[Complement[Range[Length[partition]], Flatten[connectedcomponents]]]]
							]
				]
			] /@
				Subsets[Subsets[Range[Length[partition]], {2}]]
			]
			/
			Times @@ (Last /@ Tally[partition]!)
		] /@ IntegerPartitions[gamma, {Max[0, gamma - (alpha - 1)], Infinity}, Range[2, gamma]]],
		{gamma, 0, 2 (alpha - 1)}
	]]
]
SyntaxInformation[BinomialResidueCount] = {"ArgumentsPattern" -> {_, _, _.}}

SetAttributes[Carries, Listable]
Carries[m_Integer?NonNegative, r_Integer?NonNegative, b_Integer /; b >= 2, j : (_Integer?NonNegative) : 0] :=
Module[{i = 0, count = 0, carry = 0},
	(
		If[# + carry < b,
			carry = 0
			,
			carry = 1;
			If[i >= j, count++]
		];
		i++
	) & /@ Reverse[Prepend[Total[PadLeft[IntegerDigits[{m, r}, b]]], 0]];
	count
]
SyntaxInformation[Carries] = {"ArgumentsPattern" -> {_, _, _, _.}}

SetAttributes[CoprimeFactorial, Listable]
CoprimeFactorial[n_Integer?NonNegative, p_?PositivePrimeQ] :=
	n!/(Floor[n/p]! p^Floor[n/p])
SyntaxInformation[CoprimeFactorial] = {"ArgumentsPattern" -> {_, _}}


End[]

Protect["BinomialCoefficients`*"]

EndPackage[]