(* Atul Pokharel, bbseqe code. Modified from code due to SL *)

(*10/13: There are problems of syntax*)

(* 10/12: Fixed bbstep3 so that it handles bz=0, b3=0 or both correctly. This solution involved
explicitly setting certain values to Infinity instead of allowing a division by zero. The goal
was to make bbseqe handle this case consistently without changing its behavior at any other
points.*)

(* 10/4: Rewrote code, changed variable names to be more readable, and consistent between
functions. A good learning exercise. *)


bbseqMemo={revised, version, of, bbseq, program, of, SL, to, work, when, 
bz, b3, or, both, initially, equal, to, zero}

(*startdata, psb00, eps00, std00, acc00*)

bbseqe[m1_, m2_, m3_, m4_, m5_] := 
 Module[ {ofilename,r},
	Clear[timeStart, timeFinish, psb00, eps00, std00, acc00, systart];
	timeStart=Date[];
	psb00=m2; eps00=m3; std00=m4; acc00=m5; systart=m1;
	r=bbseq[m1, m2, m3, m4, m5];
	bboutput=r;
	timeFinish=Date[];

(*	ofilename="out.bbseqe."<>ToString[timeFinish[[1]]]<>
"."<>ToString[timeFinish[[2]]]<>"."<>ToString[timeFinish[[3]]]<>
"."<>ToString[timeFinish[[4]]]<>"."<>ToString[timeFinish[[5]]]<>
"."<>ToString[timeFinish[[6]]];
*)
	ofilename="bbseqe.out";
	Save[ofilename,psb00, eps00, std00, acc00, timeStart, timeFinish, result];
	Return[r];
	]

bbseq[systart_,psb00_, eps00_,stepsdata_, acc00_Integer] :=
 Module[ {bza, b3a, b, k1, k2, j, st, st2},
	{bza,b3a}=psb00;
	{b,k1,k2,j,st,st2}=stepsdata;
	Return[bbseq4[eps00, bza, b3a, b, k1, k2, j, st, st2, systart, acc00]];
	]

bbseq4[eps00_, bza_, b3a_, b_, k1_, k2_, j_, st_Integer, st2_Integer, systart_, acc00_Integer]:=
 Module[ {ffav, ct, ct2, syva, bzv, b3v, r},
	ffav=ff0/.eps00;
	ct=0; ct2=0;
	syva=systart; 
	bzv=bza;
	b3v =b3a;
	r={bza, b3a, systart};
	Label[o1];
	{bzv, b3v} =bbstep4[b,bzv,b3v, k1, k2, j];
	syva=systep[syva, bzv,b3v,ffav,acc00];
	ct=ct+1;
	ct2=ct2+1;
	If[ct2<st2, Goto[o1]]; ct2=0;
	r=r+u^ct*{bzv, b3v, syva}; 
	If[ct<st, Goto[o1]]; Return[r];
	]

bbstep4[b_, bzv_, b3v_, k1_, k2_, j_] := b*bbstep3[bzv/b, b3v/b, k1, k2, m];

bbstep3[a1_, a2_, k1_, k2_, j_] := 
 Module [{b1, b2, s1, s2, r},
	b1=If[a1!=0, 1/a1,0];
	b2=If[a2!=0, 1/a2, 0];	
	s1=IntegerPart[k1*b1^j];
	s2 = IntegerPart[k2*b2^j]; b1=b1+s1; b2=b2+s2;
	r={0,0};
	r[[1]]=If[b1!=0, 1/b1, 0];
	r[[2]]=If[b2!=0, 1/b2, 0];
	Return[r];
	]

systep[insy_, bza, b3a_, ffa_, acc00_] :=
 Module [ {syv},
	ff=ffa/. {bz->bza, b3->b3a}; 
	syv = FRTB[ff,insy, acc];
	syv = N[syv,acc-3]; 
	Return[syv]; 
	]

end = lastline