Sine code for ANSI C

[Christian Bau]

Re: Sine code for ANSI C

In article <40A9381C.500B5 038@yahoo.com>,
CBFalconer <cbfalconer@yah oo.com> wrote:
[color=blue]
> Christian Bau wrote:[color=green]
> >[/color]
> ... snip ...[color=green]
> >
> > Somehow your mentioning of sin (DBL_MAX) sparked a little idea
> > that might help you calculating sin (x) for arbitrarily large
> > values quite quickly, especially useful if you use long double
> > with 15 exponent bits and an enormous range:
> >
> > Instead of storing huge numbers of bits of 1/2pi and doing huge
> > precision argument reduction, just store the values of
> > sin (65536 ^ N) and cos (65536 ^ N) with lets say 20 or 24
> > extra bits of precision. ... snip ...[/color]
>
> That isn't necessary. Assuming the argument is precise (i.e.
> extended by zero bits on the right) all we need to compute (arg
> mod PI) to the precision of PI is a couple of guard bits.[/color]

When you write PI you mean a floating point approximation of the
mathematical constant pi? For large x, reduction modulo an approximation
to pi will give a result that is very much different from reduction
modulo pi. That is why some libraries store pi or 1/pi with a precision
of several thousand bits.

[Irrwahn Grausewitz]

Re: Sine code for ANSI C

Joe Wright <joewwright@com cast.net> wrote:[color=blue]
>CBFalconer wrote:
>[color=green]
>> "P.J. Plauger" wrote:
>>
>> ... snip ...[/color]
>
>I got really interested in P.J.'s argument reduction to sin() et.al.
>and while constructing something to test with I have found something
>very strange. A very simple program..
>
>#include <stdio.h>
>#include <math.h>
>
>#define R PI*2/360 /* Radians per Degree */[/color]

<snip>
[color=blue]
>If you see the same difference, can you explain it?[/color]

If you actually define PI somewhere and parenthesize the
replacement text for R, the code not only compiles, but
all the strangeness should be gone.

Regards
--
Irrwahn Grausewitz (irrwahn33@free net.de)
welcome to clc: http://www.ungerhu.com/jxh/clc.welcome.txt
clc faq-list : http://www.faqs.org/faqs/C-faq/faq/
clc OT guide : http://benpfaff.org/writings/clc/off-topic.html

[CBFalconer]

Re: Sine code for ANSI C

Christian Bau wrote:[color=blue]
> CBFalconer <cbfalconer@yah oo.com> wrote:[color=green]
> > Christian Bau wrote:[color=darkred]
> > >[/color]
> > ... snip ...[color=darkred]
> > >
> > > Somehow your mentioning of sin (DBL_MAX) sparked a little idea
> > > that might help you calculating sin (x) for arbitrarily large
> > > values quite quickly, especially useful if you use long double
> > > with 15 exponent bits and an enormous range:
> > >
> > > Instead of storing huge numbers of bits of 1/2pi and doing huge
> > > precision argument reduction, just store the values of
> > > sin (65536 ^ N) and cos (65536 ^ N) with lets say 20 or 24
> > > extra bits of precision. ... snip ...[/color]
> >
> > That isn't necessary. Assuming the argument is precise (i.e.
> > extended by zero bits on the right) all we need to compute (arg
> > mod PI) to the precision of PI is a couple of guard bits.[/color]
>
> When you write PI you mean a floating point approximation of the
> mathematical constant pi? For large x, reduction modulo an
> approximation to pi will give a result that is very much different
> from reduction modulo pi. That is why some libraries store pi or
> 1/pi with a precision of several thousand bits.[/color]

You're right. I was thinking of the case where we start with a
number of the form "arg = 2 * PI * rotations + delta", where we
can recover delta provided we use the same value for PI and the
number has not lost precision in the first place.

--
Chuck F (cbfalconer@yah oo.com) (cbfalconer@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!

[Dan Pop]

Re: Sine code for ANSI C

In <%27qc.135575$G _.78424@nwrddc0 2.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:
[color=blue]
>"Dan Pop" <Dan.Pop@cern.c h> wrote in message
>news:c8ahag$qf i$17@sunnews.ce rn.ch...[color=green]
>> In <christian.ba u-3E2624.19493215 052004@slb-newsm1.svr.pol. co.uk>[/color]
>Christian Bau <christian.bau@ cbau.freeserve. co.uk> writes:
>[color=green][color=darkred]
>> >Well, you actually _can_ find the correct answer quite well. A value of
>> >type double represents a single real number.[/color]
>>
>> But, when used in a real number context (as opposed to an integer number
>> context -- floating point can be used in both contexts) it stands for a
>> whole subset of the real numbers set. The real value exactly represented
>> is no more relevant than any other value from that set.[/color]
>
>That's just one interpretation of floating-point representation. It is
>akin to the oft-repeated presumption in this thread that every value
>carries a +/- 1/2 ulp error with it. But that's not always the case.
>Many a floating-point calculation is *exact*. So the real value exactly[/color]

When I asked you for concrete examples, you provided exactly zilch.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dan Pop]

Re: Sine code for ANSI C

In <40A8F445.DBE0F 395@yahoo.com> CBFalconer <cbfalconer@yah oo.com> writes:
[color=blue]
>Dan Pop wrote:[color=green]
>>[/color]
>... snip ...[color=green]
>>
>> The point is not whether it can be calculated, but rather how much
>> precision should the calculation produce? Does it makes sense to
>> compute sin(DBL_MAX) with 53-bit precision, ignoring the fact that
>> DBL_MAX stands for an interval so large as to make this function
>> call completely devoid of any meaning?[/color]
>
>PJPs point is that that interval is the range of values specified,
>with the actual value being the most likely value in the range.
>The shape of the probability curve in the interval is known only
>to the end user. It may be an impulse function precisely at the
>specified value.[/color]

Or it may not. If I deserve 5 digits of precision, why should the
implementation waste precious CPU resources to give me 16 digits of
*bogus* precision?

To deserve 16 digits of precision, it is *my* job to call sin() with
an argument of the right magnitude. It's as simple as that.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[P.J. Plauger]

Re: Sine code for ANSI C

"Dan Pop" <Dan.Pop@cern.c h> wrote in message
news:c8dd2g$6dg $15@sunnews.cer n.ch...
[color=blue]
> In <40A8F445.DBE0F 395@yahoo.com> CBFalconer <cbfalconer@yah oo.com> writes:
>[color=green]
> >Dan Pop wrote:[color=darkred]
> >>[/color]
> >... snip ...[color=darkred]
> >>
> >> The point is not whether it can be calculated, but rather how much
> >> precision should the calculation produce? Does it makes sense to
> >> compute sin(DBL_MAX) with 53-bit precision, ignoring the fact that
> >> DBL_MAX stands for an interval so large as to make this function
> >> call completely devoid of any meaning?[/color]
> >
> >PJPs point is that that interval is the range of values specified,
> >with the actual value being the most likely value in the range.
> >The shape of the probability curve in the interval is known only
> >to the end user. It may be an impulse function precisely at the
> >specified value.[/color]
>
> Or it may not. If I deserve 5 digits of precision, why should the
> implementation waste precious CPU resources to give me 16 digits of
> *bogus* precision?[/color]

Because it doesn't know that it shouldn't? You could supply a useful
hint by calling sinf, in this particular case at least.
[color=blue]
> To deserve 16 digits of precision, it is *my* job to call sin() with
> an argument of the right magnitude. It's as simple as that.[/color]

Simple for you, perhaps, but the interface for sin currently fails
to provide any information about how many of the bits you're supplying
should be taken seriously. Absent such information, a high quality
sin function favors getting the best possible answer for the argument
given, over doing something fast in the hopes that *all* users are
as fuzzy and in a hurry as you always seem to be. And as I've indicated
several times in this thread, computing the sine of *any* angle
smaller in magnitude than pi/4 can be done pretty fast. Odds are
that any logic you add to reduce the precision in this range would
take about as long to execute as the time you save. So all you have
to do is reduce your own arguments, as quick and sloppy as suits
your needs, and you avoid all that wasteful accuracy demanded by
professional programmers.

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[P.J. Plauger]

Re: Sine code for ANSI C

"Dan Pop" <Dan.Pop@cern.c h> wrote in message
news:c8dckl$6dg $14@sunnews.cer n.ch...
[color=blue]
> In <%27qc.135575$G _.78424@nwrddc0 2.gnilink.net> "P.J. Plauger"[/color]
<pjp@dinkumware .com> writes:[color=blue]
>[color=green]
> >"Dan Pop" <Dan.Pop@cern.c h> wrote in message
> >news:c8ahag$qf i$17@sunnews.ce rn.ch...[color=darkred]
> >> In <christian.ba u-3E2624.19493215 052004@slb-newsm1.svr.pol. co.uk>[/color]
> >Christian Bau <christian.bau@ cbau.freeserve. co.uk> writes:
> >[color=darkred]
> >> >Well, you actually _can_ find the correct answer quite well. A value[/color][/color][/color]
of[color=blue][color=green][color=darkred]
> >> >type double represents a single real number.
> >>
> >> But, when used in a real number context (as opposed to an integer[/color][/color][/color]
number[color=blue][color=green][color=darkred]
> >> context -- floating point can be used in both contexts) it stands for a
> >> whole subset of the real numbers set. The real value exactly[/color][/color][/color]
represented[color=blue][color=green][color=darkred]
> >> is no more relevant than any other value from that set.[/color]
> >
> >That's just one interpretation of floating-point representation. It is
> >akin to the oft-repeated presumption in this thread that every value
> >carries a +/- 1/2 ulp error with it. But that's not always the case.
> >Many a floating-point calculation is *exact*. So the real value exactly[/color]
>
> When I asked you for concrete examples, you provided exactly zilch.[/color]

Just last week I replied to you as follows:

: > >> So, I was asking about *concrete* examples where it makes sense to
call
: > >> sine with integral arguments or with arguments that are provably
*exact*
: > >> representations of the intended value.
: > >
: > >And I gave one.
: >
: > Nope, you gave me nothing in the way of a *concrete* example. Or maybe
: > the term is beyond your grasp... Clue: "concrete" and "hypothetic al"
: > are not exactly synonyms in any language I'm familiar with.
:
: The term is beyond my grasp.

Can't you remember *anything*?

P.J. Plauger
Dinkumware, Ltd.

http://www.dinkumware.com

[Dik T. Winter]

Re: Sine code for ANSI C

In article <c8dckl$6dg$14@ sunnews.cern.ch > Dan.Pop@cern.ch (Dan Pop) writes:[color=blue]
> In <%27qc.135575$G _.78424@nwrddc0 2.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:[/color]
....[color=blue][color=green][color=darkred]
> >> But, when used in a real number context (as opposed to an integer number
> >> context -- floating point can be used in both contexts) it stands for a
> >> whole subset of the real numbers set. The real value exactly represented
> >> is no more relevant than any other value from that set.[/color]
> >
> > That's just one interpretation of floating-point representation. It is
> > akin to the oft-repeated presumption in this thread that every value
> > carries a +/- 1/2 ulp error with it. But that's not always the case.
> > Many a floating-point calculation is *exact*. So the real value exactly[/color]
>
> When I asked you for concrete examples, you provided exactly zilch.[/color]

I think I gave one? (Zero's of the Riemann zeta function.)

But if you really do think what is quoted here with "> >>", you should
distrust *every* solution of linear systems where the system has an
order over 40 or somesuch.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[CBFalconer]

Re: Sine code for ANSI C

Dan Pop wrote:[color=blue]
> "P.J. Plauger" <pjp@dinkumware .com> writes:
>[/color]
.... snip ...[color=blue][color=green]
>>
>> That's just one interpretation of floating-point representation.
>> It is akin to the oft-repeated presumption in this thread that
>> every value carries a +/- 1/2 ulp error with it. But that's not
>> always the case. Many a floating-point calculation is *exact*.[/color]
>
> When I asked you for concrete examples, you provided exactly zilch.[/color]

I expanded my earlier example showing the errors arising from
large arguments to sin to show possible user argument reductions.
The first one I tried was fmod, and I am pleasantly surprised by
its accuracy here. The library is DJGPP 2.03.

#include <stdio.h>
#include <limits.h> /* INT_MAX */
#include <math.h> /* fmod, sin, atan */
#include <stdlib.h> /* strtod */

/* try delta 0.001953125 (i.e. 2**-n i.e. single bit)
0.00048828125
0.0001220703125
0.0000305175781 25
0.0000076293945 3125
0.0000019073485 328125
to show the precision preserved by reduce() */

/* -------------------1 */

/* return arg modulo PI */
double reduce(double arg, double PI)
{
return fmod(arg, PI);
} /* reduce */

/* -------------------1 */

int main(int argc, char **argv)
{
double PI = 4.0 * atan(1.0);
unsigned int i;
double arg, narg, delta;

delta = 0.0;
if (argc > 1) delta = strtod(argv[1], NULL);

for (i = 1; i < (INT_MAX / 2U); i = 4 * i) {
arg = i * PI + delta;
narg = reduce(arg, PI);
printf("%11d %28.16f %20.16f\n%11c %28.16f %20.16f\n",
i, arg, sin(arg), ' ', narg, sin(narg));
}
return 0;
}

c:\c\junk>a 0.0000019073485 328125
1 3.1415945609383 260 -0.0000019073485 328
0.0000019073485 329 0.0000019073485 329
4 12.566372521707 7058 0.0000019073485 328
0.0000019073485 333 0.0000019073485 333
16 50.265484364785 2232 0.0000019073485 314
0.0000019073485 333 0.0000019073485 333
64 201.06193173709 52785 0.0000019073485 113
0.0000019073485 191 0.0000019073485 191
256 804.24772122633 55568 0.0000019073484 878
0.0000019073485 191 0.0000019073485 191
1024 3216.9908791832 967836 0.0000019073485 074
0.0000019073486 328 0.0000019073486 328
4096 12867.963511011 1412359 0.0000019073481 312
0.0000019073486 328 0.0000019073486 328
16384 51471.854038322 5190453 0.0000019073466 264
0.0000019073486 328 0.0000019073486 328
65536 205887.41614756 80302829 0.0000019073406 072
0.0000019073486 328 0.0000019073486 328
262144 823549.66458455 00752330 0.0000019073165 305
0.0000019073486 328 0.0000019073486 328
1048576 3294198.6583324 782550335 0.0000019072202 235
0.0000019073486 328 0.0000019073486 328
4194304 13176794.633324 1909742355 0.0000019068349 957
0.0000019073486 328 0.0000019073486 328
16777216 52707178.533291 0418510437 0.0000019052940 843
0.0000019073486 328 0.0000019073486 328
67108864 210828714.13315 84453582764 0.0000018991304 387
0.0000019073486 328 0.0000019073486 328
268435456 843314856.53262 80593872070 0.0000018744758 563
0.0000019073486 328 0.0000019073486 328


--
Chuck F (cbfalconer@yah oo.com) (cbfalconer@wor ldnet.att.net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home .att.net> USE worldnet address!

[glen herrmannsfeldt]

Re: Sine code for ANSI C

Dik T. Winter wrote:

[color=blue]
> But if you really do think what is quoted here with "> >>", you should
> distrust *every* solution of linear systems where the system has an
> order over 40 or somesuch.[/color]

I would certainly look carefully at them, especially
if done with matrix inversion.

-- glen

[Dik T. Winter]

Re: Sine code for ANSI C

In article <TCCqc.23014$gr .1960915@attbi_ s52> glen herrmannsfeldt <gah@ugcs.calte ch.edu> writes:[color=blue]
> Dik T. Winter wrote:
>
>[color=green]
> > But if you really do think what is quoted here with "> >>", you should
> > distrust *every* solution of linear systems where the system has an
> > order over 40 or somesuch.[/color]
>
> I would certainly look carefully at them, especially
> if done with matrix inversion.[/color]

Also with Gauss' elimination. If you really think that the input
numbers are intervals, in almost all cases the result makes no sense.
It is quite possible that within the set of matrices given, there is
at least one singular matrix. The number 40 I use above comes from
von Neuman. And indeed, before Wilkinson, numeric mathematicians could
not deal with the error analysis of large linear systems.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/

[Dan Pop]

Re: Sine code for ANSI C

In <V2tqc.138401$G _.102626@nwrddc 02.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:
[color=blue]
>"Dan Pop" <Dan.Pop@cern.c h> wrote in message
>news:c8dd2g$6d g$15@sunnews.ce rn.ch...
>[color=green]
>> In <40A8F445.DBE0F 395@yahoo.com> CBFalconer <cbfalconer@yah oo.com> writes:
>>[color=darkred]
>> >Dan Pop wrote:
>> >>
>> >... snip ...
>> >>
>> >> The point is not whether it can be calculated, but rather how much
>> >> precision should the calculation produce? Does it makes sense to
>> >> compute sin(DBL_MAX) with 53-bit precision, ignoring the fact that
>> >> DBL_MAX stands for an interval so large as to make this function
>> >> call completely devoid of any meaning?
>> >
>> >PJPs point is that that interval is the range of values specified,
>> >with the actual value being the most likely value in the range.
>> >The shape of the probability curve in the interval is known only
>> >to the end user. It may be an impulse function precisely at the
>> >specified value.[/color]
>>
>> Or it may not. If I deserve 5 digits of precision, why should the
>> implementation waste precious CPU resources to give me 16 digits of
>> *bogus* precision?[/color]
>
>Because it doesn't know that it shouldn't? You could supply a useful
>hint by calling sinf, in this particular case at least.[/color]

sinf is useless in portable C code.
[color=blue][color=green]
>> To deserve 16 digits of precision, it is *my* job to call sin() with
>> an argument of the right magnitude. It's as simple as that.[/color]
>
>Simple for you, perhaps, but the interface for sin currently fails
>to provide any information about how many of the bits you're supplying
>should be taken seriously.[/color]

ALL of them, as I have already explained you. It's the magnitude of the
value that decides how many digits of precision I deserve in the result.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dan Pop]

Re: Sine code for ANSI C

In <G6tqc.138420$G _.93237@nwrddc0 2.gnilink.net> "P.J. Plauger" <pjp@dinkumware .com> writes:
[color=blue]
>"Dan Pop" <Dan.Pop@cern.c h> wrote in message
>news:c8dckl$6d g$14@sunnews.ce rn.ch...
>[color=green]
>> In <%27qc.135575$G _.78424@nwrddc0 2.gnilink.net> "P.J. Plauger"[/color]
><pjp@dinkumwar e.com> writes:[color=green]
>>[color=darkred]
>> >"Dan Pop" <Dan.Pop@cern.c h> wrote in message
>> >news:c8ahag$qf i$17@sunnews.ce rn.ch...
>> >> In <christian.ba u-3E2624.19493215 052004@slb-newsm1.svr.pol. co.uk>
>> >Christian Bau <christian.bau@ cbau.freeserve. co.uk> writes:
>> >
>> >> >Well, you actually _can_ find the correct answer quite well. A value[/color][/color]
>of[color=green][color=darkred]
>> >> >type double represents a single real number.
>> >>
>> >> But, when used in a real number context (as opposed to an integer[/color][/color]
>number[color=green][color=darkred]
>> >> context -- floating point can be used in both contexts) it stands for a
>> >> whole subset of the real numbers set. The real value exactly[/color][/color]
>represented[color=green][color=darkred]
>> >> is no more relevant than any other value from that set.
>> >
>> >That's just one interpretation of floating-point representation. It is
>> >akin to the oft-repeated presumption in this thread that every value
>> >carries a +/- 1/2 ulp error with it. But that's not always the case.
>> >Many a floating-point calculation is *exact*. So the real value exactly[/color]
>>
>> When I asked you for concrete examples, you provided exactly zilch.[/color]
>
>Just last week I replied to you as follows:
>
>: > >> So, I was asking about *concrete* examples where it makes sense to
>call
>: > >> sine with integral arguments or with arguments that are provably
>*exact*
>: > >> representations of the intended value.
>: > >
>: > >And I gave one.
>: >
>: > Nope, you gave me nothing in the way of a *concrete* example. Or maybe
>: > the term is beyond your grasp... Clue: "concrete" and "hypothetic al"
>: > are not exactly synonyms in any language I'm familiar with.
>:
>: The term is beyond my grasp.
>
>Can't you remember *anything*?[/color]

Of course I do. Otherwise, I couldn't point out that you have still
failed to support your assertion with concrete examples, therefore I
refuse to consider it seriously, period.

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dan Pop]

Re: Sine code for ANSI C

In <HxyHxo.Hxo@cwi .nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:
[color=blue]
>In article <TCCqc.23014$gr .1960915@attbi_ s52> glen herrmannsfeldt <gah@ugcs.calte ch.edu> writes:[color=green]
> > Dik T. Winter wrote:
> >
> >[color=darkred]
> > > But if you really do think what is quoted here with "> >>", you should
> > > distrust *every* solution of linear systems where the system has an
> > > order over 40 or somesuch.[/color]
> >
> > I would certainly look carefully at them, especially
> > if done with matrix inversion.[/color]
>
>Also with Gauss' elimination. If you really think that the input
>numbers are intervals, in almost all cases the result makes no sense.[/color]

Unless the input numbers can be exactly represented in floating point,
it is sheer foolishness to consider them as exact values. Furthermore,
in real life applications, most input numbers are obtained with a
precision even lower than that provided by the floating point
representation, so they stand for even larger intervals. That's why
the results are intervals, too. Consider the difference between
1 +/- 1000 and 1 +/- 0.001. Only a patent fool will display the result
as a plain 1 in both cases.


Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Dan.Pop@ifh.de

[Dik T. Winter]

Re: Sine code for ANSI C

In article <c8fl7p$7u4$11@ sunnews.cern.ch > Dan.Pop@cern.ch (Dan Pop) writes:[color=blue]
> In <HxyHxo.Hxo@cwi .nl> "Dik T. Winter" <Dik.Winter@cwi .nl> writes:[/color]
....[color=blue][color=green][color=darkred]
> > > I would certainly look carefully at them, especially
> > > if done with matrix inversion.[/color]
> >
> >Also with Gauss' elimination. If you really think that the input
> >numbers are intervals, in almost all cases the result makes no sense.[/color]
>
> Unless the input numbers can be exactly represented in floating point,
> it is sheer foolishness to consider them as exact values. Furthermore,
> in real life applications, most input numbers are obtained with a
> precision even lower than that provided by the floating point
> representation, so they stand for even larger intervals. That's why
> the results are intervals, too. Consider the difference between
> 1 +/- 1000 and 1 +/- 0.001. Only a patent fool will display the result
> as a plain 1 in both cases.[/color]

Ah, so, what method do *you* use to solve linear systems on the computer?
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924131
home: bovenover 215, 1025 jn amsterdam, nederland; http://www.cwi.nl/~dik/