Help:=proc()
	if args=NULL then
		print(`collapse(b,n)`, `raboter(b,n)`, `bary_expansion(b,n)`, `SumOriginal(b,n)`, `L(b,n,pow)`, `LByFinalDigit(b,n,pow,final)`, `ClosedForm(rec,n)`, `CountByFinalDigit(b,final,n)`,  `RigorousCountByFinalDigit(b,final,n)`, `GuessSumSquares(b,n)`, `RigorousSumSquares(b,n)`,`SumPowers(b,n,p)`,`SumPowersByFinalDigit(b,final,n,p)`):
	elif args='RigorousSumSquares' then
		print(`RigorousSumSquares(b,n) uses GuessSumSquares to conjecture a closed form expression for the sum of r(b,k)^2 for k=b^n..b^(n+1)-1. It then checks that this expression has the correct initial term and satisfies the proper first-order, linear, non-homogenous recurrence to rigorously prove the conjecture.`):
	elif args='SumPowers' then
		print(`SumPowers(b,n,p) gives a formula in terms of n for the sum of r(b,k)^p for k=b^n..b^(n+1)-1.`):
		print(`For Example: SumPowers(3,n,3);`):
	elif args='SumPowersByFinalDigit' then
		print(`SumPowersByFinalDigit(b,final,n,p) gives a formula in terms of n for the sum of r(b,k)^p over all b^n <= k <= b^(n+1)-1 such that k = final mod b.`):
		print(`For Example: SumPowersByFinalDigit(4,1,n,3);`):
	elif args = 'L' then
		print(`L(b,n,pow) computes the sum r(b,k)^p for k=b^n..b^(n+1)-1 for numeric n in the most naive possible way. Useful for checking that SumPowers works properly.`):
		print(`For Example: L(3,6,3);`):
	elif args = 'LByFinalDigit' then
		print(`LByFinalDigit(b,n,pow,final) computes the sum r(b,k)^p over all b^n <= k <= b^(n+1)-1 such that k = final mod b. Useful for checking that SumPowersByFinalDigit works properly,`):
		print(`For Example: LByFinalDigit(4,6,3,1);`):
	fi:
end:



collapse:=proc(b,n) local ex, i, output,outnum:
	ex:=bary_expansion(b,n):
	if ex[1] <> ex[2] then	
		ex[1] := -1:
	fi:
	if ex[-1] <> ex[-2] then
		ex[-1] := -1:
	fi:
	for i from 2 to nops(ex)-1 do
		if ex[i]<>ex[i-1] and ex[i]<>ex[i+1] then
			ex[i] := -1:
		fi:
	od:
	output:=[]:
	for i from 1 to nops(ex) do
		if ex[i] <> -1 then
			output:=[op(output),ex[i]]:
		fi:
	od:
	outnum:=add(output[i]*b^(i-1),i=1..nops(output)):
	return(outnum):
end:
		
raboter:=proc(b,n) local ex,i,output,outnum:
	ex:=bary_expansion(b,n):
	for i from 1 to nops(ex) do
		if i = nops(ex) or ex[i+1] <> ex[i] then
			ex[i] := -1:
		fi:
	od:
	output:=[]:
	for i from 1 to nops(ex) do
		if ex[i] <> -1 then
			output:=[op(output),ex[i]]:
		fi:
	od:
	outnum:=add(output[-i]*b^(i-1),i=1..nops(output)):
	return(outnum):
end:

L:=proc(b,k,pow) local n:
	add(raboter(b,n)^pow,n=b^k..(b^(k+1)-1)):
end:

LByFinalDigit:=proc(b,k,pow,final) local n:
	add(raboter(b,n)^pow,n=b^k+final..(b^(k+1)-1),b):
end:
	
bary_expansion:=proc(b,n) local output, np:
	output:=[]:
	np:=n:
	while np>0 do
		output := [np mod b, op(output)]:
		np:=np - (np mod b):
		
		np:=np/b:
	od:
	return(output):
end:

ClosedForm:=proc(rec,n)	local init,recu,order,char_poly,i,x,vars,a,bases,eqns,pow,terms,b:
	init:=rec[1]:
	recu:=rec[2]:
	order :=nops(recu):
	char_poly:= x^order - add(recu[i]*x^(order-i), i=1..order):
	bases:=[solve(char_poly)]:
	terms:=[]:
	for b in bases do
		pow := 0:
		while member(b^n*n^pow,{op(terms)}) do
			pow := pow+1:
		od:
		terms:=[op(terms),b^n*n^pow]:
	od:
	vars:={seq(a[i],i=1..order)}:
	eqns:={seq(init[i] = add(a[j]*subs(n=i,terms[j]),j=1..order),i=1..order)}:
	vars:=subs(solve(eqns,vars),[op(vars)]):
	return(add(vars[i]*terms[i],i=1..order)):
end:
	
CountByFinalDigit:=proc(b,final,n) local i,m,seq_by_digit,new_elt,rec,new_rec:
	seq_by_digit:=[seq(add(raboter(b,m), m=b^i+final..b^(i+1)-1,b),i=1..2)]:
	while GuessRec(seq_by_digit) = FAIL do
		new_elt := SumOriginal(b,nops(seq_by_digit))+(b-1)*seq_by_digit[-1]+final*(b-1)*b^(nops(seq_by_digit)-1):
		seq_by_digit:=[op(seq_by_digit),new_elt]:
	od:
	rec:=GuessRec(seq_by_digit):
	while rec[2][-1] = 0 do
		rec:=[[op(2..-1,rec[1])],[op(1..-2,rec[2])]]:
	od:
	return(ClosedForm(rec,n)):
end:

RigorousCountByFinalDigit:=proc(b,final,n) local n1,guess,prev_term,hope:
	guess:=CountByFinalDigit(b,final,n1):
	prev_term := subs(n1 = n1 -1,guess):
	hope := SumOriginal(b,n1-1)+(b-1)*prev_term+final*(b-1)*b^(n1-2):
	if ((subs(n1=1,guess) = final) and simplify(normal(hope)) = simplify(normal(guess))) then
		return(subs(n1=n,guess)):
	else
		return(FAIL):
	fi:
end:	

GuessSumSquares:=proc(b,n) local seq_overall,seq_digits,n1,m,i,new_digit,rec:
	seq_overall:=[seq(add(raboter(b,n1)^2, n1=b^i..b^(i+1)-1), i = 1 .. 2)]:
	seq_digits:=[seq(CountByFinalDigit(b,i,n1),i=0..b-1)]:
	
	while GuessRec(seq_overall) = FAIL do
		m:=nops(seq_overall)+1:
		new_digit := (b^2+b-1)*seq_overall[-1] + add(2*b*(i-1)*subs(n1=m-1,seq_digits[i]),i=1..b)+(b-1)*(b)*(2*b-1)/6*(b-1)*b^(m-2):
		seq_overall:=[op(seq_overall),new_digit]:
	od:
	rec:=GuessRec(seq_overall):
	while rec[2][-1] = 0 do
		rec:=[[op(2..-1,rec[1])],[op(1..-2,rec[2])]]:
	od:
	return(ClosedForm(rec,n)):
end:

RigorousSumSquares:=proc(b,n) local guess,prev_term,hope,i:
	guess := GuessSumSquares(b,n1):
	prev_term := subs(n1 = n1-1,guess):
	hope := (b^2+b-1)*prev_term + add(2*b*(i)*subs(n1 = n1-1, RigorousCountByFinalDigit(b,i,n1)),i=0..b-1) + (b-1)*(b)*(2*b-1)/6*(b-1)*b^(n1-2):
	if (subs(n1=1,guess)=(b-1)*b*(2*b-1)/6 and simplify(normal(hope)) = simplify(normal(guess))) then
		return(subs(n1=n,guess)):
	else
		return(FAIL):
	fi:
end:
	
SumOriginal:=proc(b,n) 
	return(b*(b-1)/(2*b-1)*(2*b-1)^n-(b-1)/2*b^n):
end:

SumSquaresByFinalDigit:=proc(b,final,n) local count,sol,overall:
	overall:=RigorousSumSquares(b,n):
	count:=RigorousCountByFinalDigit(b,final,n):
	sol:=rsolve({a(n)=subs(n=n-1,overall) + (b^2-1)*a(n-1) + 2*final*b*subs(n=n-1,count) + final^2*(b-1)*b^(n-2),a(1)=final^2},a(n)):
end:

SumPowers:=proc(b,n,p) local a,rec,i,j:
	option remember:
	if p = 0 then
		return((b-1)*b^n):
	fi:
	if p = 1 then
		return(SumOriginal(b,n)):
	fi:
	rec := a(n) = (b-1)*a(n-1) + b^p*a(n-1) + add(add(b^(p-i)*j^i*binomial(p,i)*subs(n=n-1, SumPowersByFinalDigit(b,j,n,p-i)) ,i=1..p), j=0..(b-1)):
	return(rsolve({rec, a(1) = add(i^p, i=1..(b-1))},a(n))):
end:

SumPowersByFinalDigit:=proc(b,final,n,p) local a,rec,i,j:
	option remember:
	if p = 0 then
		return((b-1)*b^(n-1)):
	fi:
	rec := a(n) = (b^p-1)*a(n-1) + subs(n=n-1,SumPowers(b,n,p)) + add(final^i*b^(p-i)*binomial(p,i)*subs(n=n-1,SumPowersByFinalDigit(b,final,n,p-i)), i=1..p):
	return(rsolve({rec,a(1) = final^p},a(n))):
end:
	

NextTerm:=proc(L,rec) local i:
	if nops(rec) > nops(L) then
		print(`Not enough terms!`):
		return(FAIL):
	fi:
	add(L[-i]*rec[i],i=1..nops(rec)):
end:


#Inspired by the function from Doron Zeilberger's package Cfinite
GuessRec:=proc(L) local order, vars,eqns,guess,i,j:
	order := 1:
	while 2*order+1<nops(L) - 2 do
		vars:={seq(a[i],i=1..order)}:
		eqns:={seq(add(a[i]*L[j-i],i=1..order)=L[j],j=order+1..2*order)}:
		vars:=subs(solve(eqns,vars),[op(vars)]):
		for i from 2*order+1 to nops(L) do
			guess:=NextTerm(L[1..i-1],vars):
			if guess <> L[i] then
				break:
			fi:
		od:
		if i = nops(L)+1 then
			return([L[1..order],vars]):
		fi:
		order:=order+1:
	od:
	return(FAIL):
end: