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 declare integer (14), kind safe subroutine, prefix mpk conversion ctp (cim, integer, $, upward), ctp (cam, hollerith) service copy (=mpstr) *describe multiple declare integer (27), kind safe subroutine, prefix mp 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) service copy (str) *describe initialize declare integer (1), kind safe subroutine, prefix mpi service copy (str), initial (nit) 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 1.+1 1.234567890123456789012345678901234567890123456789 0 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 ----- 30000 continue call mpstr (mpres,mprlt) return end