Fortran - sin(x) for large x

Author

sin(x) for large x

glen herrmannsfeldt

2004-05-14, 6:31 pm

There is a discussion in comp.lang.c titled
Sine code for ANSI C

regarding the evaluation of the sin() function for large
arguments. It might be more applicable in this
group, so I thought I would ask here.

What result do people expect from the sin() function for
very large arguments, for example sin(1e100)?

Hopefully the discussion won't get too nasty, but it
does seem to be a contentious issue over there.

-- glen

Dr Chaos

2004-05-14, 6:31 pm

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> There is a discussion in comp.lang.c titled
> Sine code for ANSI C
>
> regarding the evaluation of the sin() function for large
> arguments. It might be more applicable in this
> group, so I thought I would ask here.
>
> What result do people expect from the sin() function for
> very large arguments, for example sin(1e100)?

sin(1e100) is a symptom of a mistake.

I believe that if one's program depends on evaluating that to any
degree of accuracy, it is computing something in a mathematically
mistaken way.

Generally you need some paper-and-pencil analysis and thinking to
avoid this problem.

> Hopefully the discussion won't get too nasty, but it
> does seem to be a contentious issue over there.

Sometimes computer scientists have an odd idea what it means to be
"correct", e.g. imagining that the particular rules of IEEE floating
point define "correctness" as opposed to solving the problem in a
satisfactory way.

Nevertheless, in a legalistic fashion, one could define it thus:
sin(1.0D100) = sin(x) where x is the closest IEEE representation, in
double precision floating point, to the number 10^100.

Now, every IEEE floating point number is a rational number mathematically.
That has an exact value of sin, and you can define sin(x) as the
IEEE floating point number which is closest to that exact number.

In single precision, 1e100 doesn't have a representation I believe.

Now this definition is logical, but it may make implementation
far too difficult, and waste computational time on the cases that
almost everybody cares about, i.e. small arguments.

> -- glen

Dick Hendrickson

2004-05-14, 7:31 pm

glen herrmannsfeldt wrote:
> There is a discussion in comp.lang.c titled
> Sine code for ANSI C
>
> regarding the evaluation of the sin() function for large
> arguments. It might be more applicable in this
> group, so I thought I would ask here.
>
> What result do people expect from the sin() function for
> very large arguments, for example sin(1e100)?
>
> Hopefully the discussion won't get too nasty, but it
> does seem to be a contentious issue over there.
>
> -- glen
>
Since it's getting late on a Friday and the beer isn't
quite chilled yet, I'll try the obviously correct answer.
"I wouldn't expect anything useful."

I think there's two ways to look at this.

We know mathematically
Sin(x+2*N*pi) = Sin(x)
But once N is large we can't computationally recover
the value of x, it gets lost in the truncation to finite
representation. And then Sin(x+2N*pi) is independent of x.

But, if we look at the problem in a different way,
and if I could remember my high school trig, I'm sure I
could expand Sin(2.0) into some nice function of
Sin(1.0) (Please don't tell me the function,
I've only got a 12 pack). And by doing the evaluation
in double precision and pulling it back to single
we could probably compute sin(2.0) to within one bit.
Then it's a simple matter to compute Sin(3.0) from the
pretty accurate Sin(2.0) and Sin(1.0), again I think to
one bit. By working our way up, wouldn't we eventually
compute sin(1.e100) from a few accurately known values
to a high degree of accuracy? But, this contradicts
the first argument.

Maybe the expectaion depends on where the 1.e100 comes
from. If it comes from a long series of addition, then
you could reasonably expect to work your way up the
sequence, effectively taking out the 2*pi each time you
go around the circle. But if it comes from dividing
the age of the universe by the gravitational constant,
you'll never get the 2*pis out.

When there's two more or less different ways to look at
a problem, and they give different answers, I usually
think that defining one of them to be correct isn't
useful. It'll confuse too many people some of the time.

Dick Hendrickson

glen herrmannsfeldt

2004-05-14, 7:31 pm

Dr Chaos wrote:

> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:

[color=darkred]
[color=darkred]
[color=darkred]
> sin(1e100) is a symptom of a mistake.

I suppose it should have been sin(1D100) in this newsgroup...

> I believe that if one's program depends on evaluating that to any
> degree of accuracy, it is computing something in a mathematically
> mistaken way.

> Generally you need some paper-and-pencil analysis and thinking to
> avoid this problem.

[color=darkred]
> Sometimes computer scientists have an odd idea what it means to be
> "correct", e.g. imagining that the particular rules of IEEE floating
> point define "correctness" as opposed to solving the problem in a
> satisfactory way.

It seems that some believe that it is possible to consider
even a very large floating point number as a value in radians,
and compute the sin() for that value. It may take hundreds
or thousands of digits of pi to do the computation.
(There are floating point systems with 15 or 16 bit exponents.)

None of the arguments made so far in comp.lang.c seem
to have been convincing that there is no mathematical
use for such, except possibly when one is already using
multiple precision arithmetic.

thanks for replying,

-- glen

Gerry Thomas

2004-05-14, 8:31 pm

"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:lKbpc.50243$z06.7101055@attbi_s01...
> Dr Chaos wrote:
>
>
>
>
>
>
> I suppose it should have been sin(1D100) in this newsgroup...
>
>
>
>
>
> It seems that some believe that it is possible to consider
> even a very large floating point number as a value in radians,
> and compute the sin() for that value. It may take hundreds
> or thousands of digits of pi to do the computation.
> (There are floating point systems with 15 or 16 bit exponents.)
>
> None of the arguments made so far in comp.lang.c seem
> to have been convincing that there is no mathematical
> use for such, except possibly when one is already using
> multiple precision arithmetic.
>
> thanks for replying,
>
> -- glen
>

In MATLAB: sin(1e100) = -0.38063773100503
In IVF 8: sin(1.d100) = -0.380637731005029
In CVF 6.6C:sin(1.d100) = -0.511101926654941

The odd man out uses the Microsoft math library.

--
E&OE

Ciao,
Gerry T.

Richard Maine

2004-05-14, 8:31 pm

glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

On computing things like sin(1d100).

> None of the arguments made so far in comp.lang.c seem
> to have been convincing that there is no mathematical
> use for such...

I wouldn't try to argue that there is no use for such a thing.
never say "never" and all that.

I will say that, even with the double precision, the finite
representation of 1d100 will not exactly equal the mathematical
value. There will be rounding error even before talking about
the trig function. Once you ackowledge that the input to sin()
might be in error by up to half an ulp, then the result of the
sin() could be anything between +/-1. The interval arithmetic
folk can tell you that.

In order for particular result of the sin() to have any useful
meaning, the input to the sin() first needs to be known to be exact
(or anyway to have a lot less than haf an ulp error). There
might be cases where that is so, I suppose. Seems unlikely
for whatever 1d100 evaluates to to be one of those cases. Now
if you had something like 2**300, where the input was exactly
representable, I might find it more plausible that some trickery
to get a "good" value of the sin() might have some use.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Dr Chaos

2004-05-14, 8:31 pm

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote:
> Dr Chaos wrote:
>
>
>
>
>
>
> I suppose it should have been sin(1D100) in this newsgroup...
>
>
>
>
>
> It seems that some believe that it is possible to consider
> even a very large floating point number as a value in radians,
> and compute the sin() for that value. It may take hundreds
> or thousands of digits of pi to do the computation.
> (There are floating point systems with 15 or 16 bit exponents.)
>
> None of the arguments made so far in comp.lang.c seem
> to have been convincing that there is no mathematical
> use for such, except possibly when one is already using
> multiple precision arithmetic.

I assume you mean the reverse, no arguments are convincing
that there is a mathematical use?

Mr Plauger is adamant that the library should implement it like
Standard C says and there is no 'cop out'. What he does is admittedly
quite heroic, but ultimately pointless.

I think a 'cop-out', e.g. raise(USER_CLUE_ABSENT)

is in fact desirable.

In practice, to help the user the most, the man page for sin()
and cos() ought to point to instructions about what people
should be doing instead for various common situations.

Rich Townsend

2004-05-14, 10:31 pm

Gerry Thomas wrote:

<snip>

> In MATLAB: sin(1e100) = -0.38063773100503
> In IVF 8: sin(1.d100) = -0.380637731005029
> In CVF 6.6C:sin(1.d100) = -0.511101926654941
>
> The odd man out uses the Microsoft math library.

Mathematica 4.2: sin(1e100) = -0.567657
Lahey F95 6.2b: sin(1.d100) = 1.000000000000000E+100 (!!!)
g77 3.3.2 : sin(1.d100) = -0.567656625
Compaq F95 5.5 : sin(1.d100) = -0.380637731005029
NAG F95 4.0a: sin(1.d100) = -0.5676566254385991
Sun F95 7.1: sin(1.d100) = -0.3806377310050287

Apart from the crazy Lahey result (which could be interpreted as a flag
value saying "Are you sure you want SIN(1e100)?"), there are two values
here: -0.567 and -0.380 (the latter being the Matlab/IVF values given
above by Gerry).

I'm not really sure which one to "trust" -- supposedly, Mathematica uses
arbitrary-precision arithmetic, but I'm can't be sure whether it is
"correct".

Rich

Rich Townsend

2004-05-14, 10:31 pm

Rich Townsend wrote:
> Gerry Thomas wrote:
>
> <snip>
>
>
>
> Mathematica 4.2: sin(1e100) = -0.567657
> Lahey F95 6.2b: sin(1.d100) = 1.000000000000000E+100 (!!!)
> g77 3.3.2 : sin(1.d100) = -0.567656625
> Compaq F95 5.5 : sin(1.d100) = -0.380637731005029
> NAG F95 4.0a: sin(1.d100) = -0.5676566254385991
> Sun F95 7.1: sin(1.d100) = -0.3806377310050287
>
> Apart from the crazy Lahey result (which could be interpreted as a flag
> value saying "Are you sure you want SIN(1e100)?"), there are two values
> here: -0.567 and -0.380 (the latter being the Matlab/IVF values given
> above by Gerry).
>
> I'm not really sure which one to "trust" -- supposedly, Mathematica uses
> arbitrary-precision arithmetic, but I'm can't be sure whether it is
> "correct".
>
> Rich

Oh, and Google gives sin(1e100) = 1e100. Perhaps this is evidence that
Google use a Lahey backend? :)

Gerry Thomas

2004-05-14, 10:31 pm

"Rich Townsend" <rhdt@barVOIDtol.udel.edu> wrote in message
news:c83ofq$od9$1@scrotar.nss.udel.edu...
> Gerry Thomas wrote:
>
> <snip>
>
>
> Mathematica 4.2: sin(1e100) = -0.567657
> Lahey F95 6.2b: sin(1.d100) = 1.000000000000000E+100 (!!!)
> g77 3.3.2 : sin(1.d100) = -0.567656625
> Compaq F95 5.5 : sin(1.d100) = -0.380637731005029
> NAG F95 4.0a: sin(1.d100) = -0.5676566254385991
> Sun F95 7.1: sin(1.d100) = -0.3806377310050287
>
> Apart from the crazy Lahey result (which could be interpreted as a flag
> value saying "Are you sure you want SIN(1e100)?"), there are two values
> here: -0.567 and -0.380 (the latter being the Matlab/IVF values given
> above by Gerry).
>
> I'm not really sure which one to "trust" -- supposedly, Mathematica uses
> arbitrary-precision arithmetic, but I'm can't be sure whether it is
> "correct".
>

This is where Dick Henrickson's maxim comes to the fore:

"When there's two more or less different ways to look at
a problem, and they give different answers, I usually
think that defining one of them to be correct isn't
useful. It'll confuse too many people some of the time."

--
E&OE

Ciao,
Gerry T.

James Giles

2004-05-14, 10:31 pm

Robert Corbett

2004-05-15, 1:31 am

glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote in message news:<lKbpc.50243$z06.7101055@attbi_s01>...

> It seems that some believe that it is possible to consider
> even a very large floating point number as a value in radians,
> and compute the sin() for that value. It may take hundreds
> or thousands of digits of pi to do the computation.
> (There are floating point systems with 15 or 16 bit exponents.)

It takes a lot of space, but not much time. In most cases, the
argument reduction code needs to look at about twice as many bits
in the binary fraction representation of pi as bits in the
significand of the input argument.

Sun's implementation of the trig functions does true pi argument
reduction. It is easier to explain giving an approximation to
the true mathematical answer with an error of just over 1/2 ulp
than to explain when the answer starts to deviate from the true
mathematical result and why that is usually OK.

Sincerely,
Bob Corbett

Jesper Harder

2004-05-15, 1:31 am

Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:

> Gerry Thomas wrote:
>
>
> Mathematica 4.2: sin(1e100) = -0.567657
> Lahey F95 6.2b: sin(1.d100) = 1.000000000000000E+100 (!!!)
> g77 3.3.2 : sin(1.d100) = -0.567656625
> Compaq F95 5.5 : sin(1.d100) = -0.380637731005029
> NAG F95 4.0a: sin(1.d100) = -0.5676566254385991
> Sun F95 7.1: sin(1.d100) = -0.3806377310050287
>
> Apart from the crazy Lahey result (which could be interpreted as a
> flag value saying "Are you sure you want SIN(1e100)?"), there are two
> values here: -0.567 and -0.380 (the latter being the Matlab/IVF values
> given above by Gerry).
>
> I'm not really sure which one to "trust"

I depends on what you really mean by sin(1.d100) -- do you want:

a) the sine of the rational that approximates the float to the
accuracy of the floating-point representation, i.e.

1000000000000000015902891109759918046836
0808563945281389781327557747838772170381
060813469985856815104

b) or the sine of the exact mathematical number 10^100:

1000000000000000000000000000000000000000
0000000000000000000000000000000000000000
000000000000000000000

I think -0.38 is the correct answer to a). I confirmed it by
calculating the range reduction with exact rational arithmetic and the
convergents of a continued fractions expansion of pi.

--
Jesper Harder <http://purl.org/harder/>

Rich Townsend

2004-05-15, 1:31 am

NNTP-Posting-Host: shayol.bartol.udel.edu
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: scrotar.nss.udel.edu 1084595495 29134 128.175.14.63 (15 May 2004 04:31:35 GMT)
X-Complaints-To: abuse@udel.edu
NNTP-Posting-Date: Sat, 15 May 2004 04:31:35 +0000 (UTC)
User-Agent: Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.6) Gecko/20040318
X-Accept-Language: en-us, en
In-Reply-To: <m3llju5nkp.fsf@defun.localdomain>
Xref: number1.nntp.dca.giganews.com comp.lang.fortran:118673

Jesper Harder wrote:
> Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:
>
>
>
>
> I depends on what you really mean by sin(1.d100) -- do you want:
>
> a) the sine of the rational that approximates the float to the
> accuracy of the floating-point representation, i.e.
>
> 1000000000000000015902891109759918046836
0808563945281389781327557747838772170381
060813469985856815104
>
> b) or the sine of the exact mathematical number 10^100:
>
> 1000000000000000000000000000000000000000
0000000000000000000000000000000000000000
000000000000000000000
>

Neither, really. If I ever end up needing sin(Anything Very Larg), then
I know I'm screwing up a piece of physics somewhere.

cheers,

Rich

Gerry Thomas

2004-05-15, 1:31 am

"Robert Corbett" <robert.corbett@sun.com> wrote in message
news:cb977dbc.0405142007.64c54be2@posting.google.com...
> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote in message
news:<lKbpc.50243$z06.7101055@attbi_s01>...
>
>
> It takes a lot of space, but not much time. In most cases, the
> argument reduction code needs to look at about twice as many bits
> in the binary fraction representation of pi as bits in the
> significand of the input argument.
>
> Sun's implementation of the trig functions does true pi argument
> reduction. It is easier to explain giving an approximation to
> the true mathematical answer with an error of just over 1/2 ulp
> than to explain when the answer starts to deviate from the true
> mathematical result and why that is usually OK.

Indeed, K.C. Ng's Sun report is revealing: pi is not evaluated as a double
but from from an on-chip look-up table to ca. 55D, IIRC. Others have
implemented such an argument reduction scheme, but like C99, Fortran 2003+
couldn't give a fcuk.

--
E&OE

Ciao,
Gerry T.

>
> Sincerely,
> Bob Corbett

James Giles

2004-05-15, 2:31 am

Robert Corbett wrote:
> glen herrmannsfeldt <gah@ugcs.caltech.edu> wrote in message
> news:<lKbpc.50243$z06.7101055@attbi_s01>...
>
>
> It takes a lot of space, but not much time. In most cases, the
> argument reduction code needs to look at about twice as many bits
> in the binary fraction representation of pi as bits in the
> significand of the input argument.

Hypothetically it could need pi to over a thousand bits (the
exponent range of double precision - I'm not even considering
double extended).

> Sun's implementation of the trig functions does true pi argument
> reduction. It is easier to explain giving an approximation to
> the true mathematical answer with an error of just over 1/2 ulp
> than to explain when the answer starts to deviate from the true
> mathematical result and why that is usually OK.

Do you really use a 1024+ bit representation of pi? Double
extended (over 16000 bits needed) could start to cause a
slowdown.

--
J. Giles

glen herrmannsfeldt

2004-05-15, 3:31 am

Rich Townsend wrote:
(snip)

> Neither, really. If I ever end up needing sin(Anything Very Larg), then
> I know I'm screwing up a piece of physics somewhere.

The first time it happened to me was when an undefined
variable was used. The IBM (G, H, and VS) Fortran libraries
have a fatal error when the argument is greater than
pi*2**18 (short), pi*2**50 (long), or pi*2**100 (extended)
(the IBM names for floating point precision).

So when it happened to me, someone nicely explained why
it would do that.

Not so long after that, someone I knew had a brand new HP-55
calculator. The sine (in degrees) of 9.999999999e99 was
not zero, somewhere close to -0.53. Since 9.999999999e99
is obviously a multiple of 360 (if extended with zeros),
I knew it wasn't right, yet the calculator owner was sure
that HP wouldn't lie.

-- glen

Tim Prince

2004-05-15, 3:31 am

"Jesper Harder" <harder@myrealbox.com> wrote in message
news:m3llju5nkp.fsf@defun.localdomain...
> Rich Townsend <rhdt@barVOIDtol.udel.edu> writes:
>
That's the documented behavior of the 386 instruction set. It sets an error
flag when it gives up and returns the argument, on account of excessive
magnitude. IIRC, the value of Pi used (for arguments of reasonable
magnitude) is rounded to 66 bits binary precision.[color=darkred]
g77 implements sin() in terms of the underlying C math library. This surely
depends on the platform, probably differs among the major Windows
implementations, for example, and may differ with the compile options
chosen.[color=darkred]
>
> I depends on what you really mean by sin(1.d100) -- do you want:
>
> a) the sine of the rational that approximates the float to the
> accuracy of the floating-point representation, i.e.
>
>
1000000000000000015902891109759918046836
080856394528138978132755774783877217
0381060813469985856815104

>
> b) or the sine of the exact mathematical number 10^100:
>
>
1000000000000000000000000000000000000000
000000000000000000000000000000000000
0000000000000000000000000

>
> I think -0.38 is the correct answer to a). I confirmed it by
> calculating the range reduction with exact rational arithmetic and the
> convergents of a continued fractions expansion of pi.
>
As several of the calculations were done on hardware which has the 64-bit
precision mode available, the corresponding approximation to 1e100 could
also be a candidate for use as the argument. Strict adherence to Fortran
would require a) the value which agrees with the implementation of the
declared data type.

glen herrmannsfeldt

2004-05-15, 3:31 am

James Giles wrote:

> Let's try a different but related question and see if it sheds
> any light. Suppose PI is name of the double precision value
> nearest the mathematical value of pi. The sine of an integer
> multiple of pi is zero. Should SIN(N*PI) be zero?

> Mathematically, the answer is no. For any non-zero N, N*PI
> is *not* really an integer multiple of the mathematical value
> of pi. For N greater than 2^54, the difference between N*PI
> and the nearest integer multiple of pi is likely to be more than
> a radian. So, what value do you expect? Should SIN(N*PI)
> return zero, of should it do careful argument reduction and
> return some non-zero value (could be anywhere between
> -1.0d0 and 1.0d0).

> This is just a small subcase of the question you asked. I've
> had people argue vehemently that N*PI should be treated as
> an integer multiple of pi. Yet numerically that's the wrong
> answer.

Well, there have been no real problems posted yet that use
the exact, maybe thousands of digits, version of pi, yet
I posted

for(x=0;x+60!=x;x+=60) printf("%20.15f\n",sin(a*x));

Maybe:

X=0
20 IF(X+90.EQ.X) GOTO 10
WRITE(*,'F20.25') SIN(A*X)
GOTO 20
10 CONTINUE

where it would be nice to have a value of a such that
x was in degrees. This does require the appropriate
approximation to pi, otherwise it starts diverging
from the desired values very early.

Though a sind() function with the argument in degrees
makes a better solution to that problem. (Do any Fortran
versions have one?) Then there is no question about the
approximations used in argument reduction, though large
numbers can still be a problem.

-- glen

Gerry Thomas

2004-05-15, 3:31 am

"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:hoipc.52319$iF6.4774854@attbi_s02...
> Rich Townsend wrote:
> (snip)
>
>
> The first time it happened to me was when an undefined
> variable was used. The IBM (G, H, and VS) Fortran libraries
> have a fatal error when the argument is greater than
> pi*2**18 (short), pi*2**50 (long), or pi*2**100 (extended)
> (the IBM names for floating point precision).
>
> So when it happened to me, someone nicely explained why
> it would do that.
>
> Not so long after that, someone I knew had a brand new HP-55
> calculator. The sine (in degrees) of 9.999999999e99 was
> not zero, somewhere close to -0.53. Since 9.999999999e99
> is obviously a multiple of 360 (if extended with zeros),
> I knew it wasn't right, yet the calculator owner was sure
> that HP wouldn't lie.

Such innocence!, and I thought you were a honest-to-goodness hippie, like
.... etc.,:-).

The only thing screwed about Townsend is Townsend, not UDel, physics,
numerics, yours truly, c.l.f.., etc., it's apparently typical of his ilk to
project (just as aussies 'splork', but I'm saving that for after I drop
into TRANZAC and witness now they react to an American greased pig on the
loose, wombats being in short supply at the moment), so just take a mop and
clean up after him, 'there but for the grace ...' and all that.

--
E&OE

HTH,
Gerry T.

Gerry Thomas

2004-05-15, 3:31 am

"Tim Prince" <tprince@computer.org> wrote in message
news:Rtipc.8478$gF3.3728@newssvr27.news.prodigy.com...
>
> "Jesper Harder" <harder@myrealbox.com> wrote in message
> news:m3llju5nkp.fsf@defun.localdomain...
> That's the documented behavior of the 386 instruction set. It sets an
error
> flag when it gives up and returns the argument, on account of excessive
> magnitude. IIRC, the value of Pi used (for arguments of reasonable
> magnitude) is rounded to 66 bits binary precision.
> g77 implements sin() in terms of the underlying C math library. This
surely
> depends on the platform, probably differs among the major Windows
> implementations, for example, and may differ with the compile options
> chosen.
values[color=darkred]
>
1000000000000000015902891109759918046836
08085639452813897813275577478387721
7
> 0381060813469985856815104
>
1000000000000000000000000000000000000000
00000000000000000000000000000000000
0
> 0000000000000000000000000
> As several of the calculations were done on hardware which has the 64-bit
> precision mode available, the corresponding approximation to 1e100 could
> also be a candidate for use as the argument. Strict adherence to Fortran
> would require a) the value which agrees with the implementation of the
> declared data type.
>

I've also done my calculations in range and interval arithmetic. Let's
compare, you go first. Chicken!, true to form.

--
Ciao,
Gerry T.
______
"Nobody can get the truth out of me because even I don't know what it is. I
keep myself in a constant state of utter confusion." -- Col Sam Flagg,
ICORPS dropin to the 4077th M*A*S*H

Gerry Thomas

2004-05-15, 4:31 am

"James Giles" <jamesgiles@worldnet.att.net> wrote in message
news:P8ipc.61199$Ut1.1578455@bgtnsc05-news.ops.worldnet.att.net...
> Robert Corbett wrote:
>
> Hypothetically it could need pi to over a thousand bits (the
> exponent range of double precision - I'm not even considering
> double extended).
>
>
>
> Do you really use a 1024+ bit representation of pi? Double
> extended (over 16000 bits needed) could start to cause a
> slowdown.
>

Giles, your trying to decieve yourself, but what's new? Nobody uses a
1024+ bit representation of pi and even if they do, you'll be the last to
know, so who care?, patentlly not you..

--
E&OE

Ciao,
Gerry T.
______
"Ah, Klinger, my constant reminder that Darwin was right!." -- Maj Charles
Winchester III, the 4077th M*A*S*H

Robert Corbett

2004-05-15, 6:31 am

In article <P8ipc.61199$Ut1.1578455@bgtnsc05-news.ops.worldnet.att.net>,
James Giles <jamesgiles@worldnet.att.net> wrote:
>
>
>
>Do you really use a 1024+ bit representation of pi? Double
>extended (over 16000 bits needed) could start to cause a
>slowdown.

The method Sun uses does not require as many bits as you think.
The tables for true pi argument reduction for quadruple-precision
floating-point use about 7,500 bits. That is less than one
kilobyte. For computers available before 1980, 1 KB tables would
be intolerable. For modern computers, it is likely to be less
than 1% of 1% of 1% of the physical memory available.

Argument reduction for very large arguments is slower than for
small arguments. It is not as much slower as one might at first
think. Since quadruple-precision on SPARC is implemented in
software, the cost of argument reduction for large
quadruple-precision numbers is about 25% of the cost of
evaluating the sine of the reduced argument.

Sincerely,
Bob Corbett

James Giles

2004-05-15, 5:31 pm

Gerry Thomas wrote:
> "James Giles" <jamesgiles@worldnet.att.net> wrote in message
> news:P8ipc.61199$Ut1.1578455@bgtnsc05-news.ops.worldnet.att.net...
....[color=darkred]
>
> Giles, your trying to decieve yourself, but what's new? [...]

In what way is a question an attempt to decieve? Myself or
anyone else? You are usually the one attempting to decieve.

--
J. Giles

James Giles

2004-05-15, 5:31 pm

Robert Corbett wrote:
> In article <P8ipc.61199$Ut1.1578455@bgtnsc05-news.ops.worldnet.att.net>,
> James Giles <jamesgiles@worldnet.att.net> wrote:
>
> The method Sun uses does not require as many bits as you think.
> The tables for true pi argument reduction for quadruple-precision
> floating-point use about 7,500 bits. That is less than one
> kilobyte. For computers available before 1980, 1 KB tables would
> be intolerable. For modern computers, it is likely to be less
> than 1% of 1% of 1% of the physical memory available.

Well, the storage cost of this was never my worry. But, I'd like
to know more about the argument reduction algorithm. How do
you get by without pi to sufficient precision to cover the whole
exponent range of the argument type?

--
J. Giles

Eric J. Kostelich

2004-05-16, 2:31 am

In article <KAbpc.59094$Ut1.1541457@bgtnsc05-news.ops.worldnet.att.net>,
Dick Hendrickson <dick.hendrickson@att.net> wrote:
>
>glen herrmannsfeldt wrote:
>
>We know mathematically Sin(x+2*N*pi) = Sin(x).
>But once N is large we can't computationally recover
>the value of x, it gets lost in the truncation to finite representation.

This is the only correct way to look at the problem. The most
significant n decimal digits of all real numbers between 10**n and
(10**n + pi) are identical. Consequently, there's no meaningful
floating-point result for sin(x) in any arithmetics in common use for
|x| larger than 10**32 or so; any number from -1 to 1 is an equally
plausible return value for sin(x) for such large arguments.

>But, if we look at the problem in a different way,
>and if I could remember my high school trig, I'm sure I
>could expand Sin(2.0) into some nice function of Sin(1.0)

Well, sin(2x) = 2 sin(x) cos(x), and you can continue this
process to expand sin(nx) in terms of various powers of sin(x)
and cos(x), but the resulting series alternates in sign and
would suffer catastrophic cancelation. So you don't get anything
useful from this route.

>Maybe the expectaion depends on where the 1.e100 comes from.

You can't possibly know what the programmer might have in mind
in the general case. In my opinion, it would be reasonable for
a library sine or cosine function to generate some sort of exception
when the argument is too large to return a meaningful result.

--Eric

Dr Chaos

2004-05-16, 3:32 pm

On Fri, 14 May 2004 20:46:52 -0400, Gerry Thomas <gfthomas@sympatico.ca> wrote:
> This is where Dick Henrickson's maxim comes to the fore:
>
> "When there's two more or less different ways to look at
> a problem, and they give different answers, I usually
> think that defining one of them to be correct isn't
> useful. It'll confuse too many people some of the time."

This is what computer scientists who think they're helping
numerical users don't always see.

Robert Corbett

2004-05-17, 3:31 am

In article <Btvpc.65043$Ut1.1641959@bgtnsc05-news.ops.worldnet.att.net>,
James Giles <jamesgiles@worldnet.att.net> wrote:
>Robert Corbett wrote:
>
>Well, the storage cost of this was never my worry. But, I'd like
>to know more about the argument reduction algorithm. How do
>you get by without pi to sufficient precision to cover the whole
>exponent range of the argument type?

A colleague saved me from repeating my mistake. I did some
bad arithmetic in computing the size of the tables. They
are actually 21,984 bits (2.75 kilobytes). That is far more
bits than will ever to be needed.

Sincerely,
Bob Corbett

Paul Van Delst

2004-05-17, 11:37 am

Rich Townsend wrote:[color=darkred]
> Rich Townsend wrote:
>

JFTHOI, one more data point for IDL Version 6.0.3 (linux x86 m32)

IDL> print, sin(1.0d100), format='(f20.15)'
-0.567656625438599

And, FWIW, if I received code where the argument of a sin/cos function could exceed
several multiples of 2pi, I'd probably send it back and tell them to fix it (or just not
use it). I'd also wonder what other variables in the code could also contain nonsense
values. My experience is admittedly limited, but (apart from, say, unwrapping phase
ambiguities in FTS spectra) I don't see too much need for a sin/cos function arg to go
beyond 2pi.

cheers,

paulv

keith bierman

2004-05-17, 6:31 pm

glen herrmannsfeldt wrote:
>
> What result do people expect from the sin() function for
> very large arguments, for example sin(1e100)?
>
>
> -- glen
>

From the Sun Numerical Computation Guide

> Trigonometric functions for radian arguments outside the range [-pi/4,pi/4] are
> usually computed by reducing the argument to the indicated range
by subtracting
> integral multiples of ./2.
>
> Because pi is not a machine-representable number, it must be approximated.
> The error in the final computed trigonometric function depends on
the
> rounding errors in argument reduction (with an approximate pi as well as
> the rounding), and approximation errors in computing the
trigonometric
> function of the reduced argument. Even for fairly small arguments,
> the relative error in the final result may be dominated by the
argument
> reduction error, while even for fairly large arguments, the error
due to
> argument reduction may be no worse than the other errors.
>
> There is widespread misapprehension that trigonometric functions of all
> large arguments are inherently inaccurate, and all small
arguments relatively
> accurate. This is based on the simple observation that large enough
> machine-representable numbers are separated by a distance greater
than pi.
>
> There is no inherent boundary at which computed trigonometric function
> values suddenly become bad, nor are the inaccurate function
values useless.
> Provided that the argument reduction is done consistently, the
fact that
> the argument reduction is performed with an approximation to pi is
> practically undetectable, because all essential identities and
> relationships are as well preserved for large arguments as for small.
>
> libm and libsunmath trigonometric functions use an "infinitely" precise pi
> for argument reduction. The value 2/. is computed to 916
hexadecimal digits
> and stored in a lookup table to use during argument reduction.
>
> The group of functions sinpi, cospi, tanpi (see Table 3-2 on page 34), scales the
> input argument by pi to avoid inaccuracies introduced by range
reduction.

I have an old draft of a paper which was published some years ago, but I
can't find a version online (it was in a sun developer newsletter). I
will happily send the pdf to any interested party.

keith bierman

2004-05-17, 6:31 pm

glen herrmannsfeldt wrote:
>

The classic paper by Goldberg on what Every Computer Scientist should
know about Floating Point arithmetic (originally an ACM paper) can be
found at docs.sun.com. It's a good read.[color=darkred]
>
> None of the arguments made so far in comp.lang.c seem
> to have been convincing that there is no mathematical
> use for such, except possibly when one is already using
> multiple precision arithmetic.
>

For those who have never seen the result before, it may be surprising.
But "infinite" pi only needs to be big enough for the wordsize. So it's
not too bad, and you only need to use as many digits as you need for the
particular argument. So it's a little slower than doing it sloppy, but
for "normal" sized arguments it's not excessively so.

So, I personally expect a decent implementation to give me at least as
good a result as fdlibm (it's in netlib; contributed by Sun's floating
point group). Hopefully faster :>

Since a freeware implementation is out there, there is no need to settle
for less.

keith bierman

2004-05-17, 6:31 pm

Gerry Thomas wrote:
>
> Indeed, K.C. Ng's Sun report is revealing: pi is not evaluated as a double
> but from from an on-chip
It's all software. Just where did you find something which you could
interpret differently?

Dr Chaos

2004-05-17, 6:31 pm

On Mon, 17 May 2004 14:38:36 -0700, keith bierman <keith.bierman@sun.com> wrote:
> glen herrmannsfeldt wrote:
>
> The classic paper by Goldberg on what Every Computer Scientist should
> know about Floating Point arithmetic (originally an ACM paper) can be
> found at docs.sun.com. It's a good read.
>
> For those who have never seen the result before, it may be surprising.
> But "infinite" pi only needs to be big enough for the wordsize. So it's
> not too bad, and you only need to use as many digits as you need for the
> particular argument. So it's a little slower than doing it sloppy, but
> for "normal" sized arguments it's not excessively so.
>
> So, I personally expect a decent implementation to give me at least as
> good a result as fdlibm (it's in netlib; contributed by Sun's floating
> point group). Hopefully faster :>
>
> Since a freeware implementation is out there, there is no need to settle
> for less.

again, what is the mathematical value of it?

Gerry Thomas

2004-05-17, 7:32 pm

Gerry Thomas

2004-05-17, 7:32 pm

"keith bierman" <keith.bierman@sun.com> wrote in message
news:40A9330C.9090005@sun.com...
> Gerry Thomas wrote:
double[color=darkred]
> It's all software. Just where did you find something which you could
> interpret differently?
>

If you had ever read KC Ng's SunPro report, "Argument Reduction for Huge
Arguments: Good to the Last Bit", 1992, he correctly states:
"... there is no standard choice of the precision of machine pi, which
means the value of a trigonometric function is machine dependent as well as
precision dependent. Current commercial vendors of math co-processors (like
the Intel x87 and motorola 68882) use a 66 bit pi for their trigonometric
functions, which is adequate for IEEE double (53 bits) or double extended
(64 bits) floating-point arithmetic,...."
Where does x87, motorola 68882, etc., squirrel away their values of pi?,
never mind.

--
E&OE

HTH,
Gerry T.

Gerry Thomas

2004-05-17, 8:31 pm

"keith bierman" <keith.bierman@sun.com> wrote in message
news:40A9330C.9090005@sun.com...
> Gerry Thomas wrote:
double[color=darkred]
> It's all software. Just where did you find something which you could
> interpret differently?
>

Do not send me personal email, my SysOp has already bounced your latest
missive. Kindly, confine you comments, which are presumably in your
capacity as a Sun Microsystems, Inc. employee, to the newsgroup.

Thank you,
Gerry T.

Robert Corbett

2004-05-17, 9:31 pm

Paul Van Delst <paul.vandelst@noaa.gov> wrote in message news:<c8akli$6ad$1@news.nems.noaa.gov>...

> And, FWIW, if I received code where the argument of a sin/cos function could exceed
> several multiples of 2pi, I'd probably send it back and tell them to fix it (or just not
> use it). I'd also wonder what other variables in the code could also contain nonsense
> values. My experience is admittedly limited, but (apart from, say, unwrapping phase
> ambiguities in FTS spectra) I don't see too much need for a sin/cos function arg to go
> beyond 2pi.

Whether sin(X) for large values of X causes a problem depends on the context.
For example, in the expression

X**2 + SIN(X)

the value of SIN(X) has a large effect when X is small, and a small effect
when X is large, which makes it a safe computation. I would be very annoyed
by the implementations described in earlier postings that delivered 1D100
in such a context.

Robert Corbett

Gerry Thomas

2004-05-17, 10:31 pm

"Robert Corbett" <robert.corbett@sun.com> wrote in message
news:cb977dbc.0405171647.68ebfb8e@posting.google.com...
> Paul Van Delst <paul.vandelst@noaa.gov> wrote in message
news:<c8akli$6ad$1@news.nems.noaa.gov>...
>
could exceed[color=darkred]
fix it (or just not[color=darkred]
contain nonsense[color=darkred]
unwrapping phase[color=darkred]
function arg to go[color=darkred]
>
> Whether sin(X) for large values of X causes a problem depends on the
context.
> For example, in the expression
>
> X**2 + SIN(X)
>
> the value of SIN(X) has a large effect when X is small, and a small
effect
> when X is large, which makes it a safe computation. I would be very
annoyed
> by the implementations described in earlier postings that delivered 1D100
> in such a context.

This is irresponsible baloney:

From SunPro KC Ng's report, "Argument Reduction for Huge
Arguments: Good to the Last Bit", 1992:
"It is often argued that being concerned about large arguments is
unnecessary, because

sophisticated users simply know better than to compute with large angles.
It

is our contention that this position is suboptimal, because:

1. It places an unnecessary burden on the user.

2. The consequences of producing incorrect (inaccurate) answers may be

catastrophic; many people assume that computers can do arithmetic very

well. While numerical analysts know better, not all programmers are

numerical analysts, nor should they be.

3. It is a vendors responsibility to provide answers that are as correct as

possible"

Don't you get the message or is Microsoft's 2b$ not enough? Get you're
resume updated, boy.

--
"Indians are not niggers. They're wogs." -- Major in Fawlty Towers.

One approach to argument reduction

Gib Bogle

2004-05-18, 12:31 am

James Giles wrote:
> Gerry Thomas wrote:
>
>
> ...
>
>
>
> In what way is a question an attempt to decieve? Myself or
> anyone else? You are usually the one attempting to decieve.

Gentlemen, please: "I before E, except after C". The word is "deceive".

Gib

Gerry Thomas

2004-05-18, 1:31 am

"Gib Bogle" <bogle@too.much.spam.ihug.co.nz> wrote in message
news:c8bufv$pmh$1@lust.ihug.co.nz...
> James Giles wrote:
>
> Gentlemen, please: "I before E, except after C". The word is "deceive".
>
> Gib
>
Gib, just received your message, hope I havent't deceived you. Jes's, Gib,
I'm still typing like a reprter (yeh, can't spel netiher, wat kin i do, go
Microsooft, or wat?)

--
E&OE

Ciao,
Gerry T.

Gerry Thomas

2004-05-18, 1:31 am

"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message
news:FRapc.12381$6f5.977509@attbi_s54...
> There is a discussion in comp.lang.c titled
> Sine code for ANSI C
>
> regarding the evaluation of the sin() function for large
> arguments. It might be more applicable in this
> group, so I thought I would ask here.
>
> What result do people expect from the sin() function for
> very large arguments, for example sin(1e100)?
>
> Hopefully the discussion won't get too nasty, but it
> does seem to be a contentious issue over there.
>
> -- glen
>

Glen is to be commended for bringing this issue here to c.l.f., the
newsgroup that lays claim to all things right in numerics (chortle,
chortle!).

Several years ago I suggested to the hawker Appleyard of
polyputthekettleonweallwantteaandbecchma
rks that he include a test on
trignometrics and the exponential of large argument as a QA test of the
current Fortran offerings. Hell no, this bucks Herr Maine (whose
defensiveness of polyputthekettleonweallwantteaandbecchma
rks Appleyard is
now clear, nodnod, winkwink, synomore, saynomore), uncovers vendor deciet
(sic deceit, WGF), Fortran standard technical committee bankruptcy, click
your heels!, stand erect!, and make the true Father-of-Fortran Yorkey
proud!, also read two chapters of Imperium, whip yourself but don't get
carried away, nodnod, winkwink, saynomore, saynomore, and get that wombat
out of Maine's eyes, he's got his eye on it!.

IMNHO, most practioneer of Fortran are the most naive suckers when it comes
to numerics: they know such noises as eps, tols, single, double, etc., but
they lack any deep understanding of fp numerics, particularly of its
pitfalls. The N(umerical) A(nalysis) fraternity have bent over backwards to
accommodate the uninitiated but the Fortran set would rather stick with
Numerical Recipes, Rodney Dangerfield's bible of NA, rather than face the
reality that Fortran 2++ (whatever) just doesn't give a fig about the
output of Fortran programs.

Today on c.l.f., we've had someone from .gov say more or less "well I
don't know, I'm all right Jack!" Lucky for him I'm not Jack! This is why
the Plugger's of this world can get buy with impunity on the 'the what the
idiots won't know won't hurt them' m.o. Evidently, this cavalier attitude
extends beyond Fortran. So, Mr Maine and your clients, you're exposed for
the charlatans that you are.

BTW, this from Cleve Moler: In Appendix 3 of 'Recent Progress of the Java
Grande Numerics Working Group, June 7, 1999':

1. "The exponential, sine and cosine functions in Microsoft's math library
are so inaccurate for large arguments that the library is unsuitable for
general-purpose use."

2. "(At the MathWorks, we use FDLIBM as the underlying library for MATLAB
on
all platforms. We did this primarily because we got tired of working around
the bugs, poor algorithms and inconsistent handling of exceptional cases
that we found in the libm's provided by other hardware and operating
systems
vendors.)"

So, for 1/2 decade the bimbo's of .gov have been going around telling
everyone to shutter up and pray with Fr Richard M., S.J. NOT, and his
coterie, for all is well and things turn out for the better and all that
polyanna crap. Be Jeasus.

--
E&OE

Ciao,
Gerry T.
______
"The old order changeth, yielding place to new,
And God fulfils Himself in many ways,
Lest one good custom should corrupt the world." -- Alfred Lord Tennyson, in
Morte d'Arthur.

Dr Chaos

2004-05-18, 1:32 pm

On Mon, 17 May 2004 21:04:44 -0400, Gerry Thomas <gfthomas@sympatico.ca> wrote:
>
> From SunPro KC Ng's report, "Argument Reduction for Huge
> Arguments: Good to the Last Bit", 1992:
> "It is often argued that being concerned about large arguments is
> unnecessary, because
>
> sophisticated users simply know better than to compute with large angles.
> It
>
> is our contention that this position is suboptimal, because:
>
> 1. It places an unnecessary burden on the user.

which is?

> 2. The consequences of producing incorrect (inaccurate) answers may be
>
> catastrophic; many people assume that computers can do arithmetic very
>
> well. While numerical analysts know better, not all programmers are
>
> numerical analysts, nor should they be.

That's precisely why doing infinitely precise argument reduction is
silly.

> 3. It is a vendors responsibility to provide answers that are as correct as
>
> possible"

It is a vendors responsibility to provide answers which lead to their
client's programs working the best.

I submit, that if a change in the least significant bit of the
argument (which is very frequently but a finite-precision rational
approximation to a real in itself) is larger than pi, the best thing
to return is 'NaN', and most of the time the users should be trapping
on NaN's.

The program almost certainly has a bug. Doing otherwise is
indulging silliness, like declaring that dereferencing NULL gives zero
or something like that instead an error.

James Giles

2004-05-18, 2:32 pm

Dr Chaos wrote:
....
> It is a vendors responsibility to provide answers which lead to their
> client's programs working the best.
>
> I submit, that if a change in the least significant bit of the
> argument (which is very frequently but a finite-precision rational
> approximation to a real in itself) is larger than pi, the best thing
> to return is 'NaN', and most of the time the users should be trapping
> on NaN's.

Well, that argument would have more support if we were
using intervals, or significance arithmetic, or some other
means to denote the accuracy to which values were known.
Then the transition from significant to NaN wouldn't be just
a sudden jump. The answers would be larger and larger
intervals (or fewer and fewer digits) until the whole interval
[-1.0, 1.0] was returned (or, "no significance"). Without that,
where are you going to draw the line?

NaN would be a good representation for "no significance". But
how do you propose telling the naive user that there was some,
but not much, significance? And, within ordinary ranges, it's not
entirely foolish to want a full precision answer (so, whatever
your solution is, you need to allow real experts to turn it off).

--
J. Giles

glen herrmannsfeldt

2004-05-18, 5:33 pm

James Giles wrote:

(someone wrote)

In the days before NaN's it could have just been a fatal
error. In one IBM manual on the DSIN function, (comparable
to the C sin function)...

They limit the domain to +/- PI*2**50, which I believe leaves
four bits after the binary point (or one digit after the
hexadecimal point for purists). The then quote relative and
absolute accuracy figures as follows:

Relative Absolute
Max Std.Dev. Max. Std.Dev.

abs(X)<= PI/2 3.60e-16 4.82e-17 7.74e-17 1.98e-17
pi/2<abs(X)<=10 1.64e-16 6.49e-17
10<abs(X)<=100 2.68e-15 1.03e-15

So, in the range where they expect most users to use
the function they give reasonably detailed numbers. Do any
of the libraries more modern than this, from other companies,
give this information?

They also give details on the argument reduction used:

Divide abs(x) by pi/4 and separate into integer part (q),
and fraction part (r). 0<=r<1.

If cosine is desired, add 2 to q. If sine is desired and
x is negative add 4 to q.

Using q0= q mod 8,

for q0=0 use sin((pi/4)*r)
q0=1 cos((pi/4)*(1-r))
q0=2 cos((pi/4)*r)
q0=3 sin((pi/4)*(1-r))
q0=4 -sin((pi/4)*r)
q0=5 -cos((pi/4)*(1-r))
q0=6 -cos((pi/4)*r)
q0=7 -sin((pi/4)*(1-r))

Using polynomial interpolation of degree 6 in r**2 for sin,
or degree 7 in r**2 for cosine. In either case, the interpolation
points were the roots of the Chebyshev polynomial of one higher
degree. The maximum relative error for the sin polynomial
is 2**-58, of the cosine polynomial is 2**-64.3.

Not that all the details are important here, but the amount
of detail that went into the explanation. This is from
the Application Programming: Library Reference, meant to
be read by application programmers.

[color=darkred]
> Well, that argument would have more support if we were
> using intervals, or significance arithmetic, or some other
> means to denote the accuracy to which values were known.
> Then the transition from significant to NaN wouldn't be just
> a sudden jump. The answers would be larger and larger
> intervals (or fewer and fewer digits) until the whole interval
> [-1.0, 1.0] was returned (or, "no significance"). Without that,
> where are you going to draw the line?

In the case above, no accuracy numbers were given past 100,
and one can see how fast the accuracy decays at that point.

> NaN would be a good representation for "no significance". But
> how do you propose telling the naive user that there was some,
> but not much, significance? And, within ordinary ranges, it's not
> entirely foolish to want a full precision answer (so, whatever
> your solution is, you need to allow real experts to turn it off).

Real experts can do their own range reduction in the way
appropriate to the problem they are working on.

-- glen

keith bierman

2004-05-18, 6:31 pm

Gerry Thomas wrote:
> "keith bierman" <keith.bierman@sun.com> wrote in message
>
>
>
> You're assuming that only Sun ever implemented large argument reduction
I have no idea how you draw that conclusion from what I wrote.
for
> trigonometric: BSD libm had it

I was under the impression that BSD uses fdlibm (but I could be
misremembering. Some of the same people contributed to both).

In any event, I imagine we've beaten this horse to death again.

Dr Chaos

2004-05-18, 9:31 pm

James Giles <jamesgiles@worldnet.att.net> wrote:
> Dr Chaos wrote:
> ...
>
> Well, that argument would have more support if we were
> using intervals, or significance arithmetic, or some other
> means to denote the accuracy to which values were known.

> Then the transition from significant to NaN wouldn't be just
> a sudden jump. The answers would be larger and larger
> intervals (or fewer and fewer digits) until the whole interval
> [-1.0, 1.0] was returned (or, "no significance"). Without that,
> where are you going to draw the line?

where the change in a least significant bit equals pi is as good
as any.

> NaN would be a good representation for "no significance". But
> how do you propose telling the naive user that there was some,
> but not much, significance?

Obviously you can't.

Here's the realistic situation: if you get up to taking sin of a large
number, it is fairly likely it's going to get up to a larger number
still, because the program is misconceived. (Say adaptive numerical
integration of exp(-a*x)*sin(1/x^2) from 0 to 1.)

At that point you get NaN, trap it, and recognize the conceptual flaw.

> And, within ordinary ranges, it's not
> entirely foolish to want a full precision answer (so, whatever
> your solution is, you need to allow real experts to turn it off).

what is a mathematically sensible situation for needing the sin of
the floating point number closest to 10^100 to full accuracy?

> J. Giles

James Giles

2004-05-18, 11:31 pm

Dr Chaos wrote:
> James Giles <jamesgiles@worldnet.att.net> wrote:
....
>
> what is a mathematically sensible situation for needing the sin of
> the floating point number closest to 10^100 to full accuracy?

I was thinking more in terms of the representable number greater
than the arbitrary point at which you choose to return NaN.

--
J. Giles

glen herrmannsfeldt

2004-05-19, 12:31 am

James Giles wrote:

> Dr Chaos wrote:
(snip)

[color=darkred]
> I was thinking more in terms of the representable number greater
> than the arbitrary point at which you choose to return NaN.

It seems that VS Fortran has an arbitrary point, prints
a message some number of times before it is fatal,
(the number might be one), and, if it does return, returns
sqrt(2)/2. Have there been any complaints about it
here in the last 40 years or so that they have been doing
it? None that I have seen. (HFP doesn't have NaN.)

-- glen

Dr Chaos

2004-05-19, 4:31 am

James Giles <jamesgiles@worldnet.att.net> wrote:
> Dr Chaos wrote:
> ...
>
> I was thinking more in terms of the representable number greater
> than the arbitrary point at which you choose to return NaN.

My arbitrary point is that number,

real(kind(0.0d0)), parameter :: sindumbarg = X.XXXXXXXXDXX

where a variation in the least significant bit contributes at least a
fluctuation of pi in the argument.

What is a mathematically sensible situation of needing the sin,
highly accurately, of those specific floating point numbers which
are representable and .GT. sindumbarg?

And is this need greater than the general likelihood that people who
pass arguments .GT. sindumbarg very probably have a mathematical
conceptual bug in their code, or are trying to solve a problem in the
wrong way?

> J. Giles

James Giles

2004-05-19, 5:31 am

Dr Chaos wrote:
> James Giles <jamesgiles@worldnet.att.net> wrote:
>
> My arbitrary point is that number,
>
> real(kind(0.0d0)), parameter :: sindumbarg = X.XXXXXXXXDXX
>
> where a variation in the least significant bit contributes at least a
> fluctuation of pi in the argument.

OK, your arbitrary point is already, by your own argument, far
beyond sensible. Why so large? The problems you are trying
to protect people from still arise. Indeed, it's likely that many
(most) programs that are faulty in the way you describe never
get large enough to trigger the NaN. Yet, they are not giving
sensible numbers.

The problem with an arbitrary decision in any design is that
it's arbitrary. The argument for generating an error (or a NaN
in this case) is equally applicable to values you fail to even warn
about. What compelling reason can you give for that particular
value and not for some smaller one (or larger)?

--
J. Giles

jan van oosterwijk

2004-05-19, 12:32 pm

Rich Townsend <rhdt@barVOIDtol.udel.edu> wrote in message news:<c846f7$see$1@scrotar.nss.udel.edu>...[color=darkred]
> Jesper Harder wrote:

My compilers seem to do the same thing as Lahey, even for x > 1d14
See results of this program.

! [JvO] 2004-05-19 s1d100.f95
program s1d100
double precision :: d = 0.1d0
do while(d < 1d22)
print *, d, sin(d)
d = 10 * d
end do
d = 1.0001d100
print *, d, sin(d)
end

! D:\Fortran\Test>s1d100 ! Absoft
! 0.100000000000000 9.983341664682815E-002
! 1.00000000000000 0.841470984807897
! 10.0000000000000 -0.544021110889370
! 100.000000000000 -0.506365641109759
! 1000.00000000000 0.826879540532003
! 10000.0000000000 -0.305614388888252
! 100000.000000000 3.574879797201638E-002
! 1000000.00000000 -0.349993502171292
! 10000000.0000000 0.420547793190771
! 100000000.000000 0.931639027109679
! 1000000000.00000 0.545843449449778
! 10000000000.0000 -0.487506025076270
! 100000000000.000 0.928693660544335
! 1000000000000.00 -0.611238701357970
! 10000000000000.0 -0.288885282492284
! 100000000000000. -0.209408433383500
! 1.000000000000000E+015 0.858272132476373
! 1.000000000000000E+016 0.779679945161067
! 1.000000000000000E+017 -0.464644108935764
! 1.000000000000000E+018 -0.992816104053004
! 1.000000000000000E+019 1.000000000000000E+019
! 1.000000000000000E+020 1.000000000000000E+020
! 1.000000000000000E+021 1.000000000000000E+021
! 1.000100000000000E+100 1.000100000000000E+100
!
! D:\Fortran\Test>ftn95 s1d100 -lgo
!
! NO ERRORS [<main program> FTN95/Win32 v4.0] ! Salford
! Creating executable: D:\Fortran\Test\lgotemp@.exe
! Program entered
! 0.100000000000 9.983341664683E-02
! 1.00000000000 0.841470984808
! 10.0000000000 -0.544021110889
! 100.000000000 -0.506365641110
! 1000.00000000 0.826879540532
! 10000.0000000 -0.305614388888
! 100000.000000 3.574879797202E-02
! 1000000.00000 -0.349993502171
! 10000000.0000 0.420547793191
! 100000000.000 0.931639027110
! 1000000000.00 0.545843449450
! 10000000000.0 -0.487506025076
! 100000000000. 0.928693660544
! 0.100000000000E+13 -0.611238701358
! 1.000000000000E+13 -0.288885282492
! 1.000000000000E+14 -0.209408433383
! 1.000000000000E+15 0.858272132476
! 1.000000000000E+16 0.779679945161
! 1.000000000000E+17 -0.464644108936
! 1.000000000000E+18 -0.992816104053
! 1.000000000000E+19 1.000000000000E+19
! 1.000000000000E+20 1.000000000000E+20
! 1.000000000000E+21 1.000000000000E+21
! 1.000100000000E+0100 1.000100000000E+0100

[JvO] transfer((/ 778985834, 1869504886, 1702130543, &
1785296754, 1635205227, 1868849518, 1819160175 /), (/'x'/) )

Gerry Thomas

2004-05-19, 11:32 pm

"Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
news:slrncam3v0.dvo.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...

Fact:
sin(ultralargeargument) can be written as a finite polynomial in
sin(argumentin0to360deg) without having to know anything about pi, least of
all its value to any number of bits. A reference can be supplied to the
measure-zero discerning on CLF, but the obdurate rump is too busy at the
moment and couldn't be bothered right as is my entitlement (I've got the
union to back me in grievance, so don't mess please!)

> What is a mathematically sensible situation of needing the sin,
> highly accurately, of those specific floating point numbers which
> are representable and .GT. sindumbarg?
>
> And is this need greater than the general likelihood that people who
> pass arguments .GT. sindumbarg very probably have a mathematical
> conceptual bug in their code, or are trying to solve a problem in the
> wrong way?
>

Question:

Here's a chance to specify sindumbarg for the idiot's of the world you
await in baited breathed breath, ...

If one has a system of ODE's with periodic coefficients (it could be a NASA
tin can replete with an expendable crew, and what the heck if anything goes
wrong we can blame it all on 'management' and continue our sinecure duties
in pajamas from the local Starbucks, etc), as you already know its
solutions can be written in Floquet form and they are not necessarily
periodic. In fact if the SL operator is bounded, there are a finite number
of solutions for arbitrary initial conditions, and they're 'almost
periodic' in compliance with the Hogg-Huberman theorem (cf. Liouville
theorem) . There is no _a priori_ prescription for the recurrence time of
the solutions, which in principal might be an astronomical multiple of the
fundamental period. Ergo the bankruptcy of your all embracing sindumbarg.
If one has to integrate way beyond the fundamental period, one can cascade
the solution a la Floquet once the integral and the Lyapunov exponent
matrices are known. However, assuming the Floquet matrix is normal (because
the SL is self adjoint, for instance), in finite precision it might whittle
away (cf. Highams book, reference to a SIAM paper by Higham and Knight,
IIRC) to naught before one ever reaches the terminal point in which case
the show is over, or perhaps it isn't in which case one's math library had
better be able to reliably deliver trigonometric without the imposition of
sindumbarg. The point of KC Ng's Sun Microsystems report is that it can and
ought to be done. And it's available freely in the C library (always to the
rescue for the inadequate Fortran!) fdlibm. If a compiler doesn't do so
then, IMO, Giles' suggestion is preferable provided of course that a NaN
can be issued no thanks to retrovariant Fortran. Plugger's approach for C99
doesn't issue a NaN or anything else, for C or Fortran, but rather a random
in the range -1 to +1, user silly user beware. Ditto for Borland C/C++,
take a look at their CRT source for their contemptuos attitude of their
library writers have for joe user. Some Fortran compilers are no better but
then again Fortran's sham claim to superiority in matters numeric has long
been exposed and only its defenders and their groveling sycophantic toadies
persist in their denial.

--
E&OE

Ciao,
Gerry T.
______
"Some of those writing letters to this program are slightly fictitious but
others definitely are not." -- Dick Cavett, host of The Detroit Symphony
Orchestra, PBS Radio.

Dr Chaos

2004-05-20, 4:31 am

On Wed, 19 May 2004 22:29:11 -0400, Gerry Thomas <gfthomas@sympatico.ca> wrote:
>
> Here's a chance to specify sindumbarg for the idiot's of the world you
> await in baited breathed breath, ...
>
> If one has a system of ODE's with periodic coefficients (it could be a NASA
> tin can replete with an expendable crew, and what the heck if anything goes
> wrong we can blame it all on 'management' and continue our sinecure duties
> in pajamas from the local Starbucks, etc), as you already know its
> solutions can be written in Floquet form and they are not necessarily
> periodic. In fact if the SL operator is bounded, there are a finite number
> of solutions for arbitrary initial conditions, and they're 'almost
> periodic' in compliance with the Hogg-Huberman theorem (cf. Liouville
> theorem) . There is no _a priori_ prescription for the recurrence time of
> the solutions, which in principal might be an astronomical multiple of the
> fundamental period. Ergo the bankruptcy of your all embracing sindumbarg.
> If one has to integrate way beyond the fundamental period, one can cascade
> the solution a la Floquet once the integral and the Lyapunov exponent
> matrices are known. However, assuming the Floquet matrix is normal (because
> the SL is self adjoint, for instance), in finite precision it might whittle
> away (cf. Highams book, reference to a SIAM paper by Higham and Knight,
> IIRC) to naught before one ever reaches the terminal point in which case
> the show is over, or perhaps it isn't in which case one's math library had
> better be able to reliably deliver trigonometric without the imposition of
> sindumbarg.

that's interesting, but wouldn't most useful results here also depend
on being able to represent the argument of the trigonometric with high
absolute accuracy and not relative accuracy?

And really in such a circumstance, why wouldn't people be integrating
ODEs as integrating ODE's?

For instance the simplest case of accumulating many upon many
rotations as complex rotations?

And in such a circumstance with sensitive dependence on initial conditions
(which can happen even without positive kolmogorov-sinai entropy) is
some exact method but for poorly represented arguments really valuable?

A super fancy sin() is a special case of an arbitrary precision
facility. Which is necessary sometimes.

> The point of KC Ng's Sun Microsystems report is that it can and
> ought to be done. And it's available freely in the C library (always to the
> rescue for the inadequate Fortran!) fdlibm. If a compiler doesn't do so
> then, IMO, Giles' suggestion is preferable provided of course that a NaN
> can be issued no thanks to retrovariant Fortran. Plugger's approach for C99
> doesn't issue a NaN or anything else, for C or Fortran, but rather a random
> in the range -1 to +1, user silly user beware. Ditto for Borland C/C++,
> take a look at their CRT source for their contemptuos attitude of their
> library writers have for joe user. Some Fortran compilers are no better but
> then again Fortran's sham claim to superiority in matters numeric has long
> been exposed and only its defenders and their groveling sycophantic toadies
> persist in their denial.

WTF?

>
> --
> E&OE
>
> Ciao,
> Gerry T.
> ______
> "Some of those writing letters to this program are slightly fictitious but
> others definitely are not." -- Dick Cavett, host of The Detroit Symphony
> Orchestra, PBS Radio.
>
>

glen herrmannsfeldt

2004-05-20, 6:31 am

Dr Chaos wrote:

(snip)

> A super fancy sin() is a special case of an arbitrary precision
> facility. Which is necessary sometimes.

Yes. If one has an arbitrary precision sin() then argument
reduction over larger exponent values makes more sense.

There is a claim that Riemann Zeta function zeros can
be calculated using low precision sin() with high precision
argument reduction. I haven't verified that one way or the other.

-- glen

Gerry Thomas

2004-05-20, 6:31 am

"Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
news:slrncaomvm.dm9.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...
> On Wed, 19 May 2004 22:29:11 -0400, Gerry Thomas <gfthomas@sympatico.ca>
wrote:
NASA[color=darkred]
goes[color=darkred]
duties[color=darkred]
number[color=darkred]
of[color=darkred]
the[color=darkred]
sindumbarg.[color=darkred]
cascade[color=darkred]
(because[color=darkred]
whittle[color=darkred]
case[color=darkred]
had[color=darkred]
of[color=darkred]
>
> that's interesting,

You bet buddy, so watch it!, even on c.l.f. - the soap Opera of computing
newsgroups - there are signs of intelligent life, yours truly included.

>but wouldn't most useful results here also depend
> on being able to represent the argument of the trigonometric with high
> absolute accuracy and not relative accuracy?

Well that's the point, or in your vernacular WTF?, and who defines what you
dub 'most useful results'? I'll decide that for me, thank you, so stick
sindumbarg, unless you want to underwrite my malpractise insurance.?

>
> And really in such a circumstance, why wouldn't people be integrating
> ODEs as integrating ODE's?
>

This is tautological nonsense. However, while I'm forrgiving of speelling
mistaxes and whatnot widout drawing the ilogix arbitrary interference (sic)
of one Sargent. (yes, he's been demoted from his prior and unwarrented
exalted pedestal to a lower NCO (sic unprofessional) rank, he being a card
carrying, possibly commie, union member) Maine woudl claim that you are
offe your rockker and were on some kind of drugggs (remember 'Just say no
to drugggs' as Nancy R. advised), maybe even coffee, what would Betty
Bowers say?, praise the Lordy Lord!

> For instance the simplest case of accumulating many upon many
> rotations as complex rotations?
>
> And in such a circumstance with sensitive dependence on initial
conditions
> (which can happen even without positive kolmogorov-sinai entropy) is
> some exact method but for poorly represented arguments really valuable?
>
> A super fancy sin() is a special case of an arbitrary precision
> facility. Which is necessary sometimes.
>

You're adopting a predictable barndoor parochial stance, consistent with
your handle (Chaos) and that's about it.

Consider the rapid eye movement of schizophrenics or the microwave
excitation of a molecular entity. In either case, the initial conditions
are not in doubt, chaos, as in uncertainty in the initial conditions and
its impact on subsequent evolution, doesn't enter into the discussion and
can't be offered as an excuse: you can either match long term observation
with simulation (give or take ergodicy, decoherence, or whatever you want
to call it) or you can't, and when you can't it's due to a failure in your
compiler's ability to render reliable trigonmetrics. When I first
encountered this with CVF I was essentially told that I was wrong by
someone who will now claim that with IVF I was right after all. What comes
around goes around, bullshit or truth, but it at least eliminates the
intrusion of Sgt. Maine from the discourse and his fatous contribution to
this or indeed any rational discussion beyond the insular issue of Fortran
portability, as if anyone gives a fcuk. I don't, nor do many.

>
the[color=darkred]
NaN[color=darkred]
C99[color=darkred]
random[color=darkred]
but[color=darkred]
long[color=darkred]
toadies[color=darkred]
>
> WTF?
>

You just don't get it, do you?, never mind.

--
Ciao,
Gerry T.
______
"...as I further could have told you as brisk as your D.B.C.
behaviouristically paillet, with a coat of homoid icing which is in reality
only a done by chance ridiculisation of the whoo-whoo and where's hairs
theorics of Winestain." -- James Joyce, in Finnegans Wake.

Gerry Thomas

2004-05-20, 6:31 am

"Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
news:slrncaomvm.dm9.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...

[...]

> A super fancy sin() is a special case of an arbitrary precision
> facility. Which is necessary sometimes.
>

I missed this, besides it's a nonsensical claim, but never mind.

--
E&OE

Ciao,
Gerry T.
______
"Klinger, it's my considered opinion that no one is going to believe you
are pregnant." -- Col Henry Blake, the 4077th M*A*S*H

Dr Chaos

2004-05-20, 2:32 pm

On Thu, 20 May 2004 04:53:18 -0400, Gerry Thomas <gfthomas@sympatico.ca> wrote:
>
> "Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
> news:slrncaomvm.dm9.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...
> wrote:
> NASA
> goes
> duties
> number
> of
> the
> sindumbarg.
> cascade
> (because
> whittle
> case
> had
> of
>
> You bet buddy, so watch it!, even on c.l.f. - the soap Opera of computing
> newsgroups - there are signs of intelligent life, yours truly included.
>
>
> Well that's the point, or in your vernacular WTF?, and who defines what you
> dub 'most useful results'? I'll decide that for me, thank you, so stick
> sindumbarg, unless you want to underwrite my malpractise insurance.?

If you pass a double precision floating point number in standard
representation to the sin() standard function then when its
argument is very large in magnitude

>
>
> This is tautological nonsense.

No it isn't.

I thought you were discussing quasiperiodic solutions to systems
of equations: you might be able to write down solutions as
sum of sinusoidal functions, or integrate the equations of motion
using typical (or atypical) schemes for integrating ODE's, conservative
or not.

> However, while I'm forrgiving of speelling
> mistaxes and whatnot widout drawing the ilogix arbitrary interference (sic)
> of one Sargent. (yes, he's been demoted from his prior and unwarrented
> exalted pedestal to a lower NCO (sic unprofessional) rank, he being a card
> carrying, possibly commie, union member) Maine woudl claim that you are
> offe your rockker and were on some kind of drugggs (remember 'Just say no
> to drugggs' as Nancy R. advised), maybe even coffee, what would Betty
> Bowers say?, praise the Lordy Lord!
>
> conditions
>
> You're adopting a predictable barndoor parochial stance, consistent with
> your handle (Chaos) and that's about it.
>
> Consider the rapid eye movement of schizophrenics or the microwave
> excitation of a molecular entity. In either case, the initial conditions
> are not in doubt, chaos, as in uncertainty in the initial conditions and
> its impact on subsequent evolution, doesn't enter into the discussion and
> can't be offered as an excuse: you can either match long term observation
> with simulation (give or take ergodicy, decoherence, or whatever you want
> to call it) or you can't, and when you can't it's due to a failure in your
> compiler's ability to render reliable trigonmetrics.

Yes, of course you need to match statistical quantities.

In what way would that depend intimately on being able to evaluate
trigonometrics of very large arguments, still represented in standard
finite precision floating point, large enough so that one LSB change
in the argument results in large changes in the value?

Gerry Thomas

2004-05-22, 10:31 pm

"Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
news:slrncaps75.itk.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...

[...]

>
> I thought you were discussing quasiperiodic solutions to systems
> of equations: you might be able to write down solutions as
> sum of sinusoidal functions, or integrate the equations of motion
> using typical (or atypical) schemes for integrating ODE's, conservative
> or not.
>

[...]

>
> In what way would that depend intimately on being able to evaluate
> trigonometrics of very large arguments, still represented in standard
> finite precision floating point, large enough so that one LSB change
> in the argument results in large changes in the value?
>

If pi isn't known beyond the working precision you're SOL and one sees that
with the results reported in this thread for numerous compilers. Nature
doesn't have this problem so mere mortals, and so called computer
scientists in particular (although RPF didn't believe they rated as
scientists, the object of their interest being of their own design), must
be doing something wrong.

--
E&OE

Ciao,
Gerry T.
______
"Facts are meaningless. You could use facts to prove anything that's even
remotely true." -- Homer Simpson.

Dr Chaos

2004-05-23, 9:35 pm

On Sat, 22 May 2004 21:10:37 -0400, Gerry Thomas <gfthomas@sympatico.ca> wrote:
>
> "Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
> news:slrncaps75.itk.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...
>
> [...]
>
>
> [...]
>
>
> If pi isn't known beyond the working precision you're SOL and one sees that
> with the results reported in this thread for numerous compilers.

That's an answer to a different question. I too think it's desirable
to get trigonometrics to high accuracy for reasonably small arguments.

the original issue is this:

"In what way would that depend intimately on being able to evaluate
trigonometrics of very large arguments, still represented in standard
finite precision floating point, large enough so that one LSB change
in the argument results in large changes in the value?

> Nature
> doesn't have this problem so mere mortals, and so called computer
> scientists in particular (although RPF didn't believe they rated as
> scientists, the object of their interest being of their own design), must
> be doing something wrong.

Nature doesn't have to approximate reals with floating point either.

I question the wisdom of doing a tiny part of the calculation to
arbitrary accuracy and none of the rest of it, as in when floating
point numbers of very large magnitude are passed to standard library
subroutines. The trig routine was already passed something which
is almost completely 'wrong' (if you care about sins and cosines of it).

(hence the thread name "sin(1d100)"

> Ciao,
> Gerry T.

Gerry Thomas

2004-05-24, 11:38 pm

"Dr Chaos" <mbkennelSPAMBEGONE@NOSPAMyahoo.com> wrote in message
news:slrncb2fqi.igl.mbkennelSPAMBEGONE@lyapunov.ucsd.edu...
> On Sat, 22 May 2004 21:10:37 -0400, Gerry Thomas <gfthomas@sympatico.ca>
wrote:
conservative[color=darkred]
that[color=darkred]
>
> That's an answer to a different question. I too think it's desirable
> to get trigonometrics to high accuracy for reasonably small arguments.
>

We're in agreement here but even sin(pi) /= 0, to say nothing of lim
sin(Npi), N->oo.

> the original issue is this:
>
> "In what way would that depend intimately on being able to evaluate
> trigonometrics of very large arguments, still represented in standard
> finite precision floating point, large enough so that one LSB change
> in the argument results in large changes in the value?
>
must[color=darkred]
>
> Nature doesn't have to approximate reals with floating point either.
>

True, that's why She's so successful: no ISO and it's toadying minions and
gofers to waste time on.

> I question the wisdom of doing a tiny part of the calculation to
> arbitrary accuracy and none of the rest of it, as in when floating
> point numbers of very large magnitude are passed to standard library
> subroutines.

Production quality numeric libraries are not a pretty sight: they're ugly,
ad hoc, empirical, etc., but efficient and ought to be reliable (sic
deliver trustworthy results, but don't count on Fortran to guarantee this
because it doesn't and would appear to be quite proud of its inadequacies
in this regard).

> The trig routine was already passed something which
> is almost completely 'wrong' (if you care about sins and cosines of it).

In finite fp arithmetic the majority of arguments are irrepresentably
'wrong': large arg trigs are no exception except that here reliable
libraries don't need user assistance in rendering accurate results via
perhaps the most hedious of numerical schemes.

To close, this topic would appear to detract from the usual mundanity of
such profound Fortran issues as to tab or not?, what path is this?, how to
read and write?, where can I get whatever for free or cheap, etc., but best
of all how do I translate this from Fortran to anything besides.

--
E&OE

Ciao,
Gerry T.
______
"In a world in which the price of calculation continues to decrease
rapidly, but the price of theorem proving continues to hold steady or
increase, elementary economics indicates that we ought to spend a larger
and larger fraction of our time on calculation." -- J.W. Tukey.

2004-06-03, 7:24 pm

Gerry Thomas wrote:
>
> To close, this topic would appear to detract from the usual mundanity of
> such profound Fortran issues as to tab or not?, what path is this?, how to

No, it ought to stay open until we find out how the compilers ended up
in two (-.567 & -.380) numeric camps, and some good soul tips off Lahey
that,
sin(x) ~ x

is for small args, not the other way around.

keith bierman

2004-06-03, 7:24 pm

bv wrote:
> Gerry Thomas wrote:
>
>
> No, it ought to stay open until we find out how the compilers ended up
> in two (-.567 & -.380) numeric camps, and some good soul tips off Lahey
> that,
> sin(x) ~ x
>
> is for small args, not the other way around.
>

Doesn't seem like much of a mystery. Clearly:

*) Some vendors use "extreme methods" to compute with enough Pi to get
the best answer possible (that is, assume the user had some rationale
reason for the computation, and deliver the best result based on that
assumption).

*) Some vendors assume that such input is the input; use the same
approximation of Pi for whatever input (and one which is "good enough"
for "sensible" input) and let the bits fall where they may.

*) Some vendors assume that such input is essentially confession by the
user that the input is meaningless, and therefore is free to output some
arbitrary value(s) presumably based on implementation convenience.

equally, from this discussion, it would appear that some users fall into
similar camps. Hopefully the users who think a particular method is
appropriate can match themselves up with a compiler/runtime which
matches their predisposition ;>

beliavsky@aol.com

2004-06-03, 7:24 pm

keith bierman <keith.bierman@sun.com> wrote:

>*) Some vendors assume that such input is essentially confession by the

>user that the input is meaningless, and therefore is free to output some

>arbitrary value(s) presumably based on implementation convenience.

I do not agree with this approach, except perhaps as a debugging option,
in which case it might be better to print an error message and stop. In math
abs(sin(x)) <= 1 for all real x, and it should not be otherwise in Fortran.
I understand that calculating sin(x) for very large x may not be feasible,
but at least some value between -1 and 1 should be produced. I think I have
a "right" to expect that

y = 1/(2 - sin(x))

will not trigger a division by zero error, for any possible value of x.
Otherwise a defensive programmer would have to clutter his program with error
checks for any expression with a sine or cosine in the denominator.

----== Posted via Newsfeed.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeed.com The #1 Newsgroup Service in the World! >100,000 Newsgroups
---= 19 East/West-Coast Specialized Servers - Total Privacy via Encryption =---

glen herrmannsfeldt

2004-06-03, 7:24 pm

beliavsky@aol.com wrote:

(snip)

> I do not agree with this approach, except perhaps as a debugging option,
> in which case it might be better to print an error message and stop. In math
> abs(sin(x)) <= 1 for all real x, and it should not be otherwise in Fortran.
> I understand that calculating sin(x) for very large x may not be feasible,
> but at least some value between -1 and 1 should be produced.

I believe that it is using the hardware SIN function without
any library help in between. It may not be hard to do in software,
but it is somewhat harder in hardware.

-- glen

Dr Chaos

2004-06-03, 7:24 pm

On 25 May 2004 16:38:30 -0500, beliavsky@aol.com <beliavsky@127.0.0.1> wrote:
>
> keith bierman <keith.bierman@sun.com> wrote:
>
>
>
>
> I do not agree with this approach, except perhaps as a debugging option,
> in which case it might be better to print an error message and stop. In math
> abs(sin(x)) <= 1 for all real x, and it should not be otherwise in Fortran.
> I understand that calculating sin(x) for very large x may not be feasible,
> but at least some value between -1 and 1 should be produced. I think I have

I think a NaN is the best out-of-band result, and most of the time
people should be trapping on NaNs.

> a "right" to expect that
>
> y = 1/(2 - sin(x))
>
> will not trigger a division by zero error, for any possible value of x.
> Otherwise a defensive programmer would have to clutter his program with error
> checks for any expression with a sine or cosine in the denominator.

A numerical result outside -1 to 1 is not desirable.

keith bierman

2004-06-03, 7:24 pm

glen herrmannsfeldt wrote:
> beliavsky@aol.com wrote:

>
>