(* :Title:  Egypt  *)

(* :Author:
	David Eppstein
	eppstein@ics.uci.edu
	http://www.ics.uci.edu/~eppstein/
*)

(* :Version:  1.3  *)

(* :Copyright:
	Copyright 2000, David Eppstein.

	Permission is hereby granted to modify and/or make copies of this file 
	for any purpose other than direct profit, or as part of a commercial 
	product, provided this copyright notice is left intact.  Sale, other 
	than for the cost of media, is prohibited.  Permission is hereby 
	granted to reproduce part or all of this file, provided that the source 
	is acknowledged.
      
	I hereby acknowledge that this copyright notice and other package 
	boilerplate (but no actual code) were reproduced from CleanSlate.m, by 
	Todd Gayley, included in the Mathematica 3.0 Distribution.
*)

(* :History:

	Created, 7 Nov 1997, based on a suggestion by Stan Wagon to improve the 
	user interface for the Egyptian Fraction routines I published in MiER.

	Sometime in the middle here released on MathSource.

	20 Mar 2000, v1.1, modified SmallMultiples to use dynamic programming
	instead of brute force

	25 Mar 2000, v1.2, modified greedy and lexicographic to prune based on
	same EstimatedTerms function used by new SmallMultiples code.

	27 Jun 2000, v1.3, add Method->BaseConversion thanks to a suggestion of
	Mark Armbrust <marka@ff.com>. Fix doc for "Verbose".
*)

(* :Context:  Egypt`  *)

(* :Mathematica Version:  3.0 
	You must have v.3.0 or later to use this package.
*)

(*********************** START OF PUBLIC CODE ***********************)

BeginPackage["Egypt`"]

Unprotect[EgyptianFraction]

EgyptianFraction::usage = "EgyptianFraction[x/y, opts...] converts the
rational number x/y to a sum of distinct unit fractions."

Optimize::usage = "Optimize is an option to EgyptianFraction that determines
which of multiple possible decompositions to return.  Its value should be a 
function which computes a number from a list of denominators; only those 
decompositions minimizing this function are returned. For instance 
Optimize->Max tries to minimize the maximum denominator, while 
Optimize->Length tries to find representations that are as short as 
possible. Optimize->All returns all decompositions found by the chosen 
algorithm."

Monotonic::usage = "Monotonic is an option to EgyptianFraction.  If
Monotonic is True, EgyptianFraction will assume that the value of the
Optimize function on a partial solution is no larger than the value of any
complete solution it leads to, and prune the partial
solution if better solutions have already been found."

MinTerms::usage = "MinTerms is an option to EgyptianFraction that, when set
to an integer n, causes only those decompositions having at least n terms
to be returned."

MaxTerms::usage = "MaxTerms is an option to EgyptianFraction that, when set
to an integer n, causes only those decompositions having at most n terms
to be returned."

MinDenominator::usage = "MinDenominator is an option to EgyptianFraction
that, when set to an integer n, causes only those decompositions having all
denominators at least n to be returned."

MaxDenominator::usage = "MaxDenominator is an option to EgyptianFraction
that, when set to an integer n, causes only those decompositions having all
denominators at most n to be returned."

Duplicates::usage = "Duplicates is an option to EgyptianFraction that
determines what to do when the method being used generates multiple copies
of the same unit fraction. Possible choices are Allow, Disallow, Pairing,
and Splitting."

Allow::usage = "Allow is a choice for the option Duplicates of
EgyptianFraction. Duplicates->Allow allows decompositions with repeated terms
to be returned."

Disallow::usage = "Disallow is a choice for the option Duplicates of
EgyptianFraction.  Duplicates->Disallow prohibits decompositions with
repeated terms from being returned."

Pairing::usage = "Pairing is a choice for the options Duplicates and Method
of EgyptianFraction. Duplicates->Pairing deals with repeated terms by
repeatedly replacing pairs of terms 1/x + 1/x with either 2/x (if x is even)
or 2/(x+1) + 2/x(x+1) (if x is odd). Method->Pairing applies this duplicate
removal process to an initial decomposition formed by repeated copies of
the inverse of the input denominator."

Splitting::usage = "Splitting is a choice for the options Duplicates and
Method of EgyptianFraction.  Duplicates->Pairing deals with repeated terms
by repeatedly replacing duplicated terms 1/x with 1/(x+1) + 1/x(x+1).
Method->Pairing applies this duplicate removal process to an initial
decomposition formed by repeated copies of the inverse of the input denominator."

Exclude::usage = "Exclude is an option to EgyptianFraction that, when set
to a list of integers causes only those decompositions that don't use
denominators in the list to be returned."

Filter::usage = "Filter is an option to EgyptianFraction that, when set to
a function, causes only those decompositions for which the function is true
of all denominators to be returned."

OutputFormat::usage = "OutputFormat is an option to EgyptianFraction that
determines how it presents its decompositions.  Possible choices are Plain
and Fancy."

Plain::usage = "Plain is a choice for the option OutputFormat of
EgyptianFraction.  OutputFormat->Plain causes the output to be returned
as a list of lists of unit fractions."

Fancy::usage = "Fancy is a choice for the option OutputFormat of
EgyptianFraction.  OutputFormat->Fancy causes the output to be formatted
as a column of sums of unit fractions."

(*

Already documented in Mma as a general-purpose option with this effect.

Method::usage = "Method is an option to EgyptianFraction that determines
which algorithm it will use to find unit fraction decompositions.  Possible
values are BinaryRemainder, ContinuedFraction, FactorMultiples, Greedy,
Lexicographic, MinimizeTerms, Remainder, ReverseGreedy, and SmallMultiples.
Specifying a list of methods causes the first method to be used at the top
level, the second method to be used for recursive calls from the first
method, etc."

Individual method doc strings below are modeled after Mma's Newton::usage, 
so don't tell me they're nonstandard.

*)

BaseConversion::usage = "BaseConversion is a choice for the option Method
of EgyptianFraction.  Method->BaseConversion generates a mixed base such
that the digital expansion of the fraction in that base is 0.1111...  which
leads to an Egyptian fraction representation in which each denominator is a
multiple of the previous one.  Suggested by Midhat GazalŽ in his book _Number:
From Ahmes to Cantor_."

Binary::usage = "Binary is a choice for the option Method of
EgyptianFraction.  Method->Binary uses fractions that are either a power of
two or a difference of two such powers, reducing the complexity of the
binary representation of the remaining fraction at each step."

BinaryRemainder::usage = "BinaryRemainder is a choice for the option Method
of EgyptianFraction.  Method->BinaryRemainder uses powers of two until the
Remainder method succeeds in decomposing the remaining fraction."

ContinuedFraction::usage = "ContinuedFraction is a choice for the option
Method of EgyptianFraction.  Method->ContinuedFraction produces
decompositions in which the partial sums are truncations of the input's
continued fraction expansion."

Greedy::usage = "Greedy is a choice for the option Method of
EgyptianFraction. Method->Greedy uses the largest possible unit fraction
at each step."

Lexicographic::usage = "Lexicographic is a choice for the option Method of
EgyptianFraction. Method->Lexicographic tries all decompositions in lexicographic
order: largest fraction first, then decompositions with the same largest
fraction are ordered by their second largest fraction, etc.  Another option
such as Optimize->Length or MaxTerms->4 is required to avoid infinite
loops."

MinimizeTerms::usage = "MinimizeTerms is a choice for the option Method of
EgyptianFraction. Method->MinimizeTerms tries all decompositions in order by
number of terms, stopping when some number of terms leads to a solution."

Remainder::usage = "Remainder is a choice for the option Method of
EgyptianFraction. Method->Remainder looks for decompositions in which each
term is a divisor of the input denominator.  The additional option
Multiplier->n causes the input x/y to be treated as an improper fraction
nx/ny in this method."

ReverseGreedy::usage = "ReverseGreedy is a choice for the option Method of
EgyptianFraction. Method->ReverseGreedy forms decompositions by removing a unit
fraction that causes the remaining value to have a reduced denominator, and
recursively decomposing that remaining value."

SmallMultiples::usage = "SmallMultiples is a choice for the option Method
of EgyptianFraction.  Method->SmallMultiples tries all combinations of
denominators that are multiples of the largest prime power dividing the
original denominator, until finding a combination such that the power does
not divide the remaining fraction. Best used in combination with another
method e.g. Method->{SmallMultiples,ReverseGreedy}."

SquareGreedy::usage = "SquareGreedy is a choice for the option Method of 
EgyptianFraction.  Method->SquareGreedy repeatedly removes the largest unit 
fraction the denominator of which is a perfect square.  Ron Graham proved 
in his thesis that a perfect-square unit fraction decomposition exists for 
all fractions less than Pi^2/6 - 1."

Verbose::usage = "Verbose is an option to EgyptianFraction that, when set to
True, turns on various print statements that show the fractions remaining
and other information showing how the computation is progressing."

Multiplier::usage = "Multiplier is an option to EgyptianFraction that, when 
set to an integer n, forms an improper fraction by multiplying both 
numerator and denominator by n.  This is primarily useful in the remainder 
and binary remainder methods; for those methods, n should be a number with 
many factors such as a factorial or power of two."

Options[EgyptianFraction] =
   {Optimize->All, Monotonic->True, MinTerms->Automatic, MaxTerms->Automatic,
    MinDenominator->Automatic, MaxDenominator->Automatic,
    Exclude->{}, Filter->All, Method->Greedy, Duplicates->Disallow,
    OutputFormat->Fancy, Verbose->False, Multiplier->1};

(*********************** START OF PRIVATE CODE ***********************)

Begin["`Private`"]

EgyptianFraction[Rational[p_,q_],opts___] :=
	Module[{numerator = p, denominator = q, output = {}, current = {},
			optfn, optval=Infinity,, minterms, maxterms, minden, maxden,
			exclude, filter, method, debug, options = Options[EgyptianFraction],
			methodtranslations, mult, monotonic, duplicates, allowdups,
			
			(* the rest are not variables, but subroutines that
			   depend on the module context as part of their defn *)

			good,addsolution,firstmethod,add,try,noterms,oneterm,
			brute,minimizeterms,remainder,binary,binrem,pairing,splitting,
			reversegreedy,greedy,continuedfraction,smallmult,smallmultrecurse},
		
		(* set options *)
		optfn = (Optimize /. {opts} /. options /. {All -> (1&)});
		monotonic = (Monotonic /. {opts} /. options);
		minterms = If[IntegerQ[#],#,1]& @ (MinTerms /. {opts} /. options);
		maxterms = If[IntegerQ[#],#,Infinity]& @ (MaxTerms /. {opts} /. options);
		minden = If[IntegerQ[#],#,2]& @ (MinDenominator /. {opts} /. options);
		maxden = If[IntegerQ[#],#,Infinity]& @ (MaxDenominator /. {opts} /. options);
		exclude = (Exclude /. {opts} /. options);
		filter = (Filter /. {opts} /. options /. {All->Function[x,True]});
		debug = (Verbose /. {opts} /. options);
		mult = (Multiplier /. {opts} /. options);

		duplicates = (Duplicates /. {opts} /. options);
		allowdups = Not[duplicates === Disallow];
		duplicates = duplicates /.
			{ Allow -> (#&), Disallow -> (#&),
			  Pairing -> PairFracs, Splitting -> SplitFracs };
		
		(* utility subroutines: test if denom fits in current partial soln *)
		good = (# >= minden && # <= maxden && filter[#] &&
				Not[MemberQ[exclude,#]] &&
				(allowdups || Not[MemberQ[current,1/#]]))&;
		
		(* utility subroutines: include newly found solution into output *)
		addsolution =
		  ( current = duplicates @ current;
			Module[{solnterms = Length[current],
					solnden = (Max @@ (1/current)),
					solnval = optfn @ (1/current),
					retval = False},
				solnterms >= minterms && solnval <= optval &&
				  (	debug && Print[StringForm["found soln ``",current]];
					output =
						If[ solnval < optval, {current},
							Append[output,current] ];
					(optfn === Max) && (maxden = solnden);
					(optfn === Length) && (maxterms = solnterms);
					optval = solnval;
					True )
			] )&;
				
		(* utility subroutines: add new terms to current partial solution
		   and call appropriate recursive method on augmented solution
		   caller is responsible for making sure terms are good
		   returns boolean: did any new solutions get found recursively *)

		add = Function[terms,
			Module[{savedcurrent=current,savedmethod=method,
					savednum=numerator,savedden=denominator,
					frac = (numerator/denominator - Plus@@terms), retval},
				numerator = Numerator[frac];
				denominator = Denominator[frac];
				current = Join[current,terms];
				debug && numerator != 0 &&
					Print[StringForm["trying ``, remainder ``", current, frac]];
				method = If[Length[method] > 1, Rest[method], method];
				retval = If[ monotonic && (optfn @ (1/current)) > optval,
						     False,
							 firstmethod[][] ];
				current = savedcurrent;
				method = savedmethod;
				numerator = savednum;
				denominator = savedden;
				retval] ];
				
		(* utility subroutines: try all good denominators in list *)
		try = Function[ dens,
			Fold[ Function[{retval,den}, (good[den] && add[{1/den}]) || retval ],
				  False, Select[dens,IntegerQ] ] ];
		
		(* method: base conversion *)
		baseconv = Function[ {},
			Module[ { x = 1, n=numerator, d=denominator, b, L={} },
				While[ n>0,
					b = Floor[(d+n-1)/n];
					x = x*b;
					L = Append[L,1/x];
					n = (b*n)-d ];
				add[L]
			] ];
		
		(* method: pairing and splitting *)
		pairing =
		  (	duplicates = PairFracs;
		  	add[Table[1/(mult*denominator), {mult*numerator}]] )&;
		splitting =
		  (	duplicates = SplitFracs;
		  	add[Table[1/(mult*denominator), {mult*numerator}]] )&;

		(* methods: what to do when maxterms=0 or 1 *)
		noterms = Function[ {}, False ];
		oneterm = Function[ {},
			numerator == 1 && good[denominator] && add[{1/denominator}] ];
			
		(* methods: rlg's square greedy *)
		sqgreedy = Function[ {},
			Module[{oldminden = minden, x, retval},
				x = Ceiling[Sqrt[Max[minden,denominator/numerator]]]^2;
				good(x) && (minden = x+1;
							retval = add[{1/x}];
							minden = oldminden;
							retval)]];
					
		(* methods: reverse greedy *)
		reversegreedy = Function[ {}, try[ RGChoices[numerator, denominator] ] ];

		(* methods: remainder and binary remainder *)
		remainder = Function[ {},
			Module[ { x = mult*numerator, y = mult*denominator,
					  divs, S = {{}}, retval = False,
					  L = maxterms-Length[current] },
				divs = y/Select[Divisors[y],good];
				divs = Select[divs,(# <= x)&];
				While[ Not[divs=={}],
					S = Select[	Join[S, Append[#,First[divs]]& /@ S],
								 (Length[#]<=L && Plus@@# <= x)& ];
					divs = Rest[divs] ];
				S = Select[S, (Plus@@# == x)& ];
				While[ Not[S=={}],
					   add[First[S]/y] && (retval = True);
					   S = Rest[S] ];
				retval
			] ];
			
		binrem = Function[ {},
			remainder[] || (good[#] && add[{1/#}])& @
				(2^Ceiling[Log[denominator/numerator]/Log[2]]) ];
		
		(* methods: binary *)
		binary =
			Module[ { i = Ceiling[Log[denominator/numerator]/Log[2]],
					  cl = CycleLength[denominator],
					  bc = BinaryComplexity[numerator/denominator],
					  retval = False },
				( IntegerQ[#] && good[#] && # >= denominator/numerator &&
				  BinaryComplexity[numerator/denominator - 1/#] < bc &&
				  add[{1/#}] && (retval = True) )& /@
					{ 2^i, 2^i - 2^(i-cl) };
				retval
			]&;
		
		(* methods: greedy *)
		greedy = Module[
				{ k, L, retval = False, mt=maxterms-Length[current] },
			If[ mt <= 2,
				(	L = Select[TwoTerm[numerator,denominator],IntegerQ];
					While[	Not[retval || L=={}],
							k = First[L];
							L = Rest[L];
							good[k] && add[{1/k}] && (retval = True)	] ),
				If[ EstimatedTerms[numerator/denominator,numerator/denominator] > mt, {},
				(	k = Max[minden,Ceiling[denominator/numerator]];
					While[ Not[retval] && k <= maxden &&
							k*numerator <= denominator*mt,
					   good[k] && add[{1/k}] && (retval = True);
					   k = k+1 ] )  ] ];
			retval ]&;
			
		(* methods: lexicographic *)
		brute = Module[
				{ savedminden = minden, retval = False,
				  k = Max[minden,Ceiling[denominator/numerator]],
				  mt=maxterms-Length[current] },
			While[ k <= maxden && k*numerator <= denominator*mt,
				   If[	mt <= 2,
				   		try[TwoTerm[numerator,denominator]] &&
				   			(retval = True);
				   			k = Infinity,
						If[ EstimatedTerms[numerator/denominator,numerator/denominator] > mt,
							maxden = k-1,
							good[k] && (minden = k; add[{1/k}]) &&
								(retval = True; mt = maxterms-Length[current]);
								k = k+1 ] ] ];
			minden = savedminden;
			retval ]&;
			
		(* methods: minimize terms *)
		minimizeterms = Function[ {},
			oneterm[] ||
			Module[ {savedmaxterms = maxterms, retval = False},
				maxterms = 2+Length[current];
				While[	maxterms <= savedmaxterms && Not[retval],
						brute[] && (retval = True);
						maxterms = maxterms+1	];
				maxterms = savedmaxterms;
				retval ] ];
		
		(* methods: continued fraction *)
		continuedfraction =
			try[ 1/Rest[ FoldList[ Plus, 0,
						 		   Reverse[ ContinuedFractionExpansion [
						 		   		numerator / denominator ] ] ] ] ]&;
		
		(* methods: small multiples
		   instead of brute force expansion of all combinations of small
		   multiples, we use dynamic programming to test whether some
		   combination of multiples ni has
				sum(1/ni mod div) = numerator /(denominator/div)  (mod div)
		   then use backtracking enumeration of only acceptable combinations
		   local variable "reached" has a table of how many terms are
		   needed to reach residue i (mod div). "prevreached" is a list
		   of four values: the values of ni, invni, reached, and prevreached
		   in the previous iteration. *)

		smallmultrecurse = Function[ {target,mt,data,terms,div},
			If[ data=={}, False,
			debug && Print[StringForm["smallmultrecurse, target=``, ni=``, 1/ni=``, terms=``",
										target,data[[1]],data[[2]],terms]];
			If[ target == data[[2]],
				add[1/(div*Append[terms,data[[1]]])],
				If[ data[[3]][[target]] > mt,
					False,
					Or[ smallmultrecurse[ Mod[target - data[[2]],div], mt - 1,
										  data[[4]], Append[ terms, data[[1]]], div ],
						Evaluate[ smallmultrecurse[ target, mt, data[[4]], terms, div ]]]]]]];

		smallmult =
			Module[ { div = Max[Power@@#& /@ FactorInteger[denominator]],
					  target, reached, estleft, estaddlterms = 0, harmonic = 0,
					  ni = 0,
					  invni = 0,
					  prevreached = {},
					  savedmaxden = maxden,
					  retval = False,
					  mt = maxterms-Length[current],
					  subsets = {{}}, S },
				debug && Print[StringForm["smallmult: ``/``, div = ``", numerator, denominator,div]];
				target = Mod[numerator PowerMod[denominator/div,-1,div], div];
				debug && Print[StringForm["target = ``", target]];
				reached = Table[ Infinity, {div-1} ];
				While[ Not[retval] && ni div < savedmaxden,
					ni = ni + 1;
					harmonic = harmonic + 1/ni;
					maxden = savedmaxden;
					If[ good[ni div] && Mod[ni,div] != 0,
						debug && Print[StringForm["including multiple ``",ni]];
						maxden = ni div;
						invni = PowerMod[ni,-1,div];
						reached = Table[ If[ i == invni, 1,
										 Min[ reached[[i]],
											  reached[[Mod[i-invni,div]]]+1 ] ],
										 {i, div - 1} ];
						prevreached = {ni, invni, reached, prevreached};
						retval = smallmultrecurse[target,
									mt - EstimatedTerms[numerator/denominator - harmonic/div,
														numerator/denominator],
									prevreached,{},div] ] ];
				maxden = savedmaxden;
				retval
			]&;
					  
		(* set method option *)		
		methodtranslations = {
			BaseConversion -> baseconv,
			Binary -> binary,
			BinaryRemainder -> binrem,
			ContinuedFraction -> continuedfraction,
			Greedy -> greedy,
			Lexicographic -> brute,
			MinimizeTerms -> minimizeterms,
			Pairing -> pairing,
			Remainder -> remainder,
			ReverseGreedy -> reversegreedy,
			SmallMultiples -> smallmult,
			Splitting -> splitting,
			SquareGreedy -> sqgreedy
			};

		method = (Method /. {opts} /. options /. methodtranslations);
		If[Not[ListQ[method]], method = {method}];

		firstmethod = 
			Module[{L = maxterms-Length[current]},
				If[ numerator==0, addsolution,
					If[	L<=0 || numerator<0 || maxden<denominator/numerator,
						noterms,
						If[	L==1, oneterm, First[method]]]]]&;
		
		(* GO FOR IT *)
		firstmethod[][];

		(* output results *)
		output = Union[Reverse[Sort[#]]& /@ output];
		If[	output == {} || Not[(OutputFormat /. {opts} /. options) === Fancy],
			output,
			If[	Length[output] == 1,
				Plusify[ output[[1]] ],
				ColumnForm[FormatSum /@ output] ] ]
	];

(************************* END OF MAIN ROUTINE ************************)

(* estimate number of terms needed to represent a fraction in the closed
   interval [lb,ub].  For correctness, the answer must always be an 
   underestimate.

   The trick used is that if we instead had an open
   interval [lb,1/x), the number of terms would be bounded by how many
   iterations of the greedy method are needed to get into that interval
   (any other sum a/b of equally many unit fractions would result in
   1/x-a/b having a smaller denominator, and the greedy method gives 
   numerator=1).  So we round ub up to the nearest larger unit fraction
   then apply this method. *)

EstimatedTerms[lb_,ub_] :=
 Module[{invub,nterms=0,remaininglb},
		remaininglb=ub-lb;
		If[ub>=1,nterms=1;ub=FractionalPart[ub]];
		remaininglb=ub-remaininglb;
		invub=Ceiling[1/ub]-1;
		While[remaininglb>0,
			nterms=nterms+1;
			remaininglb=remaininglb-1/(invub+1);
			invub=invub*(invub+1)];
		nterms];

(* reverse greedy method: make list of denominators to try *)
RGChoices[x_,y_]:=
	If[	x==1,
		{y},
    	Union[ (y Mod[# PowerMod[x,-1,y],y] / #)& /@
        		Select[Divisors[y y], (Mod[# PowerMod[x,-1,y],y] x > #)&] ] ];
        		
(* Find possible terms in two-term expansion.

   Suppose x/y = 1/a + 1/b, and consider d=xa-y, so a=(d+y)/x, b=y(d+y)/xd.
   Consider a possible prime divisor p of d.  If p divides y exactly c 
   times, then (for b to be an integer) p can divide d at most 2c times,
   so d must be a divisor of y^2.
   
   To avoid duplication we only look at divisors for which 1/a > x/2y;
   equivalently we need divisors that are at less than y.
   
   Finally, if x=1 or x=2 and we allow duplicates, we need to use y itself. *)

TwoTerm[x_,y_] := Prepend[(y + Select[Divisors[y^2],(#<y)&])/x, y];

(* binary method *)

CycleLength[den_] :=
	Module[ { i = 2, d = den, len = 1 },
		While [ Mod[d,2] == 0, d = d/2 ];
		While [ i != 1,
			i = Mod[2*i,d];
			len = len + 1
		];
		len
	];

BinaryComplexity[Rational[num_,den_]] :=
	Module[ { cl = CycleLength[den], i = num, j, tl = 0, bc = 0 },
		j = Nest[Mod[2*#,den]&, i, cl];
		While [ cl > 0,
			If [ 2*i >= den, bc = bc + 1 ];
			If [ tl < cl,
				 tl = tl + 1,
				 i == j && (cl = cl - 1) ];
			i = Mod[2*i,den];
			j = Mod[2*j,den]
		];
		bc
	];

(* continued fraction method *)

ContinuedFractionList[q_] :=
	Floor /@ Drop[ FixedPointList[ If[IntegerQ[#],0,1/(#-Floor[#])]&, q], -2];

MakeOdd[L_] :=
    If[ OddQ[Length[L]], L, Join[Drop[L, -1], {Last[L] - 1, 1}] ];

PrimarySub[{a_,b_},c_] := {b + a c, a};

PrimarySequence[L_] :=
    Transpose[Drop[FoldList[PrimarySub, {0, 1}, L], 1]] [[1]];

SecondarySequence[L_] := 
    If[	Length[L] < 3, L, 
		Join[ Table[ L[[1]] + i L[[2]], {i, 0, (L[[3]]-L[[1]])/L[[2]]-1} ],
			  SecondarySequence[Drop[L,2]] ] ];

ContinuedFractionExpansion[q_] :=
	(1/(Drop[#,1]Drop[#,-1]))& @
		SecondarySequence[PrimarySequence[MakeOdd[ContinuedFractionList[q]]]];

(* duplicate removal *)

DoPairing[p___, q : Rational[1, y_], q_, r___] := 
    If[	OddQ[y],
    	DoPairing[p, 2/(y(y + 1)), 2/(y + 1), r], 
		DoPairing[p, 2/y, r] ];
DoPairing[p___, q_Integer, r_Integer, s___] := DoPairing[p, q + r, s]; 
SetAttributes[DoPairing, Orderless];
PairFracs[L_] := List @@ DoPairing @@ L;

DoSplitting[p___, q : Rational[1, y_], q_, r___] := 
    DoSplitting[p, q, 1/(y(y + 1)), 1/(y + 1), r];
SetAttributes[DoSplitting, Orderless];
SplitFracs[L_] := List @@ DoSplitting @@ L;

(* output formatting: make fraction in column not look too small *)
FormatSum[L_] := StyleForm[ Plusify[L], ScriptSizeMultipliers->1 ];

(* output formatting: hack list of unit fractions to look like sum *)
Plusify[L_] :=
	If[ Length[L]==1, L[[1]],
		ReplacePart[ ReplacePart[x[L],HoldForm,0], Plus,{1,0} ] ];

End[];

Protect[EgyptianFraction];

EndPackage[]
