comment augment description deck for the multiple-precision arithmetic
package of r. p. brent.
three types of variable are defined here -
multiple (standard multiple-precision numbers),
multipak (packed multiple-precision numbers), and
initialize (used only as a device to persuade
augment to initialize the mp package).
working space should be allocated and the mp package
initialized by the declaration
initialize mp
in the main program.
this description deck assumes that each mp number requires no
more than 27 words (14 in packed format). this is sufficient
for about 43 decimal place accuracy on a 16-bit machine, and
higher accuracy on a machine with a longer wordlength.
see comments in routine mpinit for the method of changing the
precision, also regarding common declarations.
*describe multipak
conversion ctp (cim, integer, $, upward), ctp (cam, hollerith)
*describe multiple
operator + (, null unary, prv, $), - (neg, unary),
+ (add, binary3, prv, $, $, $, comm), * (mul),
- (sub,,,,,, noncomm), / (div), ** (pwr2),
+ (addi,,,, integer), * (muli), / (divi), ** (pwr),
* (imul,,, integer, $, $, noncomm),
.eq. (eq, binary2, prv, $, logical, comm), .ne. (ne),
.ge. (ge,,,,, noncomm), .gt. (gt), .le. (le), .lt. (lt),
+ (, null unary, prv, multipak),
- (neg, unary, prv, multipak),
+ (kadd, binary3, prv, multipak, multipak, $, comm),
* (kmul), - (ksub,,,,,, noncomm), / (kdiv), ** (pwr2),
* (kmli,,,, integer), / (kdvi), ** (pwr),
* (kiml,,, integer, multipak, $, noncomm),
.eq. (eq, binary2, prv, multipak, logical, comm), .ne. (ne),
.ge. (ge,,,,, noncomm), .gt. (gt), .le. (le), .lt. (lt)
test mpsiga (siga, integer)
field sgn (siga, sigb, ($), integer), expon (expa, expb),
base (basa, basb), numdig (diga, digb), maxexp (mexa, mexb),
digit (dga, dgb, ($, integer)),
param (para, parb, (hollerith)), fio (fin, fout, (integer),
$), unfio (unfr, unfw),
expon (expa, expb, (multipak), integer), numdig (diga, digb),
sgn (siga, sigb)
function abs (abs, ($), $), dabs (abs), asin (asin), dasin (asin),
arsin (asin), darsin (asin), arcsin (asin), atan (atan),
datan (atan), artan (atan), dartan (atan), arctan (atan),
cmf (cmf), ceil (ceil), floor (flor), cmim (cmim),
cos (cos), dcos (cos), cosh (cosh), dcosh (cosh),
daw (daw), ddaw (daw), ei (ei), dei (ei), erf (erf),
derf (erf), erfc (erfc), derfc (erfc), exp (exp),
dexp (exp), exp1 (exp1), frac (cmf), gam (gam),
dgam (gam), gamma (gam), dgamma (gam), aint (cmim),
dint (cmim), li (li), dli (li), ln (ln), log (ln),
alog (ln), dlog (ln), log10 (lg10), alog10 (lg10),
dlog10 (lg10), lngm (lngm), alngm (lngm), dlngm (lngm),
lngs (lngs), lns (lns), rec (rec), sin (sin), dsin (sin),
sinh (sinh), dsinh (sinh), sqrt (sqrt), dsqrt (sqrt),
tan (tan), dtan (tan), tanh (tanh), dtanh (tanh),
sngl (str), dble (str),
art1 (art1, (integer)), ln (lni), lni (lni), log (lni),
alog (lni), dlog (lni), zeta (zeta),
cam (cam), cam (cam, (hollerith)),
max (max, ($, $)), amax1 (max), dmax1 (max), min (min),
amin1 (min), dmin1 (min), dim (dim), ddim (dim),
gcd (gcda), atan2 (atn2), datan2 (atn2), artan2 (atn2),
mod (mod), amod (mod), dmod (mod),
sign (sign), dsign (sign),
besj (besj, ($, integer)), root (root),
mpinf (inf (subroutine), ($, integer, integer, hollerith),
logical), mpoutf (outf (subroutine)),
mpinf (inf (subroutine), ($, integer, integer, integer)),
mpoutf (outf (subroutine)),
integr (intg, ($), logical),
comp (comp, ($, $), integer), cmpa (cmpa),
comp (cmpi, ($, integer)), comp (cmpr, ($, real)),
comp (cmpd, ($, double precision)),
comp (cmpq, ($, integer, integer)),
addq (addq, ($, integer, integer), $), mulq (mulq),
qpwr (qpwr, (integer, integer, integer, integer)),
cqm (cqm, (integer, integer)), ctm (cqm),
gam (gamq), dgam (gamq), gamq (gamq),
bern (bern, (integer, integer), multipak),
int (cmi (subroutine), ($), integer),
idint (cmi (subroutine)), ifix (cmi (subroutine)),
abs (abs, (multipak), multipak), atan (atan,, $), cos (cos),
cosh (cosh), exp (exp), ln (ln), log (ln), alog (ln),
log10 (lg10), alog10 (lg10), sin (sin), sinh (sinh),
sqrt (sqrt), tan (tan), tanh (tanh),
exp1 (exp1), lns (lns),
max (max, (multipak, multipak)), min (min), dim (dim),
atan2 (atn2), mod (mod), sign (sign),
root (root, (multipak, integer))
conversion ctm (cdm, double precision, $, upward), ctm (cim, integer),
ctm (crm, real), ctm (unpk, multipak), ctm (cam, hollerith),
ctd (cmd (subroutine), $, double precision, downward),
cti (cmi (subroutine),, integer), ctr (cmr (subroutine),,
real), ctp (pack,, multipak)
*describe initialize
comment end of augment description deck for mp package
(machine-dependent statements to execute augment
and add description deck as data for augment)
*begin
*disable print, output
*enable source
(machine-dependent statement(s) to invoke fortran compiler)
c
c program to verify an identity of jacobi using the mp
c package and augment.
c
c the program reads a number x in free-field format acceptable to
c mpin. if x is non-positive it halts. otherwise it computes
c and prints fn(x), fn(1/x) and (fn(x)-fn(1/x))/fn(x),
c where fn(x) is the sum from n = -infinity to +infinity of
c sqrt(x)*exp(-pi*(n*x)**2).
c the identity verified is
c fn(x) = fn(1/x)
c
c declare variables and initialize mp package. on some systems blank
c common must be declared here - see comments in routine mpinit.
c
multiple x, f1, f2, fn
initialize mp
c
c set input record length to 72 (to avoid reading
c sequence numbers), default value is 80.
c
param (6hinrecl) = 72
c
c read mp x from unit 5 in free format, stop if x not positive.
c
10 x = fio (5)
if (x.le.0) stop
c
c write heading, x, fn(x), and fn(1/x) to 40s in standard format
c
write (6, 20)
20 format (//41h x, fn(x), fn(1/x), (fn(x)-fn(1/x))/fn(x)/)
fio (40) = x
f1 = fn(x)
fio (40) = f1
f2 = fn(1/x)
fio (40) = f2
c
c write (f1-f2)/f1 to 6 significant figures.
c note that an mp expression can be written.
c
fio (6) = (f1 - f2)/f1
go to 10
end
c
(machine-dependent statement(s) to invoke fortran compiler)
c
function fn(x)
c
c returns fn(x) = the sum from n = -infinity to +infinity of
c sqrt(x)*exp(-pi*(n*x)**2), assuming x positive.
c uses the obvious method, so slow if x small.
c note that x and fn are both type multiple.
c to illustrate mixed-mode arithmetic, some local
c variables are declared to be packed.
c
multiple fn, x
multipak tm, fac, pr
if (x.le.0) call mperr
fn = 0
c
c augment can deal with the following expression as it knows that x
c is type multiple, so calls mpcam to convert 2hpi to multiple.
c
tm = exp(-2hpi*x*x)
pr = tm
fac = tm**2
c
c loop to sum infinite series
c warning - number of iterations is proportional to 1/x
c
10 fn = fn + tm
pr = pr*fac
tm = tm*pr
c
c test for convergence by comparing exponents of fn and tm.
c we could also have saved the old value of fn and seen if
c statement 10 changed it.
c
if (expon(fn)-expon(tm).lt.numdig(x)) go to 10
fn = sqrt(x)*(2*fn+1)
return
end
*end
(machine-dependent statements to compile output from augment,
link to precompiled mp library, and execute with the
following data)
.5
.3
c $$ ****** jacobi (augment output) part 1 ******
c
c program to verify an identity of jacobi using the mp
c package and augment.
c
c the program reads a number x in free-field format acceptable to
c mpin. if x is non-positive it halts. otherwise it computes
c and prints fn(x), fn(1/x) and (fn(x)-fn(1/x))/fn(x),
c where fn(x) is the sum from n = -infinity to +infinity of
c sqrt(x)*exp(-pi*(n*x)**2).
c the identity verified is
c fn(x) = fn(1/x)
c
c declare variables and initialize mp package. on some systems blank
c common must be declared here - see comments in routine mpinit.
c
c ===== processed by augment, version 4n =====
c ----- temporary storage locations -----
c multiple
integer mptmp(27,1)
c ----- local variables -----
c multiple
integer f1(27), f2(27), x(27)
c initialize
integer mp(1)
c ----- supporting package functions -----
logical mple
c ===== translated program =====
c
c set input record length to 72 (to avoid reading
c sequence numbers), default value is 80.
c
c ----- begin initialization -----
call mpinit (mp)
c ----- end initialization -----
c
c param (6hinrecl) = 72
call mpparb (72, 6hinrecl)
c
c read mp x from unit 5 in free format, stop if x not positive.
c
c
c 10 x = fio (5)
10 call mpfin (5,x)
c
c if (x.le.0) stop
c ===== mixed mode operands accepted =====
call mpcim (0,mptmp(1,1))
if (mple (x,mptmp(1,1))) stop
c
c write heading, x, fn(x), and fn(1/x) to 40s in standard format
c
c
c write (6, 20)
write (6, 20)
20 format (//41h x, fn(x), fn(1/x), (fn(x)-fn(1/x))/fn(x)/)
c
c fio (40) = x
call mpfout (x, 40)
c
c f1 = fn(x)
call fn (x,f1)
c
c fio (40) = f1
call mpfout (f1, 40)
c
c f2 = fn(1/x)
c ===== mixed mode operands accepted =====
call mpcim (1,mptmp(1,1))
call mpdiv (mptmp(1,1),x,mptmp(1,1))
call fn (mptmp(1,1),f2)
c
c fio (40) = f2
call mpfout (f2, 40)
c
c write (f1-f2)/f1 to 6 significant figures.
c note that an mp expression can be written.
c
c
c fio (6) = (f1 - f2)/f1
call mpsub (f1,f2,mptmp(1,1))
call mpdiv (mptmp(1,1),f1,mptmp(1,1))
call mpfout (mptmp(1,1), 6)
c
c go to 10
go to 10
end
c
c $$ ****** jacobi (augment output) part 2 ******
c
subroutine fn (x, mprlt)
c ===== processed by augment, version 4n =====
c ----- temporary storage locations -----
c multiple
integer mptmp(27,2)
c ----- local variables -----
c multipak
integer fac(14), pr(14), tm(14)
c multiple
integer mpres(27)
c ----- global variables -----
c multiple
integer x(27), mprlt(27)
c ----- supporting package functions -----
logical mple
integer mpdiga, mpexpa
c ===== translated program =====
c
c returns fn(x) = the sum from n = -infinity to +infinity of
c sqrt(x)*exp(-pi*(n*x)**2), assuming x positive.
c uses the obvious method, so slow if x small.
c note that x and fn are both type multiple.
c to illustrate mixed-mode arithmetic, some local
c variables are declared to be packed.
c
c
c if (x.le.0) call mperr
c ===== mixed mode operands accepted =====
call mpcim (0,mptmp(1,1))
if (mple (x,mptmp(1,1))) call mperr
c
c fn = 0
call mpcim (0,mpres)
c
c augment can deal with the following expression as it knows that x
c is type multiple, so calls mpcam to convert 2hpi to multiple.
c
c
c tm = exp(-2hpi*x*x)
c ===== mixed mode operands accepted =====
call mpcam (2hpi,mptmp(1,1))
call mpmul (mptmp(1,1),x,mptmp(1,1))
call mpmul (mptmp(1,1),x,mptmp(1,1))
call mpneg (mptmp(1,1),mptmp(1,1))
call mpexp (mptmp(1,1),mptmp(1,1))
call mppack (mptmp(1,1),tm)
c
c pr = tm
call mpstr (tm,pr)
c
c fac = tm**2
call mppwr (tm,2,mptmp(1,1))
call mppack (mptmp(1,1),fac)
c
c loop to sum infinite series
c warning - number of iterations is proportional to 1/x
c
c
c 10 fn = fn + tm
c ===== mixed mode operands accepted =====
10 call mpunpk (tm,mptmp(1,1))
call mpadd (mpres,mptmp(1,1),mpres)
c
c pr = pr*fac
call mpkmul (pr,fac,mptmp(1,1))
call mppack (mptmp(1,1),pr)
c
c tm = tm*pr
call mpkmul (tm,pr,mptmp(1,1))
call mppack (mptmp(1,1),tm)
c
c test for convergence by comparing exponents of fn and tm.
c we could also have saved the old value of fn and seen if
c statement 10 changed it.
c
c
c if (expon(fn)-expon(tm).lt.numdig(x)) go to 10
if (mpexpa (mpres)-mpexpa (tm).lt.mpdiga (x)) go to 10
c
c fn = sqrt(x)*(2*fn+1)
call mpsqrt (x,mptmp(1,1))
call mpimul (2,mpres,mptmp(1,2))
call mpaddi (mptmp(1,2),1,mptmp(1,2))
call mpmul (mptmp(1,1),mptmp(1,2),mpres)
c
c return
go to 30000
c ----- return code -----
call mpstr (mpres,mprlt)
return
end
.
