subroutine mccopy (x, z) multiple precision x(2), z(2) z(1) = x(1) z(2) = x(2) return end subroutine mcalc (ix) integer ix(2) call malc (ix(1)) call malc (ix(2)) return end subroutine mcdalc (ix) integer ix(2) call mdalc (ix(2)) call mdalc (ix(1)) return end subroutine mcadd (x, y, z) multiple precision x(2), y(2), z(2) z(1) = x(1) + y(1) z(2) = x(2) + y(2) return end subroutine mcsub (x, y, z) multiple precision x(2), y(2), z(2) z(1) = x(1) - y(1) z(2) = x(2) - y(2) return end subroutine mcmul1 (intger, y, z) multiple precision y(2), z(2) z(1) = intger*y(1) z(2) = intger*y(2) return end subroutine mcmul2 (x, intger, z) multiple precision x(2), z(2) z(1) = intger*x(1) z(2) = intger*x(2) return end subroutine mcmul (x, y, z) multiple precision x(2), y(2), z(2), tmp tmp = x(1)*y(1) - x(2)*y(2) z(2) = x(1)*y(2) + x(2)*y(1) z(1) = tmp return end subroutine mcdiv (x, y, z) multiple precision x(2), y(2), z(2), tmp1, tmp2 tmp1 = rec (y(1)*y(1)+y(2)*y(2)) tmp2 = (x(1)*y(1)+x(2)*y(2))*tmp1 z(2) = (-x(1)*y(2)+x(2)*y(1))*tmp1 z(1) = tmp2 return end subroutine mcdivi (x, intger, z) multiple precision x(2), z(2) z(1) = x(1)/intger z(2) = x(2)/intger return end multiple complex function mccexi (x, n) multiple complex x, sq c n2 = n if (n2) 20, 10, 40 c 10 mccexi = 1 return c 20 n2 = -n2 if (x.ne.0) go to 60 call seterr (30hmccexi zero to negative power, 30, 1, 2) stop c 40 if (x.ne.0) go to 60 mccexi = 0 return c 60 sq = x if (n.lt.0) sq = 1/sq mccexi = 1 c 70 ns = n2 n2 = n2/2 if (2*n2.ne.ns) mccexi = sq*mccexi if (n2.le.0) return sq = sq*sq go to 70 c end multiple complex function mccexc (x, y) multiple complex x, y c if (y.ne.0.) go to 10 mccexc = 1 return c 10 if (x.ne.0.) go to 20 mccexc = 0 return c 20 mccexc = cexp (y*clog(x)) return c end subroutine mcitoc (intger, z) multiple precision z(2) z(1) = intger z(2) = 0 return end subroutine mcetoc (r, z) multiple precision z(2) z(1) = r z(2) = 0 return end subroutine mcdtoc (d, z) multiple precision z(2) double precision d z(1) = d z(2) = 0 return end subroutine mcctoc (c, z) multiple precision z(1) complex c z(1) = real(c) z(2) = aimag(c) return end logical function mceq (x, y) multiple precision x(2), y(2) logical mne c mceq = .false. if (mne(x(1),y(1))) return if (mne(x(2),y(2))) return mceq = .true. c return end logical function mcne (x, y) multiple precision x(2), y(2) logical meq c mcne = .false. if (meq(x(1),y(1))) return if (meq(x(2),y(2))) return mcne = .true. c return end multiple precision function mcreal (z) multiple precision z(2) mcreal = z(1) return end multiple precision function mcimag (z) multiple precision z(2) mcimag = z(2) return end subroutine mccmpl (x, y, z) multiple precision x, y, z(2) z(1) = x z(2) = y return end multiple complex function mcexp (z) multiple complex z multiple precision r, y c r = dexp (dreal(z)) if (r.eq.0) go to 10 c y = dimag (z) mcexp = cmplx (r*dcos(y), r*dsin(y)) return c 10 mcexp = 0 return c end multiple complex function mclog (z) multiple complex z c mclog = cmplx (dlog(cabs(z)), carg(z)) c return end multiple precision function mcabs (z) multiple complex z multiple precision x, y, r1, r2 c x = dabs(dreal(z)) y = dabs (dimag(z)) r1 = dmin1 (x, y) r2 = dmax1 (x, y) c mcabs = r2 if (r1.ne.0.d0) mcabs = r2*dsqrt(1.d0+(r1/r2)**2) c return end multiple precision function mcarg (z) multiple complex z c mcarg = 0 if (dreal(z).ne.0) go to 10 if (dimag(z).eq.0) return c 10 mcarg = datan2 (dimag(z), dreal(z)) c return end multiple complex function mcsqrt (z) multiple complex z multiple precision x, y, r, xtmp, ytmp c x = dreal (z) y = dimag (z) r = cabs (z) c if (r.ne.0) go to 10 mcsqrt = 0 return c 10 xtmp = dsqrt ((r+dabs(x))/2) ytmp = y/(2*xtmp) c if (x.ge.0) mcsqrt = cmplx (xtmp, ytmp) if (y.eq.0) y = 1 if (x.lt.0) mcsqrt = cmplx (dabs(ytmp), dsign(xtmp,y)) c return end