subroutine mp40d (n, x) c fixed-point output routine, writes mp x on unit lun with c n decimal places, assuming -10 .lt. x .lt. 100. c for rounding options see comments in subroutine mpout. common r logical err integer i2, mpparn, n, r(1), sv, x(1) c do nothing if n nonpositive if (n.le.0) return call mpsavn (sv) c allocate n+3 words of temporary space. call mpnew2 (i2, n+3) c convert to character format and print call mpout (x, r(i2), n+3, n) call mpio (r(i2), n+3, mpparn(4), \$ 37h(8x,13a1,4(1x,10a1)/(10x,5(1x,10a1))), err) if (err) call mperrm (32herror return from mpio in mp40d\$) call mpresn (sv) return end subroutine mp40f (n, x) c this routine is redundant, but is included for compatibility c with earlier versions of the mp package. integer n, x(1) call mpfout (x, n) return end subroutine mpabs (x, y) c sets y = abs(x) for mp numbers x and y. c y will be packed if x is. integer x(1), y(1) call mpstr (x, y) y(1) = iabs(y(1)) return end subroutine mpadd (x, y, z) c adds x and y, forming result in z, where x, y and z are mp numbers. c rounding defined by parameter rndrl in common /mpcom/ as follows - c rndrl = 0 - truncate towards zero if x*y .ge. 0, c away from zero if x*y .lt. 0, c in both cases using one guard digit, so result is c exact if severe cancellation occurs. c rndrl = 1, 2 or 3 - see comments in subroutine mpnzr. c sufficient guard digits are used to ensure the c correct result. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer ed, i2, j, kg, med, re, rs, s, two integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr, \$ r(1), rndrl, sptr, t, x(1), y(1), z(1) two = 2 c compute number of guard digits required kg = 1 if (rndrl.ne.0) kg = t c allocate t + kg words for accumulator. call mpnew2 (i2, t+kg) c check for x or y zero if (x(1).ne.0) go to 20 c x = 0 or negligible, so result = y 10 call mpstr(y, z) go to 120 20 if (y(1).ne.0) go to 40 c y = 0 or negligible, so result = x 30 call mpstr (x, z) go to 120 c compare signs 40 s = x(1)*y(1) if (iabs(s).le.1) go to 50 call mperrm (38hsign not 0, +1 or -1 in call to mpadd\$) go to 80 c compare exponents 50 ed = x(two) - y(two) med = iabs(ed) c can reduce kg if med small. this saves time. kg = min0 (kg, med + (s+1)/2) if (ed) 90, 60, 130 c exponents equal so compare signs, then fractions if nec. 60 if (s.gt.0) go to 100 do 70 j = 1, t if (x(j+2) - y(j+2)) 100, 70, 140 70 continue c result is zero 80 z(1) = 0 go to 120 c here exponent(y) .gt. exponent(x) 90 if ((med.gt.t).and.(rndrl.le.1)) go to 10 c check for common special case here (for sake of speed) if ((rndrl.gt.0).or.(((y(two+1).ge.(b-1)).or.(s.lt.0)).and. \$ ((y(two+1).eq.1).or.(s.gt.0)))) go to 100 c can avoid calling mpnzr here, saving time. call mpadd3 (x, y, re, z(two+1), 0) z(1) = y(1) z(two) = re go to 120 c here abs(y) .gt. abs(x) 100 rs = y(1) call mpadd3 (x, y, re, r(i2), kg) c normalize and round or truncate 110 call mpnzr (rs, re, z, r(i2), kg) c restore stack pointer and return 120 sptr = i2 return c here exponent(x) .gt. exponent(y) c code is as above with x and y interchanged. 130 if ((med.gt.t).and.(rndrl.le.1)) go to 30 c check for common special case here (for sake of speed) if ((rndrl.gt.0).or.(((x(two+1).ge.(b-1)).or.(s.lt.0)).and. \$ ((x(two+1).eq.1).or.(s.gt.0)))) go to 140 c can avoid calling mpnzr here. call mpadd3 (y, x, re, z(two+1), 0) z(1) = x(1) z(two) = re go to 120 c here abs(x) .gt. abs(y) 140 rs = x(1) call mpadd3 (y, x, re, r(i2), kg) go to 110 end subroutine mpadd3 (x, y, re, a, kg) c called by mpadd, does inner loops of addition c assumes 0 .lt. abs(x) .le. abs(y). c returns digits of x + y in a(1), ... , a(t+kg) c and exponent in re. assumes kg .ge. 0. there is a slight c cludge to ensure correct directed rounding. this may c affect a(t+1), ... , a(t+kg). common /mpcom/ b, t, dummy integer c, i, j, med, s, ted, tg, tg1, two, \$ a(1), b, dummy(21), kg, re, t, x(1), y(1) two = 2 c y(two) - x(two) .gt. t in call only for rndrl .gt. 1 med = min0 (y(two) - x(two), t) ted = t + med s = x(1)*y(1) re = y(two) tg = t + kg i = tg c = 0 c clear guard digits to right of x digits 10 if (i.le.ted) go to 20 a(i) = 0 i = i - 1 go to 10 20 if (s.lt.0) go to 130 c here do addition, exponent(y) .ge. exponent(x) if (i.le.t) go to 40 30 j = i - med a(i) = x(j+2) i = i - 1 if (i.gt.t) go to 30 40 if (i.le.med) go to 60 j = i - med c = y(i+2) + x(j+2) + c if (c.lt.b) go to 50 c carry generated here a(i) = c - b c = 1 i = i - 1 go to 40 c no carry generated here 50 a(i) = c c = 0 i = i - 1 go to 40 60 if (i.le.0) go to 90 c = y(i+2) + c if (c.lt.b) go to 70 a(i) = 0 c = 1 i = i - 1 go to 60 70 a(i) = c i = i - 1 c no carry possible here 80 if (i.le.0) return a(i) = y(i+2) i = i - 1 go to 80 90 if (c.eq.0) return c must shift right here as carry off end tg1 = tg + 1 do 100 j = 2, tg i = tg1 - j 100 a(i+1) = a(i) a(1) = 1 re = re + 1 return c here do subtraction, abs(y) .gt. abs(x) 110 j = i - med a(i) = c - x(j+2) c = 0 if (a(i).ge.0) go to 120 c borrow generated here c = -1 a(i) = a(i) + b 120 i = i - 1 130 if (i.gt.t) go to 110 140 if (i.le.med) go to 160 j = i - med c = y(i+2) + c - x(j+2) if (c.ge.0) go to 150 c borrow generated here a(i) = c + b c = -1 i = i - 1 go to 140 c no borrow generated here 150 a(i) = c c = 0 i = i - 1 go to 140 160 if (i.le.0) return c = y(i+2) + c if (c.ge.0) go to 70 a(i) = c + b c = -1 i = i - 1 go to 160 end subroutine mpaddi (x, iy, z) c adds multiple-precision x to integer iy giving multiple-precision z. c rounding is controlled by rndrl in common /mpcom/ - see c comments in subroutine mpnzr. common r integer iy, i2, r(1), x(1), z(1), sv c save t etc., allocate temporary space. call mpsavn (sv) call mpnew (i2) c convert iy to multiple-precision and add to x. call mpcim (iy, r(i2)) call mpadd (x, r(i2), z) c restore everything. call mpresn (sv) return end subroutine mpaddq (x, i, j, y) c adds the rational number i/j to mp number x, mp result in y c rounding defined by parameter rndrl in common /mpcom/ - c effect is same as converting i/j to multiple-precision using c mpcqm and then adding to x using mpadd, so see comments in c these routines. common r integer i, i2, j, r(1), sv, x(1), y(1) call mpsavn (sv) c allocate temporary space, compute i/j and add to x call mpnew (i2) call mpcqm (i, j, r(i2)) call mpadd (x, r(i2), y) call mpresn (sv) return end subroutine mpart1 (n, y) c computes mp y = arctan(1/n), assuming integer n .gt. 1. c uses series arctan(x) = x - x**3/3 + x**5/5 - ... c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, ktunfl, dummy integer i, id, i2, i3, ktu, sv, tg integer b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr integer mxsptr, n, r(1), rndrl, sptr, t, y(1) integer mpparn, mptlb c save t etc. call mpsavn (sv) if (n .le. 1) call mperrm (27hn .le. 1 in call to mpart1\$) c increase working precision. if ((n.lt.b).and.(t.gt.2)) t = t - 1 call mpgd3 (mptlb (iabs(r(sv))), tg) c save underflow counter ktu = ktunfl c use truncated arithmetic internally if rndrl .le. 1. if (rndrl.eq.1) rndrl = 0 c allocate space for temporary storage call mpnew (i2) call mpnew (i3) c set sum to 1/n call mpcqm (1, n, r(i3)) c set additive term to 1/n call mpstr (r(i3), r(i2)) i = 1 c mpparn(16) returns a large integer (see mpset). id = (mpparn(16)/n)/n - 2 c main loop. first reduce t if possible 10 t = tg + 2 + r(i2+1) - r(i3+1) if (t.lt.2) go to 40 t = min0 (t, tg) if (rndrl.ne.0) t = tg c if (i+2)*n**2 is not representable as an integer the division c following has to be performed in several steps. if (i.ge.id) go to 20 call mpmulq (r(i2), -i, (i+2)*n*n, r(i2)) go to 30 20 call mpmulq (r(i2), -i, i+2, r(i2)) call mpmuls (r(i2), 1, 1, n, n) 30 i = i + 2 c restore t to working precision t = tg c accumulate sum call mpadd (r(i3), r(i2), r(i3)) c can only fall through next statement if mp underflow occured. if (ktu.eq.ktunfl) go to 10 c restore t etc., round result and return 40 call mpres2 (sv) call mprnd (r(i3), tg, y, iabs(t), 1) call mpresn (sv) return end subroutine mpasin (x, y) c returns y = arcsin(x), assuming abs(x) .le. 1, c for mp numbers x and y. c y is in the range -pi/2 to +pi/2. c method is to use mpatan, so time is o(m(t)t). c rounding options not yet implemented, no guard digits used. common r integer i2, i3, sv integer r(1), x(1), y(1) integer mpcomp, mpget c save t etc., use truncated arithmetic, allocate space. call mpsavn (sv) call mpsetr (0) call mpnew (i2) if (x(1).eq.0) go to 20 if (mpget (x, 2) .le. 0) go to 30 c here abs(x) .ge. 1. see if x = +-1 call mpcim (x(1), r(i2)) if (mpcomp(x, r(i2)).ne.0) go to 10 c x = +-1 so return +-pi/2 call mppi (y) call mpdivi (y, 2*r(i2), y) go to 40 10 call mperrm (32habs(x) .gt. 1 in call to mpasin\$) 20 y(1) = 0 go to 40 c here abs(x) .lt. 1 so use arctan(x/sqrt(1 - x**2)) 30 call mpnew (i3) call mpcim (1, r(i3)) call mpstr (r(i3), r(i2)) c following subtraction is exact if x close to 1, addition exact if c x close to -1. thus (1 - x**2) is computed with small relative error. call mpsub (r(i3), x, r(i3)) call mpadd (r(i2), x, r(i2)) call mpmul (r(i3), r(i2), r(i2)) call mproot (r(i2), -2, r(i2)) call mpmul (x, r(i2), y) call mpatan (y, y) c restore everything and return. 40 call mpresn (sv) return end subroutine mpatan (x, y) c c returns y = arctan(x) for mp x and y, using an o(t.m(t)) method c which could easily be modified to an o(sqrt(t)m(t)) c method (as in mpexp1). y is in the range -pi/2 to +pi/2. c x may be packed or unpacked, y is unpacked. c c for an asymptotically faster method, see - fast multiple- c precision evaluation of elementary functions c (by r. p. brent), j. acm 23 (1976), 242-251, c and the comments in mppigl. c c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. c common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, ktunfl, dummy integer i, i2, i3, i4, ktu, q, sv, tg integer b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr integer mptlb, mxsptr, r(1), rndrl, sptr, t, x(1), y(1) if (x(1).ne.0) go to 10 c arctan (0) = 0. y(1) = 0 return c save t etc. 10 call mpsavn (sv) c increase t, allocate temporary storage. call mpgd3 (mptlb (iabs(t)), tg) call mpnew (i2) call mpnew (i3) if (rndrl.eq.1) rndrl = 0 c move x to temporary storage. call mpmove (x, iabs(r(sv)), r(i3), tg) q = 1 c reduce argument if necessary before using series 20 if (r(i3+1).lt.0) go to 30 if ((r(i3+1).eq.0).and.((2*(r(i3+2)+1)).le.b)) go to 30 q = 2*q call mprevr (r(i3)) call mpmul (r(i3), r(i3), r(i2)) call mpaddi (r(i2), 1, r(i2)) call mpsqrt (r(i2), r(i2)) call mpaddi (r(i2), 1, r(i2)) call mprevr (r(i3)) call mpdiv (r(i3), r(i2), r(i3)) go to 20 c use power series now argument in (-0.5, 0.5) 30 call mpnew (i4) call mpstr (r(i3), r(i4)) ktu = ktunfl call mpmul (r(i3), r(i3), r(i2)) i = 1 c series loop. reduce t if possible. 40 t = tg + 2 + r(i3+1) if (t.le.2) go to 50 t = min0 (t, tg) if (rndrl .ne. 0) t = tg call mpmul (r(i3), r(i2), r(i3)) call mpmulq (r(i3), -i, i+2, r(i3)) i = i + 2 t = tg call mpadd (r(i4), r(i3), r(i4)) c fall through end of loop if underflow occurred above. if (ktu.eq.ktunfl) go to 40 c correct for argument reduction. 50 t = tg call mpmuli (r(i4), q, r(i4)) c round result, restore t etc. and return. call mprnd (r(i4), tg, y, iabs(r(sv)), 1) call mpresn (sv) return end subroutine mpatn2 (x, y, z) c sets z = arctan (x/y) if y nonzero, c = pi/2 if y zero, c for (packed or unpacked) mp x, unpacked mp y. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer i2, i3, sv, tg, z(1) integer b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1) c save t etc., increase working precision, allocate space. call mpsavn (sv) call mpgd3 (1, tg) call mpnew (i2) c check for y zero. if (y(1).ne.0) go to 10 call mppi (r(i2)) call mpdivi (r(i2), 2, r(i2)) go to 20 c here y nonzero, increase m so x/y does not over/underflow. 10 m = 2*m + 2 call mpnew (i3) call mpmove (x, iabs(r(sv)), r(i2), tg) call mpmove (y, iabs(r(sv)), r(i3), tg) call mpdiv (r(i2), r(i3), r(i2)) c free some space. sptr = i3 c use mpatan(x/y) call mpatan (r(i2), r(i2)) c round result, restore t etc. and return. 20 call mprnd (r(i2), tg, z, iabs(r(sv)), 1) call mpresn (sv) return end integer function mpbasa (x) c returns the mp base (first word in common /mpcom/). c x is a dummy mp argument. integer mpparn, x(1) mpbasa = mpparn (1) return end subroutine mpbasb (i, x) c sets the mp base (first word of common /mpcom/) to i. c i should be an integer such that i .ge. 2 c and (8*i*i-1) is representable as a single-precision integer. c x is a dummy mp argument (augment expects one). c warning - setting the base does not automatically convert mp c numbers to the new base. this may be done by converting to c decimal (using mpfout), changing the base (using mpbasb), c then converting from decimal to the new base (using mpfin). integer i, x(1) call mpparc (1, i) return end subroutine mpbern (n, p, x) c c computes the bernoulli numbers b2 = 1/6, b4 = -1/30, c b6 = 1/42, b8 = -1/30, b10 = 5/66, b12 = -691/2730, etc., c defined by the generating function y/(exp(y)-1). c n and p are single-precision integers, with 2*p .ge. t+2. c x should be a one-dimensional integer array of dimension at c least p*n. the bernoulli numbers b2, b4, ... , b(2n) are c returned in packed format in x, with b(2j) in locations c x((j-1)*p+1), ... , x(p*j). thus, to get b(2j) in usual c mp format in y, one should call mpunpk (x(ix), y) after c calling mpbern, where ix = (j-1)*p+1. c c alternatively (simpler but nonstandard) - c x may be a two-dimensional integer array declared with c dimension (p, n1), where n1 .ge. n and 2*p .ge. t+2. c then b2, b4, ... , b(2n) are returned in packed format in c x, with b(2j) in x(1,j), ... , x(p,j). thus, to get c b(2j) in usual mp format in y one should c call mpunpk (x(1, j), y) after calling mpbern. c c the well-known recurrence is unstable (losing about 2j bits c of relative accuracy in the computed b(2j)), so we use a c different recurrence derived by equating coefficients in c (sinh(y)/y)*(sigma b(2j)*(2y)**(2j)/factorial(2j)) = cosh y, c where the summation is from j = 0 to infinity. this method c was suggested by christian reinsch, and is faster than the method c used in earlier versions of mpbern. the relation c b(2j) = -2*((-1)**j)*factorial(2j)*zeta(2j)/((2pi)**(2j)) c is used if zeta(2j) is equal to 1 to working accuracy. c a different method is given by knuth and buckholtz in c math. comp. 21 (1967), 663-688. c c time is o(t*(min(n, t)**2) + n*m(t)). c c rounding options not implemented, no guard digits used. c the relative error in b(2j) is o((j**2)*(b**(1-t))). c c if n is negative, abs(n) bernoulli numbers are returned, but c the comment above about relative error no longer applies - c instead the precision decreases linearly from the first to the c last bernoulli number (this is usually sufficient if the c bernoulli numbers are to be used as coefficients in an c euler-maclaurin expansion). c common r common /mpcom/ b, t, dummy integer i, ix, i2, i3, i4, j, nl, np, n2, sv integer b, dummy(21), mptlb, n, p, r(1), t, x(1) nl = iabs(n) if (nl.le.0) return call mpsavn (sv) if ((2*p) .lt. (r(sv)+2)) call mperrm ( \$ 30hp too small in call to mpbern\$) c allocate working space, use truncated arithmetic. call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpsetr (0) c compute upper limit for recurrence relation method. n2 = min0 (nl, mptlb(iabs(r(sv)))/2) c set all results to zero np = nl*p do 30 i = 1, np 30 x(i) = 0 call mpcqm (1, 12, r(i2)) call mpstr (r(i2), r(i3)) c main loop to generate scaled bernoulli numbers do 70 j = 1, n2 c decrease t if n negative. if (n .lt. 0) t = min0 (r(sv), ((nl+1-j)*(r(sv)-2))/nl + 4) ix = (j-1)*p + 1 call mppack (r(i3), x(ix)) if (j.ge.n2) go to 80 call mpmuls (r(i2), 1, 1, 4*j, 4*j+6) call mpstr (r(i2), r(i3)) do 60 i = 1, j c change t if n negative. if (n .lt. 0) t = min0 (r(sv), ((nl+1-i)*(r(sv)-2))/nl + 4) ix = (i-1)*p + 1 call mpunpk (x(ix), r(i4)) call mpmuls (r(i4), 1, 1, 4*(j-i)+4, 4*(j-i)+6) call mppack (r(i4), x(ix)) 60 call mpsub (r(i3), r(i4), r(i3)) 70 continue c now unscale results 80 t = r(sv) call mpcim (1, r(i2)) if (n2.le.1) go to 100 i = n2 90 call mpmuls (r(i2), 4*(n2-i)+4, 4*(n2-i)+6, 1, 1) i = i - 1 ix = (i-1)*p + 1 call mpunpk (x(ix), r(i4)) call mpmul (r(i2), r(i4), r(i4)) call mppack (r(i4), x(ix)) if (i.gt.1) go to 90 c now have b(2j)/factorial(2j) in x call mpcim (1, r(i2)) 100 do 110 i = 1, n2 call mpmuls (r(i2), 2*i-1, 2*i, 1, 1) ix = (i-1)*p + 1 call mpunpk (x(ix), r(i4)) call mpmul (r(i2), r(i4), r(i4)) 110 call mppack (r(i4), x(ix)) c return if finished if (nl.le.n2) go to 130 c else compute remaining numbers call mppi (r(i3)) call mppwr (r(i3), -2, r(i3)) call mpdivi (r(i3), -4, r(i3)) n2 = n2 + 1 do 120 i = n2, nl call mpmul (r(i4), r(i3), r(i4)) call mpmuls (r(i4), 2*i-1, 2*i, 1, 1) ix = (i-1)*p + 1 120 call mppack (r(i4), x(ix)) c restore stack pointer etc. and return. 130 call mpresn (sv) return end subroutine mpbes2 (x, nu, y) c uses the backward recurrence method to evaluate y = j(nu,x), c where x and y are mp numbers, nu (the index) is an integer, c and j is the bessel function of the first kind. assumes that c 0 .le. nu (not too large) and x .gt. 0. c for normalization the identity c j(0,x) + 2*j(2,x) + 2*j(4,x) + ... = 1 is used. c called by mpbesj and not recommended for independent use. c rounding options not implemented, uses no guard digits. common r common /mpcom/ b, t, dummy integer b, t, dummy(21), mpchgb, mpcmpi integer i, i2, i3, i3s, i4, i5, i6, nu1, sv integer nu, r(1), x(1), y(1), mpparn c check legality of nu and x if ((nu.ge.0) .and. (x(1).eq.1)) go to 20 10 call mperrm (34hnu or x illegal in call to mpbes2\$) 20 call mpsavn (sv) c allocate temporary space call mpnew (i2) call mpnew (i3) call mpnew (i4) c use truncated arithmetic call mpsetr (0) c use about 6 decimal places to compute starting point t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1)) call mpcim (max0 (1, nu), r(i2)) call mpmuli (r(i2), 2, r(i3)) call mpdiv (r(i3), x, r(i3)) call mpln (r(i3), r(i3)) call mpaddi (r(i3), -1, r(i3)) call mpcim (1, r(i4)) call mpmax (r(i3), r(i4), r(i3)) call mpmul (r(i2), r(i3), r(i3)) call mplni (iabs(b), r(i2)) call mpmulq (r(i2), iabs(r(sv)), 2, r(i2)) call mpadd (r(i2), r(i3), r(i2)) call mpdiv (r(i2), x, r(i2)) c 125/92 .lt. e/2 call mpmulq (r(i2), 92, 125, r(i2)) call mpcim (2, r(i4)) call mpmax (r(i2), r(i4), r(i2)) call mpstr (r(i2), r(i3)) c do two newton iterations do 30 i = 1, 2 call mpadd (r(i2), r(i3), r(i4)) call mpln (r(i3), r(i3)) call mpaddi (r(i3), 1, r(i3)) 30 call mpdiv (r(i4), r(i3), r(i3)) call mpmul (r(i3), x, r(i2)) c 34/25 .gt. e/2 call mpmulq (r(i2), 34, 25, r(i2)) if (mpcmpi (r(i2), mpparn(16)-2) .gt. 0) go to 10 call mpcmi (r(i2), nu1) nu1 = nu1 + 2 c now restore t t = r(sv) call mpnew (i5) call mpnew (i6) call mpcim (mod(nu1+1,2), r(i6)) call mprec (x, r(i2)) call mpmuli (r(i2), 2, r(i2)) r(i3) = 0 call mpcim (1, r(i4)) c backward recurrence loop 40 call mpmul (r(i4), r(i2), r(i5)) call mpmuli (r(i5), nu1, r(i5)) call mpsub (r(i5), r(i3), r(i5)) nu1 = nu1 - 1 c faster to interchange pointers than mp numbers i3s = i3 i3 = i4 i4 = i5 i5 = i3s if (mod(nu1,2) .ne. 0) go to 50 c nu1 even so update normalizing sum if (nu1.eq.0) call mpmuli (r(i6), 2, r(i6)) call mpadd (r(i6), r(i4), r(i6)) c save unnormalized result if nu1 .eq. nu 50 if (nu1.eq.nu) call mpstr (r(i4), y) if (nu1.gt.0) go to 40 c normalize result and return call mpdiv (y, r(i6), y) c restore stack pointer etc. and return. call mpresn (sv) return end subroutine mpbesj (x, nu, y) c returns y = j(nu,x), the first-kind bessel function of order nu, c for small integer nu, mp x and y. c the method is hankels asymptotic expansion if c abs(x) large, the power series if abs(x) small, and the c backward recurrence method otherwise. c results for negative arguments are defined by c j(-nu,x) = j(nu,-x) = ((-1)**nu)*j(nu,x). c error could be induced by o(b**(1-t)) perturbations c in x and y. time is o(t.m(t)) for fixed x and nu, increases c as x and nu increase, unless x large enough for asymptotic c series to be used. c rounding options not yet implemented, uses truncated arithmetic. common r common /mpcom/ b, t, dummy logical error integer ie, it, i2, i3, i4, k, mxint, nua, sv, tg, tm, \$ b, dummy(21), nu, r(1), t, x(1), y(1), \$ mpchgb, mpget, mpparn nua = iabs(nu) c check for x zero if (x(1).ne.0) go to 10 c j(nu,0) = 0 if nu .ne. 0, 1 if nu .eq. 0 y(1) = 0 if (nua.eq.0) call mpcim (1, y) return c save t etc., use truncated arithmetic. 10 call mpsavn (sv) call mpsetr (0) mxint = mpparn(16) c see if abs(x) so large that no accuracy possible if (mpget (x, 2).ge.t) call mperrm ( \$ 35habs(x) too large in call to mpbesj\$) c x nonzero so try hankel asymptotic series with guard digits. call mpgd3 (1, tg) call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) call mphank (r(i2), nua, r(i2), error) c go to rounding of final result if asymptotic series accurate. if (.not. error) go to 40 c asymptotic series not good enough so restore t etc. call mpresn (sv) call mpsavn (sv) call mpsetr (0) c may need to increase t later so prepare for this c max allowable t is approximately double tm = 2*tg tg = t c no appreciable cancellation in power series if abs(x) .lt. 1 if (mpget (x, 2).le.0) go to 20 c estimate number of digits required to compensate for cancellation. c first reduce t to equivalent of 6 decimal places (could use c single-precision real arithmetic here if trusted). t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1)) call mpnew (i2) call mpnew (i3) call mpabs (x, r(i2)) call mpdivi (r(i2), 2, r(i3)) call mpln (r(i3), r(i3)) call mpmulq (r(i3), 2*nua + 1, 2, r(i3)) call mpadd (r(i2), r(i3), r(i2)) call mplni (iabs(b), r(i3)) call mpdiv (r(i2), r(i3), r(i2)) call mpcmi (r(i2), it) call mpresn (sv) call mpsavn (sv) call mpsetr (0) tg = t + max0 (0, it + 1) c if need more digits than space allows for power series then c use recurrence method instead if (tg.gt.tm) go to 50 c prepare for power series loop 20 t = tg call mpnew (i2) call mpnew (i3) call mpnew (i4) t = r(sv) call mpdivi (x, 2, r(i4)) call mppwra (r(i4), nua, r(i4)) call mpgamq (nua+1, 1, r(i2)) call mpdiv (r(i4), r(i2), r(i4)) call mpmul (x, x, r(i3)) call mpdivi (r(i3), -4, r(i3)) c clear trailing digits of r(i3) and r(i4). call mpmove (r(i3), iabs(r(sv)), r(i3), tg) call mpmove (r(i4), iabs(r(sv)), r(i4), tg) call mpmove (r(i4), tg, r(i2), tg) ie = r(i2+1) k = 0 c power series loop, reduce t if possible 30 t = min0 (tg, tg + 2 + r(i4+1) - ie) if (t.lt.2) go to 40 call mpmul (r(i3), r(i4), r(i4)) k = k + 1 if ((k+nua) .gt. (mxint/2)) call mperrm ( \$ 36habs(nu) too large in call to mpbesj\$) call mpmuls (r(i4), 1, 1, k, k+nua) c restore t for addition t = tg call mpadd (r(i2), r(i4), r(i2)) if ((r(i4).ne.0).and.(r(i4+1).ge.(r(i2+1)-r(sv)))) go to 30 c restore t etc. and round final result 40 call mpres2 (sv) call mprnd (r(i2), tg, y, iabs(r(sv)), 0) c now restore stack pointer call mpresn (sv) c correct sign if nu odd and negative if ((nu.lt.0).and.(mod(nua,2).ne.0)) y(1) = -y(1) return c here use backward recurrence method with guard digits 50 call mpgd3 (1, tg) call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) r(i2) = iabs(r(i2)) call mpbes2 (r(i2), nua, r(i2)) c correct sign if nua odd if (mod (nua,2) .ne. 0) r(i2) = x(1)*r(i2) go to 40 end subroutine mpcam (a, x) c converts the hollerith string a to an mp number x. c a can be a string of digits acceptable to routine mpin c and terminated by a dollar (\$), e.g. 7h-5.367\$, c or one of the following special strings - c eps (mp machine-precision, see mpeps), c eul (eulers constant 0.5772..., see mpeul), c maxr (largest valid mp number, see mpmaxr), c minr (smallest postive mp number, see mpminr), c pi (pi = 3.14..., see mppi). c only the first two characters of these strings are checked. c worst-case space = 2*n + o(t) words, where n is the length in c characters of the string a. c the error message c *** mxr too small ... *** c is usually caused by omission of the sentinel \$ common r integer a(1), d(2), i2, i3, len, n, r(1), sv, x(1), \$ ha, he, hi, hm, hp, hu logical error, mpis data ha, he, hi, hm, hp, hu \$ /1ha, 1he, 1hi, 1hm, 1hp, 1hu/ call mpsavn (sv) call mpnew (i2) r(i2) = 0 len = r(sv) + 2 c unpack first 2 characters of a call mpupk (a, d, 2, n) if (n.ne.2) go to 10 c check for special strings if (mpis (d(1), he) .and. mpis (d(2), hp)) call mpeps (r(i2)) if (mpis (d(1), he) .and. mpis (d(2), hu)) call mpeul (r(i2)) if (mpis (d(1), hm) .and. mpis (d(2), ha)) call mpmaxr (r(i2)) if (mpis (d(1), hm) .and. mpis (d(2), hi)) call mpminr (r(i2)) if (mpis (d(1), hp) .and. mpis (d(2), hi)) call mppi (r(i2)) c if r(i2) nonzero one of above tests was successful. if (r(i2) .eq. 0) go to 10 call mpstr (r(i2), x) go to 20 c len is guess at string length - double and try again if too short. c first allocate necessary extra space. 10 call mpnew2 (i3, len) len = 2*len c unpack up to len characters of string call mpupk (a, r(i2), len, n) c if n .lt. len then len was large enough, else loop if (n.ge.len) go to 10 c convert unpacked string to mp number. call mpin (r(i2), x, n, error) if (error) call mperrm ( \$ 45herror in hollerith constant in call to mpcam\$) c restore stack pointer and return 20 call mpresn (sv) return end subroutine mpcdm (dx, z) c converts double-precision number dx to multiple-precision z. c some numbers will not convert exactly on machines c with base other than two or if b is not a power of two with c b**(t-1) sufficiently large. thus mpcdm should be used only to c obtain starting approximations etc. for accurate initialization c of mp numbers use mpcim, mpcqm, mpqpwr, or mpin. c warning - the parameter rndrl in common /mpcom/ has no effect. common /mpcom/ b, t, dummy integer b, dummy(21), i, ie, isv, j, re, rs, sv, t, \$ three, z(1) double precision db, dj, dx, dble real float call mpsavn (sv) call mpsetr (0) three = 3 c check sign if (dx) 20, 10, 30 c if dx = 0d0 return 0 10 z(1) = 0 go to 100 c dx .lt. 0d0 20 rs = -1 dj = -dx go to 40 c dx .gt. 0d0 30 rs = 1 dj = dx c we want to reduce dj to range 1/16 .le. dj .lt. 1. 40 ie = 0 50 if (dj.lt.1d0) go to 60 c increase ie and divide dj by 16. ie = ie + 1 dj = 0.0625d0*dj go to 50 60 if (dj.ge.0.0625d0) go to 70 ie = ie - 1 dj = 16d0*dj go to 60 c now dj is dy divided by suitable power of 16. c set exponent to 0 70 re = 0 c dfloat is not ansi standard, believe it or not db = dble(float(b)) c conversion loop (assume single-precision ops. exact) do 80 i = 1, t dj = db*dj isv = i j = idint(dj) c check for (unlikely) effect of strange floating-point arithmetic. if (j.ge.b) go to 110 if (dj.lt.0.0d0) go to 120 z(i+2) = j 80 dj = dj - dble(float(j)) c normalize result 90 call mpnzr (rs, re, z(1), z(three), 0) c now multiply by 16**ie call mpscal (z, 16, ie) c can return now 100 call mpresn (sv) return c here a digit is .ge. b, which should never happen 110 j = b-1 go to 130 c here a digit is negative, which should never happen 120 j = 0 c set remaining digits to j (either b-1 or 0) 130 do 140 i = isv, t 140 z(i+2) = j go to 90 end subroutine mpceil (x, y) c sets y = ceiling (x), i.e. the smallest integer not less than x. c x and y are mp numbers. c rounding is defined by the parameter rndrl in common /mpcom/ c (this is only relevant if x is large and positive) - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r integer i2, mpcomp, r(1), sv, x(1), y(1) c save t etc. and truncate x (towards zero) to an integer. call mpsavn (sv) call mpnew (i2) call mpcmim (x, r(i2)) c if x positive and not an integer need to add 1 if ((x(1).gt.0) .and. (mpcomp (x, r(i2)).ne.0)) \$ call mpaddi (r(i2), 1, r(i2)) c move result, restore everything and return. call mpstr (r(i2), y) call mpresn (sv) return end subroutine mpcheb (c, nc, n, ind) c c converts the power series coefficients c(1), ... , c(n) to c chebyshev series coefficients. (it is assumed that the constant c term in the chebyshev sum is halved.) c ind = 0 means that c(i) is the coefficient of x**(i-1) on input, c of t(i-1)(x) on output, c ind = -1 means that c(i) is the coefficient of x**(2i-1) on input, c of t(2i-1)(x) on output, c ind = +1 means that c(i) is the coefficient of x**(2i-2) on input, c of t(2i-2) on output. c c is an array of mp variables with first dimension nc (for c augment users, nc = mt2). c integer i, i1, j, jj, j1, j2, sv c integer n, nc (univac fortran v does not allow this) integer c(nc,n), ind, mpparn call mpsavn (sv) c use truncated arithmetic and no guard digits. call mpsetr (0) c check legality of nc, n and ind. if ((nc .lt. (mpparn(2)+2)) .or. (n .le. 0) .or. \$ (iabs(ind) .gt. 1)) call mperrm ( \$ 36hillegal arguments on call to mpcheb\$) do 40 jj = 1, n j = n + 1 - jj call mpmuli (c(1,j), 2, c(1,j)) j2 = j + 2 - iabs(ind) if (j2 .le. n) call mpadd (c(1,j), c(1,j2), c(1,j)) j1 = j + 1 if (j1 .gt. n) go to 20 do 10 i = j1, n i1 = i + 2 - iabs(ind) if (i1 .le. n) call mpadd (c(1,i), c(1,i1), c(1,i)) 10 call mpdivi (c(1,i), 2, c(1,i)) 20 if (((j .eq. 1) .and. (ind .eq. 1)) .or. (ind .eq. 0)) \$ go to 40 do 30 i = j, n if (i .lt. n) call mpadd (c(1,i), c(1,i+1), c(1,i)) 30 call mpdivi (c(1,i), 2, c(1,i)) 40 continue call mpresn (sv) return end subroutine mpchev (c, nc, n, ind, x, y) c returns y = chebyshev series evaluated at x. c x and y are mp variables, mp coefficients are in c, c for a description of c, nc, n and ind see mpcheb. common r integer i, ii, isv, i2, i3, i4, i5, sv c integer n, nc (univac fortran v does not allow this) integer c(nc,n), ind, mpparn, r(1), x(1), y(1) call mpsavn (sv) c use truncated arithmetic, no guard digits. call mpsetr (0) c check legality of nc, n and ind. if ((nc .lt. (mpparn(2)+2)) .or. (n .le. 0) .or. \$ (iabs(ind) .gt. 1)) call mperrm ( \$ 36hillegal arguments on call to mpchev\$) c allocate working space call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpnew (i5) r(i2) = 0 r(i3) = 0 call mpmuli (x, 2, r(i5)) if (ind .eq. 0) go to 10 call mpmul (r(i5), r(i5), r(i5)) call mpaddi (r(i5), (-2), r(i5)) 10 do 20 ii = 1, n i = n + 1 - ii c interchange pointers rather than mp numbers. isv = i4 i4 = i3 i3 = i2 i2 = isv call mpmul (r(i5), r(i3), r(i2)) call mpsub (r(i2), r(i4), r(i2)) 20 call mpadd (r(i2), c(1, i), r(i2)) if (ind .ge. 0) go to 30 call mpsub (r(i2), r(i3), r(i2)) call mpmul (x, r(i2), y) go to 40 30 call mpsub (r(i2), r(i4), r(i2)) call mpdivi (r(i2), 2, y) c restore everything and return. 40 call mpresn (sv) return end integer function mpchgb (b1, b2, n) c returns j such that b1**abs(j) .ge. b2**abs(n), c i.e. abs(j) .ge. abs(n)*log(b2)/log(b1), c and sign(j) = sign(n), c assuming b1 .gt. 1, b2 .ge. 1, b1*b2 .le. mxint . c usually the value of j returned is close to minimal. integer b1, b2, j, k, mb, mpparn, mul, mxint, n, n2, q mxint = mpparn(16) mb = mxint/b1 if ((b1 .le. 1) .or. (b2 .le. 0) .or. (b2 .gt. mb)) \$ call mperrm (36hillegal arguments in call to mpchgb\$) j = 0 if ((n .eq. 0) .or. (b2 .eq. 1)) go to 60 k = 0 q = 1 c the constant 100 is large enough to give reasonable accuracy. c if it is increased the time required will increase. mul = (iabs(n)-1)/100 + 1 n2 = (iabs(n)-1)/mul + 1 c increase j and k, keeping q .le. b1**j/b2**k 10 if (k .ge. n2) go to 40 20 if (q .gt. mb) go to 30 q = b1*q j = j + 1 go to 20 30 q = q/b2 k = k + 1 go to 10 c final stage, may have overshot 40 if (q .lt. b1) go to 50 q = q/b1 j = j - 1 go to 40 c check if result would overflow 50 if (j .gt. (mxint/mul)) call mperrm ( \$ 30hn too large in call to mpchgb\$) j = isign (mul*j, n) 60 mpchgb = j return end subroutine mpchk c checks legality of parameters which should be set in c common /mpcom/, also updates mxsptr. common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, ktunfl, mxunfl, decpl, mt2, \$ mxint, exwid, inrecl, inbase, outbas, expch, chword, onescp integer b, chword, decpl, expch, exwid, inbase, inrecl, ktunfl, \$ lun, m, mnexpn, mnsptr, mt2, mxexpn, mxint, mxr, mxsptr, mxunfl, \$ onescp, outbas, rndrl, sptr, t, ib, jsp c first check that lun in range 1 to 99, if not print error c message on logical unit 6. if ((0.lt.lun) .and. (lun.lt.100)) go to 20 write (6, 10) lun 10 format (10h *** lun =, i10, 29h illegal in call to mpchk ***) lun = 6 call mperr c now check legality of b, t and m 20 if (b.gt.1) go to 40 write (lun, 30) b 30 format (8h *** b =, i10, 29h illegal in call to mpchk ***) call mperr c t may exceed mt2-2 in call from other mp routines 40 if ((t.gt.1). and. ((t.le.(mt2-2)) .or. (sptr.gt.mnsptr))) \$ go to 60 write (lun, 50) t 50 format (8h *** t =, i10, 29h illegal in call to mpchk ***) call mperr 60 if (m.gt.(4*t)) go to 80 c here m is not greater than 4*t. since many mp routines increase c t internally, it is safest to set m much larger than 4*t c initially. note, though, that 4*m must not overflow. write (lun, 70) 70 format (36h *** m .le. 4*t in call to mpchk ***) call mperr c 8*b*b-1 should be representable, if not will overflow c and may become negative, so check for this 80 ib = 4*b*b - 1 if ((ib.gt.0) .and. ((2*ib+1).gt.0)) go to 100 write (lun, 90) 90 format (37h *** b too large in call to mpchk ***) call mperr c check that stack space in common is sufficient. c sptr should point to first free word in stack, c and mxr to last available word. c at least 80 words must be available (see mperrm). 100 jsp = max0 (sptr, mnsptr + 80) - 1 if (mxr.ge.jsp) go to 120 c here stack is too small. call mpstov to expand it if possible. call mpstov (jsp - mxr) c see if mpstov actually did expand enough. if (mxr.ge.jsp) go to 120 c here common is too small (assuming its size is mxr). write (lun, 110) mxr, jsp 110 format (51h *** mxr too small or not set to dim(r) before call, \$ 21h to an mp routine *** / \$ 10h *** mxr =, i6, 24h, but should be at least, i6, 5h ***) call mperr c check legality of various other parameters in common /mpcom/ 120 if ((mxsptr.ge.mnsptr) .and. (mnsptr.gt.0) .and. \$ (mnsptr.le.sptr) .and. ((mxsptr-1).le.mxr) .and. \$ (rndrl.ge.0) .and. (rndrl.le.3) .and. \$ (ktunfl.ge.0) .and. (mxunfl.ge.0)) go to 140 write (lun, 130) 130 format (38h *** one of sptr, ... , mxunfl illegal, \$ 21h in call to mpchk ***) call mperr 140 if ((decpl .gt. 0) .and. (mt2 .gt. 3) .and. (mxint .ge. 2047) \$ .and. (exwid .gt. 2) .and. (inrecl .gt. 0) \$ .and. (inrecl .le. 80) .and. (inbase .gt. 1) \$ .and. (inbase .le. 16) .and. (outbas .gt. 1) \$ .and. (outbas .le. 16) .and. (chword .gt. 0)) go to 160 write (lun, 150) 150 format (39h *** one of decpl, ... , chword illegal, \$ 21h in call to mpchk ***) call mperr c update mxsptr and return 160 mxsptr = max0 (mxsptr, sptr) return end subroutine mpcim (ix, z) c converts integer ix to multiple-precision z. c assumes that abs(ix) .le. b**t (otherwise ix can not usually be c represented exactly as a multiple-precision number). common /mpcom/ b, t, dummy integer b, dummy(21), i, ix, mpgd, n, t, two, z(1) two = 2 c check legality of parameters in common /mpcom/ call mpchk c t should be increased if abs(ix) .lt. b**t . n = ix if (mpgd(n).gt.t) call mperrm ( \$ 29ht too small in call to mpcim\$) c set z(1) to sign of ix. z(1) = 0 if (n.eq.0) return z(1) = isign (1, n) c set exponent to t z(two) = t c clear fraction do 10 i = 2, t 10 z(i+1) = 0 c insert n z(t+2) = iabs(n) c normalize by calling mpmuli call mpmuli (z, 1, z) return end subroutine mpcis (x, c, s, both) c if both = .true., returns c = cos(x) and s = sin(x). c if both = .false., returns c = cos(x) only (s unchanged). c x, c and s are mp numbers, both is logical. c x may be packed or unpacked, c and s are unpacked. c the algorithm is described in - r. p. brent, unrestricted c algorithms for elementary and special functions, in information c processing 80, north holland, 1980. c time = o(sqrt(t)m(t)). c rounding options are implemented as follows - c rndrl = 0 or 1 - absolute error less than c 0.6*b**(-t) (but the relative error c may be large if x is close to a nonzero c multiple of pi/2), c rndrl = 2 - lower bound on true result, c rndrl = 3 - upper bound on true result. common r common /mpcom/ b, t, dummy integer i2, i3, i4, i5, j, k, mpt, q, sv, tg, tm, two, \$ b, c(1), dummy(21), mpget, mptlb, r(1), s(1), t, x(1) logical both two = 2 if (x(1) .ne. 0) go to 10 c here x = 0 call mpcim (1, c) if (both) s(1) = 0 return c here x .ne. 0 10 call mpsavn (sv) c use truncated arithmetic internally call mpsetr (0) c select optimal q q = 0 mpt = mptlb (1) 20 q = q + 1 c the constant 5 was determined empirically if ((mpt*q*q) .lt. (5*t)) go to 20 c use sufficient guard digits t = t + max0 (0, x(two)) call mpgd3 (mptlb(q), tg) c allocate temporary storage and move x call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) c divide by suitable power of two k = 0 30 if (r(i2+1) .le. (-q)) go to 40 k = k + 1 call mpdivi (r(i2), 2, r(i2)) go to 30 c allocate more storage and set up inner loop 40 call mpnew (i3) if (both) call mpnew (i4) call mpnew (i5) if (.not. both) call mpmul (r(i2), r(i2), r(i2)) r(i3) = 0 if (both) r(i4) = 0 call mpcim (1, r(i5)) j = 0 c inner loop starts here 50 j = j + 2 c can reduce t to tm for multiplications tm = min0 (tg, max0 (2, r(i5+1) + tg)) t = tm call mpmul (r(i5), r(i2), r(i5)) call mpdivi (r(i5), j-1, r(i5)) t = tg if (both) call mpadd (r(i4), r(i5), r(i4)) t = tm if (both) call mpmul (r(i5), r(i2), r(i5)) call mpdivi (r(i5), (-j), r(i5)) t = tg call mpadd (r(i3), r(i5), r(i3)) c check for underflow if (r(i5) .eq. 0) go to 60 c check for convergence if (r(i5+1) .gt. (-tg)) go to 50 c now use doubling identities to get result 60 if (k .le. 0) go to 90 do 80 j = 1, k if (.not.both) go to 70 call mpmul (r(i3), r(i4), r(i2)) call mpadd (r(i4), r(i2), r(i4)) call mpmuli (r(i4), 2, r(i4)) 70 call mpaddi (r(i3), 2, r(i2)) call mpmul (r(i2), r(i3), r(i3)) 80 call mpmuli (r(i3), 2, r(i3)) c add 1 to get cos 90 call mpaddi (r(i3), 1, r(i3)) c round result(s) if (r(sv+2) .le. 1) go to 100 c fix up for directed roundings here call mpsetr (iabs(r(sv+2))) call mpcim (2*r(sv+2)-5, r(i2)) r(i2+1) = 1 - r(sv) call mpadd (r(i3), r(i2), r(i3)) if (both) call mpadd (r(i4), r(i2), r(i4)) 100 call mprnd (r(i3), tg, c, iabs(r(sv)), 0) t = r(sv) if (c(1) .eq. 0) go to 110 if (c(two) .gt. 0) call mpcim (mpget (r, i3), c) 110 if (.not. both) go to 120 call mprnd (r(i4), tg, s, iabs(r(sv)), 0) if (s(1) .eq. 0) go to 120 if (s(two) .gt. 0) call mpcim (mpget (r, i4), s) c restore t etc and return 120 call mpresn (sv) return end subroutine mpcmd (x, dz) c converts multiple-precision x to double-precision dz. c assumes x in allowable range. common /mpcom/ b, t, dummy integer b, dummy(21), i, ie, t, tm, two, x(1) double precision db, dble, dz, dz2 real float two = 2 c check legality of parameters in common /mpcom/ call mpchk dz = 0d0 c return with dz = 0.0 if x is zero. if (x(1).eq.0) return db = dble(float(b)) c loop to compute dz. do 10 i = 1, t dz = db*dz + dble(float(x(i+2))) tm = i c check if full double-precision accuracy attained dz2 = dz + 1d0 c on some machines (dz2 .le. dz) fails if ((dz2 - dz) .le. 0.0d0) go to 20 10 continue c now allow for exponent 20 ie = x(two) - tm 30 if (ie. le. 0) go to 40 ie = ie - 1 dz = dz*db go to 30 40 if (ie .eq. 0) go to 50 ie = ie + 1 dz = dz/db go to 40 c check reasonableness of result c following message indicates that x is too large or small - c try using mpcmde instead. 50 if (dz.le.0d0) call mperrm ( \$ 39hfloating-point over/underflow in mpcmd\$) c allow for sign of x and return. if (x(1).lt.0) dz = -dz return end subroutine mpcmde (x, n, dx) c returns integer n and double-precision dx such that mp c x = dx*outbas**n (approximately), where 1 .le. abs(dx) .lt. outbas c unless dx = 0. (outbas is in common /mpcom/ - default value c is 10.) c rounding options not implemented. common r integer i2, mpparn, n, r(1), sv, x(1) double precision dabs, dble, dx real float call mpsavn (sv) if (x(1).ne.0) go to 10 n = 0 dx = 0.0d0 go to 20 c use truncated arithmetic. 10 call mpsetr (0) call mpnew (i2) call mpcmef (x, n, r(i2)) call mpcmd (r(i2), dx) if (dabs(dx) .lt. dble(float(mpparn(20)))) go to 20 c here dx was rounded up to outbas n = n + 1 dx = dble(float(r(i2))) 20 call mpresn (sv) return end subroutine mpcmef (x, n, y) c c given mp x, returns integer n and mp y such that c x = (outbas**n)*y and 1 .le. abs(y) .lt. outbas c (unless x .eq. 0, when n .eq. 0 and y .eq. 0). c outbas is a parameter in common /mpcom/, default value 10. c it is assumed that x is not so large or small that n overflows. c x may be packed or unpacked, y is unpacked. c c rounding options not yet fully implemented, but the directed c rounding options (rndrl = 2 and 3) give correct results. c common r common /mpcom/ b, t, m, dummy integer iy, i2, j, n1, outbas, sv, two integer b, dummy(20), m, n, r(1), t, x(1), y(1) integer mpchgb, mpcmpi, mpparn two = 2 c save t etc. call mpsavn (sv) outbas = mpparn (20) if ((outbas .lt. 2) .or. (outbas .gt. 16)) call mperrm ( \$ 33hillegal outbas on call to mpcmef\$) c check for x zero if (x(1).ne.0) go to 10 n = 0 y(1) = 0 go to 110 c x nonzero here. 10 call mpunpk (x, y) n = 0 call mpnew (i2) c loop up to 100 times (usually one is sufficient) do 20 j = 1, 100 c estimate log (abs(y)) to base outbas n1 = mpchgb (outbas, iabs(b), y(two)-1) + \$ mpchgb (outbas, y(two+1), 1) c following avoids possibility of r(i2) overflowing below if ((j.eq.1).and.(iabs(n1).gt.(m/4))) n1 = n1/2 c leave j loop if n1 small if (iabs(n1).le.3) go to 50 c divide by outbas**n1, taking care for directed roundings. call mpscal (y, outbas, (-n1)) if (y(1) .eq. 0) go to 30 20 n = n + n1 c error if fall through loop 30 call mperrm ( \$ 42herror in mpcmef, maybe exponent too large\$) 50 if (y(1).eq.0) go to 30 c loop dividing by outbas until abs(y) .lt. 1 60 if (y(two).le.0) go to 80 n = n + 1 call mpdivi (y, outbas, y) go to 60 c loop multiplying by outbas until abs(y) .ge. 1 70 if (y(two).gt.0) go to 90 80 n = n - 1 call mpmuli (y, outbas, y) go to 70 c check for possibility that rounding up was to outbas 90 iy = y(1) y(1) = 1 if (mpcmpi (y, outbas) .lt. 0) go to 100 c it was, so set y to 1 and add one to exponent call mpcim (1, y) n = n + 1 c fix up sign of y. 100 y(1) = iy c restore stack pointer etc. and return. 110 call mpresn (sv) return end subroutine mpcmf (x, y) c for mp x and y, returns fractional part of x in y, c i.e., y = x - integer part of x (truncated towards 0). c note that rndrl in common /mpcom/ is irrelevant. common /mpcom/ b, t, dummy integer b, dummy(21), i, t, two, x(1), y(1), y1, y2 two = 2 c move x to y call mpstr (x, y) c check if y zero, when return with result zero. y1 = y(1) if (y1.eq.0) return c can also return if exponent not positive, then result x. y2 = y(two) if (y2.le.0) return c check for zero fractional part if (y2.lt.t) go to 10 c here the fractional part of x is zero. y(1) = 0 return c here 0 .lt. y2 .lt. t. clear integer part. 10 do 20 i = 1, y2 20 y(i+2) = 0 c normalize result and return. call mpnzr (y1, y2, y(1), y(two+1), 0) return end subroutine mpcmi (x, iz) c converts multiple-precision x to integer iz, c assuming that x not too large (else use mpcmim). c x is truncated towards zero, regardless of the value of rndrl. c iz is returned as zero if abs(int(x)) .gt. mxint . c the user may check for this by testing if c ((x(1).ne.0).and.(x(2).gt.0).and.(iz.eq.0)) c is true on return from mpcmi. common /mpcom/ b, t, dummy integer i, ib, mxint, xs, x2, \$ b, dummy(21), iz, mpget, mpparn, t, x(1) xs = x(1) iz = 0 c return zero if x zero or exponent(x) .le. 0. if (xs.eq.0) return x2 = mpget (x, 2) if (x2.le.0) return mxint = mpparn (16) ib = mxint/b c loop to convert integer part of x to single-precision integer. do 10 i = 1, x2 c check if b*iz would exceed mxint if (iz .gt. ib) go to 20 iz = b*iz if (i .gt. t) go to 10 c check if x(i+2) + iz would exceed mxint if (x(i+2) .gt. (mxint - iz)) go to 20 iz = iz + x(i+2) 10 continue c restore sign and return iz = xs*iz return c here abs(int(x)) is larger than mxint, so return zero. 20 iz = 0 return end subroutine mpcmim (x, y) c returns y = integer part of x (truncated towards 0), for mp x and y. c x may be packed or unpacked, y is unpacked. c use if y too large to be representable as a single-precision integer c (else could use mpcmi). rndrl is irrelevant. integer i, il, mpget, mpparn, t, x(1), y(1) c check legality of b, t etc. call mpchk call mpunpk (x, y) if (y(1).eq.0) return il = mpget (y, 2) + 1 t = mpparn(2) c if exponent large enough return y = x if (il.gt.t) return c if exponent small enough return zero if (il.gt.1) go to 10 y(1) = 0 return c set fraction to zero 10 do 20 i = il, t 20 y(i+2) = 0 return end integer function mpcmp (x, y) c compares the unpacked multiple-precision numbers x and y, c returning +1 if x .gt. y, c -1 if x .lt. y, c and 0 if x .eq. y. integer x(1), y(1), t2, i, mpparn c compare signs of x and y. if (x(1) - y(1)) 10, 30, 20 c x .lt. y 10 mpcmp = -1 return c x .gt. y 20 mpcmp = 1 return c sign(x) = sign(y), see if zero 30 if (x(1).eq.0) go to 55 c have to compare exponents and fractions t2 = mpparn(2) do 50 i = 2, t2 if (x(i) - y(i)) 60, 50, 70 50 continue c numbers equal 55 mpcmp = 0 return c abs(x) .lt. abs(y) and signs equal. 60 mpcmp = -x(1) return c abs(x) .gt. abs(y) and signs equal. 70 mpcmp = x(1) return end integer function mpcmpa (x, y) c compares abs(x) with abs(y) for mp x and y, c returning +1 if abs(x) .gt. abs(y), c -1 if abs(x) .lt. abs(y), c and 0 if abs(x) .eq. abs(y) c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1), xs, ys xs = x(1) ys = y(1) x(1) = iabs(xs) y(1) = iabs(ys) mpcmpa = mpcomp (x, y) x(1) = xs y(1) = ys return end integer function mpcmpd (x, di) c compares mp number x with double-precision number di, returning c +1 if x .gt. di, c 0 if x .eq. di, c -1 if x .lt. di. c x may be packed or unpacked. c comments regarding rounding error in subroutine mpcdm c are relevant. c x may be packed or unpacked. common r integer i2, mpcomp, r(1), sv, x(1) double precision di call mpsavn (sv) c allocate temporary space call mpnew (i2) c convert di to multiple-precision and compare call mpcdm (di, r(i2)) mpcmpd = mpcomp (x, r(i2)) call mpresn (sv) return end integer function mpcmpi (x, i) c compares mp number x with integer i, returning c +1 if x .gt. i, c 0 if x .eq. i, c -1 if x .lt. i c x may be packed or unpacked. common r integer i, i2, mpcomp, r(1), sv, x(1) call mpsavn (sv) call mpnew (i2) c convert i to multiple-precision and compare call mpcim (i, r(i2)) mpcmpi = mpcomp (x, r(i2)) call mpresn (sv) return end integer function mpcmpq (x, i, j) c returns +1 if x .gt. i/j, c 0 if x .eq. i/j, c -1 if x .lt. i/j, c for (packed or unpacked) mp x and integer i, j (j nonzero). c result is exact, so rndrl is irrelevant. common r common /mpcom/ b, t, m, dummy integer i2, k, l, sv, tg integer b, dummy(20), i, j, m, r(1), t, x(1) integer mpcmpi, mpgd call mpsavn (sv) k = i l = j c check for zero denominator. if (l) 20, 10, 30 10 call mperrm (24hj = 0 in call to mpcmpq\$) c here j was negative. 20 l = -l k = -k c now i/j = k/l and l positive 30 call mpgcd (k, l) c check for l = 1 if (l.ne.1) go to 40 c here l = 1 so can use mpcmpi. mpcmpq = mpcmpi (x, k) go to 50 c here l .gt. 1, so increase t so l*x can be formed exactly. 40 t = t + mpgd (l) tg = t c increase m so overflow can not occur m = m + t c use truncated arithmetic (actually exact). call mpsetr (0) c move x to temporary storage and multiply by l call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) call mpmuli (r(i2), l, r(i2)) c now compare l*x with k using mpmuli mpcmpq = mpcmpi (r(i2), k) c restore everything and return. 50 call mpresn (sv) return end integer function mpcmpr (x, ri) c compares mp number x with real number ri, returning c +1 if x .gt. ri, c 0 if x .eq. ri, c -1 if x .lt. ri. c x may be packed or unpacked. c comments regarding rounding error in subroutine mpcrm c are relevant. common r integer i2, mpcomp, r(1), sv, x(1) real ri call mpsavn (sv) c allocate temporary space call mpnew (i2) c convert ri to multiple-precision and compare call mpcrm (ri, r(i2)) mpcmpr = mpcomp (x, r(i2)) call mpresn (sv) return end subroutine mpcmr (x, rz) c converts multiple-precision x to single-precision rz. c assumes x in allowable range. common /mpcom/ b, t, dummy integer b, dummy(21), i, ie, t, tm, two, x(1) real float, rb, rz, rz2 two = 2 c check legality of parameters in common /mpcom/ call mpchk rz = 0e0 c return with rz = 0.0 if x is zero. if (x(1).eq.0) return rb = float(b) c loop to compute rz. do 10 i = 1, t rz = rb*rz + float(x(i+2)) tm = i c check if full single-precision accuracy attained rz2 = rz + 1e0 if (rz2.le.rz) go to 20 10 continue c now allow for exponent 20 ie = x(two) - tm 30 if (ie. le. 0) go to 40 ie = ie - 1 rz = rz*rb go to 30 40 if (ie .eq. 0) go to 50 ie = ie + 1 rz = rz/rb go to 40 c check reasonableness of result c following message indicates that x is too large or small - c try using mpcmre instead. 50 if (rz.le.0e0) call mperrm ( \$ 39hfloating-point over/underflow in mpcmr\$) c allow for sign of x and return. if (x(1).lt.0) rz = -rz return end subroutine mpcmre (x, n, rx) c returns integer n and single-precision rx such that mp c x = rx*outbas**n (approximately), where 1 .le. abs(rx) .lt. outbas c unless rx = 0. (outbas is in common /mpcom/ - default value c is 10.) c rounding options not implemented. common r integer i2, mpparn, n, r(1), sv, x(1) real float, rx call mpsavn (sv) if (x(1).ne.0) go to 10 n = 0 rx = 0e0 go to 20 c use truncated arithmetic. 10 call mpsetr (0) call mpnew (i2) call mpcmef (x, n, r(i2)) call mpcmr (r(i2), rx) if (abs(rx) .lt. float(mpparn(20))) go to 20 c here rx was rounded up to outbas n = n + 1 rx = float(r(i2)) 20 call mpresn (sv) return end subroutine mpcmul (xr, xi, yr, yi, zr, zi) c sets z = x*y where x = xr + i*xi, etc. are complex c multiple-precision numbers. c uses truncated arithmetic, no guard digits. c rounding options not implemented. common r integer i2, i3, i4, sv integer r(1), xi(1), xr(1), yi(1), yr(1), zi(1), zr(1) call mpsavn (sv) call mpsetr (0) call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpmul (xr, yr, r(i2)) call mpmul (xi, yi, r(i3)) call mpmul (xr, yi, r(i4)) call mpmul (xi, yr, zi) call mpadd (r(i4), zi, zi) call mpsub (r(i2), r(i3), zr) call mpresn (sv) return end integer function mpcomp (x, y) c compares the multiple-precision numbers x and y, c returning +1 if x .gt. y, c -1 if x .lt. y, c and 0 if x .eq. y. c x and y may be packed or unpacked. common r integer r(1), x(1), y(1), sv, mpcmp, i2, i3 c for efficiency treat unpacked numbers separately if (max0 (iabs(x(1)), iabs(y(1))) .gt. 1) go to 10 mpcomp = mpcmp (x, y) return c here x and/or y is packed 10 call mpsavn (sv) call mpnew (i2) call mpnew (i3) call mpunpk (x, r(i2)) call mpunpk (y, r(i3)) mpcomp = mpcmp (r(i2), r(i3)) call mpresn (sv) return end subroutine mpcos (x, y) c returns y = cos(x) for mp x and y, using mpcis. c time = o(sqrt(t)m(t)). c x may be packed or unpacked, y is unpacked. c rounding options are implemented as for mpcis. integer x(1), y(1), s(1) call mpcis (x, y, s, .false.) return end subroutine mpcosh (x, y) c returns y = cosh(x) for mp numbers x and y, x not too large. c x may be packed or unpacked, y is unpacked. c rounding options not yet implemented, no guard digits used. common r integer i2, r(1), sv, x(1), y(1) c check for x zero if (x(1).ne.0) go to 10 c cosh(0) = 1 call mpcim (1, y) return c save t etc. 10 call mpsavn (sv) c allocate temporary storage. call mpnew (i2) c use truncated arithmetic, move abs(x). call mpsetr (0) call mpabs (x, r(i2)) c if abs(x) too large mpexp will print error message c increase m to avoid overflow when cosh(x) representable call mpparc (3, r(sv+1)+2) call mpexp (r(i2), r(i2)) call mprec (r(i2), y) call mpadd (r(i2), y, y) c restore m etc. if result overflows or underflows, mpdivi will c act accordingly. call mpresn (sv) call mpdivi (y, 2, y) return end subroutine mpcqm (i, j, q) c converts the rational number i/j to multiple precision q. c rounding is defined by the parameter rndrl in common /mpcom/ - c see subroutine mpnzr. integer i, i1, j, j1, q(1) i1 = i j1 = j c remove any common factor of i and j. call mpgcd (i1, j1) c check sign of denominator. if (j1) 20, 10, 30 c division by zero 10 call mperrm (23hj = 0 in call to mpcqm\$) return c make denominator positive to give correct directed rounding. 20 i1 = -i1 j1 = -j1 30 call mpcim (i1, q) if (j1.ne.1) call mpdivi (q, j1, q) return end subroutine mpcrm (rx, z) c converts single-precision number rx to multiple-precision z. c some numbers will not convert exactly on machines c with base other than two or if b is not a power of two with c b**(t-1) sufficiently large. thus mpcrm should be used only to c obtain starting approximations etc. for accurate initialization c of mp numbers use mpcim, mpcqm, mpqpwr, or mpin. c warning - the parameter rndrl in common /mpcom/ has no effect. common /mpcom/ b, t, dummy integer b, dummy(21), i, ie, isv, j, re, rs, sv, t, \$ three, z(1) real float, rb, rj, rx call mpsavn (sv) call mpsetr (0) three = 3 c check sign if (rx) 20, 10, 30 c if rx = 0e0 return 0 10 z(1) = 0 go to 100 c rx .lt. 0e0 20 rs = -1 rj = -rx go to 40 c rx .gt. 0e0 30 rs = 1 rj = rx c we want to reduce rj to range 1/16 .le. rj .lt. 1. 40 ie = 0 50 if (rj.lt.1e0) go to 60 c increase ie and divide rj by 16. ie = ie + 1 rj = 0.0625e0*rj go to 50 60 if (rj.ge.0.0625e0) go to 70 ie = ie - 1 rj = 16e0*rj go to 60 c now rj is dy divided by suitable power of 16. c set exponent to 0 70 re = 0 rb = float(b) c conversion loop (assume single-precision ops. exact) do 80 i = 1, t rj = rb*rj isv = i j = int(rj) c check for (unlikely) effect of strange floating-point arithmetic. if (j.ge.b) go to 110 if (rj.lt.0.0) go to 120 z(i+2) = j 80 rj = rj - float(j) c normalize result 90 call mpnzr (rs, re, z(1), z(three), 0) c now multiply by 16**ie call mpscal (z, 16, ie) c can return now 100 call mpresn (sv) return c here a digit is .ge. b, which should never happen 110 j = b-1 go to 130 c here a digit is negative, which should never happen 120 j = 0 c set remaining digits to j (either b-1 or 0) 130 do 140 i = isv, t 140 z(i+2) = j go to 90 end subroutine mpdaw (x, y) c returns y = dawsons integral (x) c = exp(-x**2)*(integral from 0 to x of exp(u**2)du), c for packed/unpacked multiple-precision x, unpacked mp y. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, dummy logical err integer i, i2, i3, i4, sv, tg, xs integer b, dummy(20), m, mptlb, r(1), t, x(1), y(1) c save t etc. and check for x zero. call mpsavn (sv) xs = x(1) if (xs.ne.0) go to 10 c daw(0) = 0 y(1) = 0 go to 50 c increase t and allocate temporary storage. 10 call mpgd3 (mptlb (iabs(t)), tg) call mpnew (i2) c use truncated arithmetic, increase m to avoid underflow. call mpsetr (0) m = m + t c move abs(x) to temporary storage. call mpmove (x, iabs(r(sv)), r(i2), tg) r(i2) = iabs(r(i2)) c try asymptotic series call mperf3 (r(i2), r(i2), .true., err) if (err) go to 20 c here asymptotic series was successful. fix up sign. r(i2) = xs*r(i2) go to 40 c asymptotic series. not accurate enough so use power series 20 call mpnew (i3) call mpnew (i4) call mpmove (x, iabs(r(sv)), r(i4), tg) t = tg call mpmul (r(i4), r(i4), r(i4)) call mpneg (r(i4), r(i4)) call mpexp (r(i4), r(i4)) call mpmove (x, iabs(r(sv)), r(i2), tg) call mpmul (r(i2), r(i4), r(i3)) call mpmul (r(i2), r(i2), r(i4)) call mpstr (r(i3), r(i2)) i = 0 c power series loop, reduce t if possible 30 t = tg + 2 + r(i3+1) - r(i2+1) if (t.le.2) go to 40 t = min0 (t, tg) i = i + 1 call mpmul (r(i4), r(i3), r(i3)) call mpmuls (r(i3), 2*i-1, 1, i, 2*i+1) c restore t for addition t = tg call mpadd (r(i2), r(i3), r(i2)) if (r(i3).ne.0) go to 30 c round result 40 call mpres2 (sv) call mprnd (r(i2), tg, y, iabs(r(sv)), 1) c restore everything and return. 50 call mpresn (sv) return end integer function mpdga (x, n) c returns the n-th digit of the mp number x for 1 .le. n .le. t. c returns zero if x is zero or n .le. 0 or n .gt. t. integer mpparn, n, x(1) mpdga = 0 if ((x(1).ne.0).and.(n.gt.0).and.(n.le.mpparn(2))) mpdga = x(n+2) return end subroutine mpdgb (i, x, n) c sets the n-th digit of the mp number x to i. c n must be in the range 1 .le. n .le t, c i must be in the range 0 .le. i .lt. b c (and i .ne. 0 if n .eq. 1). c the sign and exponent of x are unchanged. integer i, mpparn, n, x(1) if ((n.le.0).or.(n.gt.mpparn(2))) call mperrm ( \$ 40hdigit position illegal in call to mpdgb\$) if ((i.lt.0).or.(i.ge.mpparn(1)).or.((i+n).le.1)) call \$ mperrm (37hdigit value illegal in call to mpdgb\$) x(n+2) = i return end integer function mpdiga (x) c returns the number of mp digits (second word in common /mpcom/). c x is a dummy mp argument. integer mpparn, x(1) mpdiga = mpparn (2) return end subroutine mpdigb (i, x) c sets the number of mp digits (second word of common /mpcom/) to i. c i should be an integer such that i .ge. 2 c and i+2 .le. mt2 (where mt2 is in common /mpcom/). c x is a dummy mp argument (augment expects one). integer i, x(1) call mpparc (2, i) call mpchk return end integer function mpdigs (n) c returns the number of mp digits (to base b) required for the c equivalent of at least n floating (base outbas) places, c b**(mpdigs-1) .ge. outbas**(n-1) . c if outbas has its default value of 10, result is the number c of base b digits to give the equivalent of at least n floating c decimal places. integer mpchgb, mpparn, n mpdigs = max0 (2, mpchgb (mpparn(1), mpparn(20), n-1) + 1) return end integer function mpdigv (x) c returns 0, ... , 9, 10, ... , 15 if x is the character c 0, ... , 9, a, ... , f, c and the value returned is less than inbase (default value ten), c returns -1 otherwise. integer i, inbase, j, mpdigw, mpparn, x logical mpis inbase = mpparn (19) if ((inbase .lt. 2) .or. (inbase .gt. 16)) call mpchk do 10 i = 1, inbase j = i - 1 if (mpis (x, mpdigw(j))) go to 20 10 continue mpdigv = -1 return c here digit found 20 mpdigv = j return end integer function mpdigw (n) c returns character 0, ... , 9, a, ... , f if n is c 0, ... , 9, 10, ... , 15 respectively, c returns character * otherwise. integer d(17), i, n data d(1), d(2), d(3), d(4) /1h0, 1h1, 1h2, 1h3/ data d(5), d(6), d(7), d(8) /1h4, 1h5, 1h6, 1h7/ data d(9), d(10), d(11), d(12) /1h8, 1h9, 1ha, 1hb/ data d(13), d(14), d(15), d(16) /1hc, 1hd, 1he, 1hf/ data d(17) /1h*/ i = min0 (n, 16) if (i .lt. 0) i = 16 mpdigw = d(i+1) return end subroutine mpdim (x, y, z) c sets z = x - min (x, y) = max (0, x-y) for mp x and y, c rounding as for mpsub. x and/or y may be packed or unpacked. integer x(1), y(1), z(1) c avoid possibility of overflow if result is zero. if ((x(1).le.0).and.(y(1).ge.0)) go to 10 call mpksub (x, y, z) if (z(1).ge.0) return 10 z(1) = 0 return end subroutine mpdiv (x, y, z) c c sets z = x/y, for mp x, y and z. c uses mpdivl for small t, mprec and mpmul for large t, so time c is o(m(t)). over/underflow is detected by subroutine mpnzr. c c rounding is determined by rndrl in common /mpcom/ as follows - c rndrl = 0 - error less than 1 unit in last place (ulp), so exact c if result can be exactly represented. c rndrl = 1 - see mpnzr. c rndrl = 2 or 3 - directed roundings - see subroutine mpnzr. c common r common /mpcom/ b, t, m, dummy integer cross, i2, i3, sv, tg integer b, dummy(20), m, mpgd, r(1), t, x(1), y(1), z(1) c following crossover point determined empirically data cross /80/ c save t, m, rndrl etc. call mpsavn (sv) c check for division by zero if (y(1) .eq. 0) call mperrm ( \$ 44hattempted division by zero in call to mpdiv\$) c check for x = 0 if (x(1).ne.0) go to 10 z(1) = 0 go to 30 c see if t small enough that mpdivl is faster than mprec c or must use mpdivl to ensure correct rounded result 10 tg = t + 1 + mpgd(100) if ((tg.lt.cross) .or. (r(sv+2).ne.0)) go to 20 c here use mprec and mpmul. increase t and m temporarily. t = tg m = m + t c allocate temporary space call mpnew (i2) c move y and compute reciprocal, taking care for directed rounding call mpmove (y, iabs(r(sv)), r(i2), tg) call mprevr (-x(1)) call mprec (r(i2), r(i2)) c restore rndrl and move x for multiplication call mpsetr (iabs(r(sv+2))) call mpnew (i3) call mpmove (x, iabs(r(sv)), r(i3), tg) call mpmul (r(i3), r(i2), r(i2)) c restore m (so mpnzr can detect overflow/underflow) m = r(sv+1) c round result (to nearest if rndrl = 0 or 1) call mprnd (r(i2), tg, z, iabs(r(sv)), 0) go to 30 c here faster or necessary to use mpdivl 20 call mpdivl (x, y, z) c restore everything 30 call mpresn (sv) return end subroutine mpdiv2 (x, j, kg, re, a) c called by mpdivi, not for independent use. c assumes x(1) ... x(t) represents a base b number, x(1) .gt. 0, c j .gt. 0. divides by j, result in a(1) ... a(t+kg). the result c is left shifted so a(1) .gt. 0 and re decremented accordingly. c assumes kg. ge. 0, also kg .gt. 0 if rndrl .gt. 0 . common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, mxexpn, \$ mnexpn, rndrl, ktunfl, mxunfl, decpl, mt2, mxint, dummy integer b2, c, c2, i, iq, iqj, ir, j1, j11, j2, k, kh, khh, r1, \$ a(1), b, decpl, dummy(7), j, kg, ktunfl, lun, m, mnexpn, \$ mnsptr, mt2, mxexpn, mxint, mxr, mxsptr, mxunfl, re, rndrl, \$ sptr, t, x(1) c = 0 i = 0 khh = t + kg c if j*b not representable as an integer have to simulate c long division. b2 = mxint/b if (j.ge.b2) go to 70 c look for first nonzero digit in quotient 10 i = i + 1 c = b*c if (i.le.t) c = c + x(i) r1 = c/j if (r1) 150, 10, 20 c adjust exponent and get t+kg digits in quotient 20 re = re + 1 - i a(1) = r1 c = b*(c - j*r1) kh = 2 if (i.ge.t) go to 40 kh = t + 1 - i do 30 k = 2, kh i = i + 1 c = c + x(i) a(k) = c/j 30 c = b*(c - j*a(k)) if (c.lt.0) go to 150 kh = kh + 1 40 if (kh.gt.khh) go to 60 do 50 k = kh, khh a(k) = c/j 50 c = b*(c - j*a(k)) if (c.lt.0) go to 150 c adjust last digit for directed rounding 60 if ((rndrl.gt.1).and.(c.gt.0)) a(khh) = 1 return c here need simulated double-precision division 70 c2 = 0 j1 = j/b j2 = j - j1*b j11 = j1 + 1 c look for first nonzero digit 80 i = i + 1 c = b*c + c2 c2 = 0 if (i.le.t) c2 = x(i) if (c-j1) 80, 90, 100 90 if (c2.lt.j2) go to 80 c compute t+kg quotient digits 100 re = re + 1 - i k = 1 go to 120 c main loop for large abs(iy) case 110 k = k + 1 if (k.gt.khh) go to 60 i = i + 1 c get approximate quotient first 120 ir = c/j11 c now reduce so overflow does not occur iq = c - ir*j1 if (iq.lt.b2) go to 130 c here iq*b would possibly overflow so increase ir ir = ir + 1 iq = iq - j1 130 iq = iq*b - ir*j2 if (iq.ge.0) go to 140 c here iq negative so ir was too large ir = ir - 1 iq = iq + j 140 if (i.le.t) iq = iq + x(i) iqj = iq/j c a(k) = quotient, c = remainder a(k) = iqj + ir c = iq - j*iqj if (c.ge.0) go to 110 c carry negative so overflow must have occurred 150 call mperrm (40hinteger overflow in mpdiv2, b too large\$) return end subroutine mpdiv3 (x, y, d, kg) c called by mpdivl. assumes x(2), ... x(t+1) and c y(2), ... , y(t+1) contain base b integers, x(2) .gt. 0, c y(2) .gt. 0, kg .ge. 0. also assumes kg .ge. 2 if rndrl .gt. 0. c returns first t+kg digits of quotient x/y in d(1), ... , d(t+kg). c destroys x(1), ... , x(t+1), uses but restores y(1). c arguments must not overlap. common /mpcom/ b, t, dummy integer b1, c, i, ir, j, k, q, tg, two, t1, t2, y1s, \$ b, d(1), dummy(21), kg, mpparn, t, x(1), y(1) two = 2 c save y(1) and set to zero. y1s = y(1) y(1) = 0 c set up outer loop. x(1) = 0 b1 = b*(b-1) t1 = t + 1 t2 = t1 + 1 tg = t + kg c outer loop do 70 j = 1, tg d(j) = 0 c compare x and y (usually only once for each j). 10 do 20 i = 1, t1 if (x(i) - y(i)) 50, 20, 30 20 continue c here b*y .gt. x .ge. y, compute lower bound on quotient x/y. 30 q = max0 (1, (b*x(1) + x(two))/(y(two) + 1)) c compute sharper lower bound if overflow impossible. if ((x(1)+1) .lt. (mpparn(16)/(b*b))) q = \$ max0 (q, (b*(b*x(1)+x(two))+x(two+1))/ \$ (b*y(two)+y(two+1)+1)) c increment j-th quotient digit d(j) = d(j) + q i = t2 c k is (b + signed carry), 0 .lt. k .le. b. k = b c inner loop do 40 ir = 1, t1 i = i - 1 c = b1 + k + x(i) - q*y(i) k = c/b 40 x(i) = c - b*k c there should be no carry off end. check just in case. if (k.eq.b) go to 10 call mperrm (36hinteger overflow occurred in mpdiv3\$) c shift x left 50 do 60 i = 1, t 60 x(i) = x(i+1) c last statement of outer loop. 70 x(t1) = 0 c restore y(1). y(1) = y1s c return if rndrl .le. 1 if (mpparn(11).le.1) return c if remainder nonzero make sure a guard digit nonzero for c directed rounding (assumes kg .ge. 2 here) do 80 i = 1, t 80 d(tg) = max0 (d(tg), x(i)) return end subroutine mpdivi (x, iy, z) c divides mp x by the single-precision integer iy giving mp z. c this is much faster than division by an mp number. c rounding is defined by parameter rndrl in common /mpcom/ - c see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer i2, j, kg, mpgd, re, rs, two, \$ b, dummy(12), iy, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, r(1), rndrl, sptr, t, x(1), z(1) two = 2 c compute number of guard digits required. kg = min0 (1, rndrl) if (rndrl.eq.1) kg = 1 + mpgd (iy) c allocate t+kg words for accumulator if necessary. i2 = sptr if (kg.gt.0) call mpnew2 (i2, t + kg) c get sign of x rs = x(1) j = iabs(iy) c check sign of divisor if (iy) 20, 10, 30 10 call mperrm (45hattempted division by zero in call to mpdivi\$) c allow for negative divisor 20 rs = -rs c check for zero dividend 30 if (rs.ne.0) go to 40 z(1) = 0 go to 100 40 re = x(two) c check for division by +-b (this is common - see mpout2) if (j.ne.b) go to 70 re = re - 1 c move x to z. 50 call mpstr (x, z) c fix up sign and exponent of result. 60 z(1) = rs z(two) = re c check for underflow if (re.le.(-m)) go to 90 c update mnexpn mnexpn = min0 (re, mnexpn) go to 100 c check for division by +-1 70 if (j.eq.1) go to 50 c genuine division here. check number of guard digits. if (kg.gt.0) go to 80 c here can use z(3), ... , z(t+2) as accumulator. call mpdiv2 (x(two+1), j, 0, re, z(two+1)) go to 60 c here use r(i2), ... , r(sptr-1) as accumulator 80 call mpdiv2 (x(two+1), j, kg, re, r(i2)) call mpnzr (rs, re, z(1), r(i2), kg) go to 100 c underflow here 90 call mpunfl (z) c restore stack pointer and return 100 sptr = i2 return end subroutine mpdivl (x, y, z) c divides x by y, giving result z, for mp x, y and z. c uses long division method, time o(t**2). an alternative method c (implemented in mpdiv) with time o(m(t)) is to use mprec and c mpmul. with the present implementation of mpmul, this is faster c (by up to a factor of about 2) for large t, but mpdivl is faster c for small t (see data cross /.../ in mpdiv and mprec). c rounding is determined by rndrl in common /mpcom/ - c see subroutine mpnzr. common r integer i2, i3, kg, r(1), sv, two, x(1), y(1), z(1) two = 2 c check for division by zero. if (y(1) .eq. 0) call mperrm ( \$ 35hdivision by zero in call to mpdivl\$) c check for x zero. if (x(1).ne.0) go to 10 z(1) = 0 return c here x and y are nonzero. determine number of guard digits. 10 call mpsavn (sv) kg = min0 (2, r(sv+2)+1) if (r(sv+2).eq.1) kg = r(sv) + 2 c allocate temporary storage call mpnew2 (i2, r(sv) + kg) call mpnew (i3) c move x to temporary storage as mpdiv3 overwrites first argument. call mpstr (x, r(i3)) c do long division call mpdiv3 (r(i3+1), y(two), r(i2), kg) c normalize and round result call mpnzr (x(1)*y(1), x(two)+1-y(two), z, r(i2), kg) c restore stack pointer and return. call mpresn (sv) return end subroutine mpdump (x) c dumps out the mp number x (sign, exponent, fraction digits), c useful for debugging purposes. c embedded blanks should be interpreted as zeros. (they could be c avoided by using j instead of i format, but this is nonstandard.) c x may be packed or unpacked. common r integer b, i2, lun, mpparn, r(1), sv, t2, x(1) logical err c check legality of parameters in common /mpcom/ call mpchk lun = mpparn (4) if (x(1).ne.0) go to 10 c if x = 0 just write sign as remainder undefined call mpio (x, 1, lun, 7h(1x,i2), err) return c here x is nonzero, so unpack if necessary 10 call mpsavn (sv) b = mpparn (1) t2 = mpparn (2) + 2 call mpnew (i2) call mpunpk (x, r(i2)) if (b.le.10) call mpio (r(i2), t2, lun, \$ 30h(1x,i2,i12,4x,50i1/(19x,50i1)), err) if ((b.gt.10).and.(b.le.100)) call mpio (r(i2), t2, lun, \$ 30h(1x,i2,i12,4x,25i2/(19x,25i2)), err) if ((b.gt.100).and.(b.le.1000)) call mpio (r(i2), t2, lun, \$ 30h(1x,i2,i12,4x,17i3/(19x,17i3)), err) if ((b.gt.1000).and.(b.le.10000)) call mpio (r(i2), t2, lun, \$ 30h(1x,i2,i15,4x,12i4/(22x,12i4)), err) if (b.gt.10000) call mpio (r(i2), t2, lun, \$ 30h(1x,i2,i23,4x,4i10/(30x,4i10)), err) if (err) call mperrm (33herror return from mpio in mpdump\$) call mpresn (sv) return end subroutine mpei (x, y) c returns y = ei(x) = -e1(-x) c = (principal value integral from -infinity to x of c exp(u)/u du), c for mp numbers x and y, c using the power series for small abs(x), the asymptotic series for c large abs(x), and the continued fraction for intermediate negative c x. relative error in y is small except if x is very close to the c zero 0.37250741078136663446... of ei(x), c and then the absolute error in y is o(b**(1-t)). c in any case the error in y could be induced by an c o(b**(1-t)) relative perturbation in x. c time is o(t.m(t)). c rounding options not yet implemented, guard digits not always used. common r common /mpcom/ b, t, dummy integer i, i2, i3, i4, j, k, sv, td, tm, tm2, ts1, ts2, xs integer b, dummy(21), mpt, r(1), t, x(1), y(1) integer mpchgb, mpcmpa, mpcmpq, mpget, mptlb c save t etc. call mpsavn (sv) xs = x(1) c ei(0) is undefined. if (xs .eq. 0) call mperrm (23hx zero in call to mpei\$) c prepare to increase t. tm2 = (11*t+19)/10 tm = (6*t+9)/5 i = 0 c increase t and allocate temporary space. t = tm call mpnew (i2) call mpnew (i3) call mpnew (i4) c use truncated arithmetic call mpsetr (0) c move x and restore t. call mpmove (x, iabs(r(sv)), r(i2), tm) call mpabs (r(i2), r(i3)) t = r(sv) mpt = mptlb(iabs(t)) c see if abs(x) large enough to use asymptotic series c 7/10 .gt. ln(2) if (mpcmpq (r(i3), 7*mpt, 10) .gt. 0) go to 40 c see if x negative and continued fraction usable c the constant 1/15 was determined empirically and may be c decreased (but not increased) if desired. if ((xs.lt.0) .and. (mpcmpq (r(i3), mpt, 15) .gt. 0)) go to 60 c use power series here, but need to increase t if x negative c to compensate for cancellation. t = t + 1 ts1 = t ts2 = t if (xs.gt.0) go to 10 c if x negative result about b**(-td) and terms about b**td so c need up to 2*td extra digits to compensate for cancellation c mpchgb(...)/3 underestimates ln(b). call mpmulq (r(i3), 3, mpchgb(2,max0(2,b/2),2), r(i4)) call mpcmi (r(i4), td) td = td + 1 ts2 = t + td ts1 = min0 (ts2 + td, tm) ts2 = min0 (ts2, tm2) c use ts2 digits for ln and eulers constant computation t = ts2 c now prepare to sum power series 10 call mpln (r(i3), r(i4)) c mpei could be speeded up if eulers constant were c precomputed and saved call mpeul (r(i2)) call mpadd (r(i2), r(i4), r(i2)) c now use ts1 digits for summing power series t = ts1 c restore sign of r(i3) r(i3) = xs call mpadd (r(i2), r(i3), r(i2)) call mpstr (r(i3), r(i4)) c loop to sum power series, reducing t if possible 20 if (xs.ge.0) t = ts1 + 2 + r(i4+1) - r(i2+1) if ((xs.lt.0).and.(r(i4+1).le.0)) t = ts2 + 2 + r(i4+1) t = min0 (t, ts1) if (t.le.2) go to 30 call mpmul (r(i3), r(i4), r(i4)) i = i + 1 call mpmuls (r(i4), i, 1, i+1, i+1) c restore t for addition t = ts1 call mpadd (r(i2), r(i4), r(i2)) if (r(i4).ne.0) go to 20 c restore t, move result and return 30 t = r(sv) call mpstr (r(i2), y) go to 100 c here we can use asymptotic series, and no need to increase t 40 call mprec (x, r(i3)) c mpexp gives error message if x too large here call mpexp (x, y) if (y(1).eq.0) go to 100 call mpmul (y, r(i3), y) call mpstr (y, r(i2)) c loop to sum asymptotic series, reducing t if possible 50 t = r(sv) + 2 + r(i2+1) - mpget (y, 2) c return if terms small enough to be negligible if (t.le.2) go to 100 t = min0 (t, r(sv)) call mpstr (r(i2), r(i4)) call mpmul (r(i2), r(i3), r(i2)) i = i + 1 call mpmuli (r(i2), i, r(i2)) c return if terms increasing if (mpcmpa (r(i2), r(i4)) .ge. 0) go to 100 c restore t for addition t = r(sv) call mpadd (y, r(i2), y) if (r(i2).ne.0) go to 50 go to 100 c here 0.1*t*ln(b) .lt. -x .le t*ln(b), so use continued fraction. 60 k = t j = 0 c use equivalent of 6 decimal places for forward recurrence. c we could use single-precision real here if trusted. t = min0 (t, max0 (2, mpchgb (iabs(b), 10, 5) + 1)) call mpneg (x, r(i2)) call mpcim (1, r(i3)) c use forward recurrence to find out how many terms are c needed in backward recurrence 70 j = j + 1 call mpdivi (r(i2), j, r(i4)) call mpadd (r(i3), r(i4), r(i3)) call mpmul (x, r(i3), r(i4)) call mpsub (r(i2), r(i4), r(i2)) if (r(i3) .eq. 0) go to 70 c scaling here 80 if (r(i3+1) .le. 1) go to 70 r(i3+1) = r(i3+1) - 1 r(i2+1) = r(i2+1) - 1 k = k - 2 if (k .gt. 0) go to 80 c restore t for backward recurrence t = r(sv) call mpabs (x, r(i3)) c now use backward recurrence with mp arithmetic call mpcim (1, r(i2)) 90 call mpdivi (r(i3), j, r(i4)) call mpadd (r(i2), r(i4), r(i2)) call mpmul (x, r(i2), r(i4)) call mpsub (r(i3), r(i4), r(i3)) c scale to avoid overflow r(i2+1) = r(i2+1) - r(i3+1) r(i3+1) = 0 j = j - 1 if (j.gt.0) go to 90 call mpdiv (r(i2), r(i3), r(i2)) call mpexp (x, y) call mpmul (y, r(i2), y) y(1) = -y(1) c restore everything and return. 100 call mpresn (sv) return end subroutine mpeps (x) c c sets mp x to the (multiple-precision) machine precision, c that is c x = 1.01*(b**(1-t)) (rounded up) if rndrl = 0, c 0.5*(b**(1-t)) (rounded up) if rndrl = 1, c b**(1-t) if rndrl = 2 or 3. c c x is an upper bound on the smallest positive representable c number such that the relative error in the basic mp operations c (addition, subtraction, multiplication and division) is at c most x (unless the result underflows). c common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr, \$ rndrl, rsv, sptr, t, two, x(1) two = 2 c check b, t etc. call mpchk c set x to b**(1-t) call mpcim (1, x) x(two) = 2 - t call mpupdt (x(two)) if (rndrl.gt.1) return c here rndrl = 0 or 1. rsv = rndrl c set rndrl for rounding up rndrl = 3 if (rsv .eq. 0) call mpmulq (x, 101, 100, x) if (rsv .ne. 0) call mpdivi (x, 2, x) rndrl = rsv return end logical function mpeq (x, y) c returns logical value of (x .eq. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mpeq = (mpcomp(x,y) .eq. 0) return end subroutine mperf (x, y) c returns y = erf(x) = sqrt(4/pi)*(integral from 0 to x of c exp(-u**2) du) for mp x and y. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, dummy logical err integer irx, i2, i3, sv, tg, two, xs, \$ b, dummy(20), m, r(1), t, x(1), y(1), \$ mpchgb, mpcmpq, mpgd, mptlb two = 2 c save t etc. call mpsavn (sv) c check for zero argument. xs = x(1) if (xs.ne.0) go to 10 c erf(0) = 0 y(1) = 0 go to 60 c increase t, allocate temporary storage, and restore t. 10 call mpgd3 (mptlb(iabs(t)), tg) call mpnew (i2) call mpnew (i3) c use truncated arithmetic internally. call mpsetr (0) m = m + t call mpmove (x, iabs(r(sv)), r(i2), tg) r(i2) = iabs (r(i2)) call mpmul (r(i2), r(i2), r(i3)) c 7/10 .gt. ln(2) if (mpcmpq (r(i3), 7*mptlb(iabs(t)), 10) .lt. 0) go to 20 c here abs(x) so large that erf(x) = +-1 to full accuracy call mpcim (xs, r(i3)) go to 50 c try using asymptotic series. c can possibly reduce t temporarily. 20 call mpmulq (r(i3), 10, mpchgb(2,iabs(b),7), r(i3)) call mpcmi (r(i3), irx) t = min0(t, max0(4*mpgd(100), t-irx)) call mperf3 (r(i2), r(i2), .false., err) if (.not. err) go to 30 c asymptotic series insufficient, so use power series t = tg call mpmove (x, iabs(r(sv)), r(i2), tg) call mperf2 (r(i2), r(i2)) c in both cases multiply by sqrt(4/pi)*exp(-x**2) 30 call mpmove (x, iabs(r(sv)), r(i3), iabs(t)) call mpmul (r(i3), r(i3), r(i3)) r(i3) = -r(i3) call mpexp(r(i3), r(i3)) call mpmul (r(i3), r(i2), r(i2)) call mppi (r(i3)) call mproot (r(i3), -2, r(i3)) call mpmul (r(i3), r(i2), r(i2)) if (.not. err) go to 40 c used power series so can return call mpmuli (r(i2), 2, r(i3)) go to 50 c used power series. 40 call mpmuli (r(i2), -2, r(i2)) call mpmove (r(i2), iabs(t), r(i2), tg) t = tg call mpaddi (r(i2), 1, r(i3)) r(i3) = xs*r(i3) c round result. 50 call mpres2 (sv) call mprnd (r(i3), iabs(t), y, iabs(r(sv)), 1) c ensure that result in -1 to +1. if ((y(1).ne.0) .and. (y(two).gt.0)) call mpcim (xs, y) c restore everything and return 60 call mpresn (sv) return end subroutine mperf2 (x, y) c returns y = exp(x**2)*(integral from 0 to x of exp(-u*u) du) c for mp numbers x and y, using the power series for small x. c called by mperf, not recommended for independent use. c rounding options not implemented, uses no guard digits. common r common /mpcom/ b, t, dummy integer i, i2, i3, mpcmpq, mptlb, sv, two, xs, \$ b, dummy(21), r(1), t, x(1), y(1) two = 2 c save t etc., check for zero argument. call mpsavn (sv) if (x(1).ne.0) go to 10 c return 0 if x .eq. 0 y(1) = 0 go to 40 c allocate temporary space. 10 call mpnew (i2) c use truncation internally call mpsetr (0) call mpmul (x, x, r(i2)) c 7/10 .gt. ln(2) if (mpcmpq (r(i2), 7*mptlb(iabs(t)), 10) .gt. 0) go to 30 c use the power series here call mpstr (x, y) call mpmuli (r(i2), 2, r(i2)) call mpnew (i3) call mpstr (x, r(i3)) i = 1 c loop to sum series, reducing t if possible 20 t = r(sv) + 2 + r(i3+1) - y(two) if (t.le.2) go to 40 t = min0 (t, r(sv)) call mpmul (r(i2), r(i3), r(i3)) i = i + 2 call mpdivi (r(i3), i, r(i3)) c restore t for addition t = r(sv) call mpadd (y, r(i3), y) if (r(i3).ne.0) go to 20 go to 40 c here abs(x) large, so integral is +-sqrt(pi/4) c if abs(x) too large mpexp gives error message 30 call mpexp (r(i2), r(i2)) xs = x(1) call mppi (y) call mpsqrt (y, y) call mpdivi (y, 2*xs, y) call mpmul (y, r(i2), y) c restore everything and return 40 call mpresn (sv) return end subroutine mperf3 (x, y, ind, error) c if ind = .false., returns y = exp(x**2)*(integral from x to c infinity of exp(-u**2) du), c if ind = .true., returns y = exp(-x**2)*(integral from 0 to c x of exp(u**2) du), c in both cases using the asymptotic series. c x and y are mp numbers, ind and error are logical. c error is returned as .false. if x is large enough for c the asymptotic series to give full accuracy, c otherwise error is returned as .true. and y as zero. c the condition on x for error .eq. .false. is approximately that c x .gt. sqrt(t*log(b)). c called by mperf, mperfc and mpdaw, not recommended for c independent use. c rounding options not implemented, uses no guard digits. common r common /mpcom/ b, t, dummy integer i, ie, i2, i3, sv, \$ b, dummy(21), r(1), t, x(1), y(1), \$ mpcmpi, mpcmpq, mpget, mptlb logical ind, error c save t etc, use truncation internally. call mpsavn (sv) call mpsetr (0) c get some temporary space call mpnew (i2) error = .false. c check that can get at least t-2 digits accuracy if (mpcmpi (x, 1) .le. 0) go to 10 if (mpget (x, 2) .ge. t) go to 30 call mpmul (x, x, r(i2)) if (mpcmpq (r(i2), 7*mptlb(t-2), 10) .gt. 0) go to 30 c here x is too small for asymptotic series to give c full accuracy, so return with y = 0 and error = .true.. 10 y(1) = 0 error = .true. c restore t etc. and return. 20 call mpresn (sv) return c here worth trying asymptotic series. 30 call mprec (x, y) c allocate more temporary space call mpnew (i3) call mpmul (y, y, r(i2)) call mpdivi (r(i2), 2, r(i2)) if (.not. ind) r(i2) = -r(i2) call mpdivi (y, 2, y) call mpstr (y, r(i3)) i = 1 c loop to sum series, reducing t if possible 40 ie = r(i3+1) t = r(sv) + 2 + ie - mpget (y, 2) if (t.le.2) go to 20 t = min0 (t, r(sv)) call mpmul (r(i2), r(i3), r(i3)) call mpmuli (r(i3), i, r(i3)) i = i + 2 c restore t for addition t = r(sv) c check if terms are getting larger - if so x is too c small for asymptotic series to be accurate if (r(i3+1).gt.ie) go to 10 call mpadd (y, r(i3), y) if (r(i3).ne.0) go to 40 go to 20 end subroutine mperfc (x, y) c returns y = erfc(x) = 1 - erf(x) for mp numbers x and y, c using mperf and mperf3. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy logical err integer it, i2, i3, mpchgb, mptlb, sv, tg, \$ b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1) c save t etc., check sign of argument x. call mpsavn (sv) if (x(1).gt.0) go to 10 c for x .le. 0 no significant loss of accuracy in using erf(x) call mprevr (1) call mperf (x, y) call mpres2 (sv) y(1) = -y(1) call mpaddi (y, 1, y) go to 40 c increase t, allocate some space, use truncation internally. 10 call mpgd3 (mptlb (iabs(t)), tg) call mpnew (i2) call mpsetr (0) c try asymptotic series call mpmove (x, iabs(r(sv)), r(i2), tg) call mperf3 (r(i2), r(i2), .false., err) if (err) go to 20 c asymptotic series worked, so multiply by c sqrt(4/pi)*exp(-x**2) call mpnew (i3) call mpmove (x, iabs(r(sv)), r(i3), tg) call mpmul (r(i3), r(i3), r(i3)) r(i3) = -r(i3) call mpexp (r(i3), r(i3)) call mpmul (r(i3), r(i2), r(i2)) call mppi (r(i3)) call mproot (r(i3), -2, r(i3)) call mpmul (r(i3), r(i2), r(i2)) call mpmuli (r(i2), 2, r(i2)) go to 30 c here asymptotic series inaccurate so have to c use mperf, increasing precision to compensate for c cancellation. an alternative method (possibly faster) is c to use the continued fraction for exp(x**2)*erfc(x). 20 t = r(sv) call mpmul (x, x, r(i2)) c ln(b) .gt. mpchgb (2, b, 80) / 120 call mpmulq (r(i2), 120, mpchgb(2,iabs(b),80), r(i2)) call mpcmi (r(i2), it) c compute new t for mperf computation t = t + it sptr = i2 call mpgd3 (1, tg) call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) call mperf (r(i2), r(i2)) r(i2) = -r(i2) call mpaddi (r(i2), 1, r(i2)) c round result. 30 call mpres2 (sv) call mprnd (r(i2), tg, y, iabs(r(sv)), 1) c restore everything and return. 40 call mpresn (sv) return end subroutine mperr c this routine is called when a fatal error condition is encountered, c and after a message has been written on logical unit lun. common /mpcom/ some, lun, rest integer some(3), lun, rest(19) write (lun, 10) some, rest 10 format (42h *** execution terminated by call to mperr, \$ 25h in mp version 800207 ***/ \$ 14h *** b, t, m =, i22, 2i20/ \$ 24h *** mxr, sptr, mxsptr =, i12, 2i20/ \$ 29h *** mnsptr, mxexpn, mnexpn =, i7, 2i20/ \$ 28h *** rndrl, ktunfl, mxunfl =, i8, 2i20/ \$ 24h *** decpl, mt2, mxint =, i12, 2i20/ \$ 36h *** exwid, inrecl, inbase, outbas =, 4i4/ \$ 28h *** expch, chword, onescp =, 11x, a1, 2i4//) c ansi version uses stop, univac 1108 version uses c return 0 in order to give a trace-back. call istkpr call abort stop end subroutine mperrm (a) c the argument a is a hollerith string terminated by \$. c mperrm writes the string out on unit lun, then calls mperr. c only the first 71 characters of a are significant, and c are preceeded and followed by three asterisks. common r common /mpcom/ b, t, m, lun, dummy integer a(1), i, i2, j, len, r(1), sv integer b, t, m, lun, dummy (19) logical err call mpchk call mpsavn (sv) c get space for unpacking a call mpnew2 (i2, 80) c unpack a to find its length and allow printing in a1 format. call mpupk (a, r(i2+5), 71, len) c insert leading and trailing asterisks. call mpupk (5h *** , r(i2), 5, j) i = i2 + len call mpupk (4h ***, r(i+5), 4, j) c write using mpio. call mpio (r(i2), len+9, lun, 6h(80a1), err) c and terminate execution by calling mperr. c (do not call mpresn because possible recursion.) call mperr return end subroutine mpeul (g) c returns mp g = eulers constant (gamma = 0.57721566...) c to almost full multiple-precision accuracy. c the method is based on bessel function identities and was c discovered by edwin mc millan and r. brent. it is faster than the c method of sweeney (math. comp. 17, 1963, 170) used in earlier c versions of mpeul. time o(t**2). c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer i2, i3, i4, i5, k, n, sv, tg integer b, dummy(17), g(1), lun, m, mxr, r(1), sptr, t integer mptlb c save t etc. call mpsavn (sv) c increase t as necessary to get accurate result n = mptlb (iabs(t)) call mpgd3 (n, tg) c use truncation internally call mpsetr (0) c compute n so truncation error (bounded by pi*exp(-4n)) is less c than 0.01*(b**(-t)). c 1/8 + 1/20 = 7/40 .gt. ln(2)/4 n = n/8 + n/20 + 4 c compute ln(n). call mpnew (i2) call mplni (n, r(i2)) c allocate more temporary space call mpnew (i3) call mpnew (i4) call mpnew (i5) c set up main loop r(i2) = -r(i2) call mpstr (r(i2), r(i5)) call mpcim (1, r(i4)) call mpstr (r(i4), r(i3)) k = 0 c main loop starts here 10 k = k + 1 c reduce t here if possible if (k.gt.n) t = min0 (t, tg + r(i4+1) - r(i3+1)) c test for convergence if (t.lt.2) go to 20 call mpmuls (r(i4), n, n, k, k) call mpmuls (r(i5), n, n, k, 1) call mpadd (r(i5), r(i4), r(i5)) call mpdivi (r(i5), k, r(i5)) c increase t here t = tg call mpadd (r(i3), r(i4), r(i3)) call mpadd (r(i2), r(i5), r(i2)) c end of main loop if (r(i4).ne.0) go to 10 c compute final quotient, releasing some space first 20 t = tg sptr = i4 call mpdiv (r(i2), r(i3), r(i2)) c restore rndrl etc and round result call mpres2 (sv) call mprnd (r(i2), tg, g, iabs(r(sv)), 1) c restore sptr etc. call mpresn (sv) return end subroutine mpexp (x, y) c returns y = exp(x) for packed or unpacked mp x, unpacked mp y. c time is o(sqrt(t)m(t)). c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, dummy integer i2, k, rp, sv, tg integer b, dummy(19), lun, m, r(1), t, x(1), y(1) integer mpgd, mpget c save t etc. call mpsavn (sv) c check for x = 0 if (x(1).ne.0) go to 10 call mpcim (1, y) go to 50 c check if abs(x) .lt. 1 10 if (mpget (x, 2) .gt. 0) go to 20 c use mpexp1 here t = t + mpgd (100) tg = t call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) call mpexp1 (r(i2), r(i2)) call mpaddi (r(i2), 1, r(i2)) go to 40 c increase t as necessary. 20 rp = max0 (mpget (x, 2), 0) t = t + 1 + rp + mpgd (100) tg = t c use truncation rather than rounding internally if (r(sv+2) .eq. 1) call mpsetr (0) call mpnew (i2) c move x and divide by b**rp call mpmove (x, iabs(r(sv)), r(i2), tg) r(i2+1) = r(i2+1) - rp c use exp1 to compute exp(...) - 1 call mpexp1 (r(i2), r(i2)) call mpaddi (r(i2), 1, r(i2)) c see if result already obtained (rp .le. 0) if (rp.le.0) go to 40 c here have to raise result to b-th power rp times c underflow/overflow will occur here if x out of allowable c range. do 30 k = 1, rp 30 call mppwra (r(i2), iabs(b), r(i2)) c round result to ts digits, restore t etc. and return 40 call mprnd (r(i2), tg, y, iabs(r(sv)), 0) 50 call mpresn (sv) return end subroutine mpexp1 (x, y) c assumes that x and y are mp numbers, -1 .lt. x .lt. 1. c x may be packed or unpacked, y is unpacked. c returns y = exp(x) - 1 using an o(sqrt(t).m(t)) algorithm c described in - r. p. brent, the complexity of multiple- c precision arithmetic (in complexity of computational problem c solving, univ. of queensland press, brisbane, 1976, 126-165). c asymptotically faster methods exist, but are not useful c unless t is very large. see comments to mpatan and mppigl. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, ktunfl, dummy integer i, i2, i3, i4, ktu, la, mpget, mptlb, q, sv, tg, \$ b, dummy(11), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, r(1), rndrl, sptr, t, x(1), y(1), \$ lim, mpparn, mul c check for x = 0 y(1) = x(1) if (x(1) .eq. 0) return c save t, m etc. call mpsavn (sv) c check that abs(x) .lt. 1 if (mpget (x, 2) .gt. 0) call mperrm ( \$ 41habs(x) not less than 1 in call to mpexp1\$) c compute approximately optimal q. q = 0 if (mpget (x, 2) .lt. (-t)) go to 30 la = mptlb (t/3) c integer square root approximation here 20 q = q + 1 if ((q*q) .le. la) go to 20 q = max0 (0, q + mpget (x, 2)*mptlb(1)) c increase t as necessary. 30 call mpgd3 (q, tg) c allocate temporary space 40 call mpnew (i2) call mpnew (i3) call mpnew (i4) c use truncation rather than rounding internally. if (rndrl.eq.1) rndrl = 0 c increase m (to avoid underflow problems) and move x m = m + 2*t call mpmove (x, iabs(r(sv)), r(i2), tg) ktu = ktunfl c halve q times c can not use mpscal here as potentially recursive if (q .le. 0) go to 60 lim = mpparn(16)/2 mul = 1 do 50 i = 1, q mul = 2*mul if ((mul .lt. lim) .and. (mul .ne. b) .and. \$ (i .ne. q)) go to 50 call mpdivi (r(i2), mul, r(i2)) mul = 1 50 continue c if underflow occured set q = 0 and try again with new t etc. if (ktu .eq. ktunfl) go to 60 q = 0 call mpresn (sv) call mpsavn (sv) call mpgd3 (mptlb(iabs(t)), tg) go to 40 60 call mpstr (r(i2), r(i4)) call mpstr (r(i2), r(i3)) i = 1 tg = t c sum series, reducing t where possible 70 t = tg + 2 + r(i3+1) - r(i4+1) if (t.le.2) go to 80 t = min0 (t, tg) if (rndrl.ne.0) t = tg call mpmul (r(i2), r(i3), r(i3)) i = i + 1 call mpdivi (r(i3), i, r(i3)) t = tg call mpadd (r(i3), r(i4), r(i4)) c can only fall through next statement if underflow occurred. if (ktu.eq.ktunfl) go to 70 c restore working precision 80 t = tg c check if rounding up required if (rndrl.ne.3) go to 90 c rounding up here, so add abs(last term) r(i3) = iabs(r(i3)) call mpadd (r(i3), r(i4), r(i4)) c finished if q .le. 0 90 if (q.le.0) go to 110 c apply (x+1)**2 - 1 = x(2 + x) for q iterations do 100 i = 1, q call mpaddi (r(i4), 2, r(i2)) 100 call mpmul (r(i2), r(i4), r(i4)) c round result 110 call mprnd (r(i4), tg, y, iabs(r(sv)), 0) c restore stack pointer etc and return call mpresn (sv) return end integer function mpexpa (x) c returns the exponent of the packed or unpacked mp number x c (or large negative exponent -m if x is zero). integer x(1), mpget, mpparn c return -m if x zero, x(2) otherwise mpexpa = - mpparn (3) if (x(1) .ne. 0) mpexpa = mpget (x, 2) return end subroutine mpexpb (i, x) c sets exponent of packed or unpacked mp number x to i c unless x is zero (when exponent is unchanged). common /mpcom/ b, t, m, dummy integer b, dummy(20), i, m, t, two, x(1) two = 2 c return if x is zero if (x(1).eq.0) return c set exponent of x to i x(two) = i c update mxexpn and mnexpn. call mpupdt (i) c check for overflow and underflow if (i.gt.m) call mpovfl (x) if (i.le.(-m)) call mpunfl (x) return end subroutine mpfin (lunit, x) c reads mp x from unit lunit in free format. common r logical err integer inrecl, i2, lunit, mpparn, r(1), sv, x(1) call mpsavn (sv) inrecl = mpparn (18) call mpnew2 (i2, inrecl) c format (80a1) assumes inrecl .le. 80 call mpio (r(i2), inrecl, (-lunit), 6h(80a1), err) if (err) call mperr call mpin (r(i2), x, inrecl, err) if (err) call mperrm ( \$ 42herror return from mpin, called from mpfin\$) call mpresn (sv) return end subroutine mpflor (x, y) c sets y = floor (x), i.e. the largest integer not exceeding x. c x and y are mp numbers. c rounding is defined by the parameter rndrl in common /mpcom/ c (this is only relevant if x is large and negative) - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r integer i2, mpcomp, r(1), sv, x(1), y(1) c save t etc. and truncate x (towards zero) to an integer. call mpsavn (sv) call mpnew (i2) call mpcmim (x, r(i2)) c if x negative and not an integer need to subtract 1 if ((x(1).lt.0) .and. (mpcomp (x, r(i2)).ne.0)) \$ call mpaddi (r(i2), -1, r(i2)) c move result, restore everything and return. call mpstr (r(i2), y) call mpresn (sv) return end subroutine mpfout (x, n) c c floating-point output routine. writes mp number x on c default locical unit (lun) in floating-point format with n c significant decimal digits, n .ge. 2. c x may be packed or unpacked if mpfout is called directly. c c exponent field width is determined by exwid - default is c 6 (including exponent character and sign). c c exponent character is determined by expch (default is e). c c output base defaults to 10, but may be any integer in range c 2, ... , 16 (see comment on outbas in mpparn). c c for rounding options see comments in subroutine mpoute. c common r logical err, mpis integer e, fwidth, i, i2, i3, j, k, las, outbas, \$ sv, n, r(1), x(1), mpdigw, mpparn, \$ dollar, minus, plus, star data dollar, minus, plus, star /1h\$, 1h-, 1h+, 1h*/ c do nothing if n.lt.2 if (n.lt.2) return call mpsavn (sv) c field width is exwid + n + 2 fwidth = mpparn(17) + n + 2 outbas = mpparn (20) e = mpparn (21) c change exponent character if it is a digit. do 10 i = 1, outbas 10 if (mpis (e, mpdigw(i-1))) e = dollar c allocate fwidth words for working space, check b, t etc. call mpnew2 (i2, fwidth) c convert to printable form. call mpoute (x, r(i2), j, n+2) c set up exponent in character form. i3 = i2 + n + 2 r(i3) = e r(i3+1) = plus if (j .lt. 0) r(i3+1) = minus j = iabs (j) las = i2 + fwidth - 3 c loop for each digit of exponent 20 k = mod (j, outbas) j = j/outbas r(las+2) = mpdigw(k) las = las - 1 if (las .ge. i3) go to 20 if (j .eq. 0) go to 40 c here exponent is too large to fit in output field. las = i2 + fwidth - 3 do 30 i = i3, las 30 r(i+2) = star c call mpio to write character form of x. 40 call mpio (r(i2), fwidth, mpparn(4), 9h(1x,70a1), err) if (err) call mperrm ( \$ 43herror return from mpio, called from mpfout\$) c restore stack pointer and return. call mpresn (sv) return end subroutine mpgam (x, y) c computes mp y = gamma(x) for mp argument x, using c mpgamq if abs(x) .le. mxint/240-1 and 240*x is an integer, c otherwise using mplngm. space required is o(t**2), see comments c in subroutine mplngm. time = o(t**3). c rounding options not implemented, uses no guard digits. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer ix, i2, i3, kt, n, sv, \$ b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1), \$ mpcmpi, mpcmpq, mpget, mpparn, mptlb c save t etc., allocate temporary space, use truncation. call mpsavn (sv) call mpnew (i2) call mpnew (i3) call mpsetr (0) call mpabs (x, r(i3)) if (mpcmpi (r(i3), mpparn(16)/240 - 1) .gt. 0) go to 20 c see if 240*x is almost an integer. call mpmuli (x, 240, r(i3)) call mpcmi (r(i3), ix) c compare with ix and ix+1 because r(i3) could be just c below an integer. do 10 kt = 1, 2 call mpaddi (r(i3), -ix, r(i2)) if ((r(i2).eq.0).or. \$ (((r(i3+1)-r(i2+1)).ge.(t-1)).and. \$ (r(i3+2).ge.r(i2+2)))) go to 30 10 ix = ix + 1 c here x is large or not simple rational, c check if abs(x) very small. if (mpget (x, 2) .le. (-t)) go to 90 c now check sign of x 20 if (x(1).lt.0) go to 40 c x is positive so use mplngm directly after releasing some space. sptr = i2 call mplngm (x, y) c see if mpexp will give overflow call mpnew (i2) call mpdivi (y, iabs(m), r(i2)) c ln(2) .lt. 7/10 if (mpcmpq (r(i2), mptlb(7), 10) .ge. 0) go to 70 c usually safe to call mpexp here. call mpexp (y, y) go to 100 c x = ix/240 so use mpgamq unless x zero or negative integer. 30 if ((ix.le.0).and.(mod(iabs(ix), 240) .eq. 0)) go to 50 sptr = i2 call mpgamq (ix, 240, y) go to 100 c here x is negative, so use reflection formula 40 call mpabs (x, y) c subtract even integer to avoid errors near poles call mpdivi (y, 2, r(i3)) call mpcmf (r(i3), r(i3)) call mpmuli (r(i3), 2, r(i3)) call mpaddq (r(i3), 1, 2, r(i2)) call mpcmi (r(i2), n) c check for integer overflow in mpcmi if ((r(i3).ne.0).and.(r(i3+1).gt.0).and.(n.eq.0)) go to 70 call mpaddi (r(i3), -n, r(i3)) c now abs(r(i3)) .le. 1/2 and sign determined by n if (r(i3).ne.0) go to 60 50 call mperrm (44hx zero or negative integer in call to mpgam\$) 60 call mppi (r(i2)) call mpmul (r(i3), r(i2), r(i3)) call mpsin (r(i3), r(i3)) call mpmul (r(i3), y, r(i3)) if (r(i3).eq.0) go to 70 call mpdiv (r(i2), r(i3), r(i2)) r(i2) = -((-1)**n)*r(i2) c release some space before call to mplngm. sptr = i3 call mplngm (y, y) y(1) = -y(1) call mpexp (y, y) call mpmul (y, r(i2), y) go to 100 c here x was too large or too close to a pole 70 write (lun, 80) 80 format (26h *** overflow in mpgam ***) call mpovfl (y) go to 100 c here abs(x) is very small 90 call mprec (x, y) c restore everything and return. 100 call mpresn (sv) return end subroutine mpgamq (i, j, x) c returns x = gamma (i/j) where x is multiple-precision and c i, j are small integers. the method used is reduction of c the argument to (0, 1) and then a direct c expansion of the defining integral truncated at a c sufficiently high limit, using 2t digits to c compensate for cancellation. c time is o(t**2) if i/j is not too large. c if i/j .gt. 100 (approximately) it is faster to use c mpgam (if enough space is available). c mpgamq is very slow if i/j is very large, because c the relation gamma(x+1) = x*gamma(x) is used repeatedly. c mpgamq could be speeded up by using the asymptotic series or c continued fraction for (integral from n to infinity of c u**(i/j-1)*exp(-u)du). c rounding options not yet implemented. common r common /mpcom/ b, t, dummy integer ibtn, id, ij, il, in, is, is2, i2, i3, js, mxint, n, sv integer b, dummy(21), i, j, r(1), t, ts2, ts3, x(1) integer mpchgb, mpparn, mptlb c save t etc. call mpsavn (sv) c use truncated arithmetic. call mpsetr (0) mxint = mpparn (16) if (j .eq. 0) call mperrm (24hj = 0 in call to mpgamq\$) is = i*isign (1, j) js = iabs (j) c now js is positive. reduce to lowest terms. call mpgcd (is, js) ij = mxint/js c see if js = 1, 2, or .gt. 2 if (js - 2) 20, 10, 40 c js = 2 here, for speed treat as special case 10 call mppi (x) call mpsqrt (x, x) go to 30 c js = 1 here, check that is is positive 20 if (is .le. 0) call mperrm ( \$ 39hi/j zero or negative in call to mpgamq\$) c i/j = positive integer here call mpcim (1, x) 30 is2 = 1 go to 60 c js .gt. 2 here so reduce to (0, 1) 40 is2 = mod (is, js) if (is2.lt.0) is2 = is2 + js c now 0 .lt. is2 .lt. js. compute upper limit of integral c n .ge. t*ln(b) n = mptlb ((7*t + 9)/10) ibtn = mxint/n ts3 = t + 2 c increase t to compensate for cancellation c increase .ge. n/ln(b) t = t + mpchgb (iabs(b), 2, (3*n+1)/2) ts2 = t c allocate temporary space call mpnew (i2) call mpnew (i3) call mpcim (n, r(i2)) call mpstr (r(i2), r(i3)) il = 0 in = js - is2 id = is2 c main loop 50 il = il + 1 c if terms decreasing may decrease t if (il.ge.n) t = r(i3+1) + ts3 t = max0 (2, min0 (t, ts2)) c check if in-js or id+js might overflow. if (max0 (iabs(in), id) .gt. (mxint-js)) go to 70 in = in - js id = id + js call mpmuls (r(i3), n, in, il, id) t = max0 (t, ts3) call mpadd (r(i2), r(i3), r(i2)) c loop until exponent small if ((r(i3).ne.0).and.(r(i3+1).ge.(-r(sv)))) go to 50 c restore t t = r(sv) call mpmulq (r(i2), js, is2, x) call mpqpwr (n, 1, is2-js, js, r(i3)) call mpmul (x, r(i3), x) c now x is gamma (is2/js), so use the recurrence relation c repeatedly to get gamma (i/j) (slow if i/j is large). 60 in = 1 id = 1 if (is - is2) 100, 110, 80 c j must have been too large. 70 call mperrm (30hj too large in call to mpgamq\$) c restore t etc. and return. 110 call mpresn (sv) return 80 in = in*is2 id = id*js is2 = is2 + js if ((id.le.ij).and.(iabs(in).le.(mxint/iabs(is2))) \$ .and.(is.ne.is2)) go to 80 90 call mpmulq (x, in, id, x) go to 60 100 in = in*js id = id*is is = is + js if ((in.le.ij).and.(iabs(id).le.(mxint/iabs(is))) \$ .and.(is.ne.is2)) go to 100 go to 90 end subroutine mpgcd (k, l) c returns k = k/gcd and l = l/gcd, where gcd is the c greatest common divisor of k and l. c save input parameters in local variables integer i, is, j, js, k, l i = k j = l is = iabs(i) js = iabs(j) if (js.eq.0) go to 30 c euclidean algorithm loop 10 is = mod (is, js) if (is.eq.0) go to 20 js = mod (js, is) if (js.ne.0) go to 10 js = is c here js is the gcd of i and j 20 k = i/js l = j/js return c if j = 0 return (1, 0) unless i = 0, then (0, 0) 30 k = 1 if (is.eq.0) k = 0 l = 0 return end subroutine mpgcda (x, y, z) c returns z = greatest common divisor of x and y. c gcd (x, 0) = gcd (0, x) = abs(x), gcd (x, y) .ge. 0. c x, y and z are integers represented as mp numbers, c and must satisfy abs(x) .lt. b**t, abs(y) .lt. b**t c x and/or y may be packed or unpacked, z is unpacked. c time o(t**2). the parameter rndrl is irrelevant as result exact. common r common /mpcom/ b, t, dummy integer i, iq, is, i2, i3, i4, j, sv integer b, dummy(21), mpcmp, r(1), t, x(1), y(1), z(1) logical mpintg c save t etc., allocate working space, use truncated arithmetic. call mpsavn (sv) call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpsetr (0) call mpunpk (x, r(i2)) call mpunpk (y, r(i3)) r(i2) = iabs (r(i2)) r(i3) = iabs (r(i3)) c check that x and y are exact integers if (.not. (mpintg (r(i2)) .and. mpintg (r(i3)))) \$ call mperrm (37hx or y non-integer in call to mpgcda\$) c check for x or y zero if (x(1).ne.0) go to 20 10 t = r(sv) call mpstr (r(i3), z) go to 140 20 if (y(1).ne.0) go to 30 call mpstr (r(i2), z) go to 140 c check that abs(x), abs(y) .lt. b**t 30 if (max0 (r(i2+1), r(i3+1)) .gt. t) call mperrm ( \$ 35hx or y too large in call to mpgcda\$) c start of main euclidean algorithm loop 60 if (r(i2).eq.0) go to 10 if (mpcmp (r(i2), r(i3))) 70, 10, 80 c exchange pointers only 70 is = i2 i2 = i3 i3 = is c check for small exponent 80 if (r(i2+1).le.2) go to 110 c reduce t (trailing digits must be zero) t = r(i2+1) call mpstr (r(i3), r(i4)) c force exponents to be equal r(i4+1) = r(i2+1) c get first two digits iq = b*r(i2+2) + r(i2+3) if (mpcmp (r(i2), r(i4)) .ge. 0) go to 90 c reduce exponent by one r(i4+1) = r(i4+1) - 1 c underestimate quotient iq = iq/(r(i4+2)+1) go to 100 c lehmers method would save some mp operations but not very c many unless we could use double-precision safely. c see american math. monthly 45 (1938), 227 - 233. 90 iq = max0 (1, iq/(b*r(i4+2) + r(i4+3) + 1)) 100 call mpmuli (r(i4), iq, r(i4)) call mpsub (r(i2), r(i4), r(i2)) go to 60 c here safe to use integer arithmetic 110 call mpcmi (r(i2), i) call mpcmi (r(i3), j) t = r(sv) 120 i = mod (i, j) if (i.eq.0) go to 130 j = mod (j, i) if (j.ne.0) go to 120 j = i 130 call mpcim (j, z) c restore everything and return. 140 call mpresn (sv) return end subroutine mpgcdb (x, y) c returns (x, y) as (x/z, y/z) where z is the gcd of x and y. c x and y are integers represented as mp numbers, c and must satisfy abs(x) .lt. b**t, abs(y) .lt. b**t c time o(t**2). the parameter rndrl is irrelevant as result is exact. common r integer is, iz, i2, sv integer r(1), x(1), y(1) integer mpcmpi, mpcomp, mpparn c save t etc., use truncated arith., allocate temporary space. call mpsavn (sv) call mpsetr (0) call mpnew (i2) c find gcd of x and y using mpgcda call mpgcda (x, y, r(i2)) c check for x and y equal (when may coincide) if (mpcomp (x, y) .ne. 0) go to 10 is = x(1) call mpcim (is, x) call mpstr (x, y) go to 30 c check if gcd is small. 10 if (mpcmpi (r(i2), mpparn (16)) .gt. 0) go to 20 c here it is safe to convert gcd to single-precision integer. call mpcmi (r(i2), iz) if (iz.eq.1) go to 30 call mpdivi (x, iz, x) call mpdivi (y, iz, y) go to 30 c here gcd is large. 20 call mpdiv (x, r(i2), x) call mpdiv (y, r(i2), y) c restore everything and return. 30 call mpresn (sv) return end integer function mpgd (n) c returns ceiling (ln (max (1, abs(n))) / ln(b)), c i.e. the minimum value of j .ge. 0 such that c b**j .ge. abs(n) . c this function is useful for computing the number of guard digits c required for various multiple-precision calculations. common /mpcom/ b, dummy integer b, dummy(22), j, k, n c check that b is legal. if (b.lt.2) call mpchk c set up loop j = 0 k = iabs(n) - 1 go to 20 c loop to compute result 10 j = j + 1 k = k/b 20 if (k.gt.0) go to 10 c return result j mpgd = j return end subroutine mpgd3 (n, tg) c sets t = t + 1 + mpgd (100*n) if safe to call c mpgd, effectively almost the same if not. c also sets tg = new value of t. common /mpcom/ b, t, dummy integer b, dummy(21), i, mpgd, mpparn, n, t, tg c check whether safe to form 100*n if (iabs(n) .ge. (mpparn(16)/100)) go to 10 c here it is i = mpgd (100*n) go to 20 c here it is not, so call mpgd twice 10 i = mpgd (n) + mpgd (100) c set t and tg 20 t = t + i + 1 tg = t return end logical function mpge (x, y) c returns logical value of (x .ge. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mpge = (mpcomp(x,y) .ge. 0) return end integer function mpget (x, n) c returns x(n). necessary to avoid some compiler diagnostics. c also useful to avoid pfort verifier unsafe reference warnings. integer n, x(1) mpget = x(n) return end logical function mpgt (x, y) c returns logical value of (x .gt. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mpgt = (mpcomp(x,y) .gt. 0) return end subroutine mphank (x, nu, y, error) c tries to compute the bessel function j(nu,x) using hankels c asymptotic series. nu is a nonnegative integer not too large. c error is logical, x and y are mp numbers. c returns error = .false. if successful (result in y), c error = .true. if unsuccessful (y unchanged) c time = o(t**3). c called by mpbesj, not recommended for independent use. c rounding options not implemented, uses no guard digits. common r integer ie, i2, i3, i4, i5, i6, i7, k, lg, sv integer nu, r(1), x(1), y(1) integer mpcmpi, mpparn, mptlb logical error c save t etc., use truncation internally. call mpsavn (sv) call mpsetr (0) error = .true. c give error return if nu is negative. if (nu.lt.0) go to 20 c allocate temporary space. call mpnew (i2) c work with abs(x) call mpabs (x, r(i2)) c check if abs(x) clearly too small for asymptotic series if (mpcmpi (r(i2), mptlb(iabs(r(sv)))/3) .le. 0) go to 20 c allocate more temporary space. call mpnew (i3) call mpnew (i4) call mpnew (i5) call mpnew (i6) call mpnew (i7) call mppwr (x, -2, r(i3)) call mpdivi (r(i3), -64, r(i3)) call mpcim (1, r(i4)) r(i5) = 0 call mpstr (r(i4), r(i6)) ie = 1 k = 0 lg = mpparn(16)/2 c loop to sum two asymptotic series 10 k = k + 2 c return with error = .true. if k too large. if ((nu+k) .ge. lg) go to 20 c error return if terms increasing if (r(i6+1).gt.ie) go to 20 ie = r(i6+1) call mpmuls (r(i6), 2*(nu+k)-3, 2*(nu-k)+3, k-1, 1) call mpadd (r(i5), r(i6), r(i5)) call mpmuls (r(i6), 2*(nu+k)-1, 2*(nu-k)+1, k, 1) call mpmul (r(i6), r(i3), r(i6)) call mpadd (r(i4), r(i6), r(i4)) c loop if terms not sufficiently small yet if ((r(i6).ne.0).and.(r(i6+1).gt.(-r(sv)))) go to 10 c end of asymptotic series, now compute result call mpdiv (r(i5), r(i2), r(i5)) call mpdivi (r(i5), 8, r(i5)) c compute pi/4 (slightly more accurate than calling c mppi and dividing by four) call mpart1 (5, r(i6)) call mpmuli (r(i6), 4, r(i6)) call mpart1 (239, r(i3)) call mpsub (r(i6), r(i3), r(i3)) c avoid too much cancellation in subtracting multiple of pi call mpmuli (r(i3), mod (2*nu+1, 8), r(i6)) call mpsub (r(i2), r(i6), r(i6)) call mpcis (r(i6), r(i6), r(i7), .true.) call mpmul (r(i4), r(i6), r(i4)) call mpmul (r(i5), r(i7), r(i5)) call mpsub (r(i4), r(i5), r(i4)) call mpmul (r(i3), r(i2), r(i3)) call mpmuli (r(i3), 2, r(i3)) call mproot (r(i3), -2, r(i3)) call mpmul (r(i3), r(i4), r(i3)) c correct sign of result if (mod (nu, 2) .ne. 0) r(i3) = r(i3)*x(1) error = .false. call mpstr (r(i3), y) c restore everything and return 20 call mpresn (sv) return end subroutine mpimul (ix, y, z) c multiplies single-precision integer ix by mp y giving mp z. c rounding defined by parameter rndrl in common /mpcom/ - c see subroutine mpnzr. integer ix, y(1), z(1) call mpmuli (y, ix, z) return end subroutine mpin (c, x, n, error) c c converts the decimal number (read under na1 format) c in c(1) ... c(n) to a multiple-precision number c in x. if c represents a valid number, error is returned c as .false. if c does not represent a valid number, error c is returned as .true. and x as zero. c leading and trailing blanks are allowed, embedded blanks c (except between the number and its sign) are forbidden. c if there is no decimal point one is assumed to lie just to c the right of the last decimal digit. c c numbers may have decimal exponents, provided twice the c exponent is representable as a single-precision integer. c the exponent may be preceeded by any character except a c digit, blank or period. the exponent has an optional sign, c and the special character preceeding it may be omitted if a c sign is present. c c examples - c c valid numbers invalid numbers c c - 123456789 12 345 c 3.14159 123.456e -67 c -44. 1.2.3 c .0001234 e123 c 123.456d789 64.4e+ c -.1234566-789 ++12.3 c +999+88 11e3. c c for efficiency choose b a power of 10. c x is an mp number, c an integer array, n integer, error logical. c c rounding determined by rndrl in common /mpcom/ as follows - c c rndrl = 0 or 1 - round to (approximately) nearest representable c base b, t-digit number. error is less than c 0.6 units in the last (base b) place. c rndrl = 2 - round down (towards -infinity) to a base b, c t-digit number. c rndrl = 3 - round up (towards +infinity) to a base b, c t-digit number. c c the default input base is 10 (i.e. decimal numbers), but c any base in the range 2, ... , 16 may be used. the base is c determined by the parameter inbase in common /mpcom/ - c see mpparn and mpset. c common r common /mpcom/ b, t, dummy integer ex, i, ib, inbase, ip, i2, j, k, mxint, s, sex, sv, tg, \$ b, c(1), dummy(21), n, r(1), t, x(1), \$ mpdigv, mpgd, mpparn, blank, minus, period, plus logical de, df, error, first, mpis, point, second data blank, minus, period, plus /1h , 1h-, 1h., 1h+/ c save t etc. call mpsavn (sv) mxint = mpparn (16) c get input base (default set by mpset is 10). inbase = mpparn (19) c check that 2 .le. inbase .le. 16 and n .gt. 0. if ((2 .gt. inbase) .or. (inbase .gt. 16) .or. (n.le.0)) go to 30 c increase t as necessary, allocate temporary space. t = t + mpgd (mxint) call mpgd3 (n, tg) call mpnew (i2) c use truncation rather than rounding internally. if (r(sv+2) .eq. 1) call mpsetr (0) first = .true. second = .false. r(i2) = 0 s = 1 error = .false. de = .false. df = .true. point = .false. ip = 0 k = 0 ex = 0 sex = 1 ib = mxint/(2*inbase) c scan c from left, skipping blanks 10 k = k + 1 if (.not. mpis (c(k), blank)) go to 60 20 if (k.lt.n) go to 10 c error here. 30 x(1) = 0 error = .true. c restore everything and return. 40 call mpresn (sv) return c sign here. 50 if (.not. first) go to 130 if (mpis (c(k), minus)) s = -1 c take care for directed roundings. call mprevr (-s) first = .false. go to 20 c scan number before exponent field. 60 if (mpis (c(k), minus) .or. mpis (c(k), plus)) go to 50 first = .false. j = mpdigv (c(k)) if (j .ge. 0) go to 80 c not sign or digit here. if (.not. mpis (c(k), period)) go to 70 if (point) go to 30 point = .true. go to 90 70 if (mpis (c(k), blank)) go to 100 c must be exponent indicator go to 140 c digit, so continue forming number 80 call mpmuli (r(i2), inbase, r(i2)) call mpaddi (r(i2), j, r(i2)) df = .false. if (point) ip = ip + 1 90 k = k + 1 if (k .le. n) go to 60 c check that at least one digit before exponent, and at least c one digit in exponent (if present). 100 if (de .or. df) go to 30 if (k .gt. n) go to 120 c check that all characters to right are blanks do 110 i = k, n if (.not. mpis (c(i), blank)) go to 30 110 continue c end of input field, scale appropriately. 120 call mpscal (r(i2), inbase, ex*sex - ip) c round result and fix up sign. call mprnd (r(i2), tg, x, iabs(r(sv)), 0) x(1) = s*x(1) go to 40 c deal with sign of exponent 130 if (second) go to 30 second = .true. if (mpis (c(k), minus)) sex = -1 c exponent follows. 140 de = .true. 150 k = k + 1 if (k.gt.n) go to 100 if (mpis (c(k), plus) .or. mpis (c(k), minus)) go to 130 second = .true. j = mpdigv (c(k)) if (j .lt. 0) go to 100 c check if exponent too large if (ex.ge.ib) go to 30 c otherwise incorporate next digit ex = inbase*ex + j de = .false. go to 150 end subroutine mpine (c, x, n, j, error) c same as mpin except that x is multiplied by inbase**j. c the default value of inbase is ten, but this may be changed c (see comments in mpparm). integer c(1), j, mpparn, n, x(1) logical error call mpin (c, x, n, error) call mpscal (x, mpparn (19), j) return end subroutine mpinf (x, n, unit, iform, err) c reads n words from logical unit iabs(unit) using format in iform, c then converts to mp number x using routine mpin. c iform should contain a format which allows for reading n words c in a1 format, e.g. 6h(80a1) c err returned as true if mpin could not interpret input as c an mp number or if n not positive, otherwise false. c (see also comments in mpio.) c if err is true then x is returned as zero. c space required = n + o(t) words. c for rounding options see comments in mpin. common r integer iform(1), i2, n, r(1), sv, unit, x(1) logical err call mpsavn (sv) c allocate n words on stack. call mpnew2 (i2, n) c read n words under format iform. call mpio (r(i2), n, (-unit), iform, err) x(1) = 0 c convert to mp number unless read error if (.not. err) call mpin (r(i2), x, n, err) c restore stack pointer and return. call mpresn (sv) return end subroutine mpinit (x) c c declares blank common and common /mpcom/ (used by mp package) and c calls mpset2 to initialize parameters. x is a dummy integer argument. c the augment declaration c initialize mp c causes a call to mpinit(mp) to be generated. c c *** assumes output unit 6, 64 decimal places, maximum of c *** 16 mp digits, working space 4000 words. if the augment c *** description deck is changed this routine should c *** be changed accordingly, and vice versa. c *** if wordlength is greater than 16 bits, space can be saved as c *** less than 25 mp digits will be required for 43 decimal place c *** accuracy. common r common /mpcom/ mppars integer x(1), mppars(23) c c the statements integer r(4000) call mpset2 (6, 64, 18, 1, 4000) c are a special case of c integer r(mxr) c call mpset2 (lun, decpl, t+2, mnsptr, mxr) c where c lun is the logical unit for output, c decpl is the equivalent number of decimal places required, c t is the number of mp digits, and c the working area free for use by mp is words c mnsptr, mnsptr+1, ... , mxr of blank common. c c to change the precision, modify the dimensions in the c declare statements in the augment description deck - c the dimension (mt2) for type multiple should be at least t+2, and c for type multipak should be exactly int ((mt2+1)/2). c see function mpdigs for the number of mp digits required to give c the equivalent of any desired number of floating decimal places. c c *** on some systems a declaration of blank common and common /mpcom/ c *** in the main program may be necessary. if so, declare c *** common mpwork c *** common /mpcom/ mppars c *** integer mpwork (500), mppars (23) c *** (or, more generally, c *** integer mpwork (mxr), mppars (23) ) c *** in the main program. return end logical function mpintg (x) c returns .true. if the mp number x is an exact integer, c .false. otherwise. common r integer i2, r(1), sv, x(1) call mpsavn (sv) call mpnew (i2) call mpcmf (x, r(i2)) mpintg = (r(i2) .eq. 0) call mpresn (sv) return end subroutine mpio (c, n, unit, iform, err) c c if unit .gt. 0 writes c(1), ... , c(n) in format iform c if unit .le. 0 reads c(1), ... , c(n) in format iform c in both cases uses logical unit iabs(unit). c if iform = 1hu then unformatted i/o is performed. c c err is returned as true if n non-positive, otherwise false. c we would like to return err as true if read/write error detected, c but this can not be done with ansi standard fortran (1966). c c univac ascii fortran (ftn 5r1ae) does not work if iform c is declared with dimension 1. most fortrans do though. c c integer n (univac fortran v does not allow) integer c(n), iform(n), iu, u, unit logical err, mpis data u /1hu/ err = (n.le.0) if (err) return iu = iabs(unit) if (mpis (iform(1), u)) go to 10 if (unit .gt. 0) write (iu, iform) c if (unit .le. 0) read (iu, iform) c return c unformatted i/o here 10 if (unit .gt. 0) write (iu) c if (unit .le. 0) read (iu) c return end logical function mpis (j, k) c c assumes that j and k are words containing one character in left-most c position with blank fill (read under a1 format or set by 1hx data c statement or unpacked by mpupk). returns .true. if and only if the c characters are equal. on some machines it may be necessary to c mask off other characters or shift right before comparison. c c on some machines change integer to real in next line. integer j, k c c on burroughs b6700 replace .eq. by .is. in next line. mpis = (j.eq.k) c return end subroutine mpk3v (mpxx, x, y, z) c calls mpxx (x, y, z) after unpacking x and y. c x and y are packed or unpacked mp numbers, c z is unpacked mp number, mpxx is a subroutine c taking three (unpacked) mp arguments. common r integer i2, i3, r(1), sv, x(1), y(1), z(1) call mpsavn (sv) call mpnew (i2) call mpnew (i3) call mpunpk (x, r(i2)) call mpunpk (y, r(i3)) call mpxx (r(i2), r(i3), z) call mpresn (sv) return end subroutine mpk3v2 (mpxx, x, n, y) c calls mpxx (x, n, y) after unpacking x. c x is packed or unpacked mp number, n integer, c y unpacked mp number, mpxx a subroutine c taking three arguments (mp, integer, mp). common r integer i2, n, r(1), sv, x(1), y(1) call mpsavn (sv) call mpnew (i2) call mpunpk (x, r(i2)) call mpxx (r(i2), n, y) call mpresn (sv) return end subroutine mpkadd (x, y, z) c adds x and y, forming result in z, where x and y are packed or c unpacked mp numbers, z is an unpacked mp number. c for further details see mpadd. integer x(1), y(1), z(1) external mpadd call mpk3v (mpadd, x, y, z) return end subroutine mpkdiv (x, y, z) c divides x and y, forming result in z, where x and y are c packed or unpacked mp numbers, z is unpacked. c for further details see mpdiv. integer x(1), y(1), z(1) external mpdiv call mpk3v (mpdiv, x, y, z) return end subroutine mpkdvi (x, n, y) c returns y = x/n, for packed or unpacked mp x, c integer n, unpacked mp y. c for further details see mpdivi. integer n, x(1), y(1) external mpdivi call mpk3v2 (mpdivi, x, n, y) return end subroutine mpkiml (n, x, y) c returns y = x*n, for packed or unpacked mp x, c integer n, unpacked mp y. c for further details see mpmuli. integer n, x(1), y(1) external mpmuli call mpk3v2 (mpmuli, x, n, y) return end subroutine mpkmli (x, n, y) c returns y = x*n, for packed or unpacked mp x, c integer n, unpacked mp y. c for further details see mpmuli. integer n, x(1), y(1) external mpmuli call mpk3v2 (mpmuli, x, n, y) return end subroutine mpkmul (x, y, z) c multiplies x and y, forming result in z, where x and y are c packed or unpacked mp numbers, z is unpacked. c for further details see mpmul. integer x(1), y(1), z(1) external mpmul call mpk3v (mpmul, x, y, z) return end subroutine mpksub (x, y, z) c subtracts x and y, forming result in z, where x and y are packed c or unpacked mp numbers, z is an unpacked mp number. c for further details see mpsub. integer x(1), y(1), z(1) external mpsub call mpk3v (mpsub, x, y, z) return end subroutine mpl235 (i, j, k, x) c called by mplni, not recommended for independent use. c returns mp x = ln((2**i)*(3**j)*(5**k)), for integer i, j and k. c the method requires time o(t**2). ln(81/80), ln(25/24) and c ln(16/15) are calculated first. mpl235 could be speeded c up if these constants were precomputed and saved. c assumed that i, j and k not too large. c rounding options not implemented and no guard digits used. common r common /mpcom/ b, t, dummy integer c(3), d(3), i2, i3, n, q, sv, two integer b, dummy(21), i, j, k, mpparn, r(1), t, x(1) data d(1), d(2), d(3) /161, 49, 31/ two = 2 c save t etc. and allocate temporary space call mpsavn (sv) call mpnew (i2) call mpnew (i3) c use truncation internally. call mpsetr (0) x(1) = 0 if (max0(iabs(i), iabs(j), iabs(k)) .ge. (mpparn(16)/68)) \$ call mperrm (37hi, j or k too large in call to mplni\$) c(1) = 3*i + 5*j + 7*k c(2) = 5*i + 8*j + 12*k c(3) = 7*i + 11*j + 16*k do 20 q = 1, 3 t = r(sv) call mpcqm (2*c(q), d(q), r(i2)) call mpstr (r(i2), r(i3)) call mpadd (x, r(i3), x) n = 1 10 n = n + 2 if (r(i2).eq.0) go to 20 c reduce t if possible t = r(sv) + r(i2+1) + 2 - x(two) if (t.le.2) go to 20 t = min0 (t, r(sv)) call mpmuls (r(i2), 1, 1, d(q), d(q)) call mpdivi (r(i2), n, r(i3)) t = r(sv) call mpadd (r(i3), x, x) go to 10 20 continue c restore stack pointer etc. and return call mpresn (sv) return end subroutine mplarg (mxint) c c ************************* (see the mp c *** machine dependent *** users guide for c ************************* conversion hints) c c returns mxint .le. the maximum representable integer of c the form 2**j - 1 . c c integer arithmetic must be performed exactly on integers c in the range -mxint ... +mxint, so on some machines c mxint must be less than the largest representable integer, c e.g. this applies on burroughs b6700 and cyber 76 series. c logical cyber, b6700, pdp11 integer i, k, mx, mxint, wdlen c c wdlen may be set to the wordlength (in bits) in the following c data statement if the code below for computing mxint does c not work. c data wdlen /32/ c c on cyber 76, burroughs 6700, or pdp 11 machines set c the corresponding logical variable to .true. c data cyber /.false./ data b6700 /.false./ data pdp11 /.false./ c c the following assumes at least 2 characters per word, c characters 0 and 1 have consecutive binary codes, and c leftmost character is in most significant position. c if these assumptions fail, wdlen (or one of cyber, b6700, c pdp11) must be set above. c c for a less restrictive algorithm see - algorithms to reveal c the representation of characters, integers and floating-point c numbers (by j. e. george), acm trans. math. software 1(1975), c 210-216. unlike george, we avoid the use of an external unit. c k = wdlen if (cyber) k = 49 if (b6700) k = 40 if (pdp11) k = 16 if (k .gt. 0) go to 10 c here try to compute lower bound on mxint. call mpsuba (2h10, 2h00, mxint) call mpsuba (2h01, 2h00, mx) mxint = 2*(mxint*(mxint/(4*mx)) - 1) + 1 return c here k is lower bound on wordlength, compute 2**(k-1) - 1 10 mxint = 0 do 20 i = 2, k 20 mxint = 2*mxint + 1 return end logical function mple (x, y) c returns logical value of (x .le. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mple = (mpcomp(x,y) .le. 0) return end subroutine mplg10 (x, y) c returns y = log10(x), for mp x and y, using mpln and mplni. c x may be packed or unpacked, y is unpacked. c for comments on time, see subroutine mpln. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r integer i2, i3, r(1), sv, tg, x(1), y(1) c save t etc. call mpsavn (sv) c use truncation rather than rounding. if (r(sv+2).eq.1) call mpsetr (0) c increase t and allocate working space. call mpgd3 (1, tg) call mpnew (i2) c compute ln(x)/ln(10), taking care for directed roundings. call mpmove (x, iabs(r(sv)), r(i2), tg) call mpln (r(i2), r(i2)) call mpnew (i3) call mprevr (r(i2)) call mplni (10, r(i3)) call mprevr (r(i2)) call mpdiv (r(i2), r(i3), r(i2)) c round result, restore t etc. and return. call mprnd (r(i2), tg, y, iabs(r(sv)), 1) call mpresn (sv) return end subroutine mpli (x, y) c returns y = li(x) = logarithmic integral of x c = (principal value integral from 0 to x of c du/log(u)), c using mpei. x and y are mp numbers, x .ge. 0, x .ne. 1. c x may be packed or unpacked, y is unpacked. c error in y could be induced by an o(b**(1-t)) relative c perturbation in x followed by similar perturbation in y. c thus relative error in y is small unless x is close to c 1 or to the zero 1.45136923488338105028... of li(x). c time is o(t.m(t)). c rounding options not yet implemented. integer mpcmpi, sv, x(1), y(1) if (x(1)) 10, 20, 30 c here x negative, give error message 10 call mperrm (25hx .lt. 0 in call to mpli\$) c li(0) = 0 20 y(1)= 0 return c here x is positive, see if equal to 1 30 if (mpcmpi (x, 1) .eq. 0) call mperrm ( \$ 25hx .eq. 1 in call to mpli\$) c here x positive and .ne. 1, so use ei(ln(x)), truncated arithmetic. call mpsavn (sv) call mpsetr (0) call mpln (x, y) call mpei (y, y) call mpresn (sv) return end subroutine mpln (x, y) c returns y = ln(x), for mp x and y, using mplni and mplns. c x may be packed or unpacked, y is unpacked. c restriction - integer part of ln(x) must be representable as a c single-precision integer. time is o(sqrt(t).m(t) + t**2). c for small integer x, mplni is faster. c asymptotically faster methods exist (eg the gauss-salamin c method, see mplngs), but are not useful unless t is large. c see comments to mpatan, mpexp1 and mppigl. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, dummy integer i2, i3, j, k, sv, tg integer b, dummy(22), r(1), x(1), y(1) c save t etc. call mpsavn (sv) c increase working precision. call mpgd3 (1, tg) c allocate temporary space call mpnew (i2) call mpnew (i3) c use truncated arithmetic. call mpsetr (0) c check that x is positive if (x(1) .le. 0) call mperrm ( \$ 30hx nonpositive in call to mpln\$) c move x to local storage call mpmove (x, iabs(r(sv)), r(i2), tg) c check if x is close to 1 call mpaddi (r(i2), -1, r(i3)) if (r(i3) .eq. 0) go to 10 if (r(i3+1) .ge. 0) go to 20 c here x close to 1, so use mplns directly 10 call mplns (r(i3), r(i2)) go to 30 c take off exponent and first two digits 20 j = b*r(i2+2) + r(i2+3) k = r(i2+1) - 2 r(i2+1) = 2 call mpdivi (r(i2), j, r(i2)) call mpaddi (r(i2), -1, r(i2)) c now should be safe to call mplns call mplns (r(i2), r(i2)) c now use mplni to correct result call mplni (iabs(b), r(i3)) call mpmuli (r(i3), k, r(i3)) call mpadd (r(i2), r(i3), r(i2)) call mplni (j, r(i3)) call mpadd (r(i2), r(i3), r(i2)) c round result 30 call mpsetr (iabs(r(sv+2))) call mprnd (r(i2), tg, y, iabs(r(sv)), 1) c restore t etc. and return. call mpresn (sv) return end subroutine mplngm (x, y) c returns mp y = ln(gamma(x)) for positive mp x, using stirlings c asymptotic expansion. slower than mpgamq (unless x large) c and uses more space, so use mpgamq and mpln if x is rational and c not too large, say x .le. 200. time is o(t**3). c x may be packed or unpacked, y is unpacked. c space required is o(t**2) words (more precisely, space is c nl*t/2 + o(t) words, where c nl = t*log2(b)/12 + o(1) (see below)). c rounding options not implemented, uses no guard digits. common r common /mpcom/ b, dummy integer i, ip, i2, i3, i4, i5, mxint, nl, nl2, nt, p, q, sv, xl integer b, dummy(22), r(1), x(1), y(1) integer mpchgb, mpcmpi, mpparn, mptlb c save t etc. call mpsavn (sv) c estimate c nl = number of terms required in stirlings c approximation, and c xl = lower bound on x for sufficient accuracy. c the constant 12 was determined empirically to minimize time. nl = mptlb(iabs(r(sv)))/12 + 2 nt = r(sv) + 1 mxint = mpparn(16) nl2 = 2*nl q = 1 do 20 i = 2, nl2 10 if (q .le. (mxint/i)) go to 20 q = q/b + 1 nt = nt + 1 go to 10 20 q = i*q xl = 1 30 xl = xl + 1 c 6*xl + xl/4 .le. 2*pi*xl if (mpchgb (6*xl + xl/4, iabs(b), nt) .gt. nl2) go to 30 c check that x is positive if (x(1) .le. 0) call mperrm ( \$ 27hx .le. 0 in call to mplngm\$) c allocate temporary space, use truncated arithemtic. call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpsetr (0) c move x and set y = 0 call mpmove (x, iabs(r(sv)), r(i3), iabs(r(sv))) y(1) = 0 c see if x large enough to use asymptotic series if (mpcmpi (r(i3), xl) .ge. 0) go to 50 c here x not large enough, so increase using the c identity gamma(x+1) = x*gamma(x) to correct result call mpcim (1, y) 40 call mpmul (y, r(i3), y) call mpaddi (r(i3), 1, r(i3)) if (mpcmpi (r(i3), xl) .lt. 0) go to 40 call mpln (y, y) y(1) = -y(1) c compute first terms in stirlings approximation 50 call mpln (r(i3), r(i4)) call mpaddq (r(i3), -1, 2, r(i2)) call mpmul (r(i2), r(i4), r(i4)) call mpsub (r(i4), r(i3), r(i4)) call mpadd (y, r(i4), y) call mppi (r(i4)) call mpmuli (r(i4), 2, r(i4)) call mpln (r(i4), r(i4)) call mpdivi (r(i4), 2, r(i4)) call mpadd (y, r(i4), y) c if x very large can return here if (r(i3+1).ge.r(sv)) go to 70 call mppwr (r(i3), -2, r(i4)) p = (r(sv)+3)/2 c allocate space for bernoulli numbers. call mpnew2 (i5, nl*p) c compute bernoulli numbers required (much time could be c saved if these were precomputed) call mpbern ((-nl), p, r(i5)) c sum asymptotic series do 60 i = 1, nl ip = i5 + (i-1)*p call mpunpk (r(ip), r(i2)) call mpmuls (r(i2), 1, 1, 2*i, 2*i-1) call mpmul (r(i3), r(i4), r(i3)) call mpmul (r(i3), r(i2), r(i2)) if ((r(i2).eq.0).or.(r(i2+1).le.(-r(sv)))) go to 70 60 call mpadd (y, r(i2), y) c restore everything and return 70 call mpresn (sv) return end subroutine mplngs (x, y) c returns y = ln(x) for mp x and y, using the gauss-salamin c algorithm based on the arithmetic-geometric mean iteration c (see analytic computational complexity (ed. by j. f. traub), c academic press, 1976, 151-176) unless x is close to 1. c time = o(log(t)m(t)) + o(t**2) if abs(x-1) .ge. 1/b, c and as for mplns otherwise. c slower than mpln unless t is large (.ge. about 500) so c mainly useful for testing purposes. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer i, it, i2, i3, i4, n, sv, tg, t2 integer b, dummy(17), lun, m, mxr, r(1), sptr, t, x(1), y(1) integer mpchgb, mpget, mptlb c save t etc. call mpsavn (sv) c check for zero argument. if (x(1) .le. 0) call mperrm ( \$ 27hx .le. 0 in call to mplngs\$) c see if x close to 1. call mpnew (i2) call mpaddi (x, -1, r(i2)) if ((r(i2).ne.0).and.(r(i2+1).ge.0)) go to 10 c here abs(x-1) .lt. 1/b so gauss-salamin algorithm could be c inaccurate because of cancellation. the precision could be c increased to compensate for this, but simpler to use mplns. call mplns (r(i2), y) go to 30 c increase working precision, allocate space etc. 10 sptr = i2 call mpgd3 (mptlb (iabs(t)), tg) call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpsetr (0) call mpmove (x, iabs(r(sv)), r(i2), tg) t2 = (t+1)/2 c modify exponent r(i2+1) = -t2 call mpmuli (r(i2), 4, r(i2)) call mpcim (1, r(i3)) c compute number of iterations required. n = max0 (1, mpchgb (2, mptlb (t2+1) * (mptlb (iabs(t)) \$ + 3), 1) - 3) c arithmetic-geometric mean loop do 20 i = 1, n call mpadd (r(i2), r(i3), r(i4)) call mpdivi (r(i4), 2, r(i4)) call mpmul (r(i2), r(i3), r(i3)) call mpsqrt (r(i3), r(i2)) c faster to exchange pointers than mp numbers it = i3 i3 = i4 20 i4 = it c check that convergence occurred call mpsub (r(i2), r(i3), r(i4)) if ((r(i4) .ne. 0) .and. ((r(i2+1) - r(i4+1)) .lt. (t-3))) \$ call mperrm (39hiteration failed to converge in mplngs\$) c could save some time by precomputing pi and ln(b) call mppi (r(i4)) call mpdiv (r(i4), r(i3), r(i3)) call mpdivi (r(i3), 2, r(i3)) call mplni (iabs(b), r(i4)) call mpmuli (r(i4), mpget(x, 2) + t2, r(i4)) c allow for modified exponent call mpsub (r(i4), r(i3), r(i3)) c round result call mpsetr (iabs(r(sv+2))) call mprnd (r(i3), tg, y, iabs(r(sv)), 1) c restore t etc. and return. 30 call mpresn (sv) return end subroutine mplni (n, x) c returns multiple-precision x = ln(n) for small positive c integer n (2*n must be representable). time is o(t**2). c method is to use a rapidly converging series and mpl235. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, dummy integer i, ia, ik, im, ip, ip2, iq, iq2, i2, i3, i4, j, mxint, n1, \$ n2, sv, tg, b, dummy(21), mpparn, mptlb, n, r(1), t, x(1) c check for n = 1 and n .lt. 1 if (n-1) 10, 20, 30 10 call mperrm (32hn not positive in call to mplni\$) c ln(1) = 0 20 x(1) = 0 return c here n .ge. 2, increase t as necessary. 30 call mpsavn (sv) call mpgd3 (mptlb (iabs(t)), tg) c use truncated arithmetic internally call mpsetr (0) mxint = mpparn (16) c allocate working space call mpnew (i2) if (n.gt.2) go to 40 c n = 2 is a special case call mpl235 (1, 0, 0, r(i2)) go to 150 c here n .ge. 3 40 if (n .gt. (mxint/2)) go to 80 j = 3 ia = 0 n2 = n/2 50 if (j.gt.n2) go to 60 ia = ia + 1 j = 2*j go to 50 c now j = 3*(2**ia) .le. n .lt. 6*(2**ia) 60 j = j/3 im = n ik = 0 do 70 i = 3, 6 n1 = i*j if (iabs(n1-n).gt.im) go to 70 im = iabs(n1-n) ik = i 70 continue n1 = ik*j c now n is close to n1 = ik*(2**ia) c and ik = 3, 4, 5 or 6, so mpl235 gives ln(n1). if (ik.eq.3) call mpl235 (ia, 1, 0, r(i2)) if (ik.eq.4) call mpl235 (ia+2, 0, 0, r(i2)) if (ik.eq.5) call mpl235 (ia, 0, 1, r(i2)) if (ik.eq.6) call mpl235 (ia+1, 1, 0, r(i2)) if (n.eq.n1) go to 150 c now need ln(n/n1). n2 = n call mpgcd (n2, n1) ip = n2 - n1 iq = n2 + n1 c check for possible integer overflow if (iq.gt.14) go to 100 80 call mperrm (29hn too large in call to mplni\$) c reduce to lowest terms 100 call mpgcd (ip, iq) call mpnew (i3) call mpnew (i4) call mpcqm (2*ip, iq, r(i4)) call mpstr (r(i4), r(i3)) call mpadd (r(i2), r(i3), r(i2)) i = 1 if (iq.gt.(mxint/iq)) go to 110 iq2 = iq**2 ip2 = ip**2 c loop to sum series for ln(n2/n1) 110 i = i + 2 if (r(i4).eq.0) go to 140 c reduce t if possible, done if can reduce below 2 t = tg + r(i4+1) + 2 if (t.le.2) go to 140 t = min0 (t, tg) c split up call to mpmulq if iq too large if (iq.gt.(mxint/iq)) go to 120 call mpmulq (r(i4), ip2, iq2, r(i4)) go to 130 c here iq too large for one call to mpmulq 120 call mpmulq (r(i4), ip, iq, r(i4)) call mpmulq (r(i4), ip, iq, r(i4)) 130 call mpdivi (r(i4), i, r(i3)) c restore t and accumulate sum t = tg call mpadd (r(i3), r(i2), r(i2)) go to 110 140 t = tg c restore rndrl and round result 150 call mpsetr (iabs(r(sv+2))) call mprnd (r(i2), tg, x, iabs(r(sv)), 1) c restore t etc. call mpresn (sv) return end subroutine mplns (x, y) c returns mp y = ln(1+x) if x is an mp number satisfying the c condition abs(x) .lt. 1/b, error otherwise. c x may be packed or unpacked, y is unpacked. c uses newtons method to solve the equation c exp1(-y) = x, then reverses sign of y. c (here exp1(y) = exp(y) - 1 is computed using mpexp1). c time is o(sqrt(t).m(t)) as for mpexp1. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, dummy integer it0, i2, i3, i4, i5, sv, tg, ts2, ts3 integer b, dummy(21), mpgd, mpget, r(1), t, x(1), y(1) c save t etc. then increase, use truncation internally. call mpsavn (sv) t = max0 (5, t + 1 + mpgd (100)) tg = t call mpsetr (0) c allocate temporary space. call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpnew (i5) c check for x = 0 exactly if (x(1).ne.0) go to 10 y(1) = 0 go to 50 c check that abs(x) .lt. 1/b 10 if (mpget (x, 2) .ge. 0) call mperrm ( \$ 33habs(x) .ge. 1/b in call to mplns\$) c get starting approximation to -ln(1+x). call mpmove (x, iabs(r(sv)), r(i3), tg) call mpstr (r(i3), r(i5)) call mpdivi (r(i5), 4, r(i2)) call mpaddq (r(i2), -1, 3, r(i2)) call mpmul (r(i5), r(i2), r(i2)) call mpaddq (r(i2), 1, 2, r(i2)) call mpmul (r(i5), r(i2), r(i2)) call mpaddi (r(i2), -1, r(i2)) call mpmul (r(i5), r(i2), r(i5)) c start newton iteration using small t, later increase t = max0 (5, 13 - 2*b) it0 = (t + 5)/2 20 call mpexp1 (r(i5), r(i4)) call mpmul (r(i3), r(i4), r(i2)) call mpadd (r(i4), r(i2), r(i4)) call mpadd (r(i3), r(i4), r(i4)) call mpsub (r(i5), r(i4), r(i5)) if (t.ge.tg) go to 40 c following loop computes next value of t to use. c because newtons method has 2nd order convergence, c we can almost double t each time. ts3 = t t = tg 30 ts2 = t t = (t+it0)/2 if (t.gt.ts3) go to 30 t = ts2 go to 20 c check that newton iteration was converging as expected 40 if ((r(i4).ne.0).and.((2*r(i4+1)).gt.(it0-t))) \$ call mperrm (41hnewton iteration not converging in mplns\$) c reverse sign and round result. r(i5) = -r(i5) call mpsetr (iabs(r(sv+2))) call mprnd (r(i5), tg, y, iabs(r(sv)), 1) c restore stack pointer etc. and return 50 call mpresn (sv) return end logical function mplt (x, y) c returns logical value of (x .lt. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mplt = (mpcomp(x,y) .lt. 0) return end subroutine mpmax (x, y, z) c sets z = max (x, y) where x, y and z are mp variables. c x and/or y may be packed or unpacked, z is unpacked. integer mpcomp, x(1), y(1), z(1) if (mpcomp (x, y) .ge. 0) go to 10 c here x .lt. y call mpunpk (y, z) return c here x .ge. y 10 call mpunpk (x, z) return end subroutine mpmaxr (x) c sets x to the largest possible positive mp number common /mpcom/ b, t, m, dummy integer b, dummy(20), i, it, m, t, two, x(1) two = 2 c check legality of parameters in common /mpcom/ call mpchk it = b - 1 c set fraction digits to b-1 do 10 i = 1, t 10 x(i+2) = it c set sign and exponent, and update mxexpn. x(1) = 1 x(two) = m call mpupdt (x(two)) return end integer function mpmexa (x) c returns the maximum allowable exponent of mp numbers (the third c word of common /mpcom/). x is a dummy mp argument. integer mpparn, x(1) mpmexa = mpparn (3) return end subroutine mpmexb (i, x) c sets the maximum allowable exponent of mp numbers (i.e. the c third word of common /mpcom/) to i. c i should be greater than 4*t, and 4*i should be representable c as a single-precision integer. c x is a dummy mp argument (augment expects one). integer i, x(1) call mpparc (3, i) return end subroutine mpmin (x, y, z) c sets z = min (x, y) where x, y and z are mp variables. c x and/or y may be packed or unpacked, z is unpacked. integer mpcomp, x(1), y(1), z(1) if (mpcomp (x, y) .ge. 0) go to 10 c here x .lt. y call mpunpk (x, z) return c here x .ge. y 10 call mpunpk (y, z) return end subroutine mpminr (x) c sets x to the smallest positive normalized mp number b**(-m) common /mpcom/ b, t, m, dummy integer b, dummy(20), i, il, m, t, two, x(1) two = 2 c clear x il = t + 2 do 10 i = 1, il 10 x(i) = 0 c set sign, exponent and first digit x(1) = 1 x(two) = 1 - m x(two+1) = 1 c update mnexpn and return call mpupdt (x(two)) return end subroutine mpmlp (u, v, w, j) c performs inner multiplication loop for mpmul c note that carries are not propagated in inner loop, c which saves time at the expense of space. integer i, j, u(1), v(1), w do 10 i = 1, j 10 u(i) = u(i) + w*v(i) return end subroutine mpmod (x, y, z) c c returns z = x - int(x/y)*y for mp x, y and z, y nonzero. c here int means integer part, truncated towards zero. c returns zero if y is zero. c x and/or y may be packed or unpacked, z is unpacked. c c warning - time is o(max (1, ln(abs(x/y))) * (t**2)), which may be c ******* large if abs(x) is much larger than abs(y). c common r common /mpcom/ b, t, m, dummy integer i2, i3, i4, i5, mpcmpa, sv, tg integer b, dummy(20), m, r(1), t, x(1), y(1), z(1) c check for y zero. if (y(1).ne.0) go to 10 z(1) = 0 return c save t etc. 10 call mpsavn (sv) c use truncation (actually exact arithmetic here) call mpsetr (0) c double t so multiplication is exact t = 2*t tg = t c allocate working space call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpnew (i5) c move x and y, increase m to avoid over/underflow call mpmove (x, iabs(r(sv)), r(i2), tg) call mpmove (y, iabs(r(sv)), r(i4), tg) call mpstr (r(i4), r(i5)) m = 2*m + 2 c jump to end of loop (simulate while loop). go to 30 c loop here until z result less (in absolute value) than y 20 t = r(sv) c use mpdivl to be sure result truncated call mpdivl (r(i2), r(i5), r(i3)) call mpcmim (r(i3), r(i3)) call mpmove (r(i3), iabs(r(sv)), r(i3), tg) c increase t so following multiplication and subtraction exact t = tg call mpmul (r(i3), r(i4), r(i3)) call mpsub (r(i2), r(i3), r(i2)) c loop termination check here 30 if (mpcmpa (r(i2), r(i4)) .ge. 0) go to 20 c move result (no digits lost), restore everything and return. call mpmove (r(i2), tg, z, iabs(r(sv))) call mpresn (sv) return end subroutine mpmove (x, tx, y, ty) c assumes mp number x of tx .ge. 2 digits, c mp number y of ty .ge. 2 digits (y not necessarily initialized). c moves x to y, either padding with trailing zeros or rounding c (as specified by rndrl - see subroutine mpnzr) depending on c whether tx is less or greater than ty. the value of t in c common /mpcom/ is restored on return. c x may be in packed or unpacked format, y is unpacked. common r common /mpcom/ b, t, dummy integer ex, i, ih, il, iy, i2, kg, sn, sv integer b, dummy(21), r(1), t, tx, ty, x(1), y(1) call mpsavn (sv) c use local variables in case t actual argument kg = tx - ty iy = ty t = tx c compare tx and ty. if (kg.gt.0) go to 20 c here tx .le. ty call mpunpk (x, y) if (kg.eq.0) go to 40 c here tx .lt. ty. il = t + 3 ih = iy + 2 c pad with trailing zeros. do 10 i = il, ih 10 y(i) = 0 go to 40 c here tx. gt. ty. check for x zero. 20 if (x(1).ne.0) go to 30 c here x is zero so set y to zero and return y(1) = 0 go to 40 c move x to temporary storage and round 30 call mpnew (i2) call mpunpk (x, r(i2)) c reduce t for call to mpnzr. t = iy sn = r(i2) ex = r(i2+1) call mpnzr (sn, ex, y, r(i2+2), kg) c restore stack pointer etc. 40 call mpresn (sv) return end subroutine mpmul (x, y, z) c c multiplies x and y, returning result in z, for mp x, y and z. c the simple o(t**2) algorithm is used, with rounding defined by the c parameter rndrl in common /mpcom/, as follows - c rndrl = 0 - truncate result towards zero, c error less than 1.01 units in the last place. c rndrl = 1, 2 or 3 - see subroutine mpnzr. c c advantage is taken of zero digits in x, but not in y. c asymptotically faster algorithms are known (see knuth, c vol. 2), but are difficult to implement in fortran in an c efficient and machine-independent manner. c in comments to other mp routines, m(t) is the time c to perform t-digit mp multiplication. thus c m(t) = o(t**2) with the present version of mpmul, c but m(t) = o(t.log(t).log(log(t))) is theoretically possible. c c mp over/underflow is detected by subroutine mpnzr. c common r common /mpcom/ b, t, dummy integer c, i, ih, ip, i2, j, j1, kg, re, ri, rs, tg, two, \$ b, dummy(21), next, r(1), sv, t, xi, x(1), y(1), z(1) two = 2 c form sign of product rs = x(1)*y(1) if (rs.ne.0) go to 10 c set result to zero z(1) = 0 return c compute number of guard digits required 10 call mpsavn (sv) kg = min0 (3, t) if ((r(sv+2) .ne. 0) .or. (t .gt. ((b*b)/100))) kg = t c allocate space for accumulator, check parameters etc. tg = t + kg call mpnew2 (i2, tg) next = i2 + tg ih = next - 1 c form exponent of product re = x(two) + y(two) c clear accumulator do 20 i = i2, ih 20 r(i) = 0 c perform multiplication c = 8 do 40 i = 1, t xi = x(i+2) c for speed, put the number with many zeros first if (xi.eq.0) go to 40 ip = i + i2 call mpmlp (r(ip), y(two+1), xi, min0 (t, tg - i)) c = c - 1 if (c.gt.0) go to 40 c check for legal base b digit if ((xi.lt.0).or.(xi.ge.b)) go to 80 c propagate carries at end and every eighth time, c faster than doing it every time. do 30 j = 1, tg j1 = next - j ri = r(j1) + c if (ri.lt.0) go to 70 c = ri/b 30 r(j1) = ri - b*c if (c.ne.0) go to 80 c = 8 40 continue if (c.eq.8) go to 60 if ((xi.lt.0).or.(xi.ge.b)) go to 80 c = 0 do 50 j = 1, tg j1 = next - j ri = r(j1) + c if (ri.lt.0) go to 70 c = ri/b 50 r(j1) = ri - b*c if (c.ne.0) go to 80 c normalize and round result 60 call mpnzr (rs, re, z, r(i2), kg) c restore stack pointer and return call mpresn (sv) return 70 call mperrm (39hinteger overflow in mpmul, b too large\$) 80 call mperrm (38hillegal base b digit in call to mpmul\$) return end subroutine mpmuli (x, iy, z) c multiplies mp x by single-precision integer iy giving mp z. c rounding defined by parameter rndrl in common /mpcom/ - c see subroutine mpnzr. c multiplication by 1 may be used to normalize a number c even if last digit is greater than b-1. common r common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer c, c1, c2, i, ij, ikg, ip, is, ix, i2, j, j1, j2, kg, re integer ri, rs, tg, two, t1 integer b, dummy(17), iy, lun, m, mxr, r(1), sptr, t, x(1), z(1) integer mpgd, mpparn two = 2 c compute number of guard digits required to give exact product. kg = mpgd (iy) c if iy is 1 last digit may be large (see e.g. mpcim). if (iy.eq.1) kg = 1 c allocate t + kg words for accumulator tg = t + kg call mpnew2 (i2, tg) c check for x zero rs = x(1) if (rs.eq.0) go to 110 c check for multiplier zero. j = iy if (j) 10, 110, 20 c here adjust sign for negative multiplier. 10 j = -j rs = -rs c check for multiplication by +-b 20 if (j.ne.b) go to 30 c check for overflow. if (x(two) .ge. m) call mpovfl (z) c here multiplication by +-b, result does not overflow. call mpstr (x, z) z(1) = rs z(two) = z(two) + 1 call mpupdt (z(two)) go to 120 c have to do genuine multiplication here. c set exponent to exponent(x) + number of guard digits. 30 re = x(two) + kg c form product in accumulator c = 0 t1 = t + 1 c if j*b not representable as an integer we have to simulate c double-precision multiplication. if (j.gt.(mpparn(16)/b)) go to 80 c inner loop for small abs(iy) case. do 40 ij = 1, t i = t1 - ij ri = j*x(i+2) + c c = ri/b ip = sptr - ij 40 r(ip) = ri - b*c c check for integer overflow if (ri.lt.0) go to 100 c have to treat first kg words of accumulator separately ikg = i2 + kg if (kg.le.0) go to 60 do 50 ij = 1, kg i = ikg - ij ri = c c = ri/b 50 r(i) = ri - b*c 60 if (c.ne.0) go to 100 c normalize and round or truncate result 70 call mpnzr (rs, re, z, r(i2), kg) go to 120 c here j may be too large for single-precision multiplication 80 j1 = j/b j2 = j - j1*b c form product do 90 ij = 1, tg c1 = c/b c2 = c - b*c1 i = t1 - ij ix = 0 if (i.gt.0) ix = x(i+2) ri = j2*ix + c2 is = ri/b c = j1*ix + c1 + is ip = sptr - ij 90 r(ip) = ri - b*is if (c.eq.0) go to 70 c can only get here if integer overflow occurred or b**t too small. 100 call mperrm (25herror occurred in mpmuli\$) c set result to zero 110 z(1) = 0 c restore stack pointer 120 sptr = i2 return end subroutine mpmulq (x, i, j, y) c multiplies mp x by i/j, giving mp result y. c rounding defined by parameter rndrl in common /mpcom/ as follows - c rndrl = 0 - as for mpmuli followed by mpdivi with rndrl = 0. c see comments in mpdivi, mpmuli and mpnzr. error is c less than b units in the last place. c rndrl = 1, 2, or 3 - see comments in subroutine mpnzr. c overflow/underflow may occur in mpmuli/mpdivi respectively. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer is, i2, js, mpgd, tg, ts integer b, dummy(12), i, j, lun, m, mnexpn, mnsptr, mxexpn, mxr integer mxsptr, r(1), rndrl, sptr, t, x(1), y(1) c check for division by zero. if (j .eq. 0) call mperrm ( \$ 37hattempted division by zero in mpmulq\$) c check for multiplication by zero. if (i.ne.0) go to 10 y(1) = 0 return c reduce to lowest terms 10 is = i js = j call mpgcd (is, js) if (iabs(is).eq.1) go to 30 c check rounding options if (rndrl.gt.0) go to 20 c here rndrl = 0, so use simple method. call mpmuli (x, is, y) call mpdivi (y, js, y) return c here rndrl. gt. 0, so use enough guard digits that multiplication c is exact. 20 ts = t t = t + mpgd (is) tg = t call mpnew (i2) call mpmove (x, ts, r(i2), tg) call mpmuli (r(i2), is, r(i2)) call mpdivi (r(i2), js, r(i2)) c round result call mpmove (r(i2), tg, y, ts) c restore t and stack pointer t = ts sptr = i2 return c here is = +-1 30 call mpdivi (x, is*js, y) return end subroutine mpmuls (x, i, j, k, l) c sets x = x*i*j/(k*l) for mp x, integer i, j, k, l, c k*l nonzero. calls mpmulq once if i*j and k*l not too large, c otherwise mpmulq called twice. c rounding not best possible, but directed rounding c is in the correct direction. integer ij, is, kl, ks, mpparn, mxint integer i, j, k, l, x(1) if ((j.eq.0).or.(l.eq.0)) go to 30 c see if i*j or k*l might overflow. mxint = mpparn(16) if ((iabs(i) .gt. (mxint/iabs(j))) .or. \$ (iabs(k) .gt. (mxint/iabs(l)))) go to 30 c here safe to form i*j and k*l. ij = i*j kl = k*l c worthwhile to check for special cases here. if (ij.ne.1) go to 10 call mpdivi (x, kl, x) return 10 if (kl.ne.1) go to 20 call mpmuli (x, ij, x) return 20 call mpmulq (x, ij, kl, x) return c here split up into two calls, taking care for directed rounding. 30 is = i if (j.lt.0) is = -is ks = k if (l.lt.0) ks = -ks call mpmulq (x, is, ks, x) call mpmulq (x, iabs(j), iabs(l), x) return end logical function mpne (x, y) c returns logical value of (x .ne. y) for mp x and y. c x and/or y may be packed or unpacked. integer mpcomp, x(1), y(1) mpne = (mpcomp(x,y) .ne. 0) return end subroutine mpneg (x, y) c sets y = -x for mp numbers x and y c y will be packed if x is. integer x(1), y(1) call mpstr (x, y) y(1) = -y(1) return end subroutine mpnew (i) c returns index i such that words i, ... , i+t+1 of blank common c are available for use, updates stack pointer sptr etc. common /mpcom/ b, t, dummy integer b, dummy(21), i, t call mpnew2 (i, t+2) return end subroutine mpnew2 (i, j) c returns i such that words i, ... , i+abs(j)-1 of blank common c are available for use, updates stack pointer sptr etc. common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, dummy integer b, dummy(15), i, j, lun, m, mxr, mxsptr, mnsptr, sptr, t c allocate space on stack i = sptr sptr = i + iabs(j) c check for stack overflow and update mxsptr. c useful to check validity of parameters in common /mpcom/ c if top level call, too. if ((sptr.gt.mxsptr) .or. (mnsptr.ge.i)) call mpchk return end subroutine mpnzr (rs, expn, z, a, kg) c c assumes (t+kg)-digit fraction in a(1), ... , a(t+kg), c kg .ge. 0, sign = rs, exponent = expn. c a(1) and z(3) may have the same address. c normalizes and returns mp result in z. a is destroyed. c c rounding depends on value of rndrl (in common /mpcom/) as follows - c c 0 - round towards zero (i.e. chop or truncate). the resulting c error is less than 1 ulp (unit in the last place). c this is the default (if mpset used). it is the fastest option c and is satisfactory for most purposes. c note - rndrl = 0 is assumed if kg is zero. c c 1 - round to nearest, preferring even last digit if tie. c must have kg.gt.0 to be meaningful. error .le. 1/2 ulp if c kg (the number of guard digits) is large enough. c this option gives worst-case error bounds half as large as c for rndrl = 0, and average-case error considerably better c than for rndrl = 0. however, if rndrl = 0 does not give c sufficient accuracy is is usually preferable to increase t c rather than use rndrl = 1. c c 2 - round down (towards -infinity). (must have kg .gt. 0) c c 3 - round up (towards +infinity). (must have kg .gt. 0) c these last two options are useful for interval arithmetic as c they give rigorous lower and upper bounds. c c note - comments in other routines often refer to the above c definitions of rounding options. the comments regarding c kg above refer only to this subroutine. c common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b2, i, inm, inp, is, j, k, re, tg, two, t1 integer a(1), b, dummy(12), expn, kg, lun, m, mnexpn, mnsptr integer mxexpn, mxr, mxsptr, rndrl, rs, sptr, t, z(1) two = 2 c check that sign = +-1 c the following message is usually caused by exceeding c array bounds. check that arrays used for mp numbers have c dimension at least t+2. if (iabs(rs) .gt. 1) call mperrm ( \$ 38hsign not 0, +1 or -1 in call to mpnzr\$) c check for zero if (rs .ne. 0) go to 20 c return zero. 10 z(1) = 0 return c save exponent locally 20 re = expn t1 = t + 1 tg = t + kg c look for first nonzero digit do 30 i = 1, tg is = i - 1 if (a(i).gt.0) go to 40 30 continue c fraction zero go to 10 c see if number is normalized. 40 if (is.eq.0) go to 70 c normalize by left shifting and padding on right with zeros. re = re - is inm = tg - is do 50 j = 1, inm k = j + is 50 a(j) = a(k) inp = inm + 1 do 60 j = inp, tg 60 a(j) = 0 c check various rounding options described above. 70 if ((rndrl.eq.0).or.(kg.eq.0).or. \$ ((rndrl.eq.2).and.(rs.gt.0)).or. \$ ((rndrl.eq.3).and.(rs.lt.0))) go to 160 if (rndrl.eq.1) go to 90 c here round away from zero if any of a(t+1), ... , a(t+kg) nonzero. do 80 i = t1, tg if (a(i).ne.0) go to 120 80 continue c here a(t+1) = ... = a(t+kg) = 0 so can truncate. go to 160 c round to nearest here. need to treat odd and even bases differently. 90 b2 = b/2 if ((2*b2).ne.b) go to 140 c b even. round away from zero if a(t+1).ge.b2 unless c a(t) even, a(t+1) = b2, and a(t+2) = ... = a(t+kg) = 0. if (a(t1) - b2) 160, 100, 120 100 if (mod(a(t),2).ne.0) go to 120 i = t1 110 i = i + 1 if (i.gt.tg) go to 160 if (a(i).eq.0) go to 110 c round away from zero 120 do 130 j = 1, t i = t1 - j a(i) = a(i) + 1 if (a(i).lt.b) go to 160 130 a(i) = 0 c exceptional case, rounded up to .10000... re = re + 1 a(1) = 1 go to 160 c odd base, round away from zero if .a(t+1) ... a(t+kg) .gt. 1/2 140 do 150 i = t1, tg if (a(i) - b2) 160, 150, 120 150 continue c check for overflow 160 if (re .gt. m) call mpovfl (z) c check for underflow if (re.le.(-m)) go to 180 c store result in z z(1) = rs z(two) = re do 170 i = 1, t 170 z(i+2) = a(i) c update maximum and minimum exponent indicators mxexpn = max0 (mxexpn, re) mnexpn = min0 (mnexpn, re) return c underflow here 180 call mpunfl (z) return end subroutine mpout (x, c, p, n) c c converts multiple-precision x to fp.n format in c, c which may be printed under pa1 format. note that c n = -1 is allowed, and effectively gives ip format. c digits after the decimal point are blanked out if c they could not be significant. c efficiency is higher if b is a power of 10 than if not. c x may be packed or unpacked. c dimension of c must be at least p. c c is an integer array, p and n are integers. c c rounding of decimal output is determined by the parameter c rndrl in common /mpcom/, as follows - c c call mpsetr (0) or 1 - round to approximatlely nearest decimal number c representable in the specified format. c error in conversion to decimal is less c than 0.6 units in the last decimal place. c rndrl = 2 - round down (towards -infinity) to a decimal c number representable in the specified format. c rndrl = 3 - round up (towards +infinity) to a decimal number c representable in the specified format. c c output base is usually 10 (default set by mpset), but c may be changed to any of 2, ... , 16 by changing the c parameter outbas in common /mpcom/ (see mpparn). c common r common /mpcom/ b, t, dummy integer i, ib, ip, ir, is, isz, itp, iz, i1, i2, i3, j, jd, \$ jp, nmax, np, outbas, sv, tg, tp, \$ b, c(1), dummy(21), n, p, r(1), t, x(1), \$ mpchgb, mpdigw, mpgd, mpparn, \$ blank, period, minus, star data blank, period, minus, star /1h , 1h., 1h-, 1h*/ c check b, t, outbas etc. for legality. call mpchk c save t etc. call mpsavn (sv) c get output base. outbas = mpparn (20) c check legality of parameters n and p. if ((n .le. (-1)) .or. (p .le. n) .or. (p .le. 0)) \$ call mperrm (36hn and/or p illegal in call to mpout\$) c compute displacements. np = p - n ip = np - 1 c compute power of outbas which we can safely multiply and c divide by (this saves time). ib = mpparn(16)/outbas tp = outbas itp = 1 10 if ((tp .gt. ib).or.(tp .eq. b)) go to 20 tp = tp*outbas itp = itp + 1 go to 10 c compute guard digits required, move x etc. 20 t = t + 1 + mpgd (max0 (50, tp)) tg = t call mpnew (i2) call mpnew (i3) call mpmove (x, iabs(r(sv)), r(i2), tg) c use truncation (usually) for internal arithmetic call mpsetr (0) c put formatted zero in c do 30 i = 1, p c(i) = blank if (i .ge. ip) c(i) = mpdigw(0) if (i .eq. np) c(i) = period 30 continue c get sign of x, check for zero is = r(i2) if (is .eq. 0) go to 160 r(i2) = 1 c compute maximum number of nonzero digits which we can c meaningfully give after decimal point. nmax = min0 (n, mpchgb (outbas, iabs(b), min0 (4*n, \$ max0 (0, r(sv) - r(i2+1))))) c compute rounding constant for decimal output ir = 1 if (r(sv+2) .ge. 2) ir = 2*(r(sv+2)-2) if (is.lt.0) ir = 2 - ir c change rounding option if necessary for directed rounding if (ir .eq. 2) call mpsetr (3) call mpcqm (ir, 2, r(i3)) if (nmax .gt. 0) call mpscal (r(i3), outbas, (-nmax)) c add rounding constant to abs(x). call mpadd (r(i2), r(i3), r(i2)) c ip places before point, so divide by outbas**ip if (ip .gt. 0) call mpscal (r(i2), outbas, (-ip)) iz = 0 c round result to r(sv)+1 digits for directed roundings (so subsequent c arithmetic exact). if (r(sv+2) .le. 1) go to 40 call mpmove (r(i2), tg, r(i2), r(sv)+1) call mpmove (r(i2), r(sv)+1, r(i2), tg) c now use chopped arithmetic 40 call mpsetr (0) c check that number is less than one if (r(i2+1) .gt. 0) go to 120 if (ip .le. 0) go to 90 c put digits before point in jd = 1 do 80 i = 1, ip if (jd .gt. 1) go to 70 if ((i+itp) .le. (ip+1)) go to 50 c multiply by outbas, truncating result call mpmuli (r(i2), outbas, r(i2)) jd = outbas go to 60 c here we can multiply by a power of outbas to save time 50 call mpmuli (r(i2), tp, r(i2)) jd = tp c get integer part 60 call mpcmi (r(i2), jp) c and fractional part call mpcmf (r(i2), r(i2)) 70 jd = jd/outbas c get next decimal digit j = jp/jd jp = jp - j*jd isz = iz if ((j .gt. 0).or.(i .eq. ip)) iz = 1 if (iz .gt. 0) c(i) = mpdigw(j) if ((iz .eq. isz).or.(is .gt. 0)) go to 80 if (i .eq. 1) go to 120 c(i-1) = minus 80 continue 90 if (nmax .le. 0) go to 140 c put in digits after decimal point jd = 1 do 110 i = 1, nmax if (jd .gt. 1) go to 100 call mpmuli (r(i2), tp, r(i2)) call mpcmi (r(i2), jp) call mpcmf (r(i2), r(i2)) jd = tp 100 jd = jd/outbas j = jp/jd jp = jp - j*jd i1 = np + i 110 c(i1) = mpdigw(j) go to 140 c error occurred, return asterisks. 120 do 130 i = 1, p 130 c(i) = star go to 160 c blank out any nonsignificant trailing zeros 140 if (nmax .ge. n) go to 160 i1 = np + nmax + 1 do 150 i = i1, p 150 c(i) = blank c restore stack pointer etc. and return. 160 call mpresn (sv) return end subroutine mpout2 (x, c, p, n, outbas) c this routine is redundant, but is included for compatibility c with earlier versions of the mp package. integer x(1), c(1), p, n, outbas, outsav, mpparn c save and change output base before calling mpout. outsav = mpparn (20) call mpparc (20, outbas) call mpout (x, c, p, n) c restore output base call mpparc (20, outsav) return end subroutine mpoute (x, c, j, p) c c assumes x is an mp number and c an integer array of dimension at c least p .ge. 4. on return j is the exponent (to base outbas) of x c and the fraction is in c, ready to be printed in a1 format. c for example, we could print j and c in (i10, 1x, pa1) format. c the fraction has one place before decimal point and p-3 after. c j and p are integers. x may be packed or unpacked. c c rounding not best possible, but directed roundings (rndrl = 2 c or 3) give correct bounds - see comments in mpout. c c default output base is 10, but this may be changed - see c comments on outbas in mpparn. c common r integer c(1), isn, i2, j, p, r(1), sv, x(1), \$ blank, minus logical mpis data blank, minus /1h , 1h-/ call mpsavn (sv) if (p .lt. 4) call mperrm (27hp .lt. 4 in call to mpoute\$) c allocate temporary space. call mpnew (i2) call mpcmef (x, j, r(i2)) call mpout (r(i2), c, p, p-3) c see if output of mpout was rounded up to outbas if (mpis (c(1), blank) .or. mpis (c(1), minus)) go to 10 c it was, so add 1 to j and convert sign to mp j = j + 1 isn = r(i2) call mpcim (isn, r(i2)) call mpout (r(i2), c, p, p-3) c restore stack pointer etc and return. 10 call mpresn (sv) return end subroutine mpoutf (x, p, n, iform, err) c c writes mp number x on logical unit lun (fourth word of c common /mpcom/) in format iform, after converting to fp.n c representation (i.e. p characters including n after c the decimal point). the default output base is 10 (decimal), c but this may be changed by changing outbas in common /mpcom/. c iform should contain a format which allows for output of p c words in a1 format, plus any desired headings, spacing etc. c (hollerith constants in iform are nonstandard, but usually ok.) c err is returned as true if p not positive or if mpio detects c an error when writing on unit lun, otherwise false. c x may be packed or unpacked of mpoutf is called directly. c space required = p + o(t) words. c for rounding options see comments in subroutine mpout. c common r integer iform(1), i2, mpparn, n, p, r(1), sv, x(1) logical err err = .true. c return with error flag set if output field width p not positive if (p.le.0) return call mpsavn (sv) c allocate p words of stack space. call mpnew2 (i2, p) c convert x to base outbas call mpout (x, r(i2), p, n) c and write on unit lun with format iform call mpio (r(i2), p, mpparn(4), iform, err) c restore stack pointer and return. call mpresn (sv) return end subroutine mpovfl (x) c called on multiple-precision overflow, i.e. when the c exponent of mp number x would exceed m. c execution is terminated with an error message c after calling mpmaxr(x) and setting mxexpn = m + 1. common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, dummy integer b, dummy(14), lun, m, mnsptr, mxexpn, mxr, mxsptr, sptr, \$ t, x(1) c m may have been overwritten, so check b, t, m etc. call mpchk c set x to largest possible positive number call mpmaxr (x) c set mxexpn to m + 1 (largest valid exponent is m) mxexpn = m + 1 call mperrm (37hcall to mpovfl, mp overflow occurred\$) c mperrm does not return, so never get here. return end subroutine mppack (x, y) c assumes that x is an mp number stored as usual in an integer c array of dimension at least t+2, and y is an integer array c of dimension at least int((t+3)/2). c x is stored in a compact format in y, and may be retrieved c by calling mpunpk (y, x). c if x is already in compact (packed) format, effect is c the same as for mpstr (x, y). c mppack and mpunpk are useful if space is critical, for example c when working with large arrays of mp numbers. common /mpcom/ b, t, dummy integer b, t, x(1), y(1), dummy(21), i, j, x2, two c check whether x is in packed format or not. if (iabs(x(1)) .eq. 1) go to 10 c here x is zero or in packed format. call mpstr (x, y) return c here x is nonzero and in unpacked format. 10 j = t/2 two = 2 x2 = x(two) c now pack two digits of x in each word of y. do 20 i = 1, j 20 y(i+1) = b*x(2*i+1) + x(2*i+2) c fix up last digit if t odd. if (mod (t, 2) .ne. 0) y(j+2) = b*x(t+2) c move second word of y, insert sign and exponent. y(1) = x(1)*y(two) y(two) = x2 return end integer function mppara (a) c returns the word of common /mpcom/ corresponding to the hollerith c string a. for details see subroutine mpparm. integer a(1), n call mpparm (a, .false., n) mppara = n return end subroutine mpparb (i, a) c sets the word of common /mpcom/ corresponding to the hollerith string c a to the integer value i. for details see subroutine mpparm. integer a(1), i, j j = i call mpparm (a, .true., j) return end subroutine mpparc (n, j) c sets the n-th word of common /mpcom/ to j if 1 .le. n .le. 23, c error otherwise. faster but less mnemonic than mpparb. c for the meanings of words 1 to 23 of common /mpcom/, see c the comments in subroutine mpparm. for their default c settings see mpset. common /mpcom/ mpc integer j, mpc(23), n if ((n .le. 0) .or. (n .gt. 23)) call mperrm ( \$ 28hillegal n in call to mpparc\$) mpc(n) = j c check legality of t etc. if called from top level, c i.e. if sptr .le. mnsptr. if (mpc(6) .le. mpc(8)) call mpchk return end subroutine mpparm (holl, set, j) c c holl is a hollerith string, set a logical value, j an integer. c if the first three characters of holl agree with one of - c c (1) base (the mp base), c (2) numdig (the number of digits), c (3) maxexp (the maximum allowable exponent), c (4) lun (logical unit for output), c (5) mxr (dimension of blank common, used for stack), c (6) sptr (pointer to 1 word above top of stack), c (7) mxsptr (maximum value of sptr which has occurred), c (8) mnsptr (pointer to bottom of stack, i.e. minimum of sptr), c (9) mxexpn (maximum exponent which has occurred), c (10) mnexpn (minimum exponent which has occurred), c (11) rndrl (indicator for type of rounding - see mpnzr), c (12) ktunfl (number of underflows which have occurred), c (13) mxunfl (maximum underflows allowed + 1, see mpunfl), c (14) decpl (equivalent number of floating decimal places - c actually the value of the second argument in the last c call to mpset2, c (15) mt2 (size of arrays for mp numbers - actually value of mt2 c in last call to mpset2), c (16) mxint (maximum positive integer of form 2**k - 1, see mplarg), c (17) exwid (width of exponent field for mpfout), c (18) inrecl (input record length for mpfin, .le. 80), c (19) inbase (input base for mpin etc., 2 .le. inbase .le. 16), c (20) outbas (output base for mpout, mpoute, mpfout etc., c 2 .le. outbas .le. 16), c (21) expch (exponent character for mpfout, default is e), c (22) chword (number of characters per word, returned by mpupw), c (23) onescp (ones complement indicator - see mpset2 and mpupw), c c then the corresponding word in common /mpcom/ is either c set to j (if set = .true.) or returned in j (if set = .false.). c c mpset2 sets the above to reasonable default values. for details c see the comments in subroutine mpset2. c c the description above of base, ... , onescp applies to c mppara, mpparb, mpparc and mpparn as well as mpparm. c integer c(3), d1(23), d2(23), d3(23), holl(1), i, j, l, n, \$ mpparn logical mpis, set data d1 (1), d2 (1), d3 (1) /1hb, 1ha, 1hs/ data d1 (2), d2 (2), d3 (2) /1hn, 1hu, 1hm/ data d1 (3), d2 (3), d3 (3) /1hm, 1ha, 1hx/ data d1 (4), d2 (4), d3 (4) /1hl, 1hu, 1hn/ data d1 (5), d2 (5), d3 (5) /1hm, 1hx, 1hr/ data d1 (6), d2 (6), d3 (6) /1hs, 1hp, 1ht/ data d1 (7), d2 (7), d3 (7) /1hm, 1hx, 1hs/ data d1 (8), d2 (8), d3 (8) /1hm, 1hn, 1hs/ data d1 (9), d2 (9), d3 (9) /1hm, 1hx, 1he/ data d1 (10), d2 (10), d3 (10) /1hm, 1hn, 1he/ data d1 (11), d2 (11), d3 (11) /1hr, 1hn, 1hd/ data d1 (12), d2 (12), d3 (12) /1hk, 1ht, 1hu/ data d1 (13), d2 (13), d3 (13) /1hm, 1hx, 1hu/ data d1 (14), d2 (14), d3 (14) /1hd, 1he, 1hc/ data d1 (15), d2 (15), d3 (15) /1hm, 1ht, 1h2/ data d1 (16), d2 (16), d3 (16) /1hm, 1hx, 1hi/ data d1 (17), d2 (17), d3 (17) /1he, 1hx, 1hw/ data d1 (18), d2 (18), d3 (18) /1hi, 1hn, 1hr/ data d1 (19), d2 (19), d3 (19) /1hi, 1hn, 1hb/ data d1 (20), d2 (20), d3 (20) /1ho, 1hu, 1ht/ data d1 (21), d2 (21), d3 (21) /1he, 1hx, 1hp/ data d1 (22), d2 (22), d3 (22) /1hc, 1hh, 1hw/ data d1 (23), d2 (23), d3 (23) /1ho, 1hn, 1he/ c unpack first three characters of a call mpupk (holl, c, 3, l) if (l.ne.3) go to 20 c see if they match any of those expected. do 10 i = 1, 23 n = i if (mpis (c(1), d1(i)) .and. \$ mpis (c(2), d2(i)) .and. \$ mpis (c(3), d3(i))) go to 30 10 continue c here a has length less than 3 or no match found. 20 call mperrm (43hillegal hollerith string in call to mpparm\$) c here a match was found at the n-th value. 30 if (set) call mpparc (n, j) if (.not. set) j = mpparn (n) return end integer function mpparn (n) c returns the n-th word of common /mpcom/ if 1 .le. n .le. 23, c zero otherwise. faster but less mnemonic than mppara. c for the meanings of words 1 to 23 of common /mpcom/, see c the comments in subroutine mpparm. for their default c settings see mpset. common /mpcom/ mpc integer mpc(23), n mpparn = 0 if ((n .gt. 0) .and. (n .le. 23)) mpparn = mpc(n) return end subroutine mppi (x) c sets mp x = pi to the available precision. c uses pi/4 = 4.arctan(1/5) - arctan(1/239), time = o(t**2). c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, dummy integer i2, i3, sv, tg integer b, dummy(21), mpgd, r(1), t, x(1) c save t etc. call mpsavn (sv) c see if any guard digits required if ((r(sv+2).ne.0).or.(b.le.100)) go to 10 c here no guard digits necessary. c allocate temporary space. call mpnew (i2) call mpart1 (5, r(i2)) call mpmuli (r(i2), 4, r(i2)) call mpart1 (239, x) call mpsub (r(i2), x, x) c use rounded multiplication by 4 call mpsetr (1) call mpmuli (x, 4, x) go to 20 c here need to use some guard digits 10 t = t + mpgd (100) tg = t call mpnew (i2) call mpnew (i3) call mpart1 (5, r(i2)) call mpmuli (r(i2), 4, r(i2)) c take care if directed rounding required. call mprevr (1) call mpart1 (239, r(i3)) call mpsetr (iabs(r(sv+2))) call mpsub (r(i2), r(i3), r(i3)) call mpmuli (r(i3), 4, r(i3)) c round result call mpmove (r(i3), tg, x, iabs(r(sv))) c restore everything 20 call mpresn (sv) return end subroutine mppigl (pi) c sets mp pi = 3.14159... to the available precision. c uses the gauss-legendre algorithm. c this method requires time o(ln(t)m(t)), so it is slower c than mppi if m(t) = o(t**2), but would be faster for c large t if a faster multiplication algorithm were used c (see comments in mpmul). c for a description of the method, see - multiple-precision c zero-finding and the complexity of elementary function c evaluation (by r. p. brent), in analytic computational c complexity (edited by j. f. traub), academic press, 1976, 151-176. c rounding options not implemented, no guard digits used. common r integer ix, i2, i3, i4, pi(1), r(1), sv c save t etc. call mpsavn (sv) c use truncated arithmetic. call mpsetr (0) c allocate working space. call mpnew (i2) call mpnew (i3) call mpnew (i4) call mpcim (1, pi(1)) call mpcqm (1, 2, r(i4)) call mpsqrt (r(i4), r(i4)) call mpcqm (1, 4, r(i3)) ix = 1 10 call mpstr (pi(1), r(i2)) call mpadd (pi(1), r(i4), pi(1)) call mpdivi (pi(1), 2, pi(1)) call mpmul (r(i2), r(i4), r(i4)) call mpsub (pi(1), r(i2), r(i2)) call mpmul (r(i2), r(i2), r(i2)) call mpmuli (r(i2), ix, r(i2)) call mpsub (r(i3), r(i2), r(i3)) call mpsqrt (r(i4), r(i4)) ix = 2*ix c check for convergence if ((r(i2).ne.0).and.(r(i2+1).ge.(-r(sv)))) go to 10 call mpmul (pi(1), r(i4), pi(1)) call mpdiv (pi(1), r(i3), pi(1)) call mpresn (sv) return end subroutine mppoly (x, y, ic, n) c sets y = ic(1) + ic(2)*x + ... + ic(n)*x**(n-1), c where x and y are multiple-precision numbers and c ic is an integer array of dimension at least n .gt. 0 c rounding not best possible, but directed rounding options c (rndrl = 2 or 3) give correct lower and upper bounds respectivley. common r integer i, ic(1), i2, n, r(1), sv, x(1), y(1) c save t etc. call mpsavn (sv) if (n .le. 0) call mperrm ( \$ 27hn .le. 0 in call to mppoly\$) c allocate temporary space. call mpnew (i2) call mpcim (ic(n), r(i2)) i = n - 1 if (i.le.0) go to 20 c loop to compute polynomial, taking care for directed roundings. 10 if (mod (i, 2) .eq. 0) call mprevr (-x(1)) call mpmul (r(i2), x, r(i2)) call mpsetr (iabs(r(sv+2))) call mpaddi (r(i2), ic(i), r(i2)) i = i - 1 if (i.gt.0) go to 10 20 call mpstr (r(i2), y) call mpresn (sv) return end subroutine mppwr (x, n, y) c returns y = x**n, for mp x and y, integer n, with 0**0 = 1. c x may be packed or unpacked, y is unpacked. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r integer i2, n, r(1), sv, tg, x(1), y(1) c save t etc. then increase t as necessary. call mpsavn (sv) call mpgd3 (n, tg) c allocate temporary space and move x. call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) if (n .ge. 0) go to 10 c n .lt. 0 if (x(1) .eq. 0) call mperrm ( \$ 38hx = 0 and n negative in call to mppwr\$) c here n .lt. 0, x .ne. 0. check for directed roundings. if (mod (iabs (n), 2) .eq. 0) call mprevr (-x(1)) call mprec (r(i2), r(i2)) call mpsetr (iabs(r(sv+2))) c now use mppwra to find r(i2)**abs(n). 10 call mppwra (r(i2), n, r(i2)) c round result call mprnd (r(i2), tg, y, iabs(r(sv)), 0) c restore t etc. and return. call mpresn (sv) return end subroutine mppwr2 (x, y, z) c returns z = x**y for mp numbers x, y and z, where x is c positive (x .eq. 0 allowed if y .gt. 0). slower than c mppwr and mpqpwr, so use them if possible. c x and/or y may be packed or unpacked, z is unpacked. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r common /mpcom/ b, t, dummy integer i2, i3, sv, tg integer b, dummy(21), r(1), t, x(1), y(1), z(1) integer mpcmpi, mpgd, mpget, mptlb c save t etc., then check sign of x. call mpsavn (sv) if (x(1)) 10, 20, 30 10 call mperrm (29hx negative in call to mppwr2\$) c here x is zero, return zero if y positive, otherwise error 20 if (y(1) .le. 0) call mperrm ( \$ 43hx zero and y nonpositive in call to mppwr2\$) c return zero here z(1) = 0 go to 50 c here x is positive. check for x = 1 or y = 0. 30 if ((y(1) .ne. 0) .and. (mpcmpi (x, 1) .ne. 0)) go to 40 c x**0 = 1, 1**y = 1. call mpcim (1, z) go to 50 c usual case here, x positive, y nonzero. allocate temporary c space and use mpln and mpexp to compute power. c increase t as necessary. 40 t = t + 1 + mpgd (100*mptlb(1)) \$ + mpgd (mpget (x, 2)) + max0 (0, mpget (y, 2)) tg = t call mpnew (i2) call mpmove (x, iabs(r(sv)), r(i2), tg) c take care for directed rounding call mprevr (-y(1)) call mpln (r(i2), r(i2)) call mpnew (i3) call mpmove (y, iabs(r(sv)), r(i3), tg) call mpmul (r(i2), r(i3), r(i2)) c restore rndrl before calling mpexp call mprevr (-y(1)) c if x**y is too large, mpexp will print error message call mpexp (r(i2), r(i2)) call mprnd (r(i2), tg, z, iabs(r(sv)), 0) 50 call mpresn (sv) return end subroutine mppwra (x, n, y) c returns y = x**abs(n) for mp x and y, integer n. (0**0 = 1) c uses no guard digits, chopped rather than rounded arithmetic. c called by mpexp, mppwr, mpqpwr, mproot, etc., c not recommended for independent use (use mppwr instead). c rounding not best possible, but directed rounding c (rndrl = 2 or 3) is correct in that lower and upper bounds c (respectively) are obtained. common r integer i2, n, n2, r(1), s, sv, x(1), y(1) c check for some special cases. n2 = iabs(n) if (n2 - 1) 10, 20, 30 c here abs(n) = 0 10 call mpcim (1, y) return c here abs(n) = 1 20 call mpstr (x, y) return c save t etc. and allocate temporary space etc. 30 call mpsavn (sv) call mpnew (i2) c compute sign of result s = x(1) if (mod(n2,2).eq.0) s = iabs(s) c move abs(x) to temporary storage call mpabs (x, r(i2)) c initialize product in y call mpcim (1, y) c use chopped arithmetic rather than rounded to save time. if (r(sv+2) .eq. 1) call mpsetr (0) c following for directed rounding call mprevr (-s) c loop, looking at bits of abs(n) from right. 40 if (mod (n2, 2) .ne. 0) call mpmul (y, r(i2), y) n2 = n2/2 if (n2.le.0) go to 50 call mpmul (r(i2), r(i2), r(i2)) go to 40 c restore everything and fix up sign of result. 50 call mpresn (sv) y(1) = s*y(1) return end subroutine mpqpwr (i, j, k, l, x) c sets multiple-precision x = (i/j)**(k/l) for integers c i, j, k and l, time = o(t**2). c rounding is determined by rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in the last place, c rndrl = 2 or 3 - directed roundings, see subroutine mpnzr. common r integer is, i2, js, ks, ls, sv, tg integer i, j, k, l, r(1), x(1) c save t etc. and increase working precision call mpsavn (sv) call mpgd3 (k, tg) c allocate temporary storage. call mpnew (i2) r(i2) = 0 is = i js = j ks = iabs(k) ls = l c for efficiency make ks positive and ls negative c (see comments in mproot and mppwr). if (k) 30, 10, 40 c (i/j)**(0/l) = 1 if j and l are nonzero. 10 call mpcim (1, r(i2)) if ((js.ne.0).and.(ls.ne.0)) go to 90 20 call mperrm (33hj = 0 or l = 0 in call to mpqpwr\$) c here k is negative 30 ls = -ls c now look at sign of ls 40 if (ls) 60, 20, 50 c ls positive so interchange is and js to make negative 50 is = j js = i ls = -ls c now ks positive, ls negative 60 if (is .eq. 0) call mperrm ( \$ 56hi=0 and k/l.lt.0, or j=0 and k/l.gt.0 in call to mpqpwr\$) c (i/0)**(negative) = 0 if i nonzero if (js.eq.0) go to 90 c to save time in mproot and mppwr, find gcd of ks and ls call mpgcd (ks, ls) c check for directed rounding if (mod (ks, 2) .eq. 0) call mprevr (is*(-(js/iabs(js)))) c check for ls = -1, treat as special case if (ls.ne.(-1)) go to 70 call mpcqm (js, is, r(i2)) go to 80 c usual case here, ls .ne. -1 70 call mpcqm (is, js, r(i2)) call mproot (r(i2), ls, r(i2)) 80 call mpsetr (iabs(r(sv+2))) call mppwra (r(i2), ks, r(i2)) c round result, restore t etc, and return 90 call mprnd (r(i2), tg, x, iabs(r(sv)), 0) call mpresn (sv) return end subroutine mprec (x, y) c returns y = 1/x, for mp x and y. c uses mpdivl if t small or rndrl .ne. 0, mproot otherwise. c time = o(m(t)). c rounding determined by rndrl in common /mpcom/ - see comments c in subroutine mpdiv. common r integer r(1), x(1), y(1), cross, i2, sv c following crossover point for use of mproot determined empirically. data cross /50/ call mpsavn (sv) c check for x zero if (x(1) .eq. 0) call mperrm ( \$ 44hattempted division by zero in call to mprec\$) c decide whether to use mpdivl or mproot. if ((r(sv) .lt. cross) .or. (r(sv+2) .ne. 0)) go to 10 call mproot (x, -1, y) go to 20 c here use mpdivl to compute 1/x 10 call mpnew (i2) call mpcim (1, r(i2)) call mpdivl (r(i2), x, y) 20 call mpresn (sv) return end subroutine mpres2 (n) c restores t, m and rndrl which were saved by c subroutine mpsavn. n must be the value returned by mpsavn. c sptr is unchanged. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr, \$ n, r(1), rndrl, sptr, t if ((n .ge. mnsptr) .and. (sptr .ge. (n+3))) go to 20 write (lun, 10) 10 format (44h *** illegal n or stpr on call to mpres2 ***) call mperr 20 t = r(n) m = r(n+1) rndrl = r(n+2) return end subroutine mpresn (n) c restores t, m and rndrl which were saved by c subroutine mpsavn. n must be the value returned by mpsavn. c sptr is restored to its value on the call to mpsavn (i.e. n). common /mpcom/ b, t, m, lun, mxr, sptr, dummy integer b, t, m, lun, mxr, sptr, n, dummy(17) call mpres2 (n) sptr = n return end subroutine mprevr (i) c if i .gt. 0 and rndrl (in common /mpcom/) is 2 or 3, c rndrl is replaced by 5 - rndrl (i.e. 2 and 3 interchanged). c this reverses the direction of directed rounding - see mpnzr. common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b, dummy(12), i, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, rndrl, sptr, t if ((i.gt.0) .and. (rndrl.gt.1)) rndrl = 5 - rndrl return end subroutine mprnd (x, tx, y, ty, k) c moves x + s*k*(b**(-ty))*abs(x) appropriately rounded to y, c where x is an mp number with tx digits, y is an mp number c with ty digits, k an integer, and c s = 0 if rndrl .le. 1, c -1 if rndrl .eq. 2, c +1 if rndrl .eq. 3. c note - rndrl = 0 has same effect as rndrl = 1, i.e. x rounded to c nearest representable number with ty digits. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b, dummy(12), k, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, r(1), rndrl, sptr, t, tx, ty, x(1), y(1), i2, sv c save t etc., set rndrl to 1 if zero. call mpsavn (sv) rndrl = max0 (1, rndrl) c see if addition necessary if ((rndrl.gt.1).and.(k.ne.0).and.(x(1).ne.0)) go to 10 call mpmove (x, tx, y, ty) go to 20 c here addition necessary. allocate temporary space etc. 10 t = tx call mpnew (i2) call mpstr (x, r(i2)) r(i2+1) = r(i2+1) - ty call mpmuli (r(i2), ((2*rndrl-5)*k*x(1)), r(i2)) call mpadd (r(i2), x, r(i2)) call mpmove (r(i2), tx, y, ty) c restore everything and return 20 call mpresn (sv) return end subroutine mproot (x, n, y) c returns y = x**(1/n) for integer n, mp x and y, using newtons c method without divisions. c x may be packed or unpacked, y is unpacked. c time is o(m(t)) unless abs(n) is large (when mppwr2 is used). c rounding determined by parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place c (so result exact if exactly representable). c rndrl = 2 or 3 - see comments in subroutine mpnzr. c abs(n) .le. m (maximum exponent) is required. common r common /mpcom/ b, t, m, dummy integer e, ep, ex, i2, i3, i4, i5, i6, j, k, ks, np, sv, tg, tg2 integer b, dummy(20), ktm, m, mptlb, n, r(1), t, x(1), y(1) c save t etc., use truncated arithmetic internally. call mpsavn (sv) call mpsetr (0) np = iabs (n) if (n .ne. 1) go to 10 c simply move x if n = 1 call mpunpk (x, y) go to 110 c check for various illegal argument combinations 10 if (n .eq. 0) call mperrm (24hn = 0 in call to mproot\$) if (x(1)) 30, 20, 40 20 y(1) = 0 if (n.lt.0) call mperrm (37hx = 0 and n .lt. 0 in call to mproot\$) go to 110 30 if (mod (np, 2) .eq. 0) call mperrm ( \$ 38hx .lt. 0 and n even in call to mproot\$) c increase t and allocate space 40 call mpgd3 (np, tg) call mpnew (i2) call mpnew (i3) c work with abs(x), fix up sign later. call mpmove (x, iabs(r(sv)), r(i2), tg) r(i2) = iabs (r(i2)) c check for large abs(n). if (np .ge. (m/4)) go to 130 call mpnew (i4) if (np .le. 2) go to 50 call mpnew (i5) call mpnew (i6) 50 e = -r(i2+1) c compute ep = floor (e/np) + 1 ep = e/np if ((np*ep) .le. e) ep = ep + 1 c scale to avoid under/overflow c scaled abs(x) is between 1 and b**np r(i2+1) = r(i2+1) + ep*np c lower and upper bounds on abs(scaled result) are 1/b and 1. call mpcqm (1, iabs(b), r(i3)) call mpstr (r(i3), r(i4)) if (np .gt. 2) call mpcim (1, r(i5)) c set ktm to maximum number of iterations allowed. ktm = t*mptlb(1) c reduce t at first. t = 2 call mpgd3 (np, tg2) ex = 0 c to speed up reciprocals and square roots, treat np = 1 or 2 as c special cases and get better starting approximations for them. if (np .eq. 1) \$ call mpcqm (iabs(b), (b*r(i2+2) + r(i2+3) + 1), r(i3)) if (np .ne. 2) go to 70 c here np = 2. compute j = upper bound on 4*abs(scaled x). call mpcmi (r(i2), j) j = 4*(j + 1) c compute k = upper bound on sqrt(j) c using integer newton iteration. k = j 60 ks = k k = (k + j/k)/2 if (k .lt. ks) go to 60 if ((k*k) .lt. j) k = k + 1 c now get lower bound on scaled abs (result). call mpcqm (2, k, r(i3)) c skip bisection phase if np .le. 2 70 if (np .le. 2) go to 80 c bisection loop here. check for infinite loop. ktm = ktm - 1 if (ktm .lt. 0) go to 120 call mpsub (r(i5), r(i3), r(i6)) call mpdivi (r(i6), 2, r(i6)) call mpadd (r(i3), r(i6), r(i6)) call mppwra (r(i6), np, r(i4)) call mpmul (r(i2), r(i4), r(i4)) call mpaddi (r(i4), -1, r(i4)) if (r(i4) .le. 0) call mpstr (r(i6), r(i3)) if (r(i4) .ge. 0) call mpstr (r(i6), r(i5)) c repeat bisection if abs (residual) .ge. 1/2 if (r(i4) .eq. 0) go to 90 if ((r(i4+1) .gt. 0) .or. \$ ((r(i4+1) .eq. 0) .and. ((2*r(i4+2)) .ge. b))) go to 70 c now newtons method should converge call mpstr (r(i6), r(i3)) go to 90 c newton loop starts here. choose good t. 80 t = min0 (tg, tg2 + 4*iabs(ex)) call mppwra (r(i3), np, r(i4)) call mpmul (r(i2), r(i4), r(i4)) call mpaddi (r(i4), -1, r(i4)) 90 ex = -t if (r(i4) .ne. 0) ex = r(i4+1) c check for infinite loop ktm = ktm - 1 if ((ktm .lt. 0) .or. (ex .gt. 0)) go to 120 call mpmul (r(i3), r(i4), r(i4)) call mpdivi (r(i4), np, r(i4)) call mpsub (r(i3), r(i4), r(i3)) c check for convergence if ((2*ex + tg) .gt. 0) go to 80 c end of loop, correct exponents r(i2+1) = r(i2+1) - ep*np r(i3+1) = r(i3+1) + ep c correct result for n positive if (n .lt. 0) go to 100 call mppwra (r(i3), n-1, r(i3)) call mpmul (r(i2), r(i3), r(i3)) 100 r(i3) = isign (r(i3), x(1)) c restore rndrl and do correct rounding call mpsetr (iabs(r(sv+2))) call mprnd (r(i3), tg, y, iabs(r(sv)), 1) 110 call mpresn (sv) return 120 call mperrm (35hiteration not converging in mproot\$) c here abs(n) is very large 130 call mpcqm (1, n, r(i3)) call mppwr2 (r(i2), r(i3), r(i3)) go to 100 end subroutine mpsavn (n) c saves t, m and rndrl on stack. c returns n = old value of sptr, c saves t in r(n), c m in r(n+1), c rndrl in r(n+2). c note - the argument n must not have same address as sptr. common r common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, dummy integer b, dummy(12), lun, m, mnexpn, mnsptr, mxexpn, mxr, mxsptr, \$ n, r(1), rndrl, sptr, t c allocate three words on stack. call mpnew2 (n, 3) c save t, m and rndrl. r(n) = t r(n+1) = m r(n+2) = rndrl return end subroutine mpscal (x, base, j) c sets x = x*(base**j), for mp x, integer j, and (small) c positive integer base (e.g. base = 10). c rounding not best possible, but directed rounding options c (rndrl = 2 or 3) give correct lower and upper bounds. common r integer b, i, ib, i2, ja, jp, sv integer base, j, mpparn, r(1), x(1) if (base.le.0) call mperrm (31hillegal base in call to mpscal\$) c return if j zero or x zero if ((j.eq.0).or.(x(1).eq.0)) return c save t etc. call mpsavn (sv) c otherwise multiply by base**j ja = iabs(j) c the numbers -500 and 100 were determined empirically for base = 10. c the optimum choice depends on base, b and t. if ((j.gt.(-500)).and.(j.lt.100)) go to 10 c here exponent large, so use mppwra to compute base**abs(j) call mpnew (i2) call mpcim (base, r(i2)) c take care for directed roundings. call mprevr ((-j)*x(1)) call mppwra (r(i2), ja, r(i2)) call mpsetr (iabs(r(sv+2))) if (j.lt.0) call mpdiv (x, r(i2), x) if (j.ge.0) call mpmul (x, r(i2), x) go to 30 c here abs(j) is small so probably faster to use mpdivi or mpmuli 10 jp = 1 ib = mpparn (16)/base b = mpparn (1) do 20 i = 1, ja jp = base*jp if ((jp.le.ib).and.(jp.ne.b).and.(i.lt.ja)) go to 20 if (j.lt.0) call mpdivi (x, jp, x) if (j.ge.0) call mpmuli (x, jp, x) jp = 1 20 continue 30 call mpresn (sv) return end subroutine mpset (lun, decpl, mt2, mxr) c redundant but included for compatibility with earlier c versions of mp. integer lun, decpl, mt2, mxr call mpset2 (lun, decpl, mt2, 1, mxr) return end subroutine mpset2 (lun, decpl, mt2, mnsptr, mxr) c c lun is the logical unit to be used for error messages. c sets base (b) and number of digits (t) to give the c equivalent of at least decpl decimal digits. c decpl should be a positive integer. c mt2 is the dimension of arrays used for mp numbers, c so an error occurs if the computed t exceeds mt2 - 2. c mxr is the size of blank common declared in the calling program, c and mnsptr is the first location in blank common available for c use by mp. (mp will use words mnsptr, mnsptr+1, ... , mxr of c blank common (counting from 1) for working storage.) c c let mxint be the largest integer of the form 2**k - 1 representable c as a signed integer on the machine (see mplarg). c the computed b and t satisfy the conditions c b**(t-1) .ge. 10**(decpl-1) and 8*b*b-1 .le. mxint . c approximately minimal t and maximal b satisfying c these conditions are chosen. c c mpset2 also sets c c m = mxint/4 (maximum allowable exponent), c sptr = mnsptr (stack pointer), c mxsptr = mnsptr (maximum value of sptr is saved here), c rndrl = 0 (see subroutine mpnzr), c ktunfl = 0 (underflow counter), c mxunfl = 0 (see subroutine mpunfl), c mnexpn = m+1 (minimum exponent which occurs saved here), c mxexpn = -m (maximum exponent which occurs saved here), c mxint = large integer (see subroutine mplarg), c exwid = 6 (exponent field width, see mpfout), c inrecl = 80 (input record length, see mpfin), c inbase = 10 (default input base for mpin), c outbas = 10 (default output base for mpout), c expch = 1he (default output exponent character), c chword = number of characters per word (returned by mpupw), c onescp = 1 if ones complement arithmetic, 0 otherwise. c c if these are not all as desired, change after the call to mpset2, c for example, one could set mxunfl = 10 if c execution were to be terminated after 10 mp underflows. c common r common /mpcom/ b, t, m, lunc, mxrc, sptr, mxsptr, mnsptc, \$ mxexpn, mnexpn, rndrl, ktunfl, mxunfl, \$ decplc, mt2c, mxint, exwid, inrecl, inbase, \$ outbas, expch, chword, onescp c increase dimension of c and modify data statement for maxchw c below if there are more than 10 characters per word. c 10 characters per word. integer c(10), maxchw, \$ i, i2, k, tab(10), \$ b, chw, chword, decpl, decplc, expch, exwid, inbase, inrecl, \$ ktunfl, lun, lunc, m, mnexpn, mnsptc, mnsptr, mt2, mt2c, \$ mxexpn, mxint, mxr, mxrc, mxsptr, mxunfl, onescp, outbas, \$ r(1), rndrl, sptr, t, mpdigs, mpget logical mpis c see comment above about maxchw. data maxchw /10/ data tab(1), tab(2), tab(3), tab(4) /1h0, 1h9, 1ha, 1hz/ data tab(5), tab(6), tab(7), tab(8) /1h+, 1h-, 1h., 1h\$/ data tab(9), tab(10) /1h , 1he/ c first set output unit, maximum size of common, etc. lunc = lun mxrc = mxr mnsptc = mnsptr sptr = mnsptr mxsptr = sptr rndrl = 0 ktunfl = 0 mxunfl = 0 decplc = decpl mt2c = mt2 c set default io parameters exwid = 6 inrecl = 80 inbase = 10 outbas = 10 expch = mpget (tab, 10) c check that lun in range 1,...,99 and write error message c on unit 6 if not. if ((lun.gt.0).and.(lun.lt.100)) go to 20 write (6, 10) lun 10 format (10h *** lun =, i10, 30h illegal in call to mpset2 ***) lunc = 6 call mperr c check for burroughs 6700 20 if (.not. mpis (tab(9), tab(6))) go to 40 25 write (lunc, 30) 30 format (35h *** replace mpis and/or mplarg ***) call mperr c set mxint to large representable integer (see above). 40 call mplarg (mxint) c if wordlength .ge. 12 bits, expect mxint .ge. 2047 if (mxint .lt. 2047) go to 25 c maximum exponent is mxint/4 m = mxint/4 c base b is largest power of 2 such that 8*b*b-1 .le. mxint b = 1 50 b = 2*b if ((4*b*b) .lt. (m+1)) go to 50 c set mnexpn and mxexpn to greater than m, less than -m resp. mnexpn = m+1 mxexpn = -m if (decpl.gt.0) go to 70 write (lunc, 60) 60 format (39h *** decpl .le. 0 in call to mpset2 ***) call mperr c ensure that b**(t-1) .ge. 10**(decpl-1). 70 t = mpdigs (decpl) c see if t too large for dimension statements i2 = t + 2 if (i2.le.mt2) go to 90 write (lunc, 80) i2 80 format (40h *** mt2 too small in call to mpset2 *** / \$ 45h *** increase mt2 and dimensions of mp arrays, \$ 12h to at least, i6, 4h ***) call mperr chword = maxchw c now choose onescp so that mpupw works for all c characters of interest (see mpin, mpout etc). 90 onescp = 0 do 120 k = 1, 2 do 100 i = 1, 9 call mpupw (tab(i), c, chw) chword = chw if ((chword .le. 1) .or. (chword .gt. maxchw)) go to 120 if (.not. mpis (tab(i), c(1))) go to 110 if (.not. mpis (tab(9), c(chword))) go to 110 100 continue go to 120 110 onescp = onescp + 1 120 continue if (onescp .le. 1) go to 140 write (lunc, 130) 130 format (22h *** replace mpupw ***) call mperr c check legality of b, t etc. 140 call mpchk return end subroutine mpsetr (n) c sets the parameter rndrl in common /mpcom/ to n = 0, 1, 2 or 3. c for effect of this see subroutine mpnzr etc. integer n if ((n .lt. 0) .or. (n .gt. 3)) call mperrm ( \$ 35hillegal argument in call to mpsetr\$) call mpparc (11, n) return end integer function mpsiga (x) c returns sign of packed or unpacked mp number x, sign (0) = 0. integer x(1) mpsiga = 0 if (x(1) .ne. 0) mpsiga = isign (1, x(1)) return end subroutine mpsigb (i, x) c sets sign of packed or unpacked mp number x to sign of i. c note that sign (0) = 0. integer i, j, x(1) c exponent and digits of x are unchanged. j = i if (j .eq. 0) x(1) = 0 if (j .ne. 0) x(1) = isign (x(1), j) return end subroutine mpsign (x, y, z) c sets z = abs(x)*sign(y) for mp x, y and z. c x and/or y may be packed or unpacked, z is unpacked. c note that z = 0 if y = 0 (i.e. sign (0) = 0). integer iy, x(1), y(1), z(1) iy = y(1) call mpunpk (x, z) z(1) = isign (z(1), iy) if (iy .eq. 0) z(1) = 0 return end subroutine mpsin (x, y) c returns y = sin(x) for (packed or unpacked) mp x, unpacked mp y. c uses mpcis, so time = o(sqrt(t)m(t)). c rounding options are implemented as for mpcis. common r integer i2, r(1), sv, x(1), y(1) call mpsavn (sv) call mpnew (i2) call mpcis (x, r(i2), y, .true.) call mpresn (sv) return end subroutine mpsinh (x, y) c returns y = sinh(x) for mp numbers x and y, x not too large. c x may be packed or unpacked, y is unpacked. c uses mpexp or mpexp1. c rounding options not yet implemented, no guard digits used. common r integer i2, i3, mpget, mpmexa, r(1), sv, x(1), xs, y(1) c save t etc. call mpsavn (sv) c save sign of x and check for zero, sinh(0) = 0 xs = x(1) if (xs.ne.0) go to 10 y(1) = 0 go to 40 c allocate temporary space. 10 call mpnew (i2) c use truncated arithmetic. call mpsetr (0) c work with abs(x) call mpabs (x, r(i2)) if (r(i2+1).le.0) go to 20 c here abs(x) .ge. 1, if too large mpexp gives error message c increase m to avoid overflow if sinh(x) representable call mpmexb (mpmexa (x) + 2, x) call mpexp (r(i2), r(i2)) call mprec (r(i2), y) call mpsub (r(i2), y, y) c restore m. if result overflows or underflows, mpdivi at c statement 30 will act accordingly. call mpmexb (mpget (r, sv+1), x) go to 30 c here abs(x) .lt. 1 so use mpexp1 to avoid cancellation 20 call mpnew (i3) call mpexp1 (r(i2), r(i3)) call mpaddi (r(i3), 2, r(i2)) call mpmul (r(i2), r(i3), y) call mpaddi (r(i3), 1, r(i2)) call mpdiv (y, r(i2), y) c divide by two and restore sign 30 call mpdivi (y, isign (2, xs), y) call mpresn (sv) 40 return end subroutine mpsqrt (x, y) c returns y = sqrt(x), using subroutine mproot if x .gt. 0. c x may be packed or unpacked, y is unpacked. c rounding is determined by rndrl in common /mpcom/ as follows - c rndrl = 0 or 1 - error less than 0.6 ulp, so exact if result can c be exactly represented. c rndrl = 2 or 3 - directed roundings - see subroutine mpnzr. integer x(1), y(1) c check legality of b, t, etc. call mpchk c check sign of x. if (x(1)) 10, 20, 30 10 call mperrm (36hnegative argument in call to mpsqrt\$) c here return zero. 20 y(1) = 0 return c here use mproot. 30 call mproot (x, 2, y) return end subroutine mpstov (n) c c called if stack space in blank common too small. c if possible, space should be expanded by at least n words c and mxr increased by the amount expanded. c c the new space must be contiguous with r(mxr), i.e. it must c include locations r(mxr+1), ... , r(mxr+n). c c at present mpstov does nothing because the method of expanding c the stack is machine and operating-system dependent. c common r common /mpcom/ b, t, m, lun, mxr, dummy integer b, dummy(18), lun, m, mxr, n, r(1), t write (lun, 10) 10 format (42h *** mpstov called but not implemented ***) return end subroutine mpstr (x, y) c sets y = x for mp x and y. c x may be in packed or unpacked format, c y will be in same format as x. common /mpcom/ b, t, dummy integer b, t, dummy(21), i, j, x(1), y(1), t2 c see if x and y have the same address (they often do) c note - this test is for efficiency only, and it does not matter c if the fact that x and y have the same address is not detected. j = x(1) y(1) = j + 1 if (j.eq.x(1)) go to 10 c here x(1) and y(1) must have the same address x(1) = j return c here x(1) and y(1) may have different addresses 10 y(1) = j c no need to move x(2), ... if x(1) = 0 if (j.eq.0) return t2 = t + 2 c reduce t2 if x is in packed format. if (iabs(j) .gt. 1) t2 = (t+3)/2 c now move remaining words of x. do 20 i = 2, t2 20 y(i) = x(i) return end subroutine mpsub (x, y, z) c subtracts y from x, forming result in z, for mp x, y and z. c rounding is controlled by the parameter rndrl in common /mpcom/ c as follows - c c rndrl = 0 - round towards zero if x*y .le. 0, c away from zero if x*y .gt. 0, c in both cases one guard digit is used, so result c is exact if significant cancellation occurs. c c rndrl = 1, 2 or 3 - see comments in subroutine mpnzr. c integer x(1), y(1), z(1), s c save sign of x s = x(1) c reverse sign of y, see if sign of x changed y(1) = -y(1) if (s.eq.x(1)) go to 10 c here x and y have the same address so result is zero y(1) = s z(1) = 0 return c here x and y have different addresses (or both zero) so use mpadd. 10 call mpadd (x, y, z) c restore sign of y, being careful in case y and z have same address s = z(1) y(1) = -y(1) z(1) = s return end subroutine mpsuba (x, y, i) c called by mplarg and mpupw. integer i, x(1), y(1) i = max0 (1, x(1) - y(1)) return end subroutine mptan (x, y) c sets y = tan(x) for mp x and y, c x may be packed or unpacked, y is unpacked. c uses mpcis, so time = o(sqrt(t)m(t)). c rounding options not implemented but some guard digits used. common r integer i2, i3, r(1), sv, tg, x(1), y(1) c tan(0) = 0 if (x(1) .ne. 0) go to 10 y(1) = 0 return c save t etc., increase t, and allocate working space. 10 call mpsavn (sv) call mpgd3 (1, tg) call mpnew (i2) call mpnew (i3) c use truncated arithmetic. call mpsetr (0) c move x and compute cos(x) and sin(x) call mpmove (x, iabs(r(sv)), r(i2), tg) call mpcis (r(i2), r(i2), r(i3), .true.) c following message indicates that x is too close to an odd c multiple of pi/2 if (r(i2) .eq. 0) call mperrm ( \$ 34htan(x) too large in call to mptan\$) call mpdiv (r(i3), r(i2), r(i2)) c round result, restore t etc. and return call mprnd (r(i2), tg, y, iabs(r(sv)), 0) call mpresn (sv) return end subroutine mptanh (x, y) c returns y = tanh(x) for mp numbers x and y, c x may be packed or unpacked, y is unpacked. c uses mpexp or mpexp1. c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. common r integer i2, i3, sv, tg, xs integer mpget, r(1), x(1), y(1) integer mpcmpq, mptlb if (x(1).ne.0) go to 10 c tanh(0) = 0 y(1) = 0 return c save t etc. 10 call mpsavn (sv) c increase t. call mpgd3 (1, tg) c allocate working space. call mpnew (i2) call mpnew (i3) c use truncated arithmetic. call mpsetr (0) c save sign and work with abs(x) call mpmove (x, iabs(r(sv)), r(i2), tg) xs = r(i2) r(i2) = 1 c see if abs(x) so large that result is +-1 if (mpcmpq (r(i2), 2*mptlb(iabs(r(sv))), 5) .le. 0) go to 20 c here abs(x) is very large call mpcim (xs, r(i3)) go to 50 c here abs(x) not so large 20 call mpmuli (r(i2), 2, r(i2)) if (r(i2+1).le.0) go to 30 c here abs(x) .ge. 1/2 so use mpexp call mpexp (r(i2), r(i2)) call mpaddi (r(i2), -1, r(i3)) call mpaddi (r(i2), 1, r(i2)) call mpdiv (r(i3), r(i2), r(i3)) go to 40 c here abs(x) .lt. 1/2, so use mpexp1 to avoid cancellation 30 call mpexp1 (r(i2), r(i2)) call mpaddi (r(i2), 2, r(i3)) call mpdiv (r(i2), r(i3), r(i3)) c restore sign 40 r(i3) = xs*r(i3) c round result 50 call mpres2 (sv) call mprnd (r(i3), tg, y, iabs(r(sv)), 1) c restore everything call mpresn (sv) c ensure that result is in closed interval (-1, +1). if (y(1).eq.0) return if (mpget (y, 2) .gt. 0) call mpcim (xs, y) return end integer function mptlb (j) c returns upper bound on abs(j)*log2(b). integer j, mpchgb, mpparn mptlb = mpchgb (2, mpparn(1), iabs(j)) return end subroutine mpunfl (x) c called on multiple-precision underflow, ie when the c exponent of mp number x would be less than 1-m. c action depends on parameters rndrl and mxunfl in common /mpcom/ - c if rndrl .le. 1 then c if mxunfl .eq. 0 (default option if mpset called), x is set c to zero and execution continues, c if mxunfl .gt. 0, x is set to zero and execution continues c unless mxunfl underflows have occurred, c when execution is terminated by a call c to mperr. c if rndrl .eq. 2, action as above except x is set to c -(smallest positive representable number). c if rndrl .eq. 3, action as above except x is set to c +(smallest positive representable number). common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, rndrl, ktunfl, mxunfl, dummy integer b, dummy(10), ktunfl, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, mxunfl, rndrl, sptr, t, x(1) c since m may have been overwritten, check b, t, m etc. call mpchk c the underflowing number is set to zero x(1) = 0 c set x to -+ smallest positive representable number if c rndrl = 2 or 3 respectively. if (rndrl.ge.2) call mpminr (x) if (rndrl.eq.2) x(1) = -x(1) c minimum exponent indicator (mnexpn) is set to -m mnexpn = -m c increment underflow counter (ktunfl) and call mperr if c mxunfl underflows have occurred. ktunfl = ktunfl + 1 if (ktunfl.ne.mxunfl) return write (lun, 10) ktunfl 10 format (4h ***, i10, 32h mp underflows have occurred ***) call mperr return end subroutine mpunfr (lunit, x) c reads mp number x from logical unit lunit without formatting, c assuming previously written by mpunfw. common r integer i2, lunit, r(1), sv, x(1) logical err call mpsavn (sv) call mpnew (i2) call mpio (r(i2), (r(sv)+3)/2, (-lunit), 1hu, err) if (err) call mperrm (33herror return from mpio in mpunfr\$) call mpunpk (r(i2), x) call mpresn (sv) return end subroutine mpunfw (x, lunit) c writes mp number x on logical unit lunit without formatting. c x should be a (packed or unpacked) mp variable. to save space, c x is written in packed format (so a record of length c int ((t+3)/2) words is written). common r integer i2, lunit, r(1), sv, x(1) logical err call mpsavn (sv) call mpnew (i2) call mppack (x, r(i2)) call mpio (r(i2), (r(sv)+3)/2, lunit, 1hu, err) if (err) call mperrm (33herror return from mpio in mpunfw\$) call mpresn (sv) return end subroutine mpunpk (y, x) c restores the mp number x which is stored in packed or c unpacked format in the integer array y. for further c details see subroutine mppack. common /mpcom/ b, t, dummy integer i, ib, j, two integer b, dummy(21), t, x(1), y(1) c check if y is packed or not. if (iabs(y(1)) .gt. 1) go to 10 c here y is zero or unpacked. call mpstr (y, x) return c here y is nonzero and packed. 10 two = 2 j = t/2 c do last digit separately if t odd. if (mod (t, 2) .ne. 0) x(t+2) = y(j+2)/b c work backwards in case x and y are the same array. if (j .lt. 2) go to 30 do 20 ib = 2, j i = j - ib x(2*i+6) = mod (y(i+3), b) 20 x(2*i+5) = y(i+3)/b c first word of y may be signed. 30 x(two+2) = mod (iabs(y(1)), b) x(two+1) = iabs(y(1))/b c move sign and exponent x(two) = y(two) x(1) = isign (1, y(1)) return end subroutine mpupdt (j) c updates minimum and maximum exponent indicators, c assuming j is the exponent of a newly formed nonzero mp number. c note - mpupdt does not check for mp overflow or underflow. common /mpcom/ b, t, m, lun, mxr, sptr, mxsptr, mnsptr, \$ mxexpn, mnexpn, dummy integer b, dummy(13), j, lun, m, mnexpn, mnsptr, mxexpn, mxr, \$ mxsptr, sptr, t mnexpn = min0 (j, mnexpn) mxexpn = max0 (j, mxexpn) return end subroutine mpupk (source, dest, ldest, lfield) c c this subroutine unpacks a packed hollerith string (source) c placing one character per word in the array dest (as if read in c a1 format). it continues unpacking until it finds a sentinel (\$) c or until it has filled ldest words of the array dest. c the length of the unpacked string (excluding the sentinel if c any) is returned in lfield. thus 0 .le. lfield .le. ldest. gt. 0 c common r integer dest(1), i, i2, j, k, k2, ldest, lfield, mpparn, n, \$ r(1), sentnl, source(1), sv, dollar logical mpis data dollar /1h\$/ c replace \$ by compiler-generated end-of-string sentinel (if any). data sentnl /1h\$/ call mpsavn (sv) c mpparn(22) should return upper bound on the number of characters c per word (see mpset2). call mpnew2 (i2, mpparn(22)) i = 0 j = 0 c loop to unpack each word of source 10 i = i + 1 call mpupw (source (i), r(i2), n) c assume end of string if n .le. 0 if (n .le. 0) go to 30 c move each character to dest do 20 k = 1, n k2 = i2 + k - 1 c check for sentinel or maximum length of dest if (mpis (r(k2), dollar) .or. mpis (r(k2), sentnl) .or. \$ (j .ge. ldest)) go to 30 j = j + 1 20 dest(j) = r(k2) if (j .lt. ldest) go to 10 30 lfield = j call mpresn (sv) return end subroutine mpupw (w, c, n) c c ************************* (see the mp c *** machine dependent *** users guide for c ************************* conversion hints) c c when called with packed character string in w, c returns n = number of characters/word and the c characters in c(1), ... , c(n), unpacked, left justified, c blank filled. c character w character*4 c(4) c n = 4 c(1) = w(1:1) c(2) = w(2:2) c(3) = w(3:3) c(4) = w(4:4) c return end subroutine mpzeta (n, x) c c returns mp x = zeta(n) for integer n .gt. 1, where c zeta(n) is the riemann zeta function (the sum from c i = 1 to infinity of i**(-n)). c c uses the euler-maclaurin series unless n = 2, 3, 4, 6 or 8, c so space required = o(t**2). (more accurately, c space = nl*t/2 + o(t) words, where nl is the number of terms c used in asymptotic expansion, nl = o(t). ) c time = o(t**3) (could be reduced to o(t**2) if bernoulli c numbers were precomputed). c c rounding is defined by the parameter rndrl in common /mpcom/ - c rndrl = 0 or 1 - error less than 0.6 units in last place. c rndrl = 2 or 3 - see comments in subroutine mpnzr. c common r common /mpcom/ b, t, dummy integer i, id(4), ip, i2, i3, i4, i5, kt, mxint, nl, nm, nt integer ntk, n2, p, q, sv, tg integer b, dummy(21), n, r(1), t, x(1) integer mpchgb, mpcmpi, mpparn, mptlb c zeta(n) known in terms of bernoulli numbers and c pi**n if n is even. id defines zeta(2), ... , zeta(8). data id(1), id(2), id(3), id(4) /6, 90, 945, 9450/ c save t etc. call mpsavn (sv) c increase t and allocate space for result call mpgd3 (mptlb(iabs(t)), tg) call mpnew (i2) r(i2) = 0 c use truncated arithmetic. call mpsetr (0) if (n .le. 1) call mperrm ( \$ 35hillegal argument in call to mpzeta\$) c here n .gt. 1. see if n = 2, 4, 6 or 8. if ((n.gt.8).or.(mod(n,2).ne.0)) go to 10 c here zeta(n) = (pi**n)/id(n/2) call mppi (r(i2)) call mppwra (r(i2), n, r(i2)) n2 = n/2 call mpdivi (r(i2), id(n2), r(i2)) go to 90 c see if n is very large. can return with zeta(n) = 1 to c required precision if 2**n .ge. b**(t-1) . 10 if (n .ge. mptlb (t-1)) go to 80 c check for n = 3 (see comments below). if (n .eq. 3) go to 100 c here we may use euler-maclaurin series. nl = -1 mxint = mpparn (16) q = 1 kt = t nt = n - 2 c start of loop to estimate nl = number of terms needed in c asymptotic expansion, and nm = number of terms in finite sum. 20 nl = nl + 1 c constant 5 chosen empirically to minimize time. if bernoulli c numbers were precomputed, this constant could be reduced. nm = 5*nl + 2 do 40 ntk = 1, 2 nt = nt + 1 30 if (q .lt. (mxint/nt)) go to 40 q = q/b + 1 kt = kt + 1 go to 30 40 q = nt*q c 6*nm + nm/4 is lower bound on 2*pi*nm . if (kt .ge. (mpchgb (iabs(b), nm, n-1) + \$ mpchgb (iabs(b), 6*nm + nm/4, 2*nl + 2))) go to 20 p = (t+3)/2 c allocate some working space. call mpnew (i3) c if nl .le. 0 no need to compute any bernoulli numbers if (nl .le. 0) go to 60 c allocate more space for result of mpbern. call mpnew2 (i4, nl*p) call mpnew (i5) c compute required bernoulli numbers (if zeta(n) is required c for several n, it would save time to precompute these). call mpbern ((-nl), p, r(i4)) call mpcqm (n, 2*nm, r(i3)) call mpdivi (r(i3), nm, r(i3)) c sum euler-maclaurin asymptotic series first do 50 i = 1, nl ip = i4 + (i-1)*p call mpunpk (r(ip), r(i5)) call mpmul (r(i3), r(i5), r(i5)) call mpadd (r(i2), r(i5), r(i2)) call mpmuls (r(i3), n+2*i-1, n+2*i, 2*i+1, 2*i+2) 50 call mpmuls (r(i3), 1, 1, nm, nm) c add integral approximation and multiply by nm**(1-n) call mpaddq (r(i2), 1, n-1, r(i2)) call mpscal (r(i2), nm, 1-n) c add finite sum in forward direction so can reduce t c more easily than if backward direction were used. 60 r(i3+1) = 0 do 70 i = 2, nm c reduce t for i**(-n) computation if possible t = max0 (2, tg + r(i3+1)) c compute i**(-n) using mpscal. call mpcim (1, r(i3)) call mpscal (r(i3), i, -n) c now r(i3) is i**(-n). halve last term in finite sum. if (i.eq.nm) call mpdivi (r(i3), 2, r(i3)) c restore t for addition t = tg c leave finite sum loop if mp underflow occurred if (r(i3).eq.0) go to 80 70 call mpadd (r(i3), r(i2), r(i2)) c final addition 80 call mpaddi (r(i2), 1, r(i2)) c restore rndrl and round result 90 call mpsetr (iabs(r(sv+2))) call mprnd (r(i2), tg, x, iabs(r(sv)), 1) c restore everything and check result call mpresn (sv) if ((mpcmpi (x, 1) .lt. 0) .or. (mpcmpi (x, 2) .ge. 0)) \$ call mperrm (43herror occurred in mpzeta, result incorrect\$) return c here use gospers series for zeta(3). this is faster and uses c less space than the euler-maclaurin series. c series for other odd n are given by h. riesel in c bit 13 (1973), 97-113. 100 call mpnew (i3) call mpcqm (1, 4, r(i2)) call mpcqm (5, 4, r(i3)) i = 0 c loop to sum gospers series 110 i = i + 1 c reduce t if possible t = max0 (2, min0 (tg, tg + r(i3+1))) call mpmulq (r(i3), (-i), 4*i+2, r(i3)) call mpmuls (r(i3), i, i, i+1, i+1) c restore t for addition t = tg call mpadd (r(i2), r(i3), r(i2)) c loop until terms are small if ((r(i3) .ne. 0) .and. (r(i3+1) .gt. (-t))) go to 110 go to 80 end