proc(B,n,emin,emax,C,ndigits)

local epsilonmin,powerofBoverC,e,a,Plast,r,Qlast,

Q,P,NewQ,NewP,epsilon,

numbermin,expmin,l;

epsilonmin := 12345.0 ;

Digits := ndigits;

powerofBoverC := B^(emin-n)/C;

for e from emin-n+1 to emax-n+1 do

powerofBoverC := B*powerofBoverC;

a := floor(powerofBoverC);

Plast := a;

r := 1/(powerofBoverC-a);

a := floor(r);

Qlast := 1;

Q := a;

P := Plast*a+1;

while Q < B^n-1 do

r := 1/(r-a);

a := floor(r);

NewQ := Q*a+Qlast;

NewP := P*a+Plast;

Qlast := Q;

Plast := P;

Q := NewQ;

P := NewP

od;

epsilon :=

evalf(C*abs(Plast-Qlast*powerofBoverC));

print(e+n-1);

if epsilon < epsilonmin then

epsilonmin := epsilon; numbermin := Qlast;

expmin := e

fi

od;

print('mantissa',numbermin);

print('exponent',expmin);

print('epsilon',epsilonmin);

l := evalf(log(epsilonmin)/log(B),10);

print(numberofdigits,l)

end