procedure vprod(var t1,t2: scalar; var v1,v2,v3: vector);
var
  ho1,ho2,teq: boolean;
  i,j: integer;
  s,s1,s2,s3,sf,sp,ss,ts: scalar;
begin { vprod };
  write('(vprod'); flush(output);

  ho1:=not seq0(v1[vdim]);
  ho2:=not seq0(v2[vdim]);
  teq:=(addr(t1)=addr(t2));

  if teq then begin
    s:=t1;
    isdiv(2,s);
    if ho1 or ho2 then ispower(vdim,stwo,sp);
  end
  else begin
    squot(t2,t1,s1);
    ssum(sone,s1,s2);
    squot(t2,s2,s);
    ssum(t1,t2,ts);
  end;

  pprod(vdim,v1,v2,v3);                           { p*p -> p }
  szero(s3);

  sf:=sone;
  for i:=vdim to 2*(vdim-1) do begin              { p*p -> h }
    szero(ss);
    for j:=i-vdim to vdim-1 do begin
      sprod(v1[j],v2[i-j],s1);
      ssum(s1,ss,s2);
      ss:=s2;
    end;
    sprod(ss,sf,s1);
    sunion(s3,s1,s2);
    s3:=s2;
    ismult(i,sf);
    sprod(sf,s,s1);
    sf:=s1;
  end;
  sabs0(s2,s3);

  if ho1 then begin                             { h*(p+h) -> h }
    if not teq then begin
      squot(ts,t2,s1);
      ispower(vdim,s1,sp);
    end;      
    vnorm(v2,t2,s1);
    sprod(sp,v1[vdim],sf);
    sprod(s1,sf,s2);
    ssum(s3,s2,s1);
    sabs0(s1,s3);
  end;

  if ho2 then begin                                 { p*h -> h }
    if not teq then begin
      squot(ts,t1,s1);
      ispower(vdim,s1,sp);
    end;
    ss:=v1[vdim];             
    szero(v1[vdim]);
    vnorm(v1,t1,s1);
    v1[vdim]:=ss;
    sprod(sp,v2[vdim],sf);
    sprod(s1,sf,s2);
    ssum(s3,s2,s1);
    sabs0(s1,s3);
  end;

  v3[vdim]:=s3;
  write(')'); flush(output);
end { vprod };