[CONTACT]

[ABOUT]

[POLICY]

package of three types of variable

Found at: ftp.icm.edu.pl:70/packages/netlib/bmp/descr.deck

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

.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]