# 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 is 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
)
```