Implementation of the log1p (x) function (continued)

Source: Internet
Author: User

At that time, a GSL log1p (x) code implementation was provided. However, I didn't want to understand why it was like that. I suddenly understood it when I read a post on shuimu recently.


The Code on GSL is as follows:


[Cpp] view plaincopyprint? Double gsl_log1p (const double x)
{
Volatile double y, z;
Y = 1 + x;
Z = y-1;
Return log (y)-(z-x)/y;/* cancels errors with IEEE arithmetic */
}

Double gsl_log1p (const double x)
{
Volatile double y, z;
Y = 1 + x;
Z = y-1;
Return log (y)-(z-x)/y;/* cancels errors with IEEE arithmetic */
}


We know that floating point numbers on computers only have limited precision.

Therefore, the value assignment statement "y: = 1 + x" is not exactly equal to 1 + x when x is calculated in hours.

Or we can think that what we actually calculate is

Y: 1 + x'

X' is the result of the loss of the Number of valid digits of x in the decimal place when the value of 1 + x is calculated.

Therefore, directly calculating log (1 + x) actually calculates log (1 + x ').

Since x is very similar to x', the following formula is similar:

 

 
 

The above formula is changed to the following:

 

This is the code used on GSL (x' corresponds to z in the Code ).

The following statement cannot be written in a C program:

Z = 1 + X-1;

Otherwise, the compiler will intelligently optimize it to z = x, and it does not seem to be able to turn off this optimization through the compiler's compilation options.

In addition, there is a very famous article, What Every Computer Scientist shocould Know About Floating-Point Arithmetic, which also discusses this computing. Another calculation method is provided.

 


This algorithm is based on the fact that when x is very close to x,

 

The value of x (x'-x) is very small, so the error of this algorithm is very small.


 

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.