[Maxima-commits] CVS: maxima/src/numerical/slatec/fortran d9upak.f,NONE,1.1 dspenc.f,NONE,1.1

[Raymond T. <rt...@us...>]

Update of /cvsroot/maxima/maxima/src/numerical/slatec/fortran
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv28226/src/numerical/slatec/fortran

Added Files:
	d9upak.f dspenc.f 
Log Message:
Add Spence's function from SLATEC so we can evaluate dilogarithm more
accurately. 


--- NEW FILE: d9upak.f ---
      subroutine d9upak (x, y, n)
c august 1980 portable edition.  w. fullerton, los alamos scientific lab
c
c unpack floating point number x so that x = y * 2.0**n, where
c 0.5 .le. abs(y) .lt. 1.0 .
c
      double precision x, y, absx
c
      absx = dabs(x)
      n = 0
      y = 0.0d0
      if (x.eq.0.0d0) return
c
 10   if (absx.ge.0.5d0) go to 20
      n = n - 1
      absx = absx*2.0d0
      go to 10
c
 20   if (absx.lt.1.0d0) go to 30
      n = n + 1
      absx = absx*0.5d0
      go to 20
c
 30   y = dsign (absx, x)
      return
c
      end

--- NEW FILE: dspenc.f ---
      double precision function dspenc (x)
c feb 1978 edition.   w. fullerton, c3, los alamos scientific lab.
c
c evaluate a form of spence-s function defined by
c        integral from 0 to x of  -alog(abs(1-y))/y  dy.
c for abs(x).le.1, the uniformly convergent expansion
c        spenc = sum k=1,infinity  x**k / k**2     is valid.
c
c spence-s function can be used to evaluate much more general integral
c forms.  for example,
c        integral from 0 to z of  alog(a*x+b)/(c*x+d)  dx  =
c             alog(abs(b-a*d/c))*alog(abs(a*(c*z+d)/(a*d-b*c)))/c
c             - dspenc (a*(c*z+d)/(a*d-b*c)) / c.
c
c ref -- k. mitchell, philosophical magazine, 40, p.351 (1949).
c        stegun and abramowitz, ams 55, p.1004.
c
      double precision x, spencs(38), aln, pi26, xbig, d1mach, dcsevl,
     1  dlog
      external d1mach, dcsevl, dlog, initds
c
c series for spen       on the interval  0.          to  5.00000e-01
c                                        with weighted error   4.74e-32
c                                         log weighted error  31.32
c                               significant figures required  30.37
c                                    decimal places required  32.11
c
      data spencs(  1) / +.1527365598 8924058729 4668491002 8 d+0      /
      data spencs(  2) / +.8169658058 0510144035 0183818527 1 d-1      /
      data spencs(  3) / +.5814157140 7787308729 7735064118 2 d-2      /
      data spencs(  4) / +.5371619814 5415275422 4788900531 9 d-3      /
      data spencs(  5) / +.5724704675 1858262332 1060305478 2 d-4      /
      data spencs(  6) / +.6674546121 6493363436 0783543858 9 d-5      /
      data spencs(  7) / +.8276467339 7156769815 8439168901 1 d-6      /
      data spencs(  8) / +.1073315673 0306789512 7000587335 4 d-6      /
      data spencs(  9) / +.1440077294 3032394023 3459033151 3 d-7      /
      data spencs( 10) / +.1984442029 9659063678 9887713960 8 d-8      /
      data spencs( 11) / +.2794005822 1636387202 0199482161 5 d-9      /
      data spencs( 12) / +.4003991310 8833118230 7258044590 8 d-10     /
      data spencs( 13) / +.5823462892 0446384713 6813583575 7 d-11     /
      data spencs( 14) / +.8576708692 6386892780 9791477122 4 d-12     /
      data spencs( 15) / +.1276862586 2801930459 8948303343 3 d-12     /
      data spencs( 16) / +.1918826209 0425170811 6238041606 2 d-13     /
      data spencs( 17) / +.2907319206 9771381777 9579971967 3 d-14     /
      data spencs( 18) / +.4437112685 2767804625 5747364174 5 d-15     /
      data spencs( 19) / +.6815727787 4145995278 6735913560 7 d-16     /
      data spencs( 20) / +.1053017386 0155744295 4701941664 4 d-16     /
      data spencs( 21) / +.1635389806 7523771000 5182173457 0 d-17     /
      data spencs( 22) / +.2551852874 9404639323 1090164258 1 d-18     /
      data spencs( 23) / +.3999020621 9993601127 7047037951 9 d-19     /
      data spencs( 24) / +.6291501645 2168118765 1414917119 9 d-20     /
      data spencs( 25) / +.9933827435 6756776438 0388775253 3 d-21     /
      data spencs( 26) / +.1573679570 7499648167 2176380586 6 d-21     /
      data spencs( 27) / +.2500595316 8494761293 6927095466 6 d-22     /
      data spencs( 28) / +.3984740918 3838111392 1066325333 3 d-23     /
      data spencs( 29) / +.6366473210 0828438926 9132629333 3 d-24     /
      data spencs( 30) / +.1019674287 2396783670 7706197333 3 d-24     /
      data spencs( 31) / +.1636881058 9135188411 1107413333 3 d-25     /
      data spencs( 32) / +.2633310439 4176501173 4527999999 9 d-26     /
      data spencs( 33) / +.4244811560 1239768172 2436266666 6 d-27     /
      data spencs( 34) / +.6855411983 6800529168 2474666666 6 d-28     /
      data spencs( 35) / +.1109122433 4380564340 1898666666 6 d-28     /
      data spencs( 36) / +.1797431304 9998914573 6533333333 3 d-29     /
      data spencs( 37) / +.2917505845 9760951732 9066666666 6 d-30     /
      data spencs( 38) / +.4742646808 9286710613 3333333333 3 d-31     /
c
c pi26 = pi**2/6.0
      data pi26 / +1.644934066 8482264364 7241516664 6025189219 d0 /
      data nspenc, xbig / 0, 0.0d0 /
c
      if (nspenc.ne.0) go to 10
      nspenc = initds (spencs, 38, 0.1*sngl(d1mach(3)))
      xbig = 1.0d0/d1mach(3)
c
 10   if (x.gt.2.0d0) go to 60
      if (x.gt.1.0d0) go to 50
      if (x.gt.0.5d0) go to 40
      if (x.ge.0.0d0) go to 30
      if (x.gt.(-1.d0)) go to 20
c
c here if x .le. -1.0
c
      aln = dlog(1.0d0-x)
      dspenc = -pi26 - 0.5d0*aln*(2.0d0*dlog(-x)-aln)
      if (x.gt.(-xbig)) dspenc = dspenc
     1  + (1.d0 + dcsevl (4.d0/(1.d0-x)-1.d0, spencs, nspenc))/(1.d0-x)
      return
c
c -1.0 .lt. x .lt. 0.0
c
 20   dspenc = -0.5d0*dlog(1.0d0-x)**2
     1  - x*(1.d0+dcsevl(4.d0*x/(x-1.d0)-1.d0, spencs, nspenc))/(x-1.d0)
      return
c
c 0.0 .le. x .le. 0.5
c
 30   dspenc = x*(1.d0 + dcsevl (4.d0*x-1.d0, spencs, nspenc))
      return
c
c 0.5 .lt. x .le. 1.0
c
 40   dspenc = pi26
      if (x.ne.1.d0) dspenc = pi26 - dlog(x)*dlog(1.0d0-x)
     1  - (1.d0-x)*(1.d0+dcsevl(4.d0*(1.d0-x)-1.d0, spencs, nspenc))
      return
c
c 1.0 .lt. x .le. 2.0
c
 50   dspenc = pi26 - 0.5d0*dlog(x)*dlog((x-1.d0)**2/x)
     1  + (x-1.d0)*(1.d0+dcsevl(4.d0*(x-1.d0)/x-1.d0, spencs, nspenc))/x
      return
c
c x .gt. 2.0
c
 60   dspenc = 2.0d0*pi26 - 0.5d0*dlog(x)**2
      if (x.lt.xbig) dspenc = dspenc
     1  - (1.d0 + dcsevl (4.d0/x-1.d0, spencs, nspenc))/x
      return
c
      end