How Matlab computes trigonometric functions?

When I learned the CORDIC algorithm, I had a question to Matlab.

How Matlab computes trigonometric function?

for e.g. in case of cosd(27)

1. Does answer evaluate by Taylor series using floating point computation or by Look Up Table from Taylor?

If both are not, what algorithm is used for computation?

2. If it uses Taylor series, what is the last order of Taylor series? (Just N that satisfies fractional bit of floating point?)

3. What is the accuracy of trigonometric function?

(I want to know whether Matlab gives 52(53) bit fractional accuracy using (some) nearest rounding in double precision value)

[dpb]

In general, TMW doesn't comment on such internals but for library functions such as this, I'd expect Matlab simply relies on the underlying compiler standard libraries rather than reinventing the wheel...

Answer 1

[John D'Errico]

NO. MATLAB does not do a table lookup.

To compute things like a cos or sin, one first does a range reduction, to ensure the argument is in a small interval around some point where an approximation is highly accurate. That is usually around zero.

Having done that, then one typically uses a series approximation of some sort. Or you might use a Pade approximation for some functions. See the text "Computer approximations", by J.F. Hart, et al. Alternatively, do some light reading through Abramowitz & Stegun.

Do they guarantee accuracy down to the least significant bit? This is probably not doable, at least not doable in an efficient manner. You simply use an approximation that gives a slightly higher accuracy than necessary over the entire reduced domain.

It has been a few years since I wrote HPF, but it looks like I did it as a range reduction into the interval [-pi,pi]. Symmetry considerations can further reduce it to [0,pi]. Then I broke that interval down further, with break points at pi/8, 3*pi/8, 5*pi/8, 7*pi/8. The net result was always an effective range reduction into the interval [0,1/8]*pi. (I'm sure I could have used other approaches. I just chose one that worked, and was as efficient as possible.)

Within that final interval, I used a series expansion, running backwards for accuracy and efficiency. (Thus, start with the small terms. That has advantages for several reasons.) A Pade approximation would not have worked, because HPF is a variable precision tool. When you know always exactly how many digits you need, and that is a fixed number, a Pade approximation may often be a better choice, IF one is available.

As far as cosd is concerned, you would NEVER use a different scheme than cos. Just transform the problem.

cosd(x) = cos(x*pi/180)

A basic rule: NEVER write the same code twice if a direct transformation is available.

Answer 2

[Steven Lord]

You may be interested in this article from Cleve.

Answer 3

[Walter Roberson]

cosd() checks for some special cases, and if it is not one of the special cases then it transforms to radians and invokes cos()

cos() uses the underlying hardware instructions.

You might be interested in the discussion at https://members.loria.fr/PZimmermann/papers/decimalexp.pdf which gives some reasons why it is difficult to round such functions correctly.