(* :Title: Regular Sequences *)

(* :Context: RegularSequences` *)

(* :Author: Eric Rowland *)

(* :Date: {2010, 4, 11} *)

(* :Package Version: 1.05 *)

(* :Mathematica Version: 7.0 *)

(* :Discussion:
	RegularSequences is a package for identifying and analyzing automatic and regular sequences.
	Documentation can be found at  http://math.tulane.edu/~erowland/packages.html .
*)


BeginPackage["RegularSequences`"]

{
	FindRegularSequenceFunction,
	FindRegularSequenceRecurrence,
	RegularSequence,
	RegularSequenceExponent,
	RegularSequenceGeneratorTable,
	RegularSequenceMatrixForm,
	RegularSequenceRank,
	RegularSequenceRecurrence,
	ThueMorse
}


Unprotect["RegularSequences`*"]

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


FindRegularSequenceFunction::usage =
"FindRegularSequenceFunction[{" <> sub["a", {0, 1, "..."}] <> "}, " <> var["k"] <> "] attempts to find a simple RegularSequence function that yields the sequence " <> sub["a", "n"] <> " when given successive integer arguments.
FindRegularSequenceFunction[" <> var["list", "k", "n"] <> "] gives the function applied to " <> var["n"] <> "."

FindRegularSequenceRecurrence::usage =
"FindRegularSequenceRecurrence[{" <> sub["a", {0, 1, "..."}] <> "}, " <> var["k"] <> ", " <> var["a"] <> "[" <> var["n"] <> "]] attempts to find linear relations between subsequences of " <> var["a"] <> "[" <> var["n"] <> "] of the form " <> var["a"] <> "[" <> var["k"] <> "^" <> var["e"] <> " " <> var["n"] <> " + " <> var["i"] <> "]."

RegularSequence::usage =
"RegularSequence[" <> var["\[Lambda]"] <> ", {" <> sub["m", {0, 1, "...", "k-1"}] <> "}, " <> var["\[Kappa]"] <> "][" <> var["n"] <> "] gives " <> var["\[Lambda]"] <> "." <> sub["m", "n0"] <> "." <> sub["m", "n1"] <> ".\[CenterEllipsis]." <> sub["m", "nl"] <> "." <> var["\[Kappa]"] <> ", where " <> var["n"] <> " \[Equal] " <> sub["n", "l"] <> "\[CenterEllipsis]" <> sub["n", 1] <> sub["n", 0] <> " in base " <> var["k"] <> ".
RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "][{" <> sub["n", {1, 2, "..."}] <> "}] gives terms " <> sub["n", {1, 2, "..."}] <> " of a " <> var["k"] <> "\[Hyphen]regular sequence."

RegularSequenceExponent::usage =
"RegularSequenceExponent[RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "]] gives the exponent in the asymptotic growth rate of the partial sums of a " <> var["k"] <> "\[Hyphen]regular sequence."

RegularSequenceGeneratorTable::usage =
"RegularSequenceGeneratorTable[RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "], " <> var["n"] <> "] gives terms 0 through " <> var["n"] <> " of the generators of the additive group generated by the " <> var["k"] <> "\[Hyphen]kernel of a " <> var["k"] <> "\[Hyphen]regular sequence, where the first generator is the sequence itself."

RegularSequenceMatrixForm::usage =
"RegularSequenceMatrixForm[RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "]] displays " <> var["matrices"] <> " in MatrixForm."

RegularSequenceRank::usage =
"RegularSequenceRank[RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "]] gives the size of the vectors and matrices."

RegularSequenceRecurrence::usage =
"RegularSequenceRecurrence[RegularSequence[" <> var["\[Lambda]", "matrices", "\[Kappa]"] <> "], " <> var["a", "n"] <> "] gives the recurrence relations satisfied by generators " <> sub["a", "i"] <> "[" <> var["n"] <> "] of the additive group generated by the " <> var["k"] <> "\[Hyphen]kernel of a " <> var["k"] <> "\[Hyphen]regular sequence, where " <> sub["a", 1] <> "[" <> var["n"] <> "] is the sequence itself."

ThueMorse::usage =
"ThueMorse[" <> var["n"] <> "] gives the " <> var["n"] <> "th term in the Thue-Morse sequence."


DigitList[n_Integer, b : _Integer : 10] /; b >= 2 := Reverse[IntegerDigits[n, b]]

Trim[array_] :=
With[{length = Min @@ Length /@ array},
	PadRight[array, {Length[array], length}]
]

Unriffle[{}, m : (_Integer?Positive) : 2] := Table[{}, {m}]
Unriffle[list_List, m : (_Integer?Positive) : 2] :=
Module[{p},
	Transpose[Partition[list, m, m, {1, 1}, p]] /. p -> Sequence[]
]


FindRegularSequenceBasis[list : {Except[_List] ..}, k_Integer /; k >= 2] :=
Module[
	{
		basis = {}, m = {}, level = -1, newbasiselements = {{}},
		rank = 0,
		failed
	},
	failed = Catch[
		If[MatchQ[list, {} | {0}], Throw[True]];
		While[newbasiselements != {},
			level++;
			newbasiselements = Pick[
				Thread[{level, Range[0, k^level - 1]}],
				With[{newm = Take[Append[m, #], All, Length[#]]},
					Switch[MatrixRank[newm] - rank,
						0, False,
						_?Positive,
							m = RowReduce[newm];
							rank++;
							True,
						_, Throw[True]
					]
				] & /@ Unriffle[list, k^level]
			];
			basis = Join[basis, newbasiselements]
		];
		If[Length[list] < Max[(#2 + 1 + (k rank - 1) k^#1 &) @@@ basis],
			Throw[True]
		];
		False
	];
	If[failed, $Failed, basis]
]

FindRegularSequenceFunctionFromBasis[list_, k_, basis_] :=
Module[{rank, basismatrix, matrices, failed},
	failed = Catch[
		Quiet[
			Check[
				basismatrix = Transpose[Trim[(Take[list, {#2 + 1, -1, k^#1}] &) @@@ basis]];
				rank = Length[First[basismatrix]];
				matrices = Transpose[(
					(LinearSolve @@ Trim[{basismatrix, #}] &) /@ 
						Unriffle[Take[list, {#2 + 1, -1, k^#1}], k]
				&) @@@ basis],
	
				Message[FindRegularSequenceFunction::sing,
					Transpose[
						(Take[
							list,
							{#2 + 1, #2 + 1 + (rank - 1) k^#1, k^#1}
						] &) @@@ basis
					]
				];
				Throw[True],
	
				LinearSolve::nosol
			],
			{LinearSolve::nosol}
		];
		False
	];
	If[failed,
		$Failed,
		RegularSequence[
			UnitVector[rank, 1],
			matrices,
			list[[Last /@ basis + 1]]
		]
	]
]

FindRegularSequenceFunction[list : {Except[_List] ..}, k_Integer /; k >= 2, n_] :=
	FindRegularSequenceFunction[list, k][n]
FindRegularSequenceFunction[{0, 0 ..}, k_Integer /; k >= 2] :=
	RegularSequence[{1}, ConstantArray[{{0}}, k], {0}]
FindRegularSequenceFunction[list : {Except[_List] ..}, k_Integer /; k >= 2] :=
Module[{basis, regularsequence, failed},
	failed = Catch[
		basis = FindRegularSequenceBasis[list, k];
		If[basis === $Failed, Throw[True]];
		regularsequence = FindRegularSequenceFunctionFromBasis[list, k, basis];
		If[regularsequence === $Failed, Throw[True]];
		False
	];
	regularsequence /; !failed
]
SyntaxInformation[FindRegularSequenceFunction] = {"ArgumentsPattern" -> {{__}, _, _.}}

FindRegularSequenceRecurrence[{0, 0 ..}, k_Integer /; k >= 2, a_[n_]] :=
	Append[Table[a[k n + r] == a[n], {r, 0, k - 1}], a[0] == 0]
FindRegularSequenceRecurrence[list : {Except[_List] ..}, k_Integer /; k >= 2, a_[n_]] :=
Module[{basis, regularsequence, generators, failed, nn},
	failed = Catch[
		basis = FindRegularSequenceBasis[list, k];
		If[basis === $Failed, Throw[True]];
		regularsequence = FindRegularSequenceFunctionFromBasis[list, k, basis];
		If[regularsequence === $Failed, Throw[True]];
		generators = a[k^#1 n + #2] & @@@ basis;
		False
	];
	(RegularSequenceRecurrence[regularsequence, a, n] /.
		ReleaseHold[MapIndexed[
			(* Hold lets this nn receive the same unique number as the other nn. *)
			Hold[RuleDelayed][Subscript[a, First[#2]][nn_], #1] &,
			generators /. n -> nn
		]]
	) /; !failed
]
SyntaxInformation[FindRegularSequenceRecurrence] = {"ArgumentsPattern" -> {{__}, _, _}}

RegularSequenceObjectQ[rs : (List | RegularSequence)[lambda_, matrices : {__}, kappa_]] :=
	MatchQ[RegularSequenceRank[rs], _Integer]

(rs : RegularSequence[lambda_, matrices : {__}, kappa_])[n_Integer] /; RegularSequenceObjectQ[rs] :=
	Dot @@ Join[
		{lambda},
		matrices[[DigitList[n, Length[matrices]] + 1]],
		{kappa}
	]
(* Maybe MapIndexed would be cleaner than MapThread here.  (It's the usual problem of inverting TriangularNumber.) *)
(rs : RegularSequence[lambda_, matrices : {__}, kappa_])[list : {_Integer, __Integer}] /; RegularSequenceObjectQ[rs] :=
Module[{k = Length[matrices], length, zeromatrixpowers, i = -1},
	length = IntegerLength[Max[Abs[list]], k];
	zeromatrixpowers = NestList[#.First[matrices] &, IdentityMatrix[Length[First[matrices]]], length - 1];
	(#.kappa &) /@ Nest[
		Function[vectors,
			i++;
			Join[
				vectors,
				Flatten[Outer[
					#2.#1 &,
					Rest[matrices],
					MapThread[Dot, {
						vectors,
						Prepend[
							Join @@ Table[
								ConstantArray[zeromatrixpowers[[i - j]], (k - 1) k^j],
								{j, 0, i - 1}
							],
							zeromatrixpowers[[i + 1]]
						]
					}],
				1], 1]
			]
		],
		{lambda},
		length
	][[Abs[list] + 1]]
]
(rs : RegularSequence[lambda_, matrices : {__}, kappa_])[list_List] /; RegularSequenceObjectQ[rs] :=
	RegularSequence[lambda, matrices, kappa] /@ list
SyntaxInformation[RegularSequence] = {"ArgumentsPattern" -> {_, _, _}}

RegularSequenceExponent[
	rs : (List | RegularSequence)[_?VectorQ, matrices : {__?MatrixQ}, _?VectorQ] /; RegularSequenceObjectQ[rs]
] := Log[Length[matrices], RootReduce[Abs[First[Eigenvalues[Total[matrices], 1]]]]]
SyntaxInformation[RegularSequenceExponent] = {"ArgumentsPattern" -> {_}}

RegularSequenceGeneratorTable[
	rs : (List | RegularSequence)[_?VectorQ, matrices : {__?MatrixQ}, kappa_?VectorQ] /; RegularSequenceObjectQ[rs],
	n_Integer?NonNegative
] :=
Transpose[
	(Dot @@ Append[
		matrices[[DigitList[#, Length[matrices]] + 1]],
		kappa
	] &) /@ Range[0, n]
]
SyntaxInformation[RegularSequenceGeneratorTable] = {"ArgumentsPattern" -> {_, _}}

RegularSequenceMatrixForm[
	rs : (h : List | RegularSequence)[lambda_?VectorQ, matrices : {__?MatrixQ}, kappa_?VectorQ] /; RegularSequenceObjectQ[rs]
] := h[lambda, MatrixForm /@ matrices, kappa]
SyntaxInformation[RegularSequenceExponent] = {"ArgumentsPattern" -> {_}}

RegularSequenceRank[(List | RegularSequence)[lambda_, matrices : {__}, kappa_]] :=
Module[{explicitmatrices = Cases[matrices, _List], rank, ranks},
	rank = If[explicitmatrices != {}, Length[First[explicitmatrices]], Missing["Unknown"]];
	ranks = DeleteCases[{
		rank,
		If[MatchQ[lambda, _List], Length[lambda], Missing["Unknown"]],
		If[MatchQ[kappa, _List], Length[kappa], Missing["Unknown"]]
	}, _Missing];
	First[ranks] /; Length[ranks] >= 1 && Equal @@ ranks && MatchQ[Dimensions /@ explicitmatrices, {{rank, rank} ...}]
]
SyntaxInformation[RegularSequenceRank] = {"ArgumentsPattern" -> {_}}

RegularSequenceRecurrence[
	rs : (List | RegularSequence)[_?VectorQ, matrices : {__?MatrixQ}, kappa_?VectorQ] /; RegularSequenceObjectQ[rs],
	a_,
	n_
] /; FreeQ[a, n] :=
With[{k = Length[matrices], basis = (Subscript[a, #][n] &) /@ Range[Length[kappa]]},
	Transpose[Thread /@ Append[
		Table[(basis /. n -> k n + r) == matrices[[r + 1]].basis, {r, 0, k - 1}],
		(basis /. n -> 0) == kappa
	]]
]
SyntaxInformation[RegularSequenceRecurrence] = {"ArgumentsPattern" -> {_, _, _}}

SetAttributes[ThueMorse, {Listable}]
ThueMorse[n_Integer] :=
	RegularSequence[{1, 0}, {{{1, 0}, {0, 1}}, {{0, 1}, {1, 0}}}, {0, 1}][n]
SyntaxInformation[ThueMorse] = {"ArgumentsPattern" -> {_}}


End[]

Protect["RegularSequences`*"]

EndPackage[]