module ddfunmod
contains
subroutine ddabrt
implicit none
stop
end subroutine
subroutine ddadd (dda, ddb, ddc)
implicit none
real8 dda(2), ddb(2), ddc(2)
real8 e, t1, t2
t1 = dda(1) + ddb(1)
e = t1 - dda(1)
t2 = ((ddb(1) - e) + (dda(1) - (t1 - e))) + dda(2) + ddb(2)
ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine
subroutine ddang (x, y, a)
implicit none
integer i, ix, iy, k, kk, nx, ny
real8 t1, t2, t3
real8 a(2), pi(2), x(2), y(2), s0(2), s1(2), s2(2), s3(2), s4(2)
save pi
data pi/ z’400921FB54442D18’, z’3CA1A62633145C07’/
if (x(1) .eq. 0.d0 .and. y(1) .eq. 0.d0) then
write (6, 1)
1 format (‘*** DDANG: Both arguments are zero.’)
call ddabrt
return
endif
if (x(1) .eq. 0.d0) then
if (y(1) .gt. 0.d0) then
call ddmuld (pi, 0.5d0, a)
else
call ddmuld (pi, -0.5d0, a)
endif
goto 120
elseif (y(1) .eq. 0.d0) then
if (x(1) .gt. 0.d0) then
a(1) = 0.d0
a(2) = 0.d0
else
a(1) = pi(1)
a(2) = pi(2)
endif
goto 120
endif
call ddmul (x, x, s0)
call ddmul (y, y, s1)
call ddadd (s0, s1, s2)
call ddsqrt (s2, s3)
call dddiv (x, s3, s1)
call dddiv (y, s3, s2)
call ddqdc (s1, t1)
call ddqdc (s2, t2)
t3 = atan2 (t2, t1)
a(1) = t3
a(2) = 0.d0
if (abs (t1) .le. abs (t2)) then
kk = 1
s0(1) = s1(1)
s0(2) = s1(2)
else
kk = 2
s0(1) = s2(1)
s0(2) = s2(2)
endif
do k = 1, 3
call ddcssn (a, s1, s2)
if (kk .eq. 1) then
call ddsub (s0, s1, s3)
call dddiv (s3, s2, s4)
call ddsub (a, s4, s1)
else
call ddsub (s0, s2, s3)
call dddiv (s3, s1, s4)
call ddadd (a, s4, s1)
endif
a(1) = s1(1)
a(2) = s1(2)
enddo
120 continue
return
end subroutine
subroutine ddcadd (a, b, c)
implicit none
real*8 a(4), b(4), c(4)
call ddadd (a, b, c)
call ddadd (a(3), b(3), c(3))
return
end subroutine
subroutine ddcdiv (a, b, c)
implicit none
real*8 a(4), b(4), c(4), f(2), s0(2), s1(2), s2(2), s3(2), s4(2)
if (b(1) .eq. 0.d0 .and. b(3) .eq. 0.d0) then
write (6, 1)
1 format (‘*** DDCDIV: Divisor is zero.’)
call ddabrt
return
endif
f(1) = 1.d0
f(2) = 0.d0
call ddmul (a, b, s0)
call ddmul (a(3), b(3), s1)
call ddadd (s0, s1, s2)
call ddsub (s0, s1, s3)
call ddadd (a, a(3), s0)
call ddsub (b, b(3), s1)
call ddmul (s0, s1, s4)
call ddsub (s4, s3, s1)
call ddmul (b, b, s0)
call ddmul (b(3), b(3), s3)
call ddadd (s0, s3, s4)
call dddiv (f, s4, s0)
call ddmul (s2, s0, c)
call ddmul (s1, s0, c(3))
return
end subroutine
subroutine ddceq (a, b)
implicit none
real*8 a(4), b(4)
b(1) = a(1)
b(2) = a(2)
b(3) = a(3)
b(4) = a(4)
return
end subroutine
subroutine ddcmul (a, b, c)
implicit none
real*8 a(4), b(4), c(4), s0(2), s1(2), s2(2), s3(2)
call ddmul (a, b, s0)
call ddmul (a(3), b(3), s1)
call ddmul (a, b(3), s2)
call ddmul (a(3), b, s3)
call ddsub (s0, s1, c)
call ddadd (s2, s3, c(3))
return
end subroutine
subroutine ddcpr (a, b, ic)
implicit none
integer ic
real*8 a(2), b(2)
if (a(1) .lt. b(1)) then
ic = -1
elseif (a(1) .eq. b(1)) then
if (a(2) .lt. b(2)) then
ic = -1
elseif (a(2) .eq. b(2)) then
ic = 0
else
ic = 1
endif
else
ic = 1
endif
return
end subroutine
subroutine ddcpwr (a, n, b)
implicit none
integer j, kk, kn, l1, mn, n, na1, na2, nn
real8 cl2, t1
parameter (cl2 = 1.4426950408889633d0)
real8 a(4), b(4), s0(4), s1(4), s2(4), s3(4)
if (a(1) .eq. 0.d0 .and. a(3) .eq. 0.d0) then
if (n .ge. 0) then
b(1) = 0.d0
b(2) = 0.d0
b(3) = 0.d0
b(4) = 0.d0
goto 120
else
write (6, 1)
1 format (‘*** DDCPWR: Argument is zero and N is negative or zero.’)
call ddabrt
return
endif
endif
nn = abs (n)
if (nn .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
s2(3) = 0.d0
s2(4) = 0.d0
goto 120
elseif (nn .eq. 1) then
s2(1) = a(1)
s2(2) = a(2)
s2(3) = a(3)
s2(4) = a(4)
goto 110
elseif (nn .eq. 2) then
call ddcmul (a, a, s2)
goto 110
endif
t1 = nn
mn = cl2 * log (t1) + 1.d0 + 1.d-14
s0(1) = a(1)
s0(2) = a(2)
s0(3) = a(3)
s0(4) = a(4)
s2(1) = 1.d0
s2(2) = 0.d0
s2(3) = 0.d0
s2(4) = 0.d0
kn = nn
do j = 1, mn
kk = kn / 2
if (kn .ne. 2 * kk) then
call ddcmul (s2, s0, s1)
s2(1) = s1(1)
s2(2) = s1(2)
s2(3) = s1(3)
s2(4) = s2(4)
endif
kn = kk
if (j .lt. mn) then
call ddcmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
s0(3) = s1(3)
s0(4) = s1(4)
endif
enddo
110 continue
if (n .lt. 0) then
s1(1) = 1.d0
s1(2) = 0.d0
s1(3) = 0.d0
s1(4) = 0.d0
call ddcdiv (s1, s2, s0)
s2(1) = s0(1)
s2(2) = s0(2)
s2(3) = s0(3)
s2(4) = s0(4)
endif
b(1) = s2(1)
b(2) = s2(2)
b(3) = s2(3)
b(4) = s2(4)
120 continue
return
end subroutine
subroutine ddcsqrt (a, b)
implicit none
real*8 a(4), b(4), s0(2), s1(2), s2(2)
if (a(1) .eq. 0.d0 .and. a(3) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
b(3) = 0.d0
b(4) = 0.d0
goto 100
endif
call ddmul (a, a, s0)
call ddmul (a(3), a(3), s1)
call ddadd (s0, s1, s2)
call ddsqrt (s2, s0)
s1(1) = a(1)
s1(2) = a(2)
if (s1(1) .lt. 0.d0) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif
s1 = abs (s1)
call ddadd (s0, s1, s2)
call ddmuld (s2, 0.5d0, s1)
call ddsqrt (s1, s0)
call ddmuld (s0, 2.d0, s1)
if (a(1) .ge. 0.d0) then
b(1) = s0(1)
b(2) = s0(2)
call dddiv (a(3), s1, b(3))
else
call dddiv (a(3), s1, b)
if (b(1) .lt. 0.d0) then
b(1) = - b(1)
b(2) = - b(2)
endif
b(3) = s0(1)
b(4) = s0(2)
if (a(3) .lt. 0.d0) then
b(3) = - b(3)
b(4) = - b(4)
endif
endif
100 continue
return
end subroutine
subroutine ddcssh (a, x, y)
implicit none
real*8 a(2), f(2), x(2), y(2), s0(2), s1(2), s2(2)
f(1) = 1.d0
f(2) = 0.d0
call ddexp (a, s0)
call dddiv (f, s0, s1)
call ddadd (s0, s1, s2)
call ddmuld (s2, 0.5d0, x)
call ddsub (s0, s1, s2)
call ddmuld (s2, 0.5d0, y)
return
end subroutine
subroutine ddcssn (a, x, y)
implicit none
integer ia, ka, kb, kc, l1
real8 t1, t2
real8 a(2), f(2), pi(2), x(2), y(2), s0(2), s1(2), s2(2), s3(2), s4(2), &
s5(2), s6(2)
real*8 cs(2,2,4)
save cs, pi
data cs/ &
z’3FEF6297CFF75CB0’, z’3C7562172A361FD3’, &
z’3FC8F8B83C69A60B’, z’BC626D19B9FF8D82’, &
z’3FED906BCF328D46’, z’3C7457E610231AC2’, &
z’3FD87DE2A6AEA963’, z’BC672CEDD3D5A610’, &
z’3FEA9B66290EA1A3’, z’3C39F630E8B6DAC8’, &
z’3FE1C73B39AE68C8’, z’3C8B25DD267F6600’, &
z’3FE6A09E667F3BCD’, z’BC8BDD3413B26456’, &
z’3FE6A09E667F3BCD’, z’BC8BDD3413B26456’ /
data pi/ z’400921FB54442D18’, z’3CA1A62633145C07’/
if (a(1) .eq. 0.d0) then
x(1) = 1.d0
x(2) = 0.d0
y(1) = 0.d0
y(2) = 0.d0
goto 120
endif
f(1) = 1.d0
f(2) = 0.d0
call ddmuld (pi, 2.d0, s0)
call dddiv (a, s0, s1)
call ddnint (s1, s2)
call ddsub (s1, s2, s3)
t1 = s3(1)
t2 = 4.d0 * t1
ka = nint (t2)
kb = nint (8.d0 * (t2 - ka))
t1 = (8 * ka + kb) / 32.d0
s1(1) = t1
s1(2) = 0.d0
call ddsub (s3, s1, s2)
call ddmul (s0, s2, s1)
if (s1(1) .eq. 0.d0) then
s0(1) = 0.d0
s0(2) = 0.d0
goto 110
endif
s0(1) = s1(1)
s0(2) = s1(2)
call ddmul (s0, s0, s2)
l1 = 0
100 l1 = l1 + 1
if (l1 .eq. 100) then
write (6, 1)
1 format (‘*** DDCSSN: Iteration limit exceeded.’)
call ddabrt
return
endif
t2 = - (2.d0 * l1) * (2.d0 * l1 + 1.d0)
call ddmul (s2, s1, s3)
call dddivd (s3, t2, s1)
call ddadd (s1, s0, s3)
s0(1) = s3(1)
s0(2) = s3(2)
if (abs (s1(1)) .gt. 1d-33 * abs (s3(1))) goto 100
110 continue
s1(1) = s0(1)
s1(2) = s0(2)
call ddmul (s0, s0, s2)
call ddsub (f, s2, s3)
call ddsqrt (s3, s0)
kc = abs (kb)
f(1) = 2.
if (kc .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
s3(1) = 0.d0
s3(2) = 0.d0
else
s2(1) = cs(1,1,kc)
s2(2) = cs(2,1,kc)
s3(1) = cs(1,2,kc)
s3(2) = cs(2,2,kc)
endif
if (kb .lt. 0) then
s3(1) = - s3(1)
s3(2) = - s3(2)
endif
call ddmul (s0, s2, s4)
call ddmul (s1, s3, s5)
call ddsub (s4, s5, s6)
call ddmul (s1, s2, s4)
call ddmul (s0, s3, s5)
call ddadd (s4, s5, s1)
s0(1) = s6(1)
s0(2) = s6(2)
if (ka .eq. 0) then
x(1) = s0(1)
x(2) = s0(2)
y(1) = s1(1)
y(2) = s1(2)
elseif (ka .eq. 1) then
x(1) = - s1(1)
x(2) = - s1(2)
y(1) = s0(1)
y(2) = s0(2)
elseif (ka .eq. -1) then
x(1) = s1(1)
x(2) = s1(2)
y(1) = - s0(1)
y(2) = - s0(2)
elseif (ka .eq. 2 .or. ka .eq. -2) then
x(1) = - s0(1)
x(2) = - s0(2)
y(1) = - s1(1)
y(2) = - s1(2)
endif
120 continue
return
end subroutine
subroutine ddcsub (a, b, c)
implicit none
real*8 a(4), b(4), c(4)
call ddsub (a, b, c)
call ddsub (a(3), b(3), c(3))
return
end subroutine
subroutine dddiv (dda, ddb, ddc)
implicit none
real8 dda(2), ddb(2), ddc(2)
real8 a1, a2, b1, b2, cona, conb, c11, c2, c21, e, split, s1, s2, &
t1, t2, t11, t12, t21, t22
parameter (split = 134217729.d0)
s1 = dda(1) / ddb(1)
cona = s1 * split
conb = ddb(1) * split
a1 = cona - (cona - s1)
b1 = conb - (conb - ddb(1))
a2 = s1 - a1
b2 = ddb(1) - b1
c11 = s1 * ddb(1)
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2
c2 = s1 * ddb(2)
t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21
t12 = t1 + t2
t22 = t2 - (t12 - t1)
t11 = dda(1) - t12
e = t11 - dda(1)
t21 = ((-t12 - e) + (dda(1) - (t11 - e))) + dda(2) - t22
s2 = (t11 + t21) / ddb(1)
ddc(1) = s1 + s2
ddc(2) = s2 - (ddc(1) - s1)
return
end subroutine
subroutine dddivd (dda, db, ddc)
implicit none
real8 dda(2), db, ddc(2)
real8 a1, a2, b1, b2, cona, conb, e, split, t1, t2, t11, t12, t21, t22
parameter (split = 134217729.d0)
t1 = dda(1) / db
cona = t1 * split
conb = db * split
a1 = cona - (cona - t1)
b1 = conb - (conb - db)
a2 = t1 - a1
b2 = db - b1
t12 = t1 * db
t22 = (((a1 * b1 - t12) + a1 * b2) + a2 * b1) + a2 * b2
t11 = dda(1) - t12
e = t11 - dda(1)
t21 = ((-t12 - e) + (dda(1) - (t11 - e))) + dda(2) - t22
t2 = (t11 + t21) / db
ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine
subroutine dddqc (a, b)
implicit none
real*8 a, b(2)
b(1) = a
b(2) = 0.d0
return
end subroutine
subroutine ddeq (a, b)
implicit none
real*8 a(2), b(2)
b(1) = a(1)
b(2) = a(2)
end subroutine
subroutine ddeform (a, n1, n2, b)
implicit none
integer i, ln, m1, n1, n2
parameter (ln = 40)
character1 b(n1)
character40 cs
real*8 a(2)
if (n1 .lt. 0 .or. n2 .lt. 0 .or. n1 .gt. 80 .or. n2 .gt. 30 &
.or. n2 .gt. n1 - 8) then
write (6, 1) n1, n2
1 format (‘*** DDEFORM: Improper n1, n2 =’,2i6)
call ddabrt
endif
call ddoutc (a, cs)
m1 = n1 - n2 - 8
do i = 1, m1
b(i) = ’ ’
enddo
do i = 1, n2 + 3
b(i+m1) = cs(i+2:i+2)
enddo
do i = 1, 5
b(i+m1+n2+3) = cs(i+35:i+35)
enddo
return
end subroutine
subroutine ddexp (a, b)
implicit none
integer i, ia, l1, na, nq, nz, n1
real8 t1, t2
parameter (nq = 6)
real8 a(2), b(2), al2(2), f(2), s0(2), s1(2), s2(2), s3(2), tl
save al2
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.
data al2 / z’3FE62E42FEFA39EF’, z’3C7ABC9E3B39803F’/
! data al2/ 6.9314718055994529D-01, 2.3190468138462996D-17/
! Check for overflows and underflows.
if (abs (a(1)) .ge. 709.d0) then
if (a(1) .gt. 0.d0) then
write (6, 1) a(1)
1 format (‘*** DDEXP: Argument is too large’,f12.6)
call ddabrt
return
else
call dddqc (0.d0, b)
goto 130
endif
endif
f(1) = 1.d0
f(2) = 0.d0
! Compute the reduced argument A’ = A - Log(2) * Nint [A / Log(2)]. Save
! NZ = Nint [A / Log(2)] for correcting the exponent of the final result.
call dddiv (a, al2, s0)
call ddnint (s0, s1)
t1 = s1(1)
nz = t1 + sign (1.d-14, t1)
call ddmul (al2, s1, s2)
call ddsub (a, s2, s0)
! Check if the reduced argument is zero.
if (s0(1) .eq. 0.d0) then
s0(1) = 1.d0
s0(2) = 0.d0
l1 = 0
goto 120
endif
! Divide the reduced argument by 2 ^ NQ.
call dddivd (s0, 2.d0 ** nq, s1)
! Compute Exp using the usual Taylor series.
s2(1) = 1.d0
s2(2) = 0.d0
s3(1) = 1.d0
s3(2) = 0.d0
l1 = 0
100 l1 = l1 + 1
if (l1 .eq. 100) then
write (6, 2)
2 format (‘*** DDEXP: Iteration limit exceeded.’)
call ddabrt
return
endif
t2 = l1
call ddmul (s2, s1, s0)
call dddivd (s0, t2, s2)
call ddadd (s3, s2, s0)
call ddeq (s0, s3)
! Check for convergence of the series.
if (abs (s2(1)) .gt. 1d-33 * abs (s3(1))) goto 100
! Raise to the (2 ^ NQ)-th power.
do i = 1, nq
call ddmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
enddo
! Multiply by 2 ^ NZ.
120 call ddmuld (s0, 2.d0 ** nz, b)
! Restore original precision level.
130 continue
return
end subroutine
subroutine ddfform (a, n1, n2, b)
! This routine converts the DD number A to F format, i.e. Fn1.n2.
! B is the output array (type character*1 of size n1). N1 must exceed
! N2 by at least 3, and N1 must not exceed 80. N2 must not exceed 30.
implicit none
integer i, ix, kx, ln, ls, lz, mx, nx, n1, n2
parameter (ln = 40)
real8 a(2)
character1 b(n1)
character40 c
character40 chr40
if (n1 .lt. 0 .or. n2 .lt. 0 .or. n1 .gt. 80 .or. n2 .gt. 30 &
.or. n1 - n2 .lt. 3) then
write (6, 1) n1, n2
1 format (‘*** DDFFORM: Improper n1, n2 =’,2i6)
call ddabrt
endif
! Call ddoutc and extract exponent.
call ddoutc (a, c)
ix = dddigin (c(ln-3:ln), 4)
do i = 1, n1
b(i) = ’ ’
enddo
if (a(1) .ge. 0.d0) then
ls = 0
else
ls = 1
endif
if (ix .ge. 0 .and. a(1) .ne. 0.d0) then
lz = 0
else
lz = 1
endif
mx = max (ix, 0)
! Check for overflow of field length.
if (ls + lz + mx + n2 + 2 .gt. n1) then
do i = 1, n1
b(i) = ‘*’
enddo
goto 200
endif
! Check if a zero should be output.
if (a(1) .eq. 0. .or. -ix .gt. n2) then
do i = 1, n1 - n2 - 2
b(i) = ’ ’
enddo
b(n1-n2-1) = ‘0’
b(n1-n2) = ‘.’
do i = 1, n2
b(i+n1-n2) = ‘0’
enddo
goto 200
endif
! Process other cases.
if (a(1) .lt. 0.) b(n1-n2-mx-2) = ‘-’
if (ix .ge. 0) then
b(n1-n2-ix-1) = c(4:4)
kx = min (ln - 9, ix)
do i = 1, kx
b(i+n1-n2-ix-1) = c(i+5:i+5)
enddo
do i = kx + 1, ix
b(i+n1-n2-ix-1) = ‘0’
enddo
b(n1-n2) = ‘.’
kx = max (min (ln - 9 - ix, n2), 0)
do i = 1, kx
b(i+n1-n2) = c(i+ix+5:i+ix+5)
enddo
do i = kx + 1, n2
b(i+n1-n2) = ‘0’
enddo
else
nx = - ix
b(n1-n2-1) = ‘0’
b(n1-n2) = ‘.’
do i = 1, nx - 1
b(i+n1-n2) = ‘0’
enddo
b(n1-n2+nx) = c(4:4)
kx = min (ln - 8, n2 - nx)
do i = 1, kx
b(i+n1-n2+nx) = c(i+5:i+5)
enddo
do i = kx + 1, n2 - nx
b(i+n1-n2+nx) = ‘0’
enddo
endif
200 continue
return
end subroutine
subroutine ddinfr (a, b, c)
! Sets B to the integer part of the DD number A and sets C equal to the
! fractional part of A. Note that if A = -3.3, then B = -3 and C = -0.3.
implicit none
integer ic
real*8 a(2), b(2), c(2), con(2), f(2), s0(2), s1(2), t105, t52
parameter (t105 = 2.d0 ** 105, t52 = 2.d0 ** 52)
save con
data con / t105, t52/
! Check if A is zero.
if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
c(1) = 0.d0
c(2) = 0.d0
goto 120
endif
if (a(1) .ge. t105) then
write (6, 1)
1 format (‘*** DDINFR: Argument is too large.’)
call ddabrt
return
endif
f(1) = 1.d0
f(2) = 0.d0
if (a(1) .gt. 0.d0) then
call ddadd (a, con, s0)
call ddsub (s0, con, b)
call ddcpr (a, b, ic)
if (ic .ge. 0) then
call ddsub (a, b, c)
else
call ddsub (b, f, s1)
b(1) = s1(1)
b(2) = s1(2)
call ddsub (a, b, c)
endif
else
call ddsub (a, con, s0)
call ddadd (s0, con, b)
call ddcpr (a, b, ic)
if (ic .le. 0) then
call ddsub (a, b, c)
else
call ddadd (b, f, s1)
b(1) = s1(1)
b(2) = s1(2)
call ddsub (a, b, c)
endif
endif
120 continue
return
end subroutine
subroutine ddinp (iu, a)
! This routine reads the DD number A from logical unit IU. The input
! value must be placed on a single line of not more than 80 characters.
implicit none
integer iu, ln
parameter (ln = 80)
character80 cs
real8 a(2)
read (iu, ‘(a)’, end = 100) cs
call ddinpc (cs, a)
goto 110
100 write (6, 1)
1 format (‘*** DDINP: End-of-file encountered.’)
call ddabrt
110 return
end subroutine
subroutine ddinpc (a, b)
! Converts the CHARACTER*80 array A into the DD number B.
implicit none
integer i, ib, id, ie, inz, ip, is, ix, k, ln, lnn
parameter (ln = 80)
real8 bi
character80 a
character1 ai
character10 ca, dig
parameter (dig = ‘0123456789’)
real*8 b(2), f(2), s0(2), s1(2), s2(2)
id = 0
ip = -1
is = 0
inz = 0
s1(1) = 0.d0
s1(2) = 0.d0
do i = 80, 1, -1
if (a(i:i) /= ’ ') goto 90
enddo
90 continue
lnn = i
! Scan for digits, looking for the period also.
do i = 1, lnn
ai = a(i:i)
if (ai .eq. ’ ’ .and. id == 0) then
elseif (ai .eq. ‘.’) then
if (ip >= 0) goto 210
ip = id
inz = 1
elseif (ai .eq. ‘+’) then
if (id .ne. 0 .or. ip >= 0 .or. is .ne. 0) goto 210
is = 1
elseif (ai .eq. ‘-’) then
if (id .ne. 0 .or. ip >= 0 .or. is .ne. 0) goto 210
is = -1
elseif (ai .eq. ‘e’ .or. ai .eq. ‘E’ .or. ai .eq. ‘d’ .or. ai .eq. ‘D’) then
goto 100
elseif (index (dig, ai) .eq. 0) then
goto 210
else
! read (ai, ‘(f1.0)’) bi
bi = index (dig, ai) - 1
if (inz > 0 .or. bi > 0.d0) then
inz = 1
id = id + 1
call ddmuld (s1, 10.d0, s0)
f(1) = bi
f(2) = 0.d0
call dddqc (bi, f)
call ddadd (s0, f, s1)
endif
endif
enddo
100 continue
if (is .eq. -1) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif
k = i
if (ip == -1) ip = id
ie = 0
is = 0
ca = ’ ’
do i = k + 1, lnn
ai = a(i:i)
if (ai .eq. ’ ') then
elseif (ai .eq. ‘+’) then
if (ie .ne. 0 .or. is .ne. 0) goto 210
is = 1
elseif (ai .eq. ‘-’) then
if (ie .ne. 0 .or. is .ne. 0) goto 210
is = -1
elseif (index (dig, ai) .eq. 0) then
goto 210
else
ie = ie + 1
if (ie .gt. 3) goto 210
ca(ie:ie) = ai
endif
enddo
! read (ca, ‘(i4)’) ie
ie = dddigin (ca, 4)
if (is .eq. -1) ie = - ie
ie = ie + ip - id
s0(1) = 10.d0
s0(2) = 0.d0
call ddnpwr (s0, ie, s2)
call ddmul (s1, s2, b)
goto 220
210 write (6, 1)
1 format (‘*** DDINPC: Syntax error in literal string.’)
call ddabrt
220 continue
return
end subroutine
subroutine ddlog (a, b)
implicit none
integer k
real8 t1, t2
real8 a(2), al2(2), b(2), s0(2), s1(2), s2(2)
save al2
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.
data al2 / z’3FE62E42FEFA39EF’, z’3C7ABC9E3B39803F’/
! data al2/ 6.9314718055994529D-01, 2.3190468138462996D-17/
if (a(1) .le. 0.d0) then
write (6, 1)
1 format (‘*** DDLOG: Argument is less than or equal to zero.’)
call ddabrt
return
endif
! Compute initial approximation of Log (A).
t1 = a(1)
t2 = log (t1)
b(1) = t2
b(2) = 0.d0
! Perform the Newton-Raphson iteration described above.
do k = 1, 3
call ddexp (b, s0)
call ddsub (a, s0, s1)
call dddiv (s1, s0, s2)
call ddadd (b, s2, s1)
b(1) = s1(1)
b(2) = s1(2)
enddo
120 continue
return
end subroutine
subroutine ddqdc (a, b)
! This converts the DD number A to DP.
implicit none
real*8 a(2), b
b = a(1)
return
end subroutine
subroutine ddqqc (a, b, c)
! This converts DD numbers A and B to DDC form in C, i.e. C = A + B i.
implicit none
real*8 a(2), b(2), c(4)
c(1) = a(1)
c(2) = a(2)
c(3) = b(1)
c(4) = b(2)
return
end subroutine
subroutine ddmul (dda, ddb, ddc)
! This routine multiplies DD numbers DDA and DDB to yield the DD product DDC.
implicit none
real8 dda(2), ddb(2), ddc(2)
real8 a1, a2, b1, b2, cona, conb, c11, c21, c2, e, split, t1, t2
parameter (split = 134217729.d0)
cona = dda(1) * split
conb = ddb(1) * split
a1 = cona - (cona - dda(1))
b1 = conb - (conb - ddb(1))
a2 = dda(1) - a1
b2 = ddb(1) - b1
! Multilply dda(1) * ddb(1) using Dekker’s method.
c11 = dda(1) * ddb(1)
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2
!>
! Compute dda(1) * ddb(2) + dda(2) * ddb(1) (only high-order word is needed).
c2 = dda(1) * ddb(2) + dda(2) * ddb(1)
! Compute (c11, c21) + c2 using Knuth’s trick, also adding low-order product.
t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21 + dda(2) * ddb(2)
! The result is t1 + t2, after normalization.
ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine
subroutine ddmuld (dda, db, ddc)
! This routine multiplies the DD number DDA by the DP number DB to yield
! the DD product DDC.
implicit none
real8 dda(2), db, ddc(2)
real8 a1, a2, b1, b2, cona, conb, c11, c21, c2, e, split, t1, t2
parameter (split = 134217729.d0)
cona = dda(1) * split
conb = db * split
a1 = cona - (cona - dda(1))
b1 = conb - (conb - db)
a2 = dda(1) - a1
b2 = db - b1
! Multilply dda(1) * db using Dekker’s method.
c11 = dda(1) * db
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2
!>
! Compute dda(2) * db (only high-order word is needed).
c2 = dda(2) * db
! Compute (c11, c21) + c2 using Knuth’s trick.
t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21
! The result is t1 + t2, after normalization.
ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine
subroutine ddmuldd (da, db, ddc)
! This subroutine computes ddc = da x db.
implicit none
real*8 a1, a2, b1, b2, cona, conb, da, db, ddc(2), split, s1, s2
parameter (split = 134217729.d0)
cona = da * split
conb = db * split
a1 = cona - (cona - da)
b1 = conb - (conb - db)
a2 = da - a1
b2 = db - b1
! Multiply da * db using Dekker’s method.
s1 = da * db
s2 = (((a1 * b1 - s1) + a1 * b2) + a2 * b1) + a2 * b2
!>
ddc(1) = s1
ddc(2) = s2
return
end subroutine
subroutine ddnint (a, b)
! This sets B equal to the integer nearest to the DD number A.
implicit none
real*8 a(2), b(2), con(2), s0(2), t105, t52
parameter (t105 = 2.d0 ** 105, t52 = 2.d0 ** 52)
save con
data con / t105, t52/
! Check if A is zero.
if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 120
endif
if (a(1) .ge. t105) then
write (6, 1)
1 format (‘*** DDINFR: Argument is too large.’)
call ddabrt
return
endif
if (a(1) .gt. 0.d0) then
call ddadd (a, con, s0)
call ddsub (s0, con, b)
else
call ddsub (a, con, s0)
call ddadd (s0, con, b)
endif
120 continue
return
end subroutine
subroutine ddnpwr (a, n, b)
implicit none
integer j, kk, kn, l1, mn, n, na1, na2, nn
real8 cl2, t1
parameter (cl2 = 1.4426950408889633d0)
real8 a(2), b(2), s0(2), s1(2), s2(2), s3(2)
if (a(1) .eq. 0.d0) then
if (n .ge. 0) then
s2(1) = 0.d0
s2(2) = 0.d0
goto 120
else
write (6, 1)
1 format (‘*** DDCPWR: Argument is zero and N is negative or zero.’)
call ddabrt
return
endif
endif
nn = abs (n)
if (nn .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
goto 120
elseif (nn .eq. 1) then
s2(1) = a(1)
s2(2) = a(2)
goto 110
elseif (nn .eq. 2) then
call ddmul (a, a, s2)
goto 110
endif
! Determine the least integer MN such that 2 ^ MN .GT. NN.
t1 = nn
mn = cl2 * log (t1) + 1.d0 + 1.d-14
s0(1) = a(1)
s0(2) = a(2)
s2(1) = 1.d0
s2(2) = 0.d0
kn = nn
! Compute B ^ N using the binary rule for exponentiation.
do j = 1, mn
kk = kn / 2
if (kn .ne. 2 * kk) then
call ddmul (s2, s0, s1)
s2(1) = s1(1)
s2(2) = s1(2)
endif
kn = kk
if (j .lt. mn) then
call ddmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
endif
enddo
! Compute reciprocal if N is negative.
110 continue
if (n .lt. 0) then
s1(1) = 1.d0
s1(2) = 0.d0
call dddiv (s1, s2, s0)
s2(1) = s0(1)
s2(2) = s0(2)
endif
120 continue
b(1) = s2(1)
b(2) = s2(2)
return
end subroutine
subroutine ddnrt (a, n, b)
implicit none
integer i, k, n
real8 t1, t2, tn
real8 a(2), b(2), f1(2), s0(2), s1(2)
if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 140
elseif (a(1) .lt. 0.d0) then
write (6, 1)
1 format (‘*** DDNRT: Argument is negative.’)
call ddabrt
return
endif
if (n .le. 0) then
write (6, 2) n
2 format (‘*** DDNRT: Improper value of N’,i10)
call ddabrt
return
endif
! Handle cases N = 1 and 2.
if (n .eq. 1) then
b(1) = a(1)
b(2) = a(1)
goto 140
elseif (n .eq. 2) then
call ddsqrt (a, b)
goto 140
endif
f1(1) = 1.d0
f1(2) = 0.d0
! Compute the initial approximation of A ^ (-1/N).
tn = n
t1 = a(1)
t2 = exp (- log (t1) / tn)
b(1) = t2
b(2) = 0.d0
! Perform the Newton-Raphson iteration described above.
do k = 1, 3
call ddnpwr (b, n, s0)
call ddmul (a, s0, s1)
call ddsub (f1, s1, s0)
call ddmul (b, s0, s1)
call dddivd (s1, tn, s0)
call ddadd (b, s0, s1)
b(1) = s1(1)
b(2) = s1(2)
enddo
! Take the reciprocal to give final result.
call dddiv (f1, b, s1)
b(1) = s1(1)
b(2) = s1(2)
140 continue
return
end subroutine
subroutine ddout (iu, a)
! This routine writes the DD number A on logical unit iu using a standard
! E format, with lines 40 characters long.
implicit none
integer i, iu, ln
parameter (ln = 40)
character40 cs
real8 a(2)
call ddoutc (a, cs)
write (iu, ‘(a)’) cs
return
end subroutine
subroutine ddoutc (a, b)
implicit none
integer i, ii, ix, ln, nx
parameter (ln = 40)
integer ib(ln)
real8 t1
character40 b
character10 ca, digits
parameter (digits = ‘0123456789’)
real8 a(2), f(2), s0(2), s1(2)
f(1) = 10.d0
f(2) = 0.d0
do i = 1, ln
ib(i) = 0
enddo
! Determine exact power of ten for exponent.
if (a(1) .ne. 0.d0) then
t1 = log10 (abs (a(1)))
if (t1 .ge. 0.d0) then
nx = t1
else
nx = t1 - 1.d0
endif
call ddnpwr (f, nx, s0)
call dddiv (a, s0, s1)
if (s1(1) .lt. 0.d0) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif
! If we didn’t quite get it exactly right, multiply or divide by 10 to fix.
i = 0
100 continue
i = i + 1
if (s1(1) .lt. 1.d0) then
nx = nx - 1
call ddmuld (s1, 10.d0, s0)
s1(1) = s0(1)
s1(2) = s0(2)
if (i <= 3) goto 100
elseif (s1(1) .ge. 10.d0) then
nx = nx + 1
call dddivd (s1, 10.d0, s0)
s1(1) = s0(1)
s1(2) = s0(2)
goto 100
endif
else
nx = 0
s1(1) = 0.d0
s1(2) = 0.d0
endif
! Compute digits.
do i = 1, ln - 8
ii = s1(1)
ib(i) = ii
f(1) = ii
call ddsub (s1, f, s0)
call ddmuld (s0, 10.d0, s1)
enddo
! Fix negative digits.
do i = ln - 8, 2, -1
if (ib(i) .lt. 0) then
ib(i) = ib(i) + 10
ib(i-1) = ib(i-1) - 1
endif
enddo
if (ib(1) .lt. 0) then
write (6, 1)
1 format (‘ddoutc: negative leading digit’)
call ddabrt
endif
! Round.
if (ib(ln-8) .ge. 5) then
ib(ln-9) = ib(ln-9) + 1
do i = ln - 9, 2, -1
if (ib(i) .eq. 10) then
ib(i) = 0
ib(i-1) = ib(i-1) + 1
endif
enddo
if (ib(1) .eq. 10) then
ib(1) = 1
nx = nx + 1
endif
endif
! Insert digit characters in ib.
b(1:1) = ’ ’
b(2:2) = ’ ’
if (a(1) .ge. 0.d0) then
b(3:3) = ’ ’
else
b(3:3) = ‘-’
endif
ii = ib(1)
b(4:4) = digits(ii+1:ii+1)
b(5:5) = ‘.’
b(ln:ln) = ’ ’
do i = 2, ln - 9
ii = ib(i)
b(i+4:i+4) = digits(ii+1:ii+1)
enddo
! Insert exponent.
190 continue
! write (ca, ‘(i4)’) nx
ca = dddigout (dble (nx), 4)
b(ln-4:ln-4) = ‘E’
ii = 0
do i = 1, 4
if (ca(i:i) /= ’ ') then
ii = ii + 1
b(ln-4+ii:ln-4+ii) = ca(i:i)
endif
enddo
do i = ii + 1, 4
b(ln-4+i:ln-4+i) = ’ ’
enddo
return
end subroutine
subroutine ddpic (pi)
! This returns Pi to quad precision.
implicit none
real*8 pi(2), pic(2)
save pic
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.
data pic/ z’400921FB54442D18’, z’3CA1A62633145C07’/
! data pic/ 3.1415926535897931D+00, 1.2246467991473532D-16/
pi(1) = pic(1)
pi(2) = pic(2)
return
end subroutine
subroutine ddpoly (n, a, x0, x)
implicit none
integer i, it, n
real*8 a(2,0:n), ad(2,0:n), t1(2), t2(2), t3(2), t4(2), t5(2), &
x(2), x0(2), dt1, eps
parameter (eps = 1.d-30)
do i = 0, n - 1
dt1 = i + 1
call ddmuld (a(1,i+1), dt1, ad(1,i))
enddo
ad(1,n) = 0.d0
ad(2,n) = 0.d0
x(1) = x0(1)
x(2) = x0(2)
do it = 1, 20
t1(1) = 0.d0
t1(2) = 0.d0
t2(1) = 0.d0
t2(2) = 0.d0
t3(1) = 1.d0
t3(2) = 0.d0
do i = 0, n
call ddmul (a(1,i), t3, t4)
call ddadd (t1, t4, t5)
t1(1) = t5(1)
t1(2) = t5(2)
call ddmul (ad(1,i), t3, t4)
call ddadd (t2, t4, t5)
t2(1) = t5(1)
t2(2) = t5(2)
call ddmul (t3, x, t4)
t3(1) = t4(1)
t3(2) = t4(2)
enddo
call dddiv (t1, t2, t3)
call ddsub (x, t3, t4)
x(1) = t4(1)
x(2) = t4(2)
if (abs (t3(1)) .le. eps) goto 110
enddo
write (6, 1)
1 format (‘ddroot: failed to converge.’)
call ddabrt
110 continue
return
end subroutine
subroutine ddrand (a)
! This returns a pseudo-random DD number A between 0 and 1.
implicit none
real8 f7, r28, r30, r53, r58, s0, s1, s2, sd, t1, t2, t30
parameter (f7 = 78125.d0, s0 = 314159265.d0, r30 = 0.5d0 ** 30, &
r53 = 0.5d0 ** 53, r58 = 0.5d0 ** 58, t30 = 2.d0 ** 30)
real8 a(2)
save sd
data sd /s0/
t1 = f7 * sd
t2 = aint (r30 * t1)
s1 = t1 - t30 * t2
t1 = f7 * s1
t2 = aint (r30 * t1)
s2 = t1 - t30 * t2
a(1) = r30 * s1 + r58 * s2
t1 = f7 * s2
t2 = aint (r30 * t1)
s1 = t1 - t30 * t2
t1 = f7 * s1
t2 = aint (r30 * t1)
s2 = t1 - t30 * t2
a(2) = r53 * a(1) * (r30 * s1 + r58 * s2)
sd = s2
return
end subroutine
subroutine ddsqrt (a, b)
implicit none
real8 t1, t2, t3
real8 a(2), b(2), f(2), s0(2), s1(2)
if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 100
endif
t1 = 1.d0 / sqrt (a(1))
t2 = a(1) * t1
call ddmuldd (t2, t2, s0)
call ddsub (a, s0, s1)
t3 = 0.5d0 * s1(1) * t1
s0(1) = t2
s0(2) = 0.d0
s1(1) = t3
s1(2) = 0.d0
call ddadd (s0, s1, b)
100 continue
return
end subroutine
subroutine ddsub (dda, ddb, ddc)
! This subroutine computes ddc = dda - ddb.
implicit none
real8 dda(2), ddb(2), ddc(2)
real8 e, t1, t2
! Compute dda + ddb using Knuth’s trick.
t1 = dda(1) - ddb(1)
e = t1 - dda(1)
t2 = ((-ddb(1) - e) + (dda(1) - (t1 - e))) + dda(2) - ddb(2)
! The result is t1 + t2, after normalization.
ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine
real8 function dddigin (ca, n)
implicit none
real8 d1
character*(), ca
character16 digits
integer i, k, n
parameter (digits = ‘0123456789’)
d1 = 0.d0
do i = 1, n
k = index (digits, ca(i:i)) - 1
if (k < 0) then
write (6, *) ‘dddigin: non-digit in character string’
elseif (k <= 9) then
d1 = 10.d0 * d1 + k
endif
enddo
dddigin = d1
end function
character16 function dddigout (a, n)
implicit none
real8 a, d1, d2
character*16 ca, digits
parameter (digits = ‘0123456789’)
integer i, is, k, n
ca = ’ ’
is = sign (1.d0, a)
d1 = abs (a)
do i = n, 1, -1
d2 = aint (d1 / 10.d0)
k = 1.d0 + (d1 - 10.d0 * d2)
d1 = d2
ca(i:i) = digits(k:k)
if (d1 == 0.d0) goto 100
enddo
i = 0
100 continue
if (is < 0 .and. i > 1) then
ca(i-1:i-1) = ‘-’
elseif (i == 0 .or. is < 0 .and. i == 1) then
ca = ‘****************’
endif
dddigout = ca
return
end function
end module
module dddefmod
use ddfunmod
private kdb
parameter (kdb = kind (0.d0))
type dd_real
sequence
real8 ddr(2)
end type
type dd_complex
sequence
real8 ddc(4)
end type
type (dd_real), public:: ddeps
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)
contains
subroutine ddxzc (a, b)
! This converts the DC variable A to the DDC variable B.
! This routine is not intended to be called directly by the user.
complex (kdb) a
real*8 b(4)
b(1) = a
b(2) = 0.d0
b(3) = aimag (a)
b(4) = 0.d0
return
end subroutine
subroutine ddmzc (a, b)
real*8 a(2), b(4)
b(1) = a(1)
b(2) = a(2)
b(3) = 0.d0
b(4) = 0.d0
return
end subroutine
end module
module ddrealmod
use ddfunmod
use dddefmod
private kdb
parameter (kdb = kind (0.d0))
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)
private &
dd_eqqq, dd_eqqz, dd_eqrq, dd_eqqr, dd_eqdq, dd_eqqd, dd_eqxq, dd_eqqx, &
dd_addqq, dd_addqz, dd_adddq, dd_addqd, dd_addxq, dd_addqx, &
dd_subqq, dd_subqz, dd_subdq, dd_subqd, dd_subxq, dd_subqx, dd_negq, &
dd_mulqq, dd_mulqz, dd_muldq, dd_mulqd, dd_mulxq, dd_mulqx, &
dd_divqq, dd_divqz, dd_divdq, dd_divqd, dd_divxq, dd_divqx, &
dd_expqq, dd_expqi, dd_expdq, dd_expqd, &
dd_eqtqq, dd_eqtqz, dd_eqtdq, dd_eqtqd, dd_eqtxq, dd_eqtqx, &
dd_netqq, dd_netqz, dd_netdq, dd_netqd, dd_netxq, dd_netqx, &
dd_letqq, dd_letdq, dd_letqd, dd_getqq, dd_getdq, dd_getqd, &
dd_lttqq, dd_lttdq, dd_lttqd, dd_gttqq, dd_gttdq, dd_gttqd
! DDR operator extension interface blocks.
interface assignment (=)
module procedure dd_eqqq
module procedure dd_eqqz
module procedure dd_eqdq
module procedure dd_eqqd
module procedure dd_eqxq
module procedure dd_eqqx
module procedure dd_eqqa
end interface
interface operator (+)
module procedure dd_addqq
module procedure dd_addqz
module procedure dd_adddq
module procedure dd_addqd
module procedure dd_addxq
module procedure dd_addqx
end interface
interface operator (-)
module procedure dd_subqq
module procedure dd_subqz
module procedure dd_subdq
module procedure dd_subqd
module procedure dd_subxq
module procedure dd_subqx
module procedure dd_negq
end interface
interface operator (*)
module procedure dd_mulqq
module procedure dd_mulqz
module procedure dd_muldq
module procedure dd_mulqd
module procedure dd_mulxq
module procedure dd_mulqx
end interface
interface operator (/)
module procedure dd_divqq
module procedure dd_divqz
module procedure dd_divdq
module procedure dd_divqd
module procedure dd_divxq
module procedure dd_divqx
end interface
interface operator (**)
module procedure dd_expqq
module procedure dd_expqi
module procedure dd_expdq
module procedure dd_expqd
end interface
interface operator (.eq.)
module procedure dd_eqtqq
module procedure dd_eqtqz
module procedure dd_eqtdq
module procedure dd_eqtqd
module procedure dd_eqtxq
module procedure dd_eqtqx
end interface
interface operator (.ne.)
module procedure dd_netqq
module procedure dd_netqz
module procedure dd_netdq
module procedure dd_netqd
module procedure dd_netxq
module procedure dd_netqx
end interface
interface operator (.le.)
module procedure dd_letqq
module procedure dd_letdq
module procedure dd_letqd
end interface
interface operator (.ge.)
module procedure dd_getqq
module procedure dd_getdq
module procedure dd_getqd
end interface
interface operator (.lt.)
module procedure dd_lttqq
module procedure dd_lttdq
module procedure dd_lttqd
end interface
interface operator (.gt.)
module procedure dd_gttqq
module procedure dd_gttdq
module procedure dd_gttqd
end interface
contains
! DDR assignment routines.
subroutine dd_eqqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: qb
qa%ddr(1) = qb%ddr(1)
qa%ddr(2) = qb%ddr(2)
return
end subroutine
subroutine dd_eqqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: zb
qa%ddr(1) = zb%ddc(1)
qa%ddr(2) = zb%ddc(2)
return
end subroutine
subroutine dd_eqdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: da
intent (in):: qb
da = qb%ddr(1)
return
end subroutine
subroutine dd_eqqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: db
qa%ddr(1) = db
qa%ddr(2) = 0.d0
return
end subroutine
subroutine dd_eqxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: xa
intent (in):: qb
xa = qb%ddr(1)
return
end subroutine
subroutine dd_eqqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: xb
qa%ddr(1) = xb
qa%ddr(2) = 0.d0
return
end subroutine
subroutine dd_eqqa (qa, ab)
implicit real8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
character(), intent (in):: ab
intent (out):: qa
character80 t
t = ab
call ddinpc (t, qa%ddr)
return
end subroutine
! DDR add routines.
function dd_addqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_addqq
intent (in):: qa, qb
call ddadd (qa%ddr, qb%ddr, dd_addqq%ddr)
return
end function
function dd_addqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcadd (ddt1, zb%ddc, dd_addqz%ddc)
return
end function
function dd_adddq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_adddq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddadd (ddt1, qb%ddr, dd_adddq%ddr)
return
end function
function dd_addqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_addqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddadd (qa%ddr, ddt1, dd_addqd%ddr)
return
end function
function dd_addxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcadd (ddt1, ddt2, dd_addxq%ddc)
return
end function
function dd_addqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcadd (ddt1, ddt2, dd_addqx%ddc)
return
end function
! DDR subtract routines.
function dd_subqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subqq
intent (in):: qa, qb
call ddsub (qa%ddr, qb%ddr, dd_subqq%ddr)
return
end function
function dd_subqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subqz%ddc)
return
end function
function dd_subdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddsub (ddt1, qb%ddr, dd_subdq%ddr)
return
end function
function dd_subqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddsub (qa%ddr, ddt1, dd_subqd%ddr)
return
end function
function dd_subxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcsub (ddt1, ddt2, dd_subxq%ddc)
return
end function
function dd_subqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcsub (ddt1, ddt2, dd_subqx%ddc)
return
end function
! DDR negation routine.
function dd_negq (qa)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_negq
intent (in):: qa
call ddeq (qa%ddr, dd_negq%ddr)
dd_negq%ddr(1) = - qa%ddr(1)
dd_negq%ddr(2) = - qa%ddr(2)
return
end function
! DDR multiply routines.
function dd_mulqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_mulqq
intent (in):: qa, qb
call ddmul (qa%ddr, qb%ddr, dd_mulqq%ddr)
return
end function
function dd_mulqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcmul (ddt1, zb%ddc, dd_mulqz%ddc)
return
end function
function dd_muldq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_muldq
intent (in):: da, qb
call ddmuld (qb%ddr, da, dd_muldq%ddr)
return
end function
function dd_mulqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_mulqd
intent (in):: qa, db
call ddmuld (qa%ddr, db, dd_mulqd%ddr)
return
end function
function dd_mulxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcmul (ddt1, ddt2, dd_mulxq%ddc)
return
end function
function dd_mulqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcmul (ddt1, ddt2, dd_mulqx%ddc)
return
end function
! DDR divide routines.
function dd_divqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divqq
intent (in):: qa, qb
call dddiv (qa%ddr, qb%ddr, dd_divqq%ddr)
return
end function
function dd_divqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcdiv (ddt1, zb%ddc, dd_divqz%ddc)
return
end function
function dd_divdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call dddiv (ddt1, qb%ddr, dd_divdq%ddr)
return
end function
function dd_divqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divqd
intent (in):: qa, db
call dddivd (qa%ddr, db, dd_divqd%ddr)
return
end function
function dd_divxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcdiv (ddt1, ddt2, dd_divxq%ddc)
return
end function
function dd_divqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcdiv (ddt1, ddt2, dd_divqx%ddc)
return
end function
! DDR exponentiation routines.
function dd_expqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqq
intent (in):: qa, qb
call ddlog (qa%ddr, ddt1)
call ddmul (ddt1, qb%ddr, ddt2)
call ddexp (ddt2, dd_expqq%ddr)
return
end function
function dd_expqi (qa, ib)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqi
intent (in):: qa, ib
call ddnpwr (qa%ddr, ib, dd_expqi%ddr)
return
end function
function dd_expdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddlog (ddt1, ddt2)
call ddmul (ddt2, qb%ddr, ddt3)
call ddexp (ddt3, dd_expdq%ddr)
return
end function
function dd_expqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqd
intent (in):: qa, db
call ddlog (qa%ddr, ddt1)
call ddmuld (ddt1, db, ddt2)
call ddexp (ddt2, dd_expqd%ddr)
return
end function
! DDR .EQ. routines.
function dd_eqtqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .eq. 0) then
dd_eqtqq = .true.
else
dd_eqtqq = .false.
endif
return
end function
function dd_eqtqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcpr (ddt1, zb%ddc, ic1)
call ddcpr (ddt1(3), zb%ddc(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtqz = .true.
else
dd_eqtqz = .false.
endif
return
end function
function dd_eqtdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .eq. 0) then
dd_eqtdq = .true.
else
dd_eqtdq = .false.
endif
return
end function
function dd_eqtqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .eq. 0) then
dd_eqtqd = .true.
else
dd_eqtqd = .false.
endif
return
end function
function dd_eqtxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtxq = .true.
else
dd_eqtxq = .false.
endif
return
end function
function dd_eqtqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtqx = .true.
else
dd_eqtqx = .false.
endif
return
end function
! DDR .NE. routines.
function dd_netqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .ne. 0) then
dd_netqq = .true.
else
dd_netqq = .false.
endif
return
end function
function dd_netqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcpr (ddt1, zb%ddc, ic1)
call ddcpr (ddt1(3), zb%ddc(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netqz = .true.
else
dd_netqz = .false.
endif
return
end function
function dd_netdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .ne. 0) then
dd_netdq = .true.
else
dd_netdq = .false.
endif
return
end function
function dd_netqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .ne. 0) then
dd_netqd = .true.
else
dd_netqd = .false.
endif
return
end function
function dd_netxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netxq = .true.
else
dd_netxq = .false.
endif
return
end function
function dd_netqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netqx = .true.
else
dd_netqx = .false.
endif
return
end function
! DDR .LE. routines.
function dd_letqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .le. 0) then
dd_letqq = .true.
else
dd_letqq = .false.
endif
return
end function
function dd_letdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .le. 0) then
dd_letdq = .true.
else
dd_letdq = .false.
endif
return
end function
function dd_letqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .le. 0) then
dd_letqd = .true.
else
dd_letqd = .false.
endif
return
end function
! DDR .GE. routines.
function dd_getqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .ge. 0) then
dd_getqq = .true.
else
dd_getqq = .false.
endif
return
end function
function dd_getdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .ge. 0) then
dd_getdq = .true.
else
dd_getdq = .false.
endif
return
end function
function dd_getqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .ge. 0) then
dd_getqd = .true.
else
dd_getqd = .false.
endif
return
end function
! DDR .LT. routines.
function dd_lttqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .lt. 0) then
dd_lttqq = .true.
else
dd_lttqq = .false.
endif
return
end function
function dd_lttdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .lt. 0) then
dd_lttdq = .true.
else
dd_lttdq = .false.
endif
return
end function
function dd_lttqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .lt. 0) then
dd_lttqd = .true.
else
dd_lttqd = .false.
endif
return
end function
! DDR .GT. routines.
function dd_gttqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .gt. 0) then
dd_gttqq = .true.
else
dd_gttqq = .false.
endif
return
end function
function dd_gttdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .gt. 0) then
dd_gttdq = .true.
else
dd_gttdq = .false.
endif
return
end function
function dd_gttqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .gt. 0) then
dd_gttqd = .true.
else
dd_gttqd = .false.
endif
return
end function
end module
module ddcmpmod
use ddfunmod
use dddefmod
private kdb
parameter (kdb = kind (0.d0))
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)
private &
dd_eqzq, dd_eqzz, dd_eqdz, dd_eqzd, dd_eqxz, dd_eqzx, &
dd_addzq, dd_addzz, dd_adddz, dd_addzd, dd_addxz, dd_addzx, &
dd_subzq, dd_subzz, dd_subdz, dd_subzd, dd_subxz, dd_subzx, dd_negz, &
dd_mulzq, dd_mulzz, dd_muldz, dd_mulzd, dd_mulxz, dd_mulzx, &
dd_divzq, dd_divzz, dd_divdz, dd_divzd, dd_divxz, dd_divzx, &
dd_expzi, &
dd_eqtzq, dd_eqtzz, dd_eqtdz, dd_eqtzd, dd_eqtxz, dd_eqtzx, &
dd_netzq, dd_netzz, dd_netdz, dd_netzd, dd_netxz, dd_netzx
! DDR operator extension interface blocks.
interface assignment (=)
module procedure dd_eqzq
module procedure dd_eqzz
module procedure dd_eqdz
module procedure dd_eqzd
module procedure dd_eqxz
module procedure dd_eqzx
end interface
interface operator (+)
module procedure dd_addzq
module procedure dd_addzz
module procedure dd_adddz
module procedure dd_addzd
module procedure dd_addxz
module procedure dd_addzx
end interface
interface operator (-)
module procedure dd_subzq
module procedure dd_subzz
module procedure dd_subdz
module procedure dd_subzd
module procedure dd_subxz
module procedure dd_subzx
module procedure dd_negz
end interface
interface operator (*)
module procedure dd_mulzq
module procedure dd_mulzz
module procedure dd_muldz
module procedure dd_mulzd
module procedure dd_mulxz
module procedure dd_mulzx
end interface
interface operator (/)
module procedure dd_divzq
module procedure dd_divzz
module procedure dd_divdz
module procedure dd_divzd
module procedure dd_divxz
module procedure dd_divzx
end interface
interface operator (**)
module procedure dd_expzi
end interface
interface operator (.eq.)
module procedure dd_eqtzq
module procedure dd_eqtzz
module procedure dd_eqtdz
module procedure dd_eqtzd
module procedure dd_eqtxz
module procedure dd_eqtzx
end interface
interface operator (.ne.)
module procedure dd_netzq
module procedure dd_netzz
module procedure dd_netdz
module procedure dd_netzd
module procedure dd_netxz
module procedure dd_netzx
end interface
contains
! DDC assignment routines.
subroutine dd_eqzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: qb
call ddmzc (qb%ddr, za%ddc)
return
end subroutine
subroutine dd_eqzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: zb
call ddceq (zb%ddc, za%ddc)
return
end subroutine
subroutine dd_eqdz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: da
intent (in):: zb
da = zb%ddc(1)
return
end subroutine
subroutine dd_eqzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: db
xb = db
call ddxzc (xb, za%ddc)
return
end subroutine
subroutine dd_eqxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: xa
intent (in):: zb
db = zb%ddc(1)
dc = zb%ddc(3)
xa = cmplx (db, dc, kdb)
return
end subroutine
subroutine dd_eqzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: xb
call ddxzc (xb, za%ddc)
return
end subroutine
! DDC add routines.
function dd_addzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcadd (za%ddc, ddt1, dd_addzq%ddc)
return
end function
function dd_addzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzz
intent (in):: za, zb
call ddcadd (za%ddc, zb%ddc, dd_addzz%ddc)
return
end function
function dd_adddz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_adddz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcadd (ddt1, zb%ddc, dd_adddz%ddc)
return
end function
function dd_addzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcadd (za%ddc, ddt1, dd_addzd%ddc)
return
end function
function dd_addxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addxz
intent (in):: xa, zb
call ddxzc (xa, ddt1)
call ddcadd (ddt1, zb%ddc, dd_addxz%ddc)
return
end function
function dd_addzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzx
intent (in):: za, xb
call ddxzc (xb, ddt2)
call ddcadd (za%ddc, ddt2, dd_addzx%ddc)
return
end function
! DDC subtract routines.
function dd_subzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcsub (za%ddc, ddt1, dd_subzq%ddc)
return
end function
function dd_subzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzz
intent (in):: za, zb
call ddcsub (za%ddc, zb%ddc, dd_subzz%ddc)
return
end function
function dd_subdz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subdz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subdz%ddc)
return
end function
function dd_subzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcsub (za%ddc, ddt1, dd_subzd%ddc)
return
end function
function dd_subxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subxz
intent (in):: xa, zb
call ddxzc (xa, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subxz%ddc)
return
end function
function dd_subzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzx
intent (in):: za, xb
call ddxzc (xb, ddt2)
call ddcsub (za%ddc, ddt2, dd_subzx%ddc)
return
end function
! DDC negation routine.
function dd_negz (za)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_negz
intent (in):: za
call ddceq (za%ddc, dd_negz%ddc)
dd_negz%ddc(1) = - za%ddc(1)
dd_negz%ddc(2) = - za%ddc(2)
dd_negz%ddc(3) = - za%ddc(3)
dd_negz%ddc(4) = - za%ddc(4)
return
end function
! DDC multiply routines.
function dd_mulzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcmul (za%ddc, ddt1, dd_mulzq%ddc)
return
end function
function dd_mulzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzz
intent (in):: za, zb
call ddcmul (za%ddc, zb%ddc, dd_mulzz%ddc)
return
end function
function dd_muldz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_muldz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcmul (ddt1, zb%ddc, dd_muldz%ddc)
return
end function
function dd_mulzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcmul (za%ddc,