C      ALGORITHM 443, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 16, NO. 2, February, 1973, PP.123--124.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Thu Dec 15 13:28:07 2005
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver
RESULTS= Res

Objs1= driver.o src.o
driver: $(Objs1)
	$(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS)
Res: driver 
	./driver >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c***********************************************************************
c
cc TOMS443_PRB tests TOMS443.
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS443_PRB'
      write ( *, '(a)' ) '  Test TOMS algorithm 443, to evaluate'
      write ( *, '(a)' ) '  Lambert''s W function.'

      call test01
      call test02

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS435_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'

      stop
      end
      subroutine test01

c***********************************************************************
c
cc TEST01 tests WEW_A
c
      implicit none

      real en
      integer n_data
      real w1
      real w2
      real wew_a
      real x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' ) '  Test WEW_A to evaluate'
      write ( *, '(a)' ) '  Lambert''s W function.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '      X         Exact Value      Computed'
      write ( *, '(a)' ) ' '

      n_data = 0

10    continue

      call lambert_w_values ( n_data, x, w1 )

      if ( n_data <= 0 ) then
        go to 20
      end if

      if ( x == 0.0 ) then
        w2 = 0.0
      else
        w2 = wew_a ( x, en )
      end if

      write ( *, '(2x,f12.4,2x,f8.4,2x,g16.8)' ) 
     &  x, w1, w2

      go to 10

20    continue

      return
      end
      subroutine test02

c***********************************************************************
c
cc TEST02 tests WEW_B
c
      implicit none

      real en
      integer n_data
      real w1
      real w2
      real wew_b
      real x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST02'
      write ( *, '(a)' ) '  Test WEW_B to evaluate'
      write ( *, '(a)' ) '  Lambert''s W function.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '      X         Exact Value      Computed'
      write ( *, '(a)' ) ' '

      n_data = 0

10    continue

      call lambert_w_values ( n_data, x, w1 )

      if ( n_data <= 0 ) then
        go to 20
      end if

      if ( x == 0.0 ) then
        w2 = 0.0
      else
        w2 = wew_b ( x, en )
      end if

      write ( *, '(2x,f12.4,2x,f8.4,2x,g16.8)' ) 
     &  x, w1, w2

      go to 10

20    continue

      return
      end
      subroutine lambert_w_values ( n_data, x, fx )

c*******************************************************************************
c
cc LAMBERT_W_VALUES returns some values of the Lambert W function.
c
c  Discussion:
c
c    The function W(X) is defined implicitly by:
c
c      W(X) * e^W(X) = X
c
c    The function is also known as the "Omega" function.
c
c    In Mathematica, the function can be evaluated by:
c
c      W = ProductLog [ X ]
c
c  Modified:
c
c    23 February 2005
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    R M Corless, G H Gonnet, D E Hare, D J Jeffrey, D E Knuth,
c    On the Lambert W Function,
c    Advances in Computational Mathematics,
c    Volume 5, 1996, pages 329-359.
c
c    Brian Hayes,
c    "Why W?",
c    The American Scientist,
c    Volume 93, March-April 2005, pages 104-108.
c
c    Eric Weisstein,
c    "Lambert's W-Function",
c    CRC Concise Encyclopedia of Mathematics,
c    CRC Press, 1998.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Wolfram Media / Cambridge University Press, 1999.
c
c  Parameters:
c
c    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
c    first call.  On each call, the routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    Output, real X, the argument of the function.
c
c    Output, real FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 22 )

      real fx
      real fx_vec(n_max)
      integer n_data
      real x
      real x_vec(n_max)
 
      save fx_vec
      save x_vec

      data fx_vec /
     &  0.0000000000000000D+00,  0.3517337112491958D+00, 
     &  0.5671432904097839D+00,  0.7258613577662263D+00, 
     &  0.8526055020137255D+00,  0.9585863567287029D+00, 
     &  0.1000000000000000D+01,  0.1049908894964040D+01, 
     &  0.1130289326974136D+01,  0.1202167873197043D+01, 
     &  0.1267237814307435D+01,  0.1326724665242200D+01, 
     &  0.1381545379445041D+01,  0.1432404775898300D+01, 
     &  0.1479856830173851D+01,  0.1524345204984144D+01, 
     &  0.1566230953782388D+01,  0.1605811996320178D+01, 
     &  0.1745528002740699D+01,  0.3385630140290050D+01, 
     &  0.5249602852401596D+01,  0.1138335808614005D+02 /
      data x_vec /
     &  0.0000000000000000D+00,  0.5000000000000000D+00, 
     &  0.1000000000000000D+01,  0.1500000000000000D+01, 
     &  0.2000000000000000D+01,  0.2500000000000000D+01, 
     &  0.2718281828459045D+01,  0.3000000000000000D+01, 
     &  0.3500000000000000D+01,  0.4000000000000000D+01, 
     &  0.4500000000000000D+01,  0.5000000000000000D+01, 
     &  0.5500000000000000D+01,  0.6000000000000000D+01, 
     &  0.6500000000000000D+01,  0.7000000000000000D+01, 
     &  0.7500000000000000D+01,  0.8000000000000000D+01, 
     &  0.1000000000000000D+02,  0.1000000000000000D+03, 
     &  0.1000000000000000D+04,  0.1000000000000000D+07 /

      if ( n_data < 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max < n_data ) then
        n_data = 0
        x = 0.0
        fx = 0.0
      else
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS443_PRB
  Test TOMS algorithm 443, to evaluate
  Lambert's W function.
 
TEST01
  Test WEW_A to evaluate
  Lambert's W function.
 
      X         Exact Value      Computed
 
        0.0000    0.0000     0.0000000    
        0.5000    0.3517    0.35173371    
        1.0000    0.5671    0.56714326    
        1.5000    0.7259    0.72586131    
        2.0000    0.8526    0.85260552    
        2.5000    0.9586    0.95858634    
        2.7183    1.0000     1.0000000    
        3.0000    1.0499     1.0499089    
        3.5000    1.1303     1.1302893    
        4.0000    1.2022     1.2021679    
        4.5000    1.2672     1.2672379    
        5.0000    1.3267     1.3267246    
        5.5000    1.3815     1.3815453    
        6.0000    1.4324     1.4324048    
        6.5000    1.4799     1.4798568    
        7.0000    1.5243     1.5243452    
        7.5000    1.5662     1.5662310    
        8.0000    1.6058     1.6058120    
       10.0000    1.7455     1.7455280    
      100.0000    3.3856     3.3856301    
     1000.0000    5.2496     5.2496028    
  1000000.0000   11.3834     11.383359    
 
TEST02
  Test WEW_B to evaluate
  Lambert's W function.
 
      X         Exact Value      Computed
 
        0.0000    0.0000     0.0000000    
        0.5000    0.3517    0.35173371    
        1.0000    0.5671    0.56714326    
        1.5000    0.7259    0.72586137    
        2.0000    0.8526    0.85260552    
        2.5000    0.9586    0.95858634    
        2.7183    1.0000     1.0000000    
        3.0000    1.0499     1.0499089    
        3.5000    1.1303     1.1302893    
        4.0000    1.2022     1.2021680    
        4.5000    1.2672     1.2672378    
        5.0000    1.3267     1.3267246    
        5.5000    1.3815     1.3815454    
        6.0000    1.4324     1.4324049    
        6.5000    1.4799     1.4798568    
        7.0000    1.5243     1.5243453    
        7.5000    1.5662     1.5662310    
        8.0000    1.6058     1.6058120    
       10.0000    1.7455     1.7455281    
      100.0000    3.3856     3.3856311    
     1000.0000    5.2496     5.2496014    
  1000000.0000   11.3834     11.383318    
 
TOMS435_PRB
  Normal end of execution.
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
C**** Incorporates changes given in Remark by Einarsson
C**** CACM 17(4) April 1974 p.225
      REAL FUNCTION WEW_A ( X, EN )
C
C  ITERATIVE SOLUTION OF X = W * EXP ( W ) WHERE X IS GIVEN.
C  (NOVEMBER 1970)
C  (REVISED - SEPTEMBER 1971)
C  VERSION A -- CDC 6600 MACHINE ACCURACY.
C
C  INPUT PARAMETER:
C    X  ARGUMENT OF W(X)
C
C  OUTPUT PARAMETERS:
C    WEW  THE DESIRED SOLUTION.
C    EN   THE LAST RELATIVE CORRECTION TO W(X).
C
C  SET CONSTANTS...
C     .. Scalar Arguments ..
      REAL EN,X
C     ..
C     .. Local Scalars ..
C**** Next 2 statement replaced in remark by Einarsson
C     REAL C1,C2,C3,C4,FLOGX,TEMP,TEMP2,WN,Y,ZN
C     INTEGER NEWE
      REAL FLOGX,TEMP,TEMP2,WN,Y,ZN
C**** End of replaced statements
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC LOG
C**** Next 10 statements deleted in remark by Einarsson
C     ..
C     .. Save statement ..
C      SAVE C1,C2,C3,C4,NEWE
C      DATA NEWE / 1 /
C      IF ( NEWE ) 10, 20, 10
C10    NEWE = 0
C      C1 = 4. / 3.
C      C2 = 7. / 3.
C      C3 = 5. / 6.
C      C4 = 2. / 3.
C**** End of replaced statements
C
C  COMPUTE INITIAL GUESS...
20    FLOGX = LOG ( X )
      IF ( X - 6.46 ) 30, 30, 40
C**** Next statement replaced in remark by Einarsson
C30    WN = X * ( 1. + C1 * X ) / ( 1. + X * ( C2 + C3 * X ) )
30    WN = X*(3.0 + 4.0*X)/(3.0 + X*(7.0 + 2.5*X))
C**** End of replaced statements
      ZN = FLOGX - WN - LOG ( WN )
      GO TO 50
40    WN = FLOGX
      ZN = -LOG ( WN )
50    CONTINUE
C
C  ITERATION ONE...
      TEMP = 1. + WN
C**** Next statement replaced in remark by Einarsson
C     Y = 2. * TEMP * ( TEMP + C4 * ZN ) - ZN
      Y = 2. * TEMP * ( TEMP + ZN/1.5 ) - ZN
C**** End of replaced statements
      WN = WN * ( 1. + ZN * Y / ( TEMP * ( Y - ZN ) ) )
C
C  ITERATION TWO...
      ZN = FLOGX - WN - LOG ( WN )
      TEMP = 1. + WN
**** Next statement replaced in remark by Einarsson
C     TEMP2 = TEMP + C4 * ZN
      TEMP2 = TEMP + ZN/1.5
      EN = ZN * TEMP2 / ( TEMP * TEMP2 - .5 * ZN )
      WN = WN * ( 1. + EN )
C
C  RETURN...
      WEW_A = WN
      RETURN
      END
      REAL FUNCTION WEW_B ( X, EN )
C
C  ITERATIVE SOLUTION OF X = W * EXP ( W ) WHERE X IS GIVEN.
C  (NOVEMBER 1970)
C  (REVISED - SEPTEMBER 1971)
C  VERSION B -- MAXIMUM RELATIVE ERROR 3.E-7.
C
C  INPUT PARAMETER:
C    X  ARGUMENT OF W(X)
C
C  OUTPUT PARAMETERS:
C    WEW  THE DESIRED SOLUTION.
C    EN   THE LAST RELATIVE CORRECTION TO W(X).
C
C  SET CONSTANTS...
C     .. Scalar Arguments ..
      REAL EN,X,F
C     ..
C     .. Local Scalars ..
C**** Next 2 statement replaced in remark by Einarsson
C     REAL C1,C2,C3,C4,FLOGX,TEMP,WN,Y,ZN
C     INTEGER NEWE
      REAL FLOGX,TEMP,WN,Y,ZN
C**** End of replaced statements
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC LOG
C**** Next 10 statements deleted in remark by Einarsson
C     ..
C     .. Save statement ..
C      SAVE C1,C2,C3,C4,NEWE
C      DATA NEWE / 1 /
C      IF ( NEWE ) 10, 20, 10
C10    NEWE = 0
C      C1 = 4. / 3.
C      C2 = 7. / 3.
C      C3 = 5. / 6.
C      C4 = 2. / 3.
C**** End of replaced statements
C
C  COMPUTE INITIAL GUESS...
20    FLOGX = LOG ( X )
      IF ( X - .7385 ) 30, 30, 40
C**** Next statement replaced in remark by Einarsson
C30    WN = X * ( 1. + C1 * X ) / ( 1. + X * ( C2 + C3 * X ) )
30    WN = X*(3.0 + 4.0*X)/(3.0 + X*(7.0 + 2.5*X))
C**** End of replaced statements
      GO TO 50
40    WN = FLOGX - 24. * ((FLOGX+2.) * FLOGX-3.) /
     +     (( .7 * FLOGX + 58. ) * FLOGX + 127. )
50    CONTINUE
C
C  ITERATION ONE...
      ZN = FLOGX - WN - LOG ( WN )
      TEMP = 1. + WN
C**** Next statement replaced in remark by Einarsson
C     Y = 2. * TEMP * ( TEMP + C4 * ZN ) - ZN
      Y = 2. * TEMP * ( TEMP + ZN/1.5 ) - ZN
C**** End of replaced statements
      EN = ZN * Y / ( TEMP * ( Y - ZN ) )
      WN = WN * ( 1. + EN )
C
C  RETURN...
      WEW_B = WN
      RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0