/******************************************************************/ /* */ /* ACM ALGORITHM 490: The Dilogarithm Function of a Real Argument */ /* */ /* E. S. Ginsberg, Department of Physics, University of */ /* Massachusetts at Bostion, Boston, MA., and */ /* D. Zaborowski, Information Processing Center, */ /* Massachusetts Institute of Technology, */ /* Cambridge, MA. */ /* 22 June 1973. */ /* */ /* with correction by Robert Morris, 11 July 1975. */ /* */ /******************************************************************/ /* Translated from FORTRAN to PL/I by R. A. Vowels, 20 October 2007. */ /* This algorithm computes the dilogarithm of the real argument X. */ /* The dilogarithm is real for values of X <= 1, and complex for */ /* values of X > 1. */ /* Note that the algorithm returns a complex value, the imaginary */ /* part of which is zero for X <= 1, and -pi*log(X) for X > 0. */ /* This PL/I version gives greater accuracy than the published */ /* version, typically one extra digit. */ (SUBRG, SIZE, FOFL): DILOG: PROCEDURE (X) RETURNS (FLOAT (18) COMPLEX); declare (a, b, by, c1, c2, c3, c4, dx, dy, test, w, x, x0, y, z) float (18); declare (dil, idil) float (18); /* Table of constants c(n) = (n(n+1)(n+2))**2 for n=1(1)n. */ declare c(30) float (18) static initial ( 36, 576, 3600, 14400, 44100, 112896, 254016, 518400, 980100, 1742400, 2944656, 4769856, 7452900, 11289600, 16646400, 23970816, 33802596, 46785600, 63680400, 85377600, 112911876, 147476736, 190440000, 243360000, 308002500, 386358336, 480661776, 593409600, 727380900, 885657600); declare pi float (18) static initial ( (4*atan(1q0)) ); declare pisqon6 float (18) static initial ( ((4*atan(1q0))**2/6) ); declare pisqon3 float (18) static initial ( ((4*atan(1q0))**2/3) ); /* Uses pi**2/3 = 3.289868... and */ /* pi**2/6 = 1.644394... */ declare n fixed binary (7); idil = 0; if x > 12.6q0 then go to L10; if x >= 12.59q0 then do; /* DILOG computed from Taylor series about the */ /* zero of DILOG, X0. */ idil = -pi*log(x); x0 = 12.5951703698450161e0; y = x/x0 - 1; z = 1/11.5951703698450161e0; w = y*z; c1 = (3*x0 - 2)/6; c2 = ((11*x0 - 15)*x0 + 6)/24; c3 = (((50*x0 - 104)*x0 + 84)*x0 - 24)/120; c4 = ((((274*x0 - 770)*x0 + 940)*x0 - 540)*x0 + 120)/720; dil = y*(1 - y*(0.5 - y*(1/3q0 - y*(0.25 - y*(0.2q0 - y/6)))))*log(z) - w*x0*y*(0.5 - w*(c1 - w*(c2 - w*(c3 - w*c4)))); return (complex(dil, idil)); end; if x >= 2 then go to L10; if x > 1 then go to L20; if x = 1 then go to L30; if x > 0.5 then go to L40; if x > 0.01Q0 then go to L50; if x < -1 then go to L60; if x < -0.01Q0 then go to L70; /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(1) */ dil = x*(1 + x*(0.25 + x*(1/9q0 + x*(0.0625 + x*(.04Q0 + x*(1/36q0 + x*(1/49q0 + x/64))))))); return (dil); /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(6), and description of */ /* this algorithm, Eq(4). */ L10: y = 1/x; by = -1 - y*(4 + y); dil = 3.28986813369645287E0 - log(x)**2/2 + (y*(4 + 5.75*y) + 3*(1+y)*(1-y) * log(1-y))/by; idil = -pi*log(x); if dil+4*y = dil then return (complex(dil, idil)); go to L80; /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(7), and description of */ /* this algorithm, Eq(4). */ L20: y = 1 - 1/x; dx = log(x); by = 1 + y*(4+y); dil = 1.64493406684822643E0 + dx*(dx/2 - log(x-1)) + (y*(4 +5.75*y) - 3*(1 + y)*dx/x)/by; idil = -pi*log(x); go to L80; /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(2). */ L30: dil = 1.64493406684822643E0; return (dil); /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(7), and description of */ /* this algorithm, Eq(4). */ L40: y = 1 - x; dx = log(x); by = -1 - y*(4 + y); dil = 1.64493406684822643E0 - dx*log(y) + (y*(4 + 5.75*y) + 3*(1+y)*dx*x)/by; go to L80; /* DILOG computed from description of this algorithm, Eq(4). */ L50: y = x; by = 1 + y*(4+y); dil = (y*(4+5.75*y) + 3*(1+y)*(1-y)*log(1-y))/by; go to L80; /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 245, Eq(12), with X = -X, and */ /* description of this algorithm, Eq(4). */ L60: y = 1/(1-x); dx = log(-x); dy = log(y); by = 1 + y*(4+y); dil = -1.64493406684822643E0 + dy/2*(dy + 2*dx) + (y*(4 + 5.75*y) + 3*(1+y)*(1-y)*(dx+dy))/by; if dil+4*y = dil then return (dil); go to L80; /* Formula from L. Lewin, "Dilogarithms and Associated Functions", */ /* MacDonald, London, 1958, p. 244, Eq(8), and description of */ /* of this algorithm, Eq(4). */ L70: y = x/(x-1); dx = log(1-x); by = -1 - y*(4+y); dil = (y*(4 + 5.75*y) - 3*(1+y)*(1-y)*dx)/by - dx*dx/2; L80: b = 4*y*y/by; do n = 1 to 30; b = b * y; a = b/c(n); test = dil; dil = dil + a; if dil = test then return (complex(dil, idil)); end; return (complex(dil, idil)); END DILOG;