This is an automated email from the ASF dual-hosted git repository. acassis pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/nuttx.git
commit ccbd3e4c3c2e71e6b4eafb0a84dc02a47eeb6360 Author: Marco Casaroli <[email protected]> AuthorDate: Mon Oct 19 02:30:27 2026 +0200 libs/libm: Add exp2(), log1p(), cbrtf(), cbrtl() and hypot(). math.h declares exp2(), log1p(), hypot(), their float and long double versions, cbrtf() and cbrtl(), but the built-in libm did not have them. Add them: - exp2(x) is 2^n * exp(f * ln 2), with n an integer and |f| <= 0.5. - log1p() uses the series of 2 * atanh(x / (2 + x)) for -0.5 < x < 1, and log(1 + x) for the other values. - cbrtf() and cbrtl() use the same Newton iterations as cbrt(). - hypot() calculates |a| * sqrt(1 + (b / a)^2) with |a| >= |b|, so it does not overflow before the result does. Assisted-by: Claude Code:claude-opus-5-5 Signed-off-by: Marco Casaroli <[email protected]> --- libs/libm/libm/CMakeLists.txt | 11 +++++ libs/libm/libm/Make.defs | 4 ++ libs/libm/libm/lib_cbrtf.c | 69 +++++++++++++++++++++++++++++++ libs/libm/libm/lib_cbrtl.c | 74 ++++++++++++++++++++++++++++++++++ libs/libm/libm/lib_exp2.c | 63 +++++++++++++++++++++++++++++ libs/libm/libm/lib_exp2f.c | 58 ++++++++++++++++++++++++++ libs/libm/libm/lib_exp2l.c | 63 +++++++++++++++++++++++++++++ libs/libm/libm/lib_hypot.c | 70 ++++++++++++++++++++++++++++++++ libs/libm/libm/lib_hypotf.c | 65 ++++++++++++++++++++++++++++++ libs/libm/libm/lib_hypotl.c | 70 ++++++++++++++++++++++++++++++++ libs/libm/libm/lib_log1p.c | 94 +++++++++++++++++++++++++++++++++++++++++++ libs/libm/libm/lib_log1pf.c | 89 ++++++++++++++++++++++++++++++++++++++++ libs/libm/libm/lib_log1pl.c | 94 +++++++++++++++++++++++++++++++++++++++++++ 13 files changed, 824 insertions(+) diff --git a/libs/libm/libm/CMakeLists.txt b/libs/libm/libm/CMakeLists.txt index 039a2e0faa2..7ec03aca011 100644 --- a/libs/libm/libm/CMakeLists.txt +++ b/libs/libm/libm/CMakeLists.txt @@ -151,6 +151,17 @@ if(CONFIG_LIBM) __sin.c lib_gamma.c lib_lgamma.c + lib_exp2.c + lib_exp2f.c + lib_exp2l.c + lib_log1p.c + lib_log1pf.c + lib_log1pl.c + lib_cbrtf.c + lib_cbrtl.c + lib_hypot.c + lib_hypotf.c + lib_hypotl.c lib_fma.c lib_fmal.c lib_fdim.c diff --git a/libs/libm/libm/Make.defs b/libs/libm/libm/Make.defs index 83aa81be227..312abcd9cc6 100644 --- a/libs/libm/libm/Make.defs +++ b/libs/libm/libm/Make.defs @@ -57,6 +57,10 @@ CSRCS += lib_fmax.c lib_fmaxf.c lib_fmaxl.c CSRCS += lib_fmin.c lib_fminf.c lib_fminl.c CSRCS += __cos.c __sin.c lib_gamma.c lib_lgamma.c +CSRCS += lib_exp2.c lib_exp2f.c lib_exp2l.c +CSRCS += lib_log1p.c lib_log1pf.c lib_log1pl.c +CSRCS += lib_cbrtf.c lib_cbrtl.c +CSRCS += lib_hypot.c lib_hypotf.c lib_hypotl.c CSRCS += lib_fma.c lib_fmal.c CSRCS += lib_fdim.c lib_fdimf.c lib_fdiml.c CSRCS += lib_ilogb.c lib_ilogbf.c lib_ilogbl.c diff --git a/libs/libm/libm/lib_cbrtf.c b/libs/libm/libm/lib_cbrtf.c new file mode 100644 index 00000000000..e7a78106773 --- /dev/null +++ b/libs/libm/libm/lib_cbrtf.c @@ -0,0 +1,69 @@ +/**************************************************************************** + * libs/libm/libm/lib_cbrtf.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +float cbrtf(float x) +{ + float m; + float y; + int e; + int r; + int i; + + if (x == 0.0F || isnan(x) || isinf(x)) + { + return x; + } + + /* |x| = m * 2^e: move e % 3 into m, then cbrt(|x|) = cbrt(m) * 2^(e / 3) + * with m in [0.5, 4). + */ + + m = frexpf(fabsf(x), &e); + r = e % 3; + if (r < 0) + { + r += 3; + } + + m = scalbnf(m, r); + e = (e - r) / 3; + + /* Newton iterations from a linear first guess (error < 12%) */ + + y = 0.68F + 0.227F * m; + for (i = 0; i < 4; i++) + { + y -= (y - m / (y * y)) / 3.0F; + } + + return copysignf(scalbnf(y, e), x); +} diff --git a/libs/libm/libm/lib_cbrtl.c b/libs/libm/libm/lib_cbrtl.c new file mode 100644 index 00000000000..a37857b3f82 --- /dev/null +++ b/libs/libm/libm/lib_cbrtl.c @@ -0,0 +1,74 @@ +/**************************************************************************** + * libs/libm/libm/lib_cbrtl.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_LONG_DOUBLE +long double cbrtl(long double x) +{ + long double m; + long double y; + int e; + int r; + int i; + + if (x == 0.0L || isnan(x) || isinf(x)) + { + return x; + } + + /* |x| = m * 2^e: move e % 3 into m, then cbrt(|x|) = cbrt(m) * 2^(e / 3) + * with m in [0.5, 4). + */ + + m = frexpl(fabsl(x), &e); + r = e % 3; + if (r < 0) + { + r += 3; + } + + m = scalbnl(m, r); + e = (e - r) / 3; + + /* Newton iterations from a linear first guess (error < 12%) */ + + y = 0.68L + 0.227L * m; + for (i = 0; i < 7; i++) + { + y -= (y - m / (y * y)) / 3.0L; + } + + return copysignl(scalbnl(y, e), x); +} +#endif diff --git a/libs/libm/libm/lib_exp2.c b/libs/libm/libm/lib_exp2.c new file mode 100644 index 00000000000..b55118fe35d --- /dev/null +++ b/libs/libm/libm/lib_exp2.c @@ -0,0 +1,63 @@ +/**************************************************************************** + * libs/libm/libm/lib_exp2.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_DOUBLE +double exp2(double x) +{ + const double ln2 = 0.693147180559945309417232121458176568; + double n; + + if (isnan(x)) + { + return x; + } + + /* Large |x|: scalbn() gives the overflow or the underflow */ + + if (x > 2200.0) + { + x = 2200.0; + } + else if (x < -2200.0) + { + x = -2200.0; + } + + /* 2^x = 2^n * e^(f * ln 2), with n an integer and |f| <= 0.5 */ + + n = floor(x + 0.5); + return scalbn(exp((x - n) * ln2), (int)n); +} +#endif diff --git a/libs/libm/libm/lib_exp2f.c b/libs/libm/libm/lib_exp2f.c new file mode 100644 index 00000000000..c3ad9ae2dcb --- /dev/null +++ b/libs/libm/libm/lib_exp2f.c @@ -0,0 +1,58 @@ +/**************************************************************************** + * libs/libm/libm/lib_exp2f.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +float exp2f(float x) +{ + const float ln2 = 0.693147180559945309417232121458176568F; + float n; + + if (isnan(x)) + { + return x; + } + + /* Large |x|: scalbn() gives the overflow or the underflow */ + + if (x > 300.0F) + { + x = 300.0F; + } + else if (x < -300.0F) + { + x = -300.0F; + } + + /* 2^x = 2^n * e^(f * ln 2), with n an integer and |f| <= 0.5 */ + + n = floorf(x + 0.5F); + return scalbnf(expf((x - n) * ln2), (int)n); +} diff --git a/libs/libm/libm/lib_exp2l.c b/libs/libm/libm/lib_exp2l.c new file mode 100644 index 00000000000..94a56eace28 --- /dev/null +++ b/libs/libm/libm/lib_exp2l.c @@ -0,0 +1,63 @@ +/**************************************************************************** + * libs/libm/libm/lib_exp2l.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_LONG_DOUBLE +long double exp2l(long double x) +{ + const long double ln2 = 0.693147180559945309417232121458176568L; + long double n; + + if (isnan(x)) + { + return x; + } + + /* Large |x|: scalbn() gives the overflow or the underflow */ + + if (x > 33000.0L) + { + x = 33000.0L; + } + else if (x < -33000.0L) + { + x = -33000.0L; + } + + /* 2^x = 2^n * e^(f * ln 2), with n an integer and |f| <= 0.5 */ + + n = floorl(x + 0.5L); + return scalbnl(expl((x - n) * ln2), (int)n); +} +#endif diff --git a/libs/libm/libm/lib_hypot.c b/libs/libm/libm/lib_hypot.c new file mode 100644 index 00000000000..aa905cb4851 --- /dev/null +++ b/libs/libm/libm/lib_hypot.c @@ -0,0 +1,70 @@ +/**************************************************************************** + * libs/libm/libm/lib_hypot.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_DOUBLE +double hypot(double x, double y) +{ + double a = fabs(x); + double b = fabs(y); + double t; + + if (isinf(a) || isinf(b)) + { + return INFINITY; + } + + if (isnan(a) || isnan(b)) + { + return a + b; + } + + if (a < b) + { + t = a; + a = b; + b = t; + } + + if (b == 0.0) + { + return a; + } + + /* a * sqrt(1 + (b / a)^2) does not overflow before the result does */ + + t = b / a; + return a * sqrt(1.0 + t * t); +} +#endif diff --git a/libs/libm/libm/lib_hypotf.c b/libs/libm/libm/lib_hypotf.c new file mode 100644 index 00000000000..e764ab9f618 --- /dev/null +++ b/libs/libm/libm/lib_hypotf.c @@ -0,0 +1,65 @@ +/**************************************************************************** + * libs/libm/libm/lib_hypotf.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +float hypotf(float x, float y) +{ + float a = fabsf(x); + float b = fabsf(y); + float t; + + if (isinf(a) || isinf(b)) + { + return INFINITY_F; + } + + if (isnan(a) || isnan(b)) + { + return a + b; + } + + if (a < b) + { + t = a; + a = b; + b = t; + } + + if (b == 0.0F) + { + return a; + } + + /* a * sqrt(1 + (b / a)^2) does not overflow before the result does */ + + t = b / a; + return a * sqrtf(1.0F + t * t); +} diff --git a/libs/libm/libm/lib_hypotl.c b/libs/libm/libm/lib_hypotl.c new file mode 100644 index 00000000000..1b1a0ead57f --- /dev/null +++ b/libs/libm/libm/lib_hypotl.c @@ -0,0 +1,70 @@ +/**************************************************************************** + * libs/libm/libm/lib_hypotl.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_LONG_DOUBLE +long double hypotl(long double x, long double y) +{ + long double a = fabsl(x); + long double b = fabsl(y); + long double t; + + if (isinf(a) || isinf(b)) + { + return INFINITY_L; + } + + if (isnan(a) || isnan(b)) + { + return a + b; + } + + if (a < b) + { + t = a; + a = b; + b = t; + } + + if (b == 0.0L) + { + return a; + } + + /* a * sqrt(1 + (b / a)^2) does not overflow before the result does */ + + t = b / a; + return a * sqrtl(1.0L + t * t); +} +#endif diff --git a/libs/libm/libm/lib_log1p.c b/libs/libm/libm/lib_log1p.c new file mode 100644 index 00000000000..309796037ef --- /dev/null +++ b/libs/libm/libm/lib_log1p.c @@ -0,0 +1,94 @@ +/**************************************************************************** + * libs/libm/libm/lib_log1p.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <float.h> +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_DOUBLE +double log1p(double x) +{ + double s; + double s2; + double t; + double term; + double sum; + int k; + + if (isnan(x) || x == INFINITY) + { + return x; + } + + if (x == -1.0) + { + return -HUGE_VAL; + } + + if (x < -1.0) + { + return (x - x) / (x - x); + } + + if (x <= -0.5 || x >= 1.0) + { + return log(1.0 + x); + } + + if (1.0 + x == 1.0) + { + return x; + } + + /* log(1 + x) = 2 * atanh(s) = 2 * (s + s^3 / 3 + s^5 / 5 + ...), with + * s = x / (2 + x) and |s| < 1 / 3 + */ + + s = x / (2.0 + x); + s2 = s * s; + t = s; + sum = s; + for (k = 3; ; k += 2) + { + t *= s2; + term = t / k; + if (fabs(term) <= fabs(sum) * DBL_EPSILON * 0.25) + { + break; + } + + sum += term; + } + + return 2.0 * sum; +} +#endif diff --git a/libs/libm/libm/lib_log1pf.c b/libs/libm/libm/lib_log1pf.c new file mode 100644 index 00000000000..97b1eaf9dc2 --- /dev/null +++ b/libs/libm/libm/lib_log1pf.c @@ -0,0 +1,89 @@ +/**************************************************************************** + * libs/libm/libm/lib_log1pf.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <float.h> +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +float log1pf(float x) +{ + float s; + float s2; + float t; + float term; + float sum; + int k; + + if (isnan(x) || x == INFINITY_F) + { + return x; + } + + if (x == -1.0F) + { + return -HUGE_VALF; + } + + if (x < -1.0F) + { + return (x - x) / (x - x); + } + + if (x <= -0.5F || x >= 1.0F) + { + return logf(1.0F + x); + } + + if (1.0F + x == 1.0F) + { + return x; + } + + /* log(1 + x) = 2 * atanh(s) = 2 * (s + s^3 / 3 + s^5 / 5 + ...), with + * s = x / (2 + x) and |s| < 1 / 3 + */ + + s = x / (2.0F + x); + s2 = s * s; + t = s; + sum = s; + for (k = 3; ; k += 2) + { + t *= s2; + term = t / k; + if (fabsf(term) <= fabsf(sum) * FLT_EPSILON * 0.25F) + { + break; + } + + sum += term; + } + + return 2.0F * sum; +} diff --git a/libs/libm/libm/lib_log1pl.c b/libs/libm/libm/lib_log1pl.c new file mode 100644 index 00000000000..bcd45ea1a11 --- /dev/null +++ b/libs/libm/libm/lib_log1pl.c @@ -0,0 +1,94 @@ +/**************************************************************************** + * libs/libm/libm/lib_log1pl.c + * + * SPDX-License-Identifier: Apache-2.0 + * + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. The + * ASF licenses this file to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance with the + * License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT + * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the + * License for the specific language governing permissions and limitations + * under the License. + * + ****************************************************************************/ + +/**************************************************************************** + * Included Files + ****************************************************************************/ + +#include <nuttx/config.h> +#include <nuttx/compiler.h> + +#include <float.h> +#include <math.h> + +/**************************************************************************** + * Public Functions + ****************************************************************************/ + +#ifdef CONFIG_HAVE_LONG_DOUBLE +long double log1pl(long double x) +{ + long double s; + long double s2; + long double t; + long double term; + long double sum; + int k; + + if (isnan(x) || x == INFINITY_L) + { + return x; + } + + if (x == -1.0L) + { + return -HUGE_VALL; + } + + if (x < -1.0L) + { + return (x - x) / (x - x); + } + + if (x <= -0.5L || x >= 1.0L) + { + return logl(1.0L + x); + } + + if (1.0L + x == 1.0L) + { + return x; + } + + /* log(1 + x) = 2 * atanh(s) = 2 * (s + s^3 / 3 + s^5 / 5 + ...), with + * s = x / (2 + x) and |s| < 1 / 3 + */ + + s = x / (2.0L + x); + s2 = s * s; + t = s; + sum = s; + for (k = 3; ; k += 2) + { + t *= s2; + term = t / k; + if (fabsl(term) <= fabsl(sum) * LDBL_EPSILON * 0.25L) + { + break; + } + + sum += term; + } + + return 2.0L * sum; +} +#endif
