1 #include "wright_omega.h"
2
3 /**
4 * Santiago Akle
5 * ICME Stanford University 2014
6 *
7 * Computes the value \omega(z) defined as the solution y to
8 * the equation y+log(y) = z ONLY FOR z real and z>=1.
9 *
10 * Loosely follows the recommendations by
11 * PIERS W. LAWRENCE, ROBERT M. CORLESS, and DAVID J. JEFFREY.
12 * Published in:
13 * Algorithm 917: Complex Double-Precision Evaluation of the Wright \omega Function
14 * ACM Transactions on Mathematical Software (TOMS) TOMS Homepage table of contents archive
15 * Volume 38 Issue 3, April 2012
16 * Article No. 20
17 * Publication Date 2012-04-01 (yyyy-mm-dd)
18 * Publisher ACM New York, NY, USA
19 * ISSN: 0098-3500 EISSN: 1557-7295 doi>10.1145/2168773.2168779
20 */
wrightOmega(pfloat z)21 pfloat wrightOmega(pfloat z)
22 {
23 pfloat w = 0.0;
24 pfloat r = 0.0;
25 pfloat q = 0.0;
26 pfloat zi = 0.0;
27
28 if(z<0.0) return -1; /* Fail if the input is not supported */
29
30 if(z<1.0+M_PI) /* If z is between 0 and 1+pi */
31 {
32 q = z-1;
33 r = q;
34 w = 1+0.5*r;
35 r *= q;
36 w += 1/16.0*r;
37 r *= q;
38 w -= 1/192.0*r;
39 r *= q;
40 w -= 1/3072.0*q;
41 r *= q; /* (z-1)^5 */
42 w += 13/61440.0*q;
43 /* Initialize with the taylor series */
44 }
45 else
46 {
47 r = log(z);
48 q = r;
49 zi = 1.0/z;
50 w = z-r;
51 q = r*zi;
52 w += q;
53 q = q*zi;
54 w += q*(0.5*r-1);
55 q = q*zi;
56 w += q*(1/3.0*r*r-3.0/2.0*r+1);
57 /* Initialize with w(z) = z-r+r/z^2(r/2-1)+r/z^3(1/3ln^2z-3/2r+1) */
58 }
59 /* FSC iteration */
60 /* Initialize the residual */
61 r = z-w-log(w);
62 z = (1+w);
63 q = z+2/3.0*r;
64 w *= 1+r/z*(z*q-0.5*r)/(z*q-r);
65 r = (2*w*w-8*w-1)/(72.0*(z*z*z*z*z*z))*r*r*r*r;
66 /* Check residual */
67 /* if(r<1.e-16) return w; */ /*XXX Just do two rounds */
68 z = (1+w);
69 q = z+2/3.0*r;
70 w *= 1+r/z*(z*q-0.5*r)/(z*q-r);
71 r = (2*w*w-8*w-1)/(72.0*(z*z*z*z*z*z))*r*r*r*r;
72
73 return w;
74 }
75