1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
|
/*
* lambertw.c
*
* (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
*/
#include <stdlib.h>
#include <stdio.h>
#include <getopt.h>
#include <math.h>
static double precision = 1e-16;
static FILE *logfile = NULL;
double f(double x) {
return x * exp(x);
}
double lambertW0(double y) {
double x = log(1 + y);
double delta = 1;
int count = 0;
while ((count < 20) && (delta > precision)) {
double error = f(x) - y;
if (logfile) {
fprintf(logfile, "%20.16f %20.16f\n", x, error);
}
double xnew = x - error / ((1 + x) * exp(x));
delta = fabs(xnew - x);
x = xnew;
count++;
}
if (count == 0) {
fprintf(stderr, "count exceeded\n");
return -1;
}
if (logfile) {
fprintf(logfile, "%d\n", count);
}
return x;
}
double lambertW1(double y) {
double x = log(-y);
double delta = 1;
int count = 0;
while ((count < 20) && (delta > precision)) {
double error = f(x) - y;
if (logfile) {
fprintf(logfile, "%20.16f %20.16f\n", x, error);
}
double xnew = x - error / ((1 + x) * exp(x));
delta = fabs(xnew - x);
x = xnew;
count++;
}
if (count == 0) {
fprintf(stderr, "count exceeded\n");
return -1;
}
if (logfile) {
fprintf(logfile, "%d\n", count);
}
return x;
}
struct option options[] = {
{ "precision", required_argument, NULL, 'p' },
{ "logfile", required_argument, NULL, 'l' },
{ "minus", no_argument, NULL, 'm' },
{ NULL, 0, NULL, 0 }
};
int main(int argc, char *argv[]) {
int c;
int l;
int minusone = 0;
while (EOF != (c = getopt_long(argc, argv, "p:l:m", options, &l)))
switch (c) {
case 'l':
logfile = fopen(optarg, "w");
break;
case 'p':
precision = atof(optarg);
break;
case 'm':
minusone = 1;
break;
}
for (; optind < argc; optind++) {
double y = atof(argv[optind]);
double x = (y >= 0) ? lambertW0(y) :
((minusone) ? lambertW1(y) : lambertW0(y));
printf("%20.16f -> %20.16f (%20.16f)\n", x, y, y - f(x));
}
if (logfile) {
fclose(logfile);
}
return EXIT_SUCCESS;
}
|
