lineintegralC3(f,df,R,N):= /* 線分 [0,Ri] に沿う線積分 */
  block (s:0, z1:0, dz:R*%i/N,
    for j: 1 thru N
      do (w:rectform(df(z1)/f(z1)),z2:z1+dz,
          s:rectform(s+w*dz), z1:z2), 
    return(s));
lineintegralC4(f,df,R,N,pi):= /* 4半円周 Re^{it}, pi/2<t<pi に沿う線積分 */
  block(s:0, t:pi/2, dt:pi/2/N, z1:float(R*cos(t)+R*%i*sin(t)),
    for j: 1 thru N
      do (w:rectform(df(z1)/f(z1)),t:t+dt,
        z2:float(R*cos(t)+R*%i*sin(t)),s:rectform(s+w*(z2-z1)),z1:z2),
    return(s));
lineintegralC5(f,df,R,N):= /* 線分 [-R,0] に沿う線積分 */
  block (s:0, z1:-R, dz:R/N,
    for j: 1 thru N
      do (w:rectform(df(z1)/f(z1)),z2:z1+dz,
          s:rectform(s+w*(z2-z1)), z1:z2), 
    return(s));
lineintegralS1(f,df,R,N,pi):=
  block(s1:lineintegralC3(f,df,R,N),
        s2:lineintegralC4(f,df,R,N,pi),
        s3:lineintegralC5(f,df,R,N),
        s:rectform((s1+s2+s3)/2/pi/%i),
        return(s));
f(z):=z^6+z^4+3*z^2+2*z+2$
df(_z):=''(diff(f(_z),_z))$
R:2.0;
N:100;
pi:3.1416;
lineintegralS1(f,df,R,N,pi);


