procedure vconvinit(var c,r: scalar; var vk: vector);
var
  i: integer;
  s,si,ri,sc,sq,sr,st,s1,s2: scalar;
  vt0,vt1,vt2: vector; 

begin { vconvinit };
  write('(vconvinit'); flush(output);
  s1:=c;
  ismult(2,s1);
  ssum(r,s1,s);                                    { s=r+2c }
  squot(sone,s,si);                                { si=1/s }
  squot(sone,r,ri);                                { ri=1/r }
  squot(c,s,sc);                                   { sc=c/s }
  squot(s,r,sr);                                   { sr=s/r }
  ssqrt(sr,sq);                              { sq=sqrt(s/r) }

  vt0[0]:=sone;
  for i:=1 to vdim-1 do begin
    vt0[i]:=vt0[i-1];
    isdiv(i,vt0[i]);
  end;
  szero(vt0[vdim]);

  vscale(si,vt0,vt1);                  { vt1[k]=(si^k)/(k!) }
  vconv(vdim,c,vt0,vt1,vt2);   
  vscale(ri,vt0,vt1);
  vsprod(sq,vt1,vt0);               { vt0[k]=sq*(ri^k)/(k!) }
  vdiff(vt0,vt2,vt1);
  st:=s;
  st.l.si:=0;
  ispower(vdim,st,s1);
  for i:=2 to vdim do ismult(i,s1);
  vsprod(s1,vt1,vt2);                         { first bound }

  st.u:=sone.u;
  st.l.si:=0;
  for i:=vdim-1 downto 0 do begin
    ismult(i+1,st);
    ismult(2*(i+1)-1,st);
    isdiv(vdim-i,st);
    sprod(st,c,s1);
    st:=s1;
    s1:=sc;
    if i=0 then ismult(2,s1)
    else begin
      ismult(2*vdim+1,s1);
      isdiv(vdim-i+1,s1);
    end;
    sdiff(sone,s1,s2);
    if s2.l.si>0 then begin
      squot(sone,s2,s1);
      sprod(s1,st,vt1[i]);                   { second bound }
    end
    else vt1[i]:=vt2[i];
  end;

  for i:=0 to vdim-1 do begin
    sinter(vt1[i],vt2[i],vk[i]);    { min of the two bounds }
    with vk[i] do begin
      l:=u;
      l.si:=-u.si;
    end;
  end;

  ispower(vdim,sr,s2);
  sprod(sq,s2,vk[vdim]);
  vk[vdim].l.si:=0;

  write(')'); flush(output);
end { vconvinit };