Essays/Extended Precision Functions

Some algebraic and transcendental functions are computed to an arbitrary number of decimal places. In each case, x f y computes f y to x decimal places, with a default of 40 for x . The computations exploit the rational number facility in J; therefore, intermediate calculations nearly always are performed on many more digits than x and are not as efficient as can be.

The argument y should be specified as either an extended integer or rational number, to avoid errors that can be introduced in the conversion to such when specified as an ordinary number (a machine-word integer or IEEE 64-bit floating point number) and errors inherent in the representation. See the examples under Logarithm and Sine and Cosine.

Utilities

DP=: 40

round=: DP&$: : (4 : 0)
 b %~ <.1r2+y*b=. 10x^x
)

pi=: DP&$: : (4 : 0)
 b %~ (x:!.0 y) * <.@o. b=. 10x^x+8+0>.>.10^.>./|y
)

DP specifies the default number of decimal places.  x round y rounds y to x decimal places.  x pi y computes π times y to x decimal places.

Logarithm

 ^.y satisfies the identities

(^.a*b) = (^.a) + ^.b
(^.a^b) = b * ^.a

y can be written as r * 2^m where m is an integer and r is between %%:2 and %:2 ,  whence ^.y can be computed as (^.r)+m*^.2 .

Equation 4.1.27 in Abramowitz & Stegun is a power series for ^.y :

(^.y) = 2 * +/ i %~ ((y-1)%y+1) ^ i=. 1+2*i.n

The series converge for all y , rapidly if y is between %%:2 and %:2 (between approximately 0.7 and 1.4).

The logic in the preceding discussion is implemented in the verb ln :

ln=: DP&$: : (4 : 0) " 0
 assert. 0<y
 m=. <.0.5+2^.y
 t=. (<:%>:) (x:!.0 y)%2x^m
 if. x<-:#":t do. t=. (1+x) round t end.
 ln2=. 2*+/1r3 (^%]) 1+2*i.>.0.5*(%3)^.0.5*0.1^x+>.10^.1>.m
 lnr=. 2*+/t   (^%]) 1+2*i.>.0.5*(|t)^.0.5*0.1^x
 lnr + m * ln2
)

Examples:

   0j40 ": ln ,.2 10 100 1000 0.1 0.01 0.001
 0.6931471805599453094172321214581765680755
 2.3025850929940456840179914546843642076011
 4.6051701859880913680359829093687284152022
 6.9077552789821370520539743640530926228033
_2.3025850929940456840179914546843642076011
_4.6051701859880913680359829093687284152022
_6.9077552789821370520539743640530926228033

ln can not deliver more precision than the argument possesses. For example:

   0j50 ": 50 ln 1x1
0.99999999999999999284228413360371724347638656933644

1x1 is Euler's number e represented as a 64-bit floating point number and is therefore a rational approximation to e. A more accurate representation for e leads to the expected answer:

   1e_50 (>i.1:) %!i.50
42
   %!42
7.11741e_52
   0j50 ": +/%!i.42x
2.71828182845904523536028747135266249775724709369996
   0j50 ": 50 ln +/%!i.42x
1.00000000000000000000000000000000000000000000000000

Exponential

By application of the identities

(^a+b) = (^a) * ^b
(^a-b) = (^a) % ^b
(^a*b) = (^a) ^ b

^y is seen to be equal to (2^m)*^y-m*^.2 . We specify m=.<.0.5+y%^.2 so that the magnitude of the difference d=.y-m*^.2 is bounded by 2%~^.2 (approximately 0.35). The power series +/(d^i)%!i=.i.n for ^d converges rapidly for such d .

The logic is the preceding discussion is implemented in the verb exp :

exp=: DP&$: : (4 : 0) " 0
 m=. <.0.5+y%^.2
 xm=. x+>.m*10^.2
 d=. (x:!.0 y)-m*xm ln 2
 if. xm<-:#":d do. d=. xm round d end.
 e=. 0.1^xm
 n=. e (>i.1:) a (^%!@]) i.>.a^.e [ a=. |y-m*^.2
 (2x^m) * 1++/*/\d%1+i.n
)

Examples:

   0j40 ": exp ,. _1 1 2 20
        0.3678794411714423215955237701614608674458
        2.7182818284590452353602874713526624977572
        7.3890560989306502272304274605750078131803
485165195.4097902779691068305415405586846389889448

   0j40 ": ln exp 20
20.0000000000000000000000000000000000000000
   0j40 ": exp ln 1r5
0.2000000000000000000000000000000000000000

   0j40 ": exp -: ln 2
1.4142135623730950488016887242096980785697

   0j40 ": exp pi sqrt 163
262537412640768743.9999999999992500725971981194575373308106

Square Root

sqrt=: DP&$: : (4 : 0) " 0
 assert. 0<:y
 %/ <.@%: (2 x: (2*x) round y)*10x^2*x+0>.>.10^.y
)

Examples:

   0j40 ": sqrt 2
1.4142135623730950488016887242096980785697

   0j40 ": ,. (, *:) sqrt 3 %~ 10^15x
       18257418.5835055371152323260933600711317581564999
333333333333333.3333333333333333333333333333333333333333

Sine and Cosine

Since sin y and cos y are periodic only arguments in the range 0 to 2*π need to be handled by the core routines. The arguments that need to be handled can be further restricted by application of the sum formulas

(sin a+b) = ((sin a)*cos b) + (sin b)*cos a
(cos a+b) = ((cos a)*cos b) - (sin a)*sin b

The following table shows how to compute sin y and cos y in various ranges.

range (multiples of π)		sin y	cos y
0	0.25	sin y	cos y
0.25	0.5	cos y-~π*0.5	sin y-~π*0.5
0.5	0.75	cos y- π*0.5	- sin y- π*0.5
0.75	1	sin y-~π	- cos y-~π
1	1.25	- sin y- π	- cos y- π
1.25	1.5	- cos y-~π*1.5	- sin y-~π*1.5
1.5	1.75	- cos y- π*1.5	sin y- π*1.5
1.75	2	- sin y-~π*2	cos y-~π*2

In all cases, sin y and cos y can be computed as sin t or cos t where t is between 0 and π*0.25 . (However, if y is 1 or less, it is better to just compute sin y and cos y directly rather than cos y-~π*0.5 and sin y-~π*0.5 as suggested by the table.)

The Taylor series for sin and cos are as follows. The series converge for all y , rapidly if |y is less than 1.

(sin y) = -/ (y ^ i) % ! i=. 1 + 2 * i. n
(cos y) = -/ (y ^ i) % ! i=.     2 * i. n

The logic in the preceding discussion is implemented in the following verbs:

cnt=: 4 : 0  NB. range index; # terms; standardized argument
 t=. ((x+0>.>.(10^.|y)-10^.2p1) pi 2)|x:!.0 y
 c=. (1,0.25p1*2+i.6) I. x:^:_1 t
 t=. (1+x) round +/(_1x^c+0 1)*t,(1r2*>.c%2)*(1+x) pi 1
 e=. 0.1^x
 n=. e (>i.1:) d (^%!@]) 2*i.x>.<.0.5*d^.e [ d=. x:^:_1 t
 c;n;t
)

sin=: DP&$: : (4 : 0) " 0
 'c n t'=. x cnt y
 (_1^c e. 4 5 6 7) * -/ t (^%!@]) (2x*i.n) + c e. 0 3 4 7
)

cos=: DP&$: : (4 : 0) " 0
 'c n t'=. x cnt y
 (_1^c e. 2 3 4 5) * -/ t (^%!@]) (2x*i.n) + c e. 1 2 5 6
)

Examples:

   0j50 ": 50 sin 1
0.84147098480789650665250232163029899962256306079837
   0j50 ": 50 cos 1
0.54030230586813971740093660744297660373231042061792
   0j50 ": 50 (sin +&*: cos) 1
1.00000000000000000000000000000000000000000000000000

   0j40 ": sin 0.8
0.7173560908995227616271746105813853661928

Regarding the last result, Abramowitz & Stegun says (Table 4.6, page 158) that sin 0.8 is 0.71735 60908 99522 76162 718 . The last digit in the A&S figure is incorrect.

The constant 1p1 is represented as a 64-bit floating point number and is therefore a rational number. Consequently, sin 1p1 is nonzero as sin y can not be zero for any nonzero rational number y .

   1p1 - 3.141592653589793
0

Since only a true 0 displays as 0, 1p1 is exactly 3.141592653589793 , a number slight less than π.

   y=: 3141592653589793x % 10^15x
   y
3141592653589793r1000000000000000
   0j40 ": sin y
0.0000000000000002384626433832795028841972
   0j40 ": cos y
_0.9999999999999999999999999999999715677839

Inverse Sine and Cosine

For the inverse sine, Abramowitz & Stegun 4.4.41 is used when 0.5<:|y and 4.4.40 when 0.5>|y . The following table is the estimated number of required terms for 40 decimal places for each y for 4.4.41 and 4.4.40:

   f=: 4 : '(<._0.5*x%10^.y),~<.-x%10^.(1-y)%2'
   0j2 6 0 ": (,. 40&f"0) (0.02*i.6),0.2 0.4 0.5 0.6 0.8 ,0.9+0.02*i.5
0.00   132    0
0.02   129   11
0.04   125   14
0.06   121   16
0.08   118   18
0.10   115   20
0.20   100   28
0.40    76   50
0.50    66   66
0.60    57   90
0.80    40  206
0.90    30  437
0.92    28  552
0.94    26  744
0.96    23 1128
0.98    20 2279

asin1=: 4 : 0   NB. A&S 4.4.41
 z=. 1-y
 k=. 1x + i.<.-x%10^.z%2
 s=. 1x + +/ (z^k) * (>:2*k) %~ */\ (<:2*k) % 4*k
 (x pi 1r2) - s * x sqrt 2*z
)

asin0=: 4 : 0   NB. A&S 4.4.40
 k=. 1x + 2 * i.<._0.5*x%10^.y
 +/ (y^k) * k %~ }: 1 , */\ k % 1+k
)

arcsin=: DP&$: : (4 : 0) " 0
 assert. 1>:|y
 y1=. | (1+x) round x:!.0 y
 if. 0.5<:|y do. (*y)*x asin1 y1 else. (*y)*x asin0 y1 end.
)

arccos=: DP&$: : (4 : 0) " 0
 (x pi 1r2) - x arcsin y
)

Examples:

   0j40 ": arcsin 1r2
0.5235987755982988730771072305465838140329
   0j40 ": sin arcsin 1r2
0.5000000000000000000000000000000000000000

Inverse Tangent

arctan=: DP&$: : (4 : 0) " 0
 if. 0=y do. 0 return. end.
 if. 1>:|y do. x arctan1 y else. x arctan0 y end.
)

arctan0=: DP&$: : (4 : 0) " 0    NB. A&S 4.4.42
 y1=. x:!.0 |y
 r =. %^:(1<|y)            y1
 r2=. %^:(1<|y) x round *: y1
 n=. >.-x*r2^.10
 a=. r * -/(1+2*i.n)%~*/\1,(n-1)$r2
 if. 1>|y do. (*y)*a else. (*y)*a-~(2*t)%~<.@o. t=. 10x^1+x end.
)

arctan1=: DP&$: : (4 : 0) " 0    NB. Euler
 y1=. x:!.0 y
 r=. x round (%>:) *: y1
 n=. >.-x*r^.10
 y1 %~ +/ */\ r * (1>.4x**:i.n) % */"1]1>.i.n,2
)

Examples:

   0j40 ": arctan 1
0.7853981633974483096156608458198757210493
   0j40 ": pi 1r4
0.7853981633974483096156608458198757210493

   0j40 ": arctan sqrt 3
1.0471975511965977461542144610931676280657
   0j40 ": pi 1r3
1.0471975511965977461542144610931676280657

   0j40 ": arctan 1000
1.5697963271282297525647978820048308980870

Hyperbolic Sine and Cosine

sinh=: DP&$: : (4 : 0) " 0
 if. x>10^.|y do. (1r2*_1^0>y)*x exp |y end.
 if. 1<|y do. -:-/x exp y,-y end.
 e=. 0.1^x
 n=. e (>i.1:) y (^%!@]) 1+2*i.>.0.5*y^.e
 +/(x:!.0 y) (^%!@]) 1+2x*i.n
)

cosh=: DP&$: : (4 : 0) " 0
 if. x>10^.|y do. (1r2*_1^0>y)*x exp |y end.
 if. 1<|y do. -:+/x exp y,-y end.
 e=. 0.1^x
 n=. e (>i.1:) y (^%!@]) 2*i.>.0.5*y^.e
 +/(x:!.0 y) (^%!@]) 2x*i.n
)

Examples:

   0j40 ": sinh 1r3
0.3395405572561501391012606113386035850724
   0j40 ": cosh 1r3
1.0560718678299393895268647082639832525255

Error Function

The Taylor's series for erf y is:

(erf y) = 2p_0.5 * -/ (!i.n) %~ y (^%]) 1+2*i.n

This can be implemented as follows, with (y^2*k)%!k being computed as */(k$*:y)%1+i.k to reduce the magnitude of intermediate results.

erf=: DP&$: : (4 : 0) " 0
 e=. 0.5p0.5*0.1^1+x
 y1=. |x:^:_1 y=. x: y
 if. 1>y1 do. m=. >.-:y1^.e else. m=. >.2p1**:y1 end.
 n=. (^.e) (>i.1:) (^.y1)+(i*2*^.y1)-(^.1+2*i)++/\^.1>.i=. i.m
 (2 % x sqrt x pi 1) * -/(1+2*i.n)%~*/\(y,(n-1)$*:y)%1>.i.n
)

Examples:

   0j40 ": erf 4r5
0.7421009647076604861671105865029458773177
   0j40 ": erf 6
0.9999999999999999784802632875010868834067

Normal CDF

The CDF of the N(0,1) distribution can be computed readily from the error function:

n01cdf=: DP&$: : (4 : 0) " 0
 2 %~ 1 + x erf (x: y) * x sqrt 1r2
)

Examples:

   0j40 ": n01cdf 1r2
0.6914624612740131036377046106083377398836

   0j40 ": ,. n01cdf 1 2 3
0.8413447460685429485852325456320379224779
0.9772499023342778722461244926196457873979
0.9986501019683874744324597336961989547832

Collected Definitions

DP=: 40

round=: DP&$: : (4 : 0)
 b %~ <.1r2+y*b=. 10x^x
)

pi=: DP&$: : (4 : 0)
 b %~ (x:!.0 y) * <.@o. b=. 10x^x+8+0>.>.10^.>./|y
)

ln=: DP&$: : (4 : 0) " 0
 assert. 0<y
 m=. <.0.5+2^.y
 t=. (<:%>:) (x:!.0 y)%2x^m
 if. x<-:#":t do. t=. (1+x) round t end.
 ln2=. 2*+/1r3 (^%]) 1+2*i.>.0.5*(%3)^.0.5*0.1^x+>.10^.1>.m
 lnr=. 2*+/t   (^%]) 1+2*i.>.0.5*(|t)^.0.5*0.1^x
 lnr + m * ln2
)

exp=: DP&$: : (4 : 0) " 0
 m=. <.0.5+y%^.2
 xm=. x+>.m*10^.2
 d=. (x:!.0 y)-m*xm ln 2
 if. xm<-:#":d do. d=. xm round d end.
 e=. 0.1^xm
 n=. e (>i.1:) a (^%!@]) i.>.a^.e [ a=. |y-m*^.2
 (2x^m) * 1++/*/\d%1+i.n
)

sqrt=: DP&$: : (4 : 0) " 0
 assert. 0<:y
 %/ <.@%: (2 x: (2*x) round y)*10x^2*x+0>.>.10^.y
)

cnt=: 4 : 0  NB. range index; # terms; standardized argument
 t=. ((x+0>.>.(10^.|y)-10^.2p1) pi 2)|x:!.0 y
 c=. (1,0.25p1*2+i.6) I. x:^:_1 t
 t=. (1+x) round +/(_1x^c+0 1)*t,(1r2*>.c%2)*(1+x) pi 1
 e=. 0.1^x
 n=. e (>i.1:) d (^%!@]) 2*i.x>.<.0.5*d^.e [ d=. x:^:_1 t
 c;n;t
)

sin=: DP&$: : (4 : 0) " 0
 'c n t'=. x cnt y
 (_1^c e. 4 5 6 7) * -/ t (^%!@]) (2x*i.n) + c e. 0 3 4 7
)

cos=: DP&$: : (4 : 0) " 0
 'c n t'=. x cnt y
 (_1^c e. 2 3 4 5) * -/ t (^%!@]) (2x*i.n) + c e. 1 2 5 6
)

asin1=: 4 : 0   NB. A&S 4.4.41
 z=. 1-y
 k=. 1x + i.<.-x%10^.z%2
 s=. 1x + +/ (z^k) * (>:2*k) %~ */\ (<:2*k) % 4*k
 (x pi 1r2) - s * x sqrt 2*z
)

asin0=: 4 : 0   NB. A&S 4.4.40
 k=. 1x + 2 * i.<._0.5*x%10^.y
 +/ (y^k) * k %~ }: 1 , */\ k % 1+k
)

arcsin=: DP&$: : (4 : 0) " 0
 assert. 1>:|y
 y1=. | (1+x) round x:!.0 y
 if. 0.5<:|y do. (*y)*x asin1 y1 else. (*y)*x asin0 y1 end.
)

arccos=: DP&$: : (4 : 0) " 0
 (x pi 1r2) - x arcsin y
)

arctan=: DP&$: : (4 : 0) " 0
 if. 0=y do. 0 return. end.
 if. 1>:|y do. x arctan1 y else. x arctan0 y end.
)

arctan0=: DP&$: : (4 : 0) " 0    NB. A&S 4.4.42
 y1=. x:!.0 |y
 r =. %^:(1<|y)            y1
 r2=. %^:(1<|y) x round *: y1
 n=. >.-x*r2^.10
 a=. r * -/(1+2*i.n)%~*/\1,(n-1)$r2
 if. 1>|y do. (*y)*a else. (*y)*a-~(2*t)%~<.@o. t=. 10x^1+x end.
)

arctan1=: DP&$: : (4 : 0) " 0    NB. Euler
 y1=. x:!.0 y
 r=. x round (%>:) *: y1
 n=. >.-x*r^.10
 y1 %~ +/ */\ r * (1>.4x**:i.n) % */"1]1>.i.n,2
)

sinh=: DP&$: : (4 : 0) " 0
 if. x>10^.|y do. (1r2*_1^0>y)*x exp |y end.
 if. 1<|y do. -:-/x exp y,-y end.
 e=. 0.1^x
 n=. e (>i.1:) y (^%!@]) 1+2*i.>.0.5*y^.e
 +/(x:!.0 y) (^%!@]) 1+2x*i.n
)

cosh=: DP&$: : (4 : 0) " 0
 if. x>10^.|y do. (1r2*_1^0>y)*x exp |y end.
 if. 1<|y do. -:+/x exp y,-y end.
 e=. 0.1^x
 n=. e (>i.1:) y (^%!@]) 2*i.>.0.5*y^.e
 +/(x:!.0 y) (^%!@]) 2x*i.n
)

erf=: DP&$: : (4 : 0) " 0
 e=. 0.5p0.5*0.1^1+x
 y1=. |x:^:_1 y=. x: y
 if. 1>y1 do. m=. >.-:y1^.e else. m=. >.2p1**:y1 end.
 n=. (^.e) (>i.1:) (^.y1)+(i*2*^.y1)-(^.1+2*i)++/\^.1>.i=. i.m
 (2 % x sqrt x pi 1) * -/(1+2*i.n)%~*/\(y,(n-1)$*:y)%1>.i.n
)

n01cdf=: DP&$: : (4 : 0) " 0
 2 %~ 1 + x erf (x: y) * x sqrt 1r2
)

---

See also

• Chi Squared CDF
• Normal CDF
• Pi (Chudnovsky Algorithm)
• Sine
• Square Root
• t-Distribution CDF
• Extended Precision Functions

---

Contributed by Roger Hui.