what does it mean
Printable View
what does it mean
Number crunching.
One of the most famous examples.
Another exampleQuote:
program piqp
c This program employs the recently discovered digit extraction scheme
c to produce hex digits of pi. This code is valid up to ic = 2^24 on
c systems with IEEE arithmetic.
real*8 pid, s1, s2, s3, s4, series
parameter (ic = 100000, nhx = 12)
character*1 chx(nhx)
c ic is the hex digit position -- output begins at position ic + 1.
s1 = series (1, ic)
s2 = series (4, ic)
s3 = series (5, ic)
s4 = series (6, ic)
pid = mod (4.d0 * s1 - 2.d0 * s2 - s3 - s4, 1.d0) + 1.d0
call ihex (pid, nhx, chx)
write (6, 1) ic, pid, chx
1 format ('Pi hex digit computation'/'position =',i12,' + 1'/
$ f20.15/12a1)
stop
end
subroutine ihex (x, nhx, chx)
c This returns, in chx, the first nhx hex digits of the fraction of x.
real*8 x, y
character*1 chx(nhx), hx(0:15)
data hx/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'A',
$ 'B', 'C', 'D', 'E', 'F'/
y = abs (x)
do i = 1, nhx
y = 16.d0 * (y - aint (y))
chx(i) = hx(int(y))
enddo
return
end
function series (m, ic)
c This routine evaluates the series sum_k 16^(ic-k)/(8*k+m)
c using the modular exponentiation technique.
real*8 ak, eps, expm, p, s, series, t
parameter (eps = 1d-17)
s = 0.d0
c Sum the series up to ic.
do k = 0, ic - 1
ak = 8 * k + m
p = ic - k
t = expm (p, ak)
s = mod (s + t / ak, 1.d0)
enddo
c Compute a few terms where k >= ic.
do k = ic, ic + 100
ak = 8 * k + m
t = 16.d0 ** (ic - k) / ak
if (t .lt. eps) goto 100
s = mod (s + t, 1.d0)
enddo
100 continue
series = s
return
end
function expm (p, ak)
c expm = 16^p mod ak. This routine uses the left-to-right binary
c exponentiation scheme. It is valid for ak <= 2^24.
parameter (ntp = 25)
real*8 ak, expm, p, p1, pt, r, tp(ntp)
save tp
data tp/ntp*0.d0/
c If this is the first call to expm, fill the power of two table tp.
if (tp(1) .eq. 0.d0) then
tp(1) = 1.d0
do i = 2, ntp
tp(i) = 2.d0 * tp(i-1)
enddo
endif
if (ak .eq. 1.d0) then
expm = 0.d0
return
endif
c Find the greatest power of two less than or equal to p.
do i = 1, ntp
if (tp(i) .gt. p) goto 100
enddo
100 pt = tp(i-1)
p1 = p
r = 1.d0
c Perform binary exponentiation algorithm modulo ak.
110 continue
if (p1 .ge. pt) then
r = mod (16.d0 * r, ak)
p1 = p1 - pt
endif
pt = 0.5d0 * pt
if (pt .ge. 1.d0) then
r = mod (r * r, ak)
goto 110
endif
expm = mod (r, ak)
return
end
http://www.activemind.com/Mysterious/Topic...e_equation.html
You could have just used a dictionary defenition... <_<
Feel free to add to the discussion, I believe that's what forums are for?
Nah I prefer to drop a little comment here and there and then walk off again.
Should I provide a definition for that?
:lol: