Essays/Extremal Arguments
We consider the values of some algebraic and transcendental functions on extremal arguments. The arguments and results are represented as 64-bit IEEE floating point numbers, with the letters x and y denoting real numbers while the letter z denoting complex numbers.
Utilities
sin =: 1&o. cos =: 2&o. sinh =: 5&o. cosh =: 6&o. hex =: (4+#@$) }. 2&(3!:3) unhex=: (3!:2) @ ((}: 2 (3!:3) 0.5)&,) " 1
The smallest positive IEEE number is 0000000000000001 in hex and the largest is 7fefffffffffffff , respectively 2^_1074 and (2-2^_52)*2^1023 .
unhex (15$'0'),'1' 4.94066e_324 2^_1074 4.94066e_324 unhex '7fe',13$'f' 1.79769e308 (2-2^_52)*2^1023 1.79769e308 e=: 2^_1074 E=: (2-2^_52)*2^1023 hex e,E 0000000000000001 7fefffffffffffff
sin and cos on real arguments
sin and cos on real arguments have a period of 2p1 . When y is sufficiently large that y+2p1 has the same representation as y , then 2p1|!.0 x is 0 for all x>:y and sin x and cos x are meaningless. (sin x is 0 and cos x is 1 for all such x .)
The following demonstrates that y=:2^56 is the smallest number such that y+2p1 has the same representation as y .
] y=: 2^56 7.20576e16 y+2p1 7.20576e16 y - y+2p1 0 hex y+0 2p1 4370000000000000 4370000000000000 hex (unhex '436',13$'f')+0 2p1 436fffffffffffff 4370000000000000
2p1 |!.0 x is 0 for all x>:y :
] x=. unhex ('437',12$'0'),"1 0 '0123456789abcdef'
7.20576e16 7.20576e16 7.20576e16 7.20576e16 7.20576e16 7.20576e16 ...
hex x
4370000000000000
4370000000000001
4370000000000002
4370000000000003
4370000000000004
4370000000000005
4370000000000006
4370000000000007
4370000000000008
4370000000000009
437000000000000a
437000000000000b
437000000000000c
437000000000000d
437000000000000e
437000000000000f
x - {. x
0 16 32 48 64 80 96 112 128 144 160 176 192 208 224 240
2p1 |!.0 x
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
The J interpreter stops well short of the 2^56 limit. A limit error is signalled when 2p1 |!.0 x loses half of the number of bits of precision, that is, when the magnitude of x is greater than or equal to 2p1 * 2^53%2 where 53 is the number of bits in the significant (including the implied leading 1 bit). Thus:
thmax=: 2p1 * 2^53%2 thmax 5.96314e8 sin 6e8 |limit error: sin | sin 1000000000 cos _2e9 |limit error: cos | cos _2000000000 sin _ |limit error: sin | sin _
Exponential
Mathematically, ^x is positive for all x . However, if ^x is smaller than the smallest positive IEEE number e , that is, if x < ^.e , then ^x is represented as 0 .
e 4.94066e_324 emin=: ^. e emin _744.44 _1074 * ^. 2 _744.44 ^ emin 4.94066e_324 hex emin c0874385446d71c3 x=. unhex '4' ,~ _1}. ,hex emin NB. the next smaller number hex x c0874385446d71c4 ^ x 0
At the other extreme, if ^x is larger than the largest IEEE number E , that is, if x > ^.E , then ^x is represented as _ (infinity).
E 1.79769e308 emax=: ^.E emax 709.783 1024 * ^.2 709.783 hex emax,1024*^.2 40862e42fefa39ef 40862e42fefa39ef ^ emax 1.79769e308 hex emax 40862e42fefa39ef x=. unhex 'f0' ,~ _2}. ,hex emax NB. the next larger number hex x 40862e42fefa39f0 ^ x _
Now we consider complex arguments. From the identity
(^ x j. y) = (^x) * (cos y) j. (sin y)
If x<emin , ^x is (represented as) 0 and the product is 0 . Otherwise, if thmax<:|y , cos y and sin y are not accurate and a limit error is signalled.
If thmax>|y and x>:emin, the above identity is equivalent to:
(^x j. y) = ((*cos y)*^x+^.|cos y) j. ((*sin y)*^x+^.|sin y)
The expression (*cos y)^x+^.|cos y is computationally more useful than (cos y)*^x for x>emax (similarly for the imaginary part).
^ 716j1.57 7.17696e307j_ emax 709.783 ^ 716 _ ^ 716 + ^. cos 1.57 7.17696e307 ^ 716 + ^. sin 1.57 _
The procedure can be modelled as follows:
expcheck=: 3 : 0 'x y'=. +. y if. x<emin do. 0j0 return. end. assert. thmax>|y if. x<:emax do. (^x) * (cos y) j. (sin y) else. ((*cos y)*^x+^.|cos y) j. ((*sin y)*^x+^.|sin y) end. ) expc 3j4 _13.1288j_15.2008 ^ 3j4 _13.1288j_15.2008 expc 716j1.57 7.17696e307j_ ^ 716j1.57 7.17696e307j_ expc 1j1e9 |assertion failure: exp | thmax>|y ^ 1j1e9 |limit error | ^1j1e9 expc _800j1e9 0 ^ _800j1e9 0
sinh and cosh
By definition,
(sinh y) = 2 %~ (^y) - ^-y (cosh y) = 2 %~ (^y) + ^-y
For large magnitude y, the computation comes down to 2%~^|y . From the Exponential section above it is seen that if emax<(|y)-^.2 , ^(|y)-^.2 exceeds the largest IEEE number E . That is, if (|y)>emax+^.2 , then sinh y and cosh y are _ .
emax2=: emax+^.2 emax2 710.476 E +&^. 2 710.476 sinh _1 1*emax2 _1.79769e308 1.79769e308 cosh _1 1*emax2 1.79769e308 1.79769e308 hex emax2 408633ce8fb9f87d x=. unhex 'e',~}:,hex emax2 NB. the next larger number hex x 408633ce8fb9f87e sinh _1 1 * x __ _ cosh _1 1 * x _ _
For complex arguments, the identities
(sinh z) = sin&.j. z NB. A&S 4.5.7 (cosh z) = cos& j. z NB. A&S 4.5.8
reduces sinh and cosh to the trigonometric versions. (See below.)
sin and cos
(sin x j. y) = ((sin x) * cosh y) j. (cos x) * sinh y NB. A&S 4.3.55 (cos x j. y) = ((cos x) * cosh y) j. - (sin x) * sinh y NB. A&S 4.3.56
Since cosh y and sinh y can not be both 0 , for large magnitude x either the real part or the imaginary part (or both) of sin x j. y and cos x j. y are indeterminate, and a limit error should be signaled. Otherwise, the identities above apply.
The procedures can be modeled as follows:
sincheck=: 3 : 0 'x y'=. +. y assert. thmax>|x ((sin x)*cosh y) j. ( cos x)*sinh y ) coscheck=: 3 : 0 'x y'=. +. y assert. thmax>|x ((cos x)*cosh y) j. (-sin x)*sinh y ) sincheck 3j4 3.85374j_27.0168 sin 3j4 3.85374j_27.0168 sincheck 2e9j_1 |assertion failure: sincheck | thmax>|x sin 2e9j_1 |limit error: sin | sin 2e9j_1 sincheck 5j1000 __j_ sin 5j1000 __j_ sincheck 0j900 0j_ sin 0j900 0j_ coscheck 0j900 _ cos 0j900 _
Collected Definitions
sin =: 1&o. cos =: 2&o. sinh =: 5&o. cosh =: 6&o. hex =: (4+#@$) }. 2&(3!:3) unhex=: (3!:2) @ ((}: 2 (3!:3) 0.5)&,) " 1 e =: 2^_1074 NB. smallest positive finite number E =: (2-2^_52)*2^1023 NB. largest positive finite number thmax=: 2p1 * 2^53%2 NB. upper bound on arguments to sin and cos emin =: ^.e NB. ^y is 0 if y<emin; _1074*^.2 emax =: ^.E NB. ^y is _ if y>emax; 1024*^.2 emax2=: E +&^. 2 NB. cosh y is _ if y>emax2; 1025*^.2 expcheck=: 3 : 0 'x y'=. +. y if. x<emin do. 0j0 return. end. assert. thmax>|y if. x<:emax do. (^x) * (cos y) j. (sin y) else. ((*cos y)*^x+^.|cos y) j. ((*sin y)*^x+^.|sin y) end. ) sincheck=: 3 : 0 'x y'=. +. y assert. thmax>|x ((sin x)*cosh y) j. ( cos x)*sinh y ) coscheck=: 3 : 0 'x y'=. +. y assert. thmax>|x ((cos x)*cosh y) j. (-sin x)*sinh y )
See also
Contributed by Roger Hui.Categoty:Numeric Precision