1 #include "wright.hh"
2
3 /**********************************************************************/
4 /* wrightomega is the simple routine for evaluating the wright omega */
5 /* function. */
6 /* */
7 /* Calling: */
8 /* w = wrightomega(z) */
9 /* */
10 /* Input: */
11 /* z -- double complex */
12 /* Value to evaluate Wrightomega(z) at. */
13 /* */
14 /* Output: */
15 /* w -- double complex */
16 /* Value of Wrightomega(z) */
17 /* */
18 /**********************************************************************/
19
20 /*
21 Also published as ACM TOMS 917; relicensed as BSD by the author.
22
23 Copyright (C) Piers Lawrence.
24 All rights reserved.
25
26 Redistribution and use in source and binary forms, with or without
27 modification, are permitted provided that the following conditions are met:
28
29 1. Redistributions of source code must retain the above copyright notice, this
30 list of conditions and the following disclaimer.
31
32 2. Redistributions in binary form must reproduce the above copyright notice,
33 this list of conditions and the following disclaimer in the documentation
34 and/or other materials provided with the distribution.
35
36 3. Neither the name of the copyright holder nor the names of its contributors
37 may be used to endorse or promote products derived from this software without
38 specific prior written permission.
39
40 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
41 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
42 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
43 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
44 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
45 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
46 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
47 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
48 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
49 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
50 */
51
52
53 /**********************************************************************/
54 /* wrightomega_ext is the extended routine for evaluating the wright */
55 /* omega function with the option of extracting the last update step, */
56 /* the penultimate residual and the condition number estimate. */
57 /* */
58 /* Calling: */
59 /* success = wrightomega_ext(z,w,e,r,cond); */
60 /* */
61 /* Input: */
62 /* z -- double complex */
63 /* Value to evaluate Wrightomega(z) at. */
64 /* */
65 /* w -- double complex* */
66 /* Pointer to return value of Wrightomega(z). */
67 /* */
68 /* cond -- double complex* */
69 /* Pointer to the condition number estimate. If NULL the */
70 /* condition number is not calculated. */
71 /* */
72 /* Output: returns 0 on successful exit. */
73 /**********************************************************************/
74
75 #include <Python.h>
76 extern "C" {
77 #include <numpy/npy_math.h>
78 #include "sf_error.h"
79 #include "_c99compat.h"
80 #include "_round.h"
81 }
82
83 #include <cmath>
84 #include <cfloat>
85
86 using std::complex;
87
88 #define NaN NPY_NAN
89 #define TWOITERTOL DBL_EPSILON
90
91 const complex<double> I(0.0, 1.0);
92
93
94 int
wrightomega_ext(complex<double> z,complex<double> * w,complex<double> * cond)95 wright::wrightomega_ext(complex<double> z, complex<double> *w,
96 complex<double> *cond)
97 {
98 double pi = NPY_PI, s = 1.0;
99 double x, y, ympi, yppi, near;
100 complex<double> e, r, pz, wp1, t, fac;
101
102
103 /* extract real and imaginary parts of z */
104 x=real(z);
105 y=imag(z);
106
107 /* compute if we are near the branch cuts */
108 ympi=y-pi;
109 yppi=y+pi;
110 near=0.1e-1;
111
112 /* Test for floating point exceptions */
113 /*****************************/
114 /* NaN output for NaN input */
115 /*****************************/
116 if(sc_isnan(x) || sc_isnan(y))
117 {
118 *w = complex<double>(NaN, NaN);
119 return 0;
120 }
121 /*********************************/
122 /* Signed zeros between branches */
123 /*********************************/
124 else if(sc_isinf(x) && (x < 0.0) && (-pi < y) && (y<= pi))
125 {
126 if (fabs(y) <= pi/2.0)
127 {
128 if (y >= 0)
129 {
130 *w = complex<double>(0.0, 0.0);
131 }
132 else
133 {
134 *w = complex<double>(0.0, -0.0);
135 }
136 }
137 else
138 {
139 if (y >= 0)
140 {
141 *w = complex<double>(-0.0, 0.0);
142 }
143 else
144 {
145 *w = complex<double>(-0.0, -0.0);
146 }
147 }
148 return 0;
149 }
150 /**************************/
151 /* Asymptotic for large z */
152 /**************************/
153 else if(sc_isinf(x) || sc_isinf(y))
154 {
155 *w = complex<double>(x, y);
156 return 0;
157 }
158
159 /******************************************/
160 /* Test If exactly on the singular points */
161 /******************************************/
162 if((x==-1.0) && (fabs(y)==pi))
163 {
164 *w = complex<double>(-1.0, 0.0);
165 return 0;
166 }
167
168
169 /* Choose approximation based on region */
170 /**********************************/
171 /* Region 1: upper branch point */
172 /* Series about z=-1+Pi*I */
173 /**********************************/
174 if ((-2.0<x && x<=1.0 && 1.0<y && y< 2.0*pi))
175 {
176 pz=conj(sqrt(conj(2.0*(z+1.0-I*pi))));
177 *w=-1.0+(I+(1.0/3.0+(-1.0/36.0*I+(1.0/270.0+1.0/4320.0*I*pz)*pz)*pz)*pz)*pz;
178 }
179 /**********************************/
180 /* Region 2: lower branch point */
181 /* Series about z=-1-Pi*I */
182 /**********************************/
183 else if ((-2.0<x && x<=1.0 && -2.0*pi<y && y<-1.0))
184 {
185 pz=conj(sqrt(conj(2.0*(z+1.0+I*pi))));
186 *w=-1.0+(-I+(1.0/3.0+(1.0/36.0*I+(1.0/270.0-1.0/4320.0*I*pz)*pz)*pz)*pz)*pz;
187 }
188 /*********************************/
189 /* Region 3: between branch cuts */
190 /* Series: About -infinity */
191 /*********************************/
192 else if (x <= -2.0 && -pi < y && y <= pi)
193 {
194 pz=exp(z);
195 *w=(1.0+(-1.0+(3.0/2.0+(-8.0/3.0+125.0/24.0*pz)*pz)*pz)*pz)*pz;
196 if (*w == 0.0)
197 {
198 sf_error("wrightomega", SF_ERROR_UNDERFLOW, "underflow in exponential series");
199 /* Skip the iterative scheme because it computes log(*w) */
200 if (cond != NULL)
201 {
202 *cond = z/(1.0+*w);
203 }
204 return 0;
205 }
206 }
207 /**********************/
208 /* Region 4: Mushroom */
209 /* Series about z=1 */
210 /**********************/
211 else if (((-2.0 < x) && (x<=1.0) && (-1.0 <= y) && (y <= 1.0))
212 || ((-2.0 < x) && (x-0.1e1)*(x-0.1e1)+y*y<=pi*pi))
213 {
214 pz=z-1.0;
215 *w=1.0/2.0+1.0/2.0*z+(1.0/16.0+(-1.0/192.0+(-1.0/3072.0+13.0/61440.0*pz)*pz)*pz)*pz*pz;
216 }
217 /*************************/
218 /* Region 5: Top wing */
219 /* Negative log series */
220 /*************************/
221 else if (x<=-0.105e1 && pi<y && y-pi<=-0.75e0*(x+0.1e1))
222 {
223 t=z-I*pi;
224 pz=log(-t);
225 *w = t - pz;
226 fac = pz/t;
227 *w += fac;
228 fac /= t;
229 *w += fac*(0.5*pz - 1.0);
230 fac /= t;
231 *w += fac*(pz*pz/3.0 - 3.0*pz/2.0 + 1.0);
232 if (abs(z) > 1e50)
233 /* Series is accurate and the iterative scheme could overflow */
234 {
235 if (cond != NULL)
236 {
237 *cond = z/(1.0+*w);
238 }
239 return 0;
240 }
241 }
242 /***************************/
243 /* Region 6: Bottom wing */
244 /* Negative log series */
245 /***************************/
246 else if (x<=-0.105e1 && 0.75e0*(x+0.1e1)< y+pi && y+pi<=0.0e0)
247 {
248 t=z+I*pi;
249 pz=log(-t);
250 *w = t - pz;
251 fac = pz/t;
252 *w += fac;
253 fac /= t;
254 *w += fac*(0.5*pz - 1.0);
255 fac /= t;
256 *w += fac*(pz*pz/3.0 - 3.0*pz/2.0 + 1.0);
257 if (abs(z) > 1e50)
258 /* Series is accurate and the iterative scheme could overflow */
259 {
260 if (cond != NULL)
261 {
262 *cond = z/(1.0+*w);
263 }
264 return 0;
265 }
266 }
267 /************************************/
268 /* Region 7: Everywhere else */
269 /* Series solution about infinity */
270 /************************************/
271 else
272 {
273 pz=log(z);
274 *w = z - pz;
275 fac = pz/z;
276 *w += fac;
277 fac /= z;
278 *w += fac*(0.5*pz - 1.0);
279 fac /= z;
280 *w += fac*(pz*pz/3.0 - 3.0*pz/2.0 + 1.0);
281 if (abs(z) > 1e50)
282 /* Series is accurate and the iterative scheme could overflow */
283 {
284 if (cond != NULL)
285 {
286 *cond = z/(1.0+*w);
287 }
288 return 0;
289 }
290 }
291
292 /**********************************/
293 /* Regularize if near branch cuts */
294 /**********************************/
295 if (x <= -0.1e1 + near && (fabs(ympi) <= near || fabs(yppi) <= near))
296 {
297 s = -1.0;
298 if (fabs(ympi) <= near)
299 {
300 /* Recompute ympi with directed rounding */
301 ympi = add_round_up(y, -pi);
302
303 if( ympi <= 0.0)
304 {
305 ympi = add_round_down(y, -pi);
306 }
307
308 z = x + I*ympi;
309 }
310 else
311 {
312 /* Recompute yppi with directed rounding */
313 yppi = add_round_up(y, pi);
314
315 if( yppi <= 0.0)
316 {
317 yppi = add_round_down(y, pi);
318 }
319
320 z = x + I*yppi;
321 }
322 }
323
324 /*****************/
325 /* Iteration one */
326 /*****************/
327 *w=s**w;
328 r=z-s**w-log(*w);
329 wp1=s**w+1.0;
330 e=r/wp1*(2.0*wp1*(wp1+2.0/3.0*r)-r)/(2.0*wp1*(wp1+2.0/3.0*r)-2.0*r);
331 *w=*w*(1.0+e);
332
333 /*****************/
334 /* Iteration two */
335 /*****************/
336 if(abs((2.0**w**w-8.0**w-1.0)*pow(abs(r),4.0)) >= TWOITERTOL*72.0*pow(abs(wp1),6.0) )
337 {
338 r=z-s**w-log(*w);
339 wp1=s**w+1.0;
340 e=r/wp1*(2.0*wp1*(wp1+2.0/3.0*r)-r)/(2.0*wp1*(wp1+2.0/3.0*r)-2.0*r);
341 *w=*w*(1.0+e);
342 }
343
344 /***********************/
345 /* Undo regularization */
346 /***********************/
347 *w=s**w;
348
349 /***************************************************/
350 /* Provide condition number estimate if requested */
351 /***************************************************/
352 if(cond != NULL)
353 {
354 *cond = z/(1.0+*w);
355 }
356
357 return 0;
358 }
359
360 complex<double>
wrightomega(complex<double> z)361 wright::wrightomega(complex<double> z)
362 {
363 complex<double> w;
364 wrightomega_ext(z,&w,NULL);
365 return w;
366 }
367
368
369 double
wrightomega_real(double x)370 wright::wrightomega_real(double x)
371 {
372 double w, wp1, e, r;
373
374 /* NaN output for NaN input */
375 if (sc_isnan(x))
376 {
377 return x;
378 }
379
380 /* Positive infinity is asymptotically x */
381 /* Negative infinity is zero */
382 if (sc_isinf(x))
383 {
384 if (x > 0.0)
385 {
386 return x;
387 }
388 else
389 {
390 return 0.0;
391 }
392 }
393
394 if (x < -50.0)
395 {
396 /*
397 * Skip the iterative scheme because it exp(x) is already
398 * accurate to double precision.
399 */
400 w = exp(x);
401 if (w == 0.0)
402 {
403 sf_error("wrightomega", SF_ERROR_UNDERFLOW, "underflow in exponential series");
404 }
405 return w;
406 }
407 if (x > 1e20)
408 {
409 /*
410 * Skip the iterative scheme because the result is just x to
411 * double precision
412 */
413 return x;
414 }
415
416 /* Split int three distinct intervals (-inf,-2), [-2,1), [1,inf) */
417 if (x < -2.0)
418 {
419 /* exponential is approx < 1.3e-1 accurate */
420 w = exp(x);
421 }
422 else if (x < 1)
423 {
424 /* on [-2,1) approx < 1.5e-1 accurate */
425 w = exp(2.0*(x-1.0)/3.0);
426 }
427 else
428 {
429 /* infinite series with 2 terms approx <1.7e-1 accurate */
430 w = log(x);
431 w = x - w + w/x;
432 }
433
434 /* Iteration one of Fritsch, Shafer, and Crowley (FSC) iteration */
435 r = x - w - log(w);
436 wp1 = w + 1.0;
437 e = r/wp1*(2.0*wp1*(wp1+2.0/3.0*r)-r)/(2.0*wp1*(wp1+2.0/3.0*r)-2.0*r);
438 w = w*(1.0+e);
439
440 /* Iteration two (if needed based on the condition number) */
441 if (fabs((2.0*w*w-8.0*w-1.0)*pow(fabs(r),4.0)) >= TWOITERTOL*72.0*pow(fabs(wp1),6.0))
442 {
443 r = x - w - log(w);
444 wp1 = w+1.0;
445 e = r/wp1*(2.0*wp1*(wp1+2.0/3.0*r)-r)/(2.0*wp1*(wp1+2.0/3.0*r)-2.0*r);
446 w = w*(1.0+e);
447 }
448
449 return w;
450 }
451