2
$\begingroup$

As commonly known for programming on computer, if x is a finite precision floating-point number such as double/float in C language, arcsin(sin(x)) is usually not equal to x due to the numerical issue. I'm wondering whether there is a way to design a sin function and an arcsin function such that arcsin(sin(x))=x where x is an fixed-size multi-precision floating-point number. Ideally, both sin and arcsin function would have error guarantee.

On solution would be to approximate the sin using linear piecewise function. However, it would take a lot of memory to store all functions if we want to achieve a very accurate approximation.

$\endgroup$
2
  • 2
    $\begingroup$ I suppose you are restricting yourself to a specific domain such as $[-\pi/2,\pi/2]$, because otherwise you'd have to deal with the fact that the sine is not injective. It isn't impossible, but it would surely return a wildly discontinuous arcsin() function, since values close to 0 have to be mapped back to 0, $\pi$, $2\pi$, ... $\endgroup$ Mar 21, 2023 at 13:41
  • $\begingroup$ Yes, x is only in $[-\pi/2, \pi/2]$ $\endgroup$ Mar 21, 2023 at 13:45

2 Answers 2

8
$\begingroup$

That seems impossible to do with a small numerical error, because of cardinality reasons. Consider for instance the much simpler case of the function $f(x) = x^2$ over $[0,1]$. This function maps $[0,1/2]$ into $[0,1/4]$. But the floating-point numbers in $[0,1/4]$ are much fewer than the numbers in $[0,1/2]$, and hence this is impossible to do.

To turn this into a mathematical argument: assume that your number system contains $N$ numbers in $[0,1/2)$, the number $1/2$ (which we assume to be an exact floating-point number, for simplicity) and $M$ numbers in $(1/2,1]$. Consider your implementation $g(x)$ which is supposed to approximate the map $f(x) = x^2$. Since you want it to be injective, the $N+1$ numbers in $[0,1/2]$ must be mapped to $N+1$ different floating point numbers. In particular, at least one of them must be in $[1/2,1]$, since there are only $N$ floating-point numbers inside $[0,1/2)$. So there must be a number $a\in [0,1/2]$ for which $g(a) \geq 1/2$, which is a very large error to have since $f(a) \leq 1/4$.

It should be simple to adapt this argument to the sine function.

$\endgroup$
8
  • $\begingroup$ Thanks for the answer! I'm not sure if it is the right place to ask a follow-up question here or if I should open a new question. I would appreciate if you could answer. What if x is a rational number with infinite precision, is there way to design a sin and an arcsin function such that arcsin(sin(x))=x? $\endgroup$ Mar 21, 2023 at 18:10
  • 2
    $\begingroup$ Yes, you can construct such maps. I have a construction in mind but it's not particularly pretty as it is based on set-theoretic ideas. The idea is that you can partition the domain X = [-pi/2, pi/2] into intervals $I_1, I_2, \dots, I_n$ of width $\varepsilon$, and map in any way the rationals in the interval $I_k$ into those in the interval $f(I_k)$ bijectively. This should give you a map that is close to the sine in the L_inf distance and is a bijection on the rationals. $\endgroup$ Mar 21, 2023 at 18:32
  • $\begingroup$ But note that this is more of a mathematical curiosity and it is not particularly useful for computational purposes. $\endgroup$ Mar 21, 2023 at 18:34
  • 3
    $\begingroup$ @user3677630 The pigeon hole principle this answer is based on points out the problem when the inner function is contracting while the outer function is expanding. However, it can work fine the other way around: For example, in the absence of overflow/underflow, IEEE-754 arithmetic guarantees that $\sqrt{x^2}=x$ while this cannot be guaranteed for $(\sqrt{x})^2$. For your use case, could you reverse the order of operations, that is, use $\sin^{-1} x$ as the inner function? $\endgroup$
    – njuffa
    Mar 22, 2023 at 6:49
  • 2
    $\begingroup$ @user3677630 For some $f$. Where "some" is likely "few". I would not know how to prove that it works in a particular case, but some people do: "We have proved that, using radix $2$ and a precision greater than $1$, the rounded computation of $\sqrt{a^2}$ is indeed $|a|$." You might want to try something like $\ln(\exp(x))$ first as an easier case and see what you can prove about it. Note that correctly-rounded transcendental functions fall only under the optional provisions of IEEE-754 and there are few libraries providing them (see crlibm). $\endgroup$
    – njuffa
    Mar 24, 2023 at 5:18
0
$\begingroup$

I think it is possible. You would need two ingredients: very accurate sin function and very accurate arcsin function. How accurate? Take this with grain of salt (I don't have the proof!), but if x is binary64 (aka double in C), it is very likely that triple double would enough, to be (almost) sure quad double. If x is binary32 (aka float in C), using double double should do the trick. Keep in mind that such functions are going to be very, very slow and relatively hostile to modern CPUs due huge number of data dependencies.

You can try with simpler functions like polynomials (Taylor style) in Gappa and see how increasing number of bits reduces evaluation error and in some cases completely eliminates it. It goes something like: use higher precision for everything, except final rounding to for example double and compare that to true mathematical result. Error will be much smaller than when double is used for everything, for same coefficients on same interval.

$\endgroup$
3
  • 1
    $\begingroup$ Re "Error will be much smaller than when double is used for everything, for same coefficients on same interval" It is not unusual for math libraries to deliver simple functions faithfully rounded, i.e. the maximum error is bounded by 1 ulp, without resorting to multi-precision computation. For example, for the Intel math library the maximum error in sin and asin is < 0.55 ulp. The best one can do is deliver correctly-rounded result with maximum error <= 0.5 ulp. Is that much smaller error? I would call it incrementally smaller. $\endgroup$
    – njuffa
    Jan 26 at 10:04
  • $\begingroup$ To get errors less than 1ulp at least some (one or two) double-double can do wonders. But, that assumes Horner or Horner-like evaluation scheme. With Estrin's scheme error is larger. With extended precision it becomes much like noise. Intel had done something like that with Itanium DE. It was just about evaluation. If yo don't believe me, try approximation for some function, I usually start with exponential, difference between Horner and Estrin is significant. You made implicit assumption that for sin and arcsin input value and return value are double, in this case that is not mandatory. $\endgroup$ Jan 30 at 3:03
  • $\begingroup$ Estrin's scheme is okay for higher part, in combination with some double-double it gives both fast and accurate results. Good to exploit deep FMA pipelines too. It has some other handy applications. My point was that increased precision can make rounding errors to vanish even for Estrin's scheme in it's raw form and thus it is good example. Gappa offers fine control so it is easier to observe effects. To return to the question: specially designed functions could do that, somewhat slower with multi word arithmetic, still faster than arbitrary precision and without external dependencies. $\endgroup$ Jan 30 at 14:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.