procedure vsqr(var t: scalar; var v1,v2: vector);
var
  i,j,k: integer;
  s,s1,s2,s3,sf,ss: scalar;
  u1: vector;
begin { vsqr }
  write('(vsqr'); flush(output);
  for i:=1 to vdim-1 do begin                     { p*p -> p }
    k:=(i-1) div 2;
    szero(ss);
    for j:=0 to k do begin
      sprod(v1[j],v1[i-j],s1);
      ssum(s1,ss,s2);
      ss:=s2;
    end;
    ismult(2,s2);
    if odd(i) then v2[i]:=s2
    else begin
      sprod(v1[k+1],v1[k+1],s1);
      ssum(s1,s2,v2[i]);
    end;
  end;
  sprod(v1[0],v1[0],v2[0]);

  for i:=0 to vdim-1 do sabs0(v1[i],u1[i]);
  szero(u1[vdim]);

  s:=t;
  isdiv(2,s);
  szero(s3);
  sf.u:=sone.u;
  sf.l.si:=0;
  for i:=vdim to 2*(vdim-1) do begin            { p*p -> h }
    k:=(i-1) div 2;
    szero(ss);
    for j:=i-vdim to k do begin
      sprod(u1[j],u1[i-j],s1);
      ssum(s1,ss,s2);
      ss:=s2;
    end;
    ismult(2,s2);
    if odd(i) then ss:=s2
    else begin
      sprod(u1[k+1],u1[k+1],s1);
      ssum(s1,s2,ss);
    end;
    sprod(ss,sf,s1);
    sunion(s3,s1,s2);
    s3:=s2;
    ismult(i,sf);
    sprod(sf,s,s1);
    sf:=s1;
  end;

  vnorm(u1,t,s1);                              { p*h+h*p -> h }
  ispower(vdim,stwo,s2);
  sprod(s2,v1[vdim],sf);
  sprod(s1,sf,s2);
  ismult(2,s2);
  ssum(s3,s2,s1);

  s3:=s;
  ismult(4,s3);
  ispower(vdim,s3,s2);
  isfactorial(vdim,sf);
  sprod(s2,sf,s3);
  sprod(v1[vdim],v1[vdim],s2);
  sprod(s3,s2,sf);                                { h*h -> h }

  sunion(sf,s1,s2);
  sabs0(s2,v2[vdim]);
  write(')'); flush(output);
end { vsqr };