Issue-#157 Extending the domain of Math.Sin, Cos, & Tan

[Robert]

I propose extending the effective domain of Sin, Cos, & Tan; see https://redmine.blackboxframework.org/issues/157.
Currently for angles above 1.E18 the code returns 0. for Sin & 1. for Cos. The FPU provides accurate answers upto angles of over 9.22E18. I illustrate this improved accuracy for one angle in this range: 5.E18 (which can be represented exactly as a 64-bit REAL).

With the proposed change:
Sin (5.E18) = −0.825512
Cos (5.E18) = 0.564385
Tan (5.E18) = −1.462675

The precise answers are:
Sin (5.E18) = −0.829128
Cos (5.E18) = 0.559059

Sin (5000000000000000000.006437) = −0.825512
Cos (5000000000000000000.006437) = 0.564385
Tan (5000000000000000000.006437) = −1.462675

Note that an ULP (Unit in the Last Place) is 1024, so the approximate result is in fact the exact result for some x within 6 millionths of an ULP of the argument 5.E18.

This shows that the FPU does provide a good (useful) accuracy over this range, and accessing it has no cost (in fact it simplifies and speeds the code) so I hope you will support this proposal. I will post the changed Mod files shortly; no changes to Docu or Rsrc files are anticipated since the current unnecessary limitation is not documented.

To avoid merge conflicts I will wait until issue 151 is resolved; the proposed changes are illustrated below:

Code: Select all

	PROCEDURE Sin* (x: REAL): REAL;
	BEGIN
		(* 20, ABS(x) # INF *)
		FLD(x); FSIN; WAIT;
		IF 10 IN FSWs() THEN RETURN 0. END;
		RETURN TOP()
	END Sin;

[Josef Templ]

I propose extending the effective domain of Sin, Cos, & Tan; see https://redmine.blackboxframework.org/issues/157.
Currently for angles above 1.E18 the code returns 0. for Sin & 1. for Cos. The FPU provides accurate answers upto angles of over 9.22E18. I illustrate this improved accuracy for one angle in this range: 5.E18 (which can be represented exactly as a 64-bit REAL).

Just for curiosity:
1. In which cases is it important to support those very big arguments?
2. Why is it acceptable to return 0 for those very big arguments?
3. Is this common practice in numerics? Is it the same in say GO, Java, C++?

- Josef

[Josef Templ]

Code: Select all

	PROCEDURE Sin* (x: REAL): REAL;
	BEGIN
		(* 20, ABS(x) # INF *)
		FLD(x); FSIN; WAIT;
		IF 10 IN FSWs() THEN RETURN 0. END;
		RETURN TOP()
	END Sin;

I think the trap number "20" in the comment should be changed to "143" in order to correspond with
the real runtime trap number and the entry in the System/Rsrc/Strings fle.
Same for Cos and Tan, I assume.

- Josef

[Robert]

1. In which cases is it important to support those very big arguments?

No idea - probably never!. Maybe we never need arguments greater than 1.E17; should we use that as an arbitrary cut-off point?
My concerns are two-fold: I don't like the arbitrary nature of the 1.E18 limit, and I have been nervous, for years, about the comment in Reduce "(* to be completed *)". We do require pretty large arguments (eg 1.E8), so reasonable argument reduction performance is needed. So, with the recent interest in this module, I investigated a bit closer, and found that the hardware is quite safe upto 2^63, and this distracting comment can be removed. The side-effect is smaller, simpler, faster code. But the main motivation is to remove the confidence removing comment.

2. Why is it acceptable to return 0 for those very big arguments?

I assume you mean > 2^63.
In some cases it won't be, but covering those very special cases is a new feature that would be expensive to implement, and no one is asking for it. The limited extension proposed is very cheap - in fact negative cost.

3. Is this common practice in numerics? Is it the same in say GO, Java, C++?

I don't know. I could investigate MATLAB in a few days.

I am (now) also proposing adding a new procedure SinCos. This does bring significant tangible benefits to signal processing (and other) applications: both simpler source code and significant speed advantages.
The corrected Sin suggestion is

Code: Select all

	PROCEDURE Sin* (x: REAL): REAL;
	BEGIN
		(* 20, ABS(x) # INF *)
		FLD(x); FSIN; WAIT;
		IF 10 IN FSWs() THEN FSTPst0; RETURN 0. ELSE  END;
		RETURN TOP()
	END Sin;

The suggested SinCos procedure is

Code: Select all

	PROCEDURE SinCos* (x: REAL; OUT s, c: REAL);
	BEGIN
		(* 20, ABS(x) # INF *)
		FLD(x); FSINCOS; WAIT;
		IF 10 IN FSWs() THEN FSTPst0; s := 0.; c := 1. ELSE c := TOP(); s := TOP() END
	END SinCos;

[Robert]

I think the trap number "20" in the comment should be changed to "143" in order to correspond with
the real runtime trap number and the entry in the System/Rsrc/Strings fle.
Same for Cos and Tan, I assume.

I agree the number 20 in both the code comment and the Docu has no apparent significance; and I think it should be removed.

I don't understand the 143; is it explained anywhere? It does not appear in the Trap window. I would not put it into either the comment or the Docu; I would have no number.

Some comments, and Strings entries, have ABS(x) # INF, some x # -INF or INF. Should this be rationalised?

[Robert]

3. Is this common practice in numerics? Is it the same in say GO, Java, C++?

I don't know. I could investigate MATLAB in a few days.

Just done some MATLAB tests (on a Mac).

Code: Select all

Sin(5.E18) = -0.8291276
Cos(5.E18) =  0.5590594
Sin(5.E19) =  0.3435328
Cos(5.E19) = -0.9391406

And some Maxima tests

Code: Select all

Sin(5.E18) = -0.8255117
Cos(5.E18) =  0.5643849
Sin(5.E19) = -0.8121308
Cos(5.E19) = -0.5834753

Mystery on mystery!
The MATLAB results are all correct to the 7 digits quoted. It must be using multi-precision of some kind, which must be slow?
The Maxima 5.E18 results agree with the FPU (ie proposed BlackBox) results to the 7 digits recorded.
I have no idea where the Maxima 5.E19 numbers come from. This is the native Maxima floating point format.
(If I use the Maxima "bigfloat" format I get the correct answers, presumably to many thousand digits if I ask!)

A Maxima 5.E19 theory - it is doing 64-bit software range reduction - anyway these results (ie above 2^63) are useless!

[Robert]

The Math Docu includes some definitions of additional "Elementary" transcendental functions that are not included in the library.
One of the definitions (ArcCot) is not universally agreed, and I don't think we should mention it in isolation.

A different definition has strong support (NISDT; the US National Institute of Standards and Technology, amongst others), so I think it ought to be mentioned, if not recommended.

Some wording for discussion:

ArcCot (x) = Math.Pi() / 2.0 - Math.ArcTan(x)
(* Refs: https://en.wikipedia.org/wiki/Inverse_t ... _functions, Maple *)

ArcCot (x) = Math.ArcTan(1.0 / x)
(* Refs: http://dlmf.nist.gov/search/search?q=arccot, MATLAB, Mathematica, Maxima, MuPAD, Wikipedia *)

Although a slightly different topic, this is a minor matter, and I see no harm in including it in this issue.

[Josef Templ]

I think the trap number "20" in the comment should be changed to "143" in order to correspond with
the real runtime trap number and the entry in the System/Rsrc/Strings fle.
Same for Cos and Tan, I assume.

I agree the number 20 in both the code comment and the Docu has no apparent significance; and I think it should be removed.

I don't understand the 143; is it explained anywhere? It does not appear in the Trap window. I would not put it into either the comment or the Docu; I would have no number.

You must look into the System/Rsrc/Strings file and you will immediately see the significance of the number 143.
This must not be removed but changed from 20 to 143 in order to be consistent with the Kernel's error number
and the Strings file.
FYI: the mentioned precondition ABS(x) # INF is checked implicitly by the FSIN instruction and produces Trap 143.

I agree that adding SinCos is a good idea.
But the implementation has a bug, I think.
If it produces two values on the FPU stack then there must be two pops
in order to remove both values in THEN. (This is only a theory. I did not do any tests.)
Please check.

Otherwise the implementation looks OK to me.

- Josef

[Robert]

[I agree that adding SinCos is a good idea.
But the implementation has a bug, I think.
If it produces two values on the FPU stack then there must be two pops
in order to remove both values in THEN. (This is only a theory. I did not do any tests.)
Please check.

I don't really understand this assembler sufficiently to give a convincing answer (and I mean convincing me, never mind anyone else).
The manual says:

Description
Computes both the sine and the cosine of the source operand in register ST(0),
stores the sine in ST(0), and pushes the cosine onto the top of the FPU register stack.
(This instruction is faster than executing the FSIN and FCOS instructions in succes-
sion.)

which might mean it only stacks one value, so only one pop is required?

I have tested it, but not very heavily yet. I intend to modify my FFT algorithms to use it. Since (some of) my multi-precision arithmetic uses FFTs, I find that computing PI to 1000000 places flushes out most programming errors.

[Robert]

You must look into the System/Rsrc/Strings file and you will immediately see the significance of the number 143.

Yes - I understood that the Strings file used the 143, but I didn't understand where it came from.

This must not be removed but changed from 20 to 143 in order to be consistent with the Kernel's error number
and the Strings file.
FYI: the mentioned precondition ABS(x) # INF is checked implicitly by the FSIN instruction and produces Trap 143.

So, I now think the source code comment would confuse me less if it said "(* Kernel error 143, ABS(x) # INF *)".

Regarding the Docu pre-condition I still think that there should be no number. Neither 20 or 143 is visible to the non system programmer who uses the Math.Sin function.