Essays/Extremal Arguments

From J Wiki
Jump to navigation Jump to search

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