Efficiency of math.h

[Dik T. Winter]

> Christian Bau wrote:
>
>
> So... you're a newbie eh?

Hrm. The first language I used, Algol 60, did not even *have* the
arccosine as a function. It had only sin, cos and arctan from this
group. Only when I started programming in Fortran did I see the
ACOS as a standard function, but have actually never needed it.
Done quite a bit of numerical programming actually. Like implementing
an arccosine function, but never used it ;-).

 

[Dik T. Winter]

> Christian Bau wrote:
>
>
> Please show us the formula.

You mean:
T_0(x) = 1
T_1(x) = x
T_n(x) = 2x.T_{n-1}(x) - T_{n-2}(x) for n > 1
? (This from memory...)

 

[E. Robert Tisdale]

1. Calculating Taylor coefficients on the fly requires division,

Nonsense! Multiply the Taylor coefficients through
by the least common denominator, sum using Horner's rule
then multiply the result by the inverse of the least common denominator.
All of the coefficients can be calculated using integer multiplication.

which is on most processors much slower than multiplication.

2. The number of coefficients required in a Tschebychev approximation
is less than the number of Taylor coefficients needed.

By one or, at best, two terms.

For instance, to get a good sine in double precision on [0, pi/2],
you can use only 6 terms
of the Tschebychev expansion of (sin(x)/x - 1).
(This formula to get good relative precision near x = 0.)

Let's suppose that the Taylor series requires 8 terms
to obtain the minimum precision everywhere.
The Tschebychev approximation could be, at best, 25% faster.
But it could easily require much more time to load
the high precision Tschebychev coefficients from memory.
But don't take my word for it.
Implement both and benchmark them against each other.
The winner will be determined by your computer architecture.

 

[E. Robert Tisdale]

E.Robert.Tisdale writes:

The first language I used, Algol 60,

Me too. On a CDC 6600.

did not even *have* the arccosine as a function.
It had only sin, cos and arctan from this group.
Only when I started programming in Fortran did I see the
ACOS as a standard function, but have actually never needed it.
Done quite a bit of numerical programming actually.

I only have a little experience

http://csdl.computer.org/comp/proceedings/arith/1995/7089/00/70890163abs.htm

 

[Nick Landsberg]

I think the point of Christian's posting was that
he would like more information from the [Dave Rudolf]
[about] why the acos() function was called so often.
If it is something that requires high precision
and MUST be called often
(as in a signal processing application),
then it is one thing.

If it does not require high precision
or [high precision] is used gratuitously,
there might be better ways to accomplish the same thing.
As one poster pointed out, a rough [interpolation]
between values in a static array may satisfy the requirement.
Without knowing the degree of precision needed,
we just don't know.

[You misunderstood at least one part in your paraphrase
of my post. In the second paragraph I meant ...
"If it does not require high precision or *acos()* is
used gratuitiouly, ..." not, as you have interpreted it.
But, re-reading it, I can see how the sentence
could be so interpreted. I also did not put []'s
around the word "interpolation."]

Unfortunately,
I don't know what degree of precision is required either.
We can only hope that Dave Rudolf is still following this thread
and will reply with the answer.

Then why attack various proposed algorithms? Without knowing
the problem the poster is trying to solve, why argue about
the details of the solution?

But I'll tell you what you can do without knowing any more
about Dave Rudolf's precision requirements.
You can implement you proposed interpolation scheme
and benchmark it against the acos provided with your C compiler.
I predict that you will find that your interpolation scheme
is *slower* even if you need just a few decimal digits precision.
Please post your implementation so that we can repeat
your benchmark experiments ourselves.

I did not propose the interpolation scheme, some other
poster did, I believe it was Eric Sosman. Since you are the
one questioning the suggestion, may I suggest you post your
benchmark results first to prove your contention?

There is a reason why math functions like acos
are included in the standard library.
There is *no* single portable implementation
that will deliver performance and accuracy everywhere.
Compiler developers hire experts to customize implementations
for each target platform and, usually, they are hard to beat.

All that is true, for applications which require precision.
It is not clear that this application requires that precision.
We come back to full circle regarding:
1 - what is the precision required for this application?
2 - is there a better way rather than the mathematically precise
way to get a "good enough" answer?

And this whole digression belongs in something like

comp.algorithms.philosophy

rather than a discussion of the language, IMO.

 

[Guillaume]

Only an idiot would calculate them using the acos () function.

Besides, that would be stupid to re-compute the window function
every time you have to use it: there is most likely some means
of computing it once, place the discret values in an array and
then use this. After all, we are dealing with discret algorithms.

Unless the number of points was gigantic.

 

[Dik T. Winter]

> Dik T. Winter wrote:
>
>
> Nonsense! Multiply the Taylor coefficients through
> by the least common denominator, sum using Horner's rule
> then multiply the result by the inverse of the least common denominator.
> All of the coefficients can be calculated using integer multiplication.

Pray show the algorithm you intend with Taylor coefficients 1/(2k+1)!.
I have no idea what you mean with the above.

>
> By one or, at best, two terms.

As typically the Tschebychev coefficients are about 2^(-n) times the
original coefficients, I wonder where you get your one or two from.

> > For instance, to get a good sine in double precision on [0, pi/2],
> > you can use only 6 terms
> > of the Tschebychev expansion of (sin(x)/x - 1).
> > (This formula to get good relative precision near x = 0.)

>
> Let's suppose that the Taylor series requires 8 terms
> to obtain the minimum precision everywhere.
> The Tschebychev approximation could be, at best, 25% faster.

OK, the Tschebyshev approximation in the above case:
double dc[] = {-0.16666666666666662966,
0.00833333333333094971,
-0.00019841269836759707,
0.00000275573161030613,
-0.00000002505113193827,
0.00000000015918135277};
The algorithm reads (y = x*x):
sin = ((((((dc[5] * y + dc[4]) * y + dc[3]) * y + dc[2]) * y + dc[1]) * y
+ dc[0]) * y + 1.0) * x;
(Actually this is for x in [-pi/2,pi/2] with a relative precision of about
1 ULP, or, in tangible terms, the relative error is bounded by about 2e-16.)

> But it could easily require much more time to load
> the high precision Tschebychev coefficients from memory.
> But don't take my word for it.
> Implement both and benchmark them against each other.

Now implement your algorithm, without loading high precision coefficients
from memory or using division...

 

[Mike Wahler]

So... you're a newbie eh?

I also have been using C for a long time (about 17 years).
I've never used *any* of the 'math.h' functions in my
production code (well, there might be one or two I've forgotten
about, but you get the idea). That's because my application
domains never needed them. So am I a 'newbie' too?

-Mike

 

[E. Robert Tisdale]

E.Robert.Tisdale writes:

Pray show the algorithm you intend with Taylor coefficients 1/(2k+1)!.
I have no idea what you mean with the above.

Suppose that the highest order coefficient is 1/(2n+1)!. Compute

[( ... ((x^2 + (2n+1)!/(2n-1)!)*x^2 + (2n+1)!/(2n-3)!)*x^2 + ... + 1)*x]
*[1/(2n+1)!]

The first coefficient

(2n+1)!/(2n-1)! = (2n+1)*2n

The second coefficient

(2n+1)!/(2n-3)!) = ((2n+1)!/(2n-1)!)*(2n-1)*(2n-2)

etc.

For small n, you can compute these coefficients
in an integer register on-the-fly.
All you need to do is multiply the Horner product through by

[1/(2n+1)!]

at the end.

 

[E. Robert Tisdale]

OK, the Tschebyshev approximation in the above case:
double dc[] = {-0.16666666666666662966,
0.00833333333333094971,
-0.00019841269836759707,
0.00000275573161030613,
-0.00000002505113193827,
0.00000000015918135277};
The algorithm reads (y = x*x):
sin = ((((((dc[5] * y + dc[4]) * y + dc[3]) * y + dc[2]) * y + dc[1]) * y
+ dc[0]) * y + 1.0) * x;
(Actually this is for x in [-pi/2, pi/2] with a relative precision of about
1 ULP, or, in tangible terms, the relative error is bounded by about 2e-16.)

That's only about 4 or 5 decimal digits. You will probably find that
your algorithm would benefit from a "range reduction".

Anyway, Have you tried to "benchmark" your implementation
against the sin(x) function in your standard C library.
If so, please report your results
including the architecture and operating system
where you ran the benchmarks.

 

[E. Robert Tisdale]

I also have been using C for a long time (about 17 years).
I've never used *any* of the 'math.h' functions
in my production code (well, there might be one or two
that I've forgotten about, but you get the idea).
That's because my application domains never needed them.

So am I a 'newbie' too?

To numerical computing evidently.

 

[Mac]

Nick Landsberg wrote:

[snip]

If it does not require high precision
or [high precision] is used gratuitously,
there might be better ways to accomplish the same thing.
As one poster pointed out, a rough [interpolation]
between values in a static array may satisfy the requirement.
Without knowing the degree of precision needed,
we just don't know.

[You misunderstood at least one part in your paraphrase
of my post. In the second paragraph I meant ...
"If it does not require high precision or *acos()* is
used gratuitiouly, ..." not, as you have interpreted it.
But, re-reading it, I can see how the sentence
could be so interpreted. I also did not put []'s
around the word "interpolation."]

For whatever reason, E. Robert Tisdale modifies posts he is replying to on
occasion. Tisdale has not, to my knowledge, offered any explicit
explanation for this behavior.


-Mac

 

[Mike Wahler]

To numerical computing evidently.

Well, yes, but so what? Numerical computing is only one
of a virtually limitless number of domains in which
C is used. Accounting/business/database applications have been
putting bread on my table for decades.

Here's what I was responding to:

Since I have actually not used the acos() function _once_
in over twenty years of C programming,

You replied:

So... you're a newbie eh?

The only conclusion I could draw from context was that you're
calling Christian a 'newbie' with C. IMO this newsgroup's
archives provide much evidence to the contrary.

-Mike

 

[E. Robert Tisdale]

For whatever reason,
E. Robert Tisdale modifies posts he is replying to on occasion.
Tisdale has not, to my knowledge,
offered any explicit explanation for this behavior.

What's to explain?
There is no need to quote and it is discouraged.
It is as evil as top posting.
I fix capitalization, punctuation, spelling and grammar errors
to reflect my understanding of what was said.
I am pretty sure that Nick Landsberg meant interpolation
when he said extrapolation and that he meant Dave Rudolf
when he said OP.

I also quietly fix style and programming errors in other people's code
if they aren't essential to understanding the point in question.

 

[E. Robert Tisdale]

The only conclusion I could draw from context was that
you're calling Christian a 'newbie' with C.

Yes. At 20 years, he is just getting started.

This newsgroup's archives provide much evidence to the contrary.

Get a life Mike. And a sense of humor to make it bearable

 

[CBFalconer]

What's to explain?
There is no need to quote and it is discouraged.
It is as evil as top posting.
I fix capitalization, punctuation, spelling and grammar errors
to reflect my understanding of what was said.
I am pretty sure that Nick Landsberg meant interpolation
when he said extrapolation and that he meant Dave Rudolf
when he said OP.

I also quietly fix style and programming errors in other people's code
if they aren't essential to understanding the point in question.

No, you misquote in order to impose your own misunderstandings.
If you have some different understanding, write what that is. We
can then be amused at the prattling of children, rather than
outraged.

 

[Christian Bau]

Have you every implemented and tested this "algorithm"?
If so please publish it here so that we can "benchmark" it
against our C library implementations of acos. I'll bet you $1.00 that
I can find one that beats your implementation.

Just recently I read something about the typical troll's methods: Demand
documentation, references, everything from others. Always refuse to
provide the same yourself.

So you ask me to invest considerable amount of time to implement an
algorithm for acos, then you do a web search, and all for one dollar?
Asshole.

Lets do it the other way round: You publish ANY implementation of acos
in C, and I bet you UK?20,000 that I can beat it. I have a Macintosh at
home, so we will run speed tests on a PowerPC with a Macintosh G4, using
CodeWarrior 7 for MacOS as the compiler. 64 bit double precision, error
less than 1ulp for the whole range are required, so the only thing to be
measured is speed.

Twenthy thousand pounds makes it worth my while, hurts you enough, and I
am quite confident that I can beat anything you can find.

 

[Christian Bau]

Please show us the formula.

Google for "Dolph Chebychev". Idiot.

 

[Christian Bau]

Yes. At 20 years, he is just getting started.

After 20 years, I am still learning. Your posts make you look as if you
stopped learning 30 years ago.

Get a life Mike. And a sense of humor to make it bearable

Another typical troll/bully tactics: When called to order, accuse people
of having no sense of humour.

 

[Richard Heathfield]

What's to explain?
There is no need to quote and it is discouraged.

Appropriate quotation of context is encouraged.

It is as evil as top posting.
No.

I fix capitaliation, punctuation, spelling and grammar errors[.]

I have redacted the above line to fit in with my own personal tastes. Note
that I have used standard editorial notation to mark where changes have
been made. I suggest that, if you must edit quotations, you adopt a similar
strategy.

to reflect my understanding of what was said.
I am pretty sure that Nick Landsberg meant interpolation
when he said extrapolation and that he meant Dave Rudolf
when he said OP.

I also quietly fix style and programming errors in other people's code
if they aren't essential to understanding the point in question.

Even people who are widely recognised to be competent C programmers don't do
this. Your grasp of C has proved to be a little shaky, so it might be wiser
to refrain from introducing what might turn out to be your bugs in what you
claim is other people's code.