Edit - Stack Overflow

You are not logged in. Your edit will be placed in a queue until it is peer reviewed.

We welcome edits that make the post easier to understand and more valuable for readers. Because community members review edits, please try to make the post substantially better than how you found it, for example, by fixing grammar or adding additional resources and hyperlinks.

Required fields*

Rev

Required fields*

floating-point implementation of a function with removable singularity

Consider the following function:

ƒ(x) = sin(sqrt(x)) / sqrt(x) if x > 0
ƒ(x) = 1 if x = 0
ƒ(x) = sinh(sqrt(-x)) / sqrt(-x) if x < 0

This is an entire function with a Taylor series at 0 of: 1/1!, -1/3!, 1/5!, -1/7!, ... that I want to implement in floating point (think IEEE-754 single, IEEE-754 double or arbitrary precision like with MPFR.

My current C-ish implementation with a floating-point type F looks something like:

F func (F x)
{
 F q = sqrt (fabs (x));
 if (x * x == 0)
 return 1 - x / 6;
 else
 return x > 0
 ? sin (q) / q
 : sinh (q) / q;
}

From what I can tell the implementation works reasonably well in an environment of 0, but I am wondering if there are loophole values where it does not work (well).

Answer*

**Loophole values where it does not work**

**Infinities:** When `x` is +infinity, the expected math result is 0.0, yet OP's code may return Not-a-number due to `sin(infinity)`.

When `x` is -infinity, the expected math result is _infinity_, yet OP's code may return Not-a-number due to `sinh(infinity)/infinity`.

Code could test for these non-finite `x` values.

**Overflow:** When `x` is modestly negative, `sinh(q)` may overflow into infinity, yet the extend math quotient `sinh(q)/q` is still finite.

```c
// float example
float yn = sinhf(q) / q;
// If a non-finite value resulted, try a slower, less precise, yet wider range approach.
if (!isfinite(yn)) {
 yn = expf(q - logf(q) - logf(2));
}
```

@njuffa offers [a good alternative](https://stackoverflow.com/questions/79615576/floating-point-implementation-of-a-function-with-removable-singularity/79615833?noredirect=1#comment140414360_79615833).

Another improvement:

```c
float yn = sinhf(q) / q;
if (!isfinite(yn)) {
 // I expect this to work for all `x` where the math `sinhf(q) / q <= FLT_MAX`.
 yn = expf(q/2);
 yn *= yn/(2*q);
}
```

**Modest magnitude `x`:** I found that, using `float`, `q = sqrtf(fabsf(x)); return sinhf(q) / q;` loses about 6 bits of correctness before `sinhf(q)` hits overflow. This is because the error is forming the square root `q` and then `q` applied to `sinhf(g)` incurs an error proportional to `q`.

To really reduce this function's error, a better, TBD, approach is needed as `x` grows negative (about `x < -4`).

**Trig errors**

`q = sqrtf(x)` has up to 0.5 ULP error in the result. When `q < 99.99%*π`, `sinf(q)` is well behaved. Yet at/near every following multiple of π, the ULP error of `sinf(q)` can be terrible, even if the absolution value is very reasonable. As `q` grows and its least significant significant digit nears 1.0, _sign_ and magnitude of `sinf(q)` becomes meaningless.

Additional code is needed to compensate for square root errors. Considerations start at `x > 1.57` for ULP errors. 

There are no for absolute errors concerns for finite `x`. 

----

**OP likely area of concern**.

By using `if (x * x == 0) return 1 - x / 6;` code is well behaved _near_ zero. The end-points of the region is about `x = [-sqrt(DBL_TRUE_MIN) ... sqrt(DBL_TRUE_MIN)]`. Now let us look at the end points as `x` is just within this region and just outside (e.g. near `+/-sqrt(DBL_TRUE_MIN)`)

Taylor's series for `sin(x)` is `x - x^3/3! + x^5/5! +/- ...`

Taylor's series for `sinh(x)` is `x + x^3/3! + x^5/5! + ...`

Mathematically, all looks fine as the function is continuous, the 1st,2nd derivative also.

We could look close and compute with `double x` values near this area, chart that and report _loophole values_ where the `y` steps are a tad jagged. Yet doing this examination is dubious as there is scant benefit over the simpler:

```
F func2 (F x)
{
 if (x == 0)
 return 1;
 else
 F q = sqrt (fabs (x));
 return x > 0
 ? sin (q) / q
 : sinh (q) / q;
}
```

Here, the nature of `x` near 0.0 for `sin(x), sinh(x)` is well researched and tested.

Still, to test, consider using `float` as it is feasible to test all `float`.

**Find first difference**

```c
#include <float.h>
#include <math.h>
#include <stdio.h>

int main() {
 float x = sqrtf(FLT_TRUE_MIN);
 printf("sqrtf(FLT_TRUE_MIN):%24a\n", sqrtf(FLT_TRUE_MIN));
 x = FLT_TRUE_MIN;
 for (int i = 0; 1; i++) {
 float q = sqrtf(x);
 float yp = sinf(q) / q;
 float yn = sinhf(q) / q;
 float y0 = 1.0f - x / 6.0f;
 if (y0 != yp || y0 != yn) {
 printf("x:%a, q:%a, sin(q)/q:%a, sinh(q)/q:%a, 1-x/6:%a\n", //
 x, q, yp, yn, y0);
 break;
 }
 x = nextafterf(x, 2 * x);
 }
 return 0;
}
```

Sample output:

```
sqrtf(FLT_TRUE_MIN): 0x1.6a09e6p-75
x: 0x1p-46, q:0x1p-23, sin(q)/q:0x1p+0, sinh(q)/q:0x1.fffffep-1, 1-x/6:0x1p+0
```

It is clear that OP's original code has no discernable issues in the `+/-sqrtf(FLT_TRUE_MIN)` range as all 3 sub-function result in the same `y` values for `x` until much larger than `sqrtf(FLT_TRUE_MIN)` - as expected.

Still, I do not see a benefit (unless most usage values are near 0.0) in doing the `if (x*x == 0)` test versus the simple `if (x==0)`.

**[Soap box](https://en.wikipedia.org/wiki/Soapbox)**

Certainly OP was most interesting is dealing with `x` at or near 0.0 and the "loopholes" there about `x*x == 0`. So was I when I first saw this post. Later I realized that many `x` values in the range [-8828.3 ... -7995.2] failed.

It highlights the importance of _functional_ testing across the spectrum of `x` before attempting optimization in speed or low error performance.

Draft saved

Draft discarded

Edit Summary*

Cancel

1

(+1) For "very negative" x, why not simply compute 1 / (2 * sqrt (-x) * exp (-sqrt (-x))), taking advantage of subnormals? This should produce results with smaller error than going through logarithmic space and exponentiating at the end, which incurs a substantial amount of error magnification at the extreme end of the range.

njuffa
– njuffa

2025年05月10日 21:26:12 +00:00
Commented May 10, 2025 at 21:26
@njuffa Interesting idea. Answer updated. Just went through a similar issue here with best approach. Will ponder some more.

chux
– chux

2025年05月10日 21:36:56 +00:00
Commented May 10, 2025 at 21:36
I decided to be less lazy and actually tried my idea and your proposed fix, and the results of a quick experiment seem clear. For a double implementation, the maximum error with my approach is about 6050 ulp, while it is about 1500 ulp with yours (i.e., when going through logarithmic space). My initial thought was that both approaches would lose on the order of 12 bits due to error magnification, but mine would reduce the error being magnified a bit for overall better accuracy. Clearly that is not the case.

njuffa
– njuffa

2025年05月10日 22:35:52 +00:00
Commented May 10, 2025 at 22:35
@njuffa The trick bring the computation into the fat FP zone near 0.0 to 1.0, forms the pre-answer and then scale exactly by 2^n. I'll look more into it later - real life calls.

chux
– chux

2025年05月11日 01:14:00 +00:00
Commented May 11, 2025 at 1:14
1

You latest proposal results in a maximum error of 513 ulps. I am sampling double arguments in [-514165, -504776], BTW.

njuffa
– njuffa

2025年05月11日 03:16:51 +00:00
Commented May 11, 2025 at 3:16

| Show 1 more comment

How to Edit

Correct minor typos or mistakes
Clarify meaning without changing it
Add related resources or links
Always respect the author’s intent
Don’t use edits to reply to the author

How to Format

create code fences with backticks ` or tildes ~
```
like so
```
add language identifier to highlight code
```python
def function(foo):
print(foo)
```
put returns between paragraphs
for linebreak add 2 spaces at end
_italic_ or **bold**
indent code by 4 spaces
backtick escapes `like _so_`
quote by placing > at start of line
to make links (use https whenever possible)

<https://example.com>

[example](https://example.com)

<a href="https://example.com">example</a>

formatting help »
answering help »

How to Tag

A tag is a keyword or label that categorizes your question with other, similar questions. Choose one or more (up to 5) tags that will help answerers to find and interpret your question.

complete the sentence: my question is about...
use tags that describe things or concepts that are essential, not incidental to your question
favor using existing popular tags
read the descriptions that appear below the tag

If your question is primarily about a topic for which you can't find a tag:

combine multiple words into single-words with hyphens (e.g. python-3.x), up to a maximum of 35 characters
creating new tags is a privilege; if you can't yet create a tag you need, then post this question without it, then ask the community to create it for you

popular tags »

CollectivesTM on Stack Overflow

floating-point implementation of a function with removable singularity

Answer*