J Rational Number Utilities & Verification of Knuth’s Constants
Vector 17:4, page 131
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 131
Overview¶
The utilities presented produce higher precision results than normal floating point. It is hoped that they, and their use to verify Knuth’s constants (Ref 1), are of interest.
The utility routines are applicable to other calculations requiring up to about 50 decimal places. Given further effort, the utilities could be presented as a high precision package.
Practical Calculations with Rational Numbers¶
J provides rational arithmetic giving exact rational results. However, long calculations can lead to lengthy rational number representations which may require extensive storage and are slow to process. In these circumstances the theoretical infinite precision of rational numbers can be sacrificed by rounding values to an appropriate precision. In any case, any irrational values required, have to be rounded (approximated) in order to be held as finite length rational numbers.
For example, the repeated squaring of pi (an 80 decimal place approximation to the well known mathematical constant), shows that the length of the rational representation of the results just about doubles at each squaring:
0j5": drl *: ^: (i.5) pi
drl: 79 79
drl: 157 157
drl: 314 313
drl: 628 625
drl: 1256 1249
3.14159 9.86960 97.40909 9488.53102 90032220.84293
Where the drl utility displays the numerator and denominator lengths of the rational number results. Continued squaring would eventually fill the available memory or be so slow that one is not prepared to wait for the result. Even using the simple rational, 22r7, as an approximation to pi, creates rational values with thousands of digits after squaring 10 times.
The utilities presented round [intermediate] results to keep memory requirements and execution times within practical limits to achieve high, but not infinite, precision. The programmer and user then have the added responsibility of ensuring that the rounding is appropriate.
The Utilities¶
Maths Rational Utilities¶
The Maths utilities mimic similar built-in J functions, but with a restricted argument range:
- No complex arguments or results (applies to all the utility functions)
- Limitation of arguments to the range of J floating values (say ±1.6e308)
- Result values may be restricted to those near unity
| Name | Description | J built-in function mimicked | Algorithm Type | Approximate Arg. Limits at 45dp | Arg. Type x | y |
|---|---|---|---|---|---|---|
ca |
Cosine | monad 2&o. |
as sa | as sa | [N] | I, Q, X |
ea |
Exponential | monad ^ |
Series | 70x | [N] | I, Q, X |
fa |
Factorial | monad ! |
Series | _1< to <=6 | [N] | I, Q, X |
ge |
Euler’s constant, ( | (none) | Series | not applicable | [N] | I, Q, X |
pa |
Pi | (o. 1) |
J (o.) |
not applicable | [N] | I, Q, X |
la |
Natural logarithm | monad ^. |
Series | 10^300x | [N] | I, Q, X |
ra |
(Square) Root | monad %: |
J (%:) |
10^80x | [N] | I, Q, X |
sa |
Sine | monad 1&o. |
Series | +/- pi*10^3x | [N] | I, Q, X |
pa and ra make use of J built-in o. and %: monadic functions which produce extended integer results as a special case.
Argument types are: N any numeric; I integer; Q rational and X extended integer. In all cases the result is a rational value. Where the optional x argument is present, it is an integer specifying the number of decimal places of the result of the function (the default is the value in global variable T). The right argument may not be a floating or complex type numeric.
Basic Rational Utilities¶
| Name | Valence | Description |
|---|---|---|
BN |
Monad | Bernoulli numbers (requires extended integer argument for rational result) |
cdp |
Dyad | Number of common decimal places in two rational numbers. cdp is |
Cdp |
Dyad | limited to argument differences of magnitude less than about 5e17 |
drl |
Monad | Display rational length (numerator and denominator) and return rational |
flg |
Monad | Floor of base 10 logarithm of the rational argument (currently unused) |
l2r |
Monad | Literal representation (of an integer, decimal or floating value) to rational |
rr |
Ambivalent | Round rational. Reduces the length of a rational value where it can. No attempt is made to find the best approximation possible. |
Notes:
cdp,Cdp, anddrlare mainly for diagnosticsl2renables rational numbers to be written as literal floating or decimal numbers
Constants for the Rational Utilities¶
These are rational constants required by the utilities, but some may be useful to the J programmer directly. They are given to 80 decimal places including 5 guard digits.
| Name(s) | Description |
|---|---|
e1 |
Base of natural logarithms, e |
ln2, ln10 |
Natural logarithms of 2 10 |
lnu, lnx, lny, lnz |
Natural logarithms of 33r32 1025r1024 32769r32768 1048577r1048576 |
pi |
Pi |
The constants were produced using the J script itself, with the global control variable ctr set to zero. Thus, the script is not dependant on numerical values from any earlier work except in so far as such values are included within J itself.
Calculation of Knuth’s Constants¶
Run Constant Table script¶
The call RCT KC produces a listing of most of the constants in Table 1, Appendix 1 of Knuth’s book ( Ref 1). It might be expected that these could be extended to more and more decimal places. However, 50 decimal places is about the limit with the current utilities.
The 31 constants are calculated in about 5 seconds on an IBM compatible PC. The time is the time to run RCT KC on a system with an 233MHz Pentium II processor and 64 MB RAM. The time does not include the precalculated constants in the source file accessed by RCT, and their calculation times would need to be taken into account for an objective timing.
Run Constant Table script Diagnostic¶
The RCTD function the displays number of common decimal places between pairs of results in order of increasing precision. It can be a useful diagnostic. A typical call is:
Knuth’s Constants Table Data¶
The KC table contains three fields for each of Knuth’s constants (Ref 1):
1) Line No.; 2) Short description of constant; 3) J script to calculate the constant
The relation (x^y) = ea y * la x is used in 3) to calculate the power of a number in terms of the available rational utilities, ea and la.
Alternative J source has been used to calculate the constants in some cases. As an example, consider an alternative J script for pi:
When compared with Line 12 in the calculated constants, the difference is _266e_45 within the 5 guard digits, giving increased confidence in the utility results.
Reference¶
- Knuth, Donald E: The Art of Computer Programming Vol 1, First edition 1968; Addison-Wesley.
The J Source¶
The J source code is provided ‘as is’, and consists of about 300 lines (4 A4 pages of J10pt).
If problems are encountered, a reduction in the number of decimal places will often enable a result to be obtained. This may be achieved by reducing the value in:
- the global variable
T, the default precision, or - the left arguments of a utility, where a utility precision is set explicitly
Verification of Knuth’s Constants (Ref 1)¶
The constants which have been recalculated are in Table 1, Appendix 1 of Knuth’s book ‘The Art of Computer Programming’. Of the 32 constants, one, zeta, was omitted due to problems finding a suitable algorithm. The remaining 31 results agree with Knuth’s values (40 decimal places).
Although the results agree with Knuth’s excellent work, the exercise was a very useful debugging exercise for the Rational Number utilities.
RCT KC
Nominal calculation precision, T, 45 decimal places. Results correct to 40 decimal places
LN Constant Time Memory J Code D.Places: 10 20 30 40
0) Sqrt 2 0.007 50240 ra 2 1.41421 35623 73095 04880 16887 24209 69807 85696 71875
1) Sqrt 3 0.007 63744 ra 3 1.73205 08075 68877 29352 74463 41505 87236 69428 05253
2) Sqrt 5 0.006 42688 ra 5 2.23606 79774 99789 69640 91736 68731 27623 54406 18359
3) Sqrt 10 0.005 39104 ra 10 3.16227 76601 68379 33199 88935 44432 71853 37195 55139
4) Cube root 2 0.284 184256 ea (%&3) la 2 1.25992 10498 94873 16476 72106 07278 22835 05702 51466
5) Cube root 3 0.360 153664 ea (%&3) la 3 1.44224 95703 07408 38232 16383 10780 10958 83918 69253
6) Fourth root 2 0.013 50368 ra ra 2 1.18920 71150 02721 06671 74999 70560 47591 52929 72092
7) ln 2 0.020 79680 la 2 0.69314 71805 59945 30941 72321 21458 17656 80755 00134
8) ln 3 0.129 153088 la 3 1.09861 22886 68109 69139 52452 36922 52570 46474 90557
9) ln 10 0.148 140608 la 10 2.30258 50929 94045 68401 79914 54684 36420 76011 01490
10) 1/ln 2 0.023 79680 % la 2 1.44269 50408 88963 40735 99246 81001 89213 74266 45955
11) 1/ln 10 0.149 140608 % la 10 0.43429 44819 03251 82765 11289 18916 60508 22943 97006
12) pi 0.010 50048 pa 0 3.14159 26535 89793 23846 26433 83279 50288 41971 69399
13) pi/180 0.005 67712 pi%180 0.01745 32925 19943 29576 92369 07684 88612 71344 28719
14) 1/pi 0.005 64064 % pi 0.31830 98861 83790 67153 77675 26745 02872 40689 19291
15) pi squared 0.010 167552 *: pi 9.86960 44010 89358 61883 44909 99876 15113 53136 99407
16) Sqrt pi, Gamma 1/2 0.007 43648 ra pi 1.77245 38509 05516 02729 81674 83341 14518 27975 49456
17) Gamma 1/3 1.287 171712 G1r3=.fa _2r3 2.67893 85347 07747 63365 56929 40974 67764 41286 89529
18) Gamma 2/3 0.026 147200 +:pi%G1r3*ra 3 1.35411 79394 26400 41694 52880 28154 51378 55193 27190
19) e 0.019 82240 ea 1 2.71828 18284 59045 23536 02874 71352 66249 77572 47094
20) 1/e 0.018 82368 ea _1 0.36787 94411 71442 32159 55237 70161 46086 74458 11131
21) e squared 0.029 130176 ea 2 7.38905 60989 30650 22723 04274 60575 00781 31803 15571
22) Euler's const gamma 0.687 170112 EC=.ge 1e2 0.57721 56649 01532 86060 65120 90082 40243 10421 59334
23) ln pi 0.170 156032 la pi 1.14472 98858 49400 17414 34273 51353 05871 16472 94812
24) phi, Golden Ratio 0.012 63744 -:>: ra 5 1.61803 39887 49894 84820 45868 34365 63811 77203 09180
25) e to power gamma 0.242 194816 ea EC 1.78107 24179 90197 98523 65041 03107 17954 91696 45209
26) e to power pi/4 0.249 174080 ea (%&4) pa 0 2.19328 00507 38015 45655 97696 59278 73822 34616 37641
27) sin 1 0.466 161280 s1=.sa 1 0.84147 09848 07896 50665 25023 21630 29899 96225 63063
28) cos 1 0.018 83136 ra 1-*:s1 0.54030 23058 68139 71740 09366 07442 97660 37323 10417
29) zeta 3, Sigma 1/n^3 0.000 2112 Not calculated
30) ln phi 0.259 148672 la -:>: ra 5 0.48121 18250 59603 44749 77589 13424 36842 31351 84335
Appendix: Listing of J script¶
NB. File: ratnum9.ijs Rational Utilities (J Ver 405d) Wed 15/11/2000
ctr =: 1 NB. 0: Initialise constants; 1: Normal
diag =: 0 NB. 0: Normal; 1: Display loop counts in ea etc
dsp =: 1!:2&2 NB. Display y. and result y.
T =: 45 NB. Default rounding decimal places
NB. Utilities for Rational Numbers
BN =: 3 : 0 NB. Bernoulli Numbers. Argument is an extended integer
NB. SYNTAX: BN rat_no_of_Bernoulli_nos_to_create (starting at BN0)
B=.'' NB. Initial results vector
d=.1,(<:y.)#0 NB. Initial denominator coefficients
r=.%<.@!>:i.y. NB. Initial divisor coefficients
for. i.y. do. NB. Loop for algebraic division
B=.B, q=. d %&{. r NB. Quotient is scalar q
d=.}.d-q*r NB. Update partial dividend (drop initial 0)
r=.}:r NB. Drop first/last items of d & r
end.
B * <.@!i.y. NB. Bernoulli Numbers
)
NB. Bernoulli Numbers for k=0 1 2 (2) 50 ie. omitting bn having zero values
bnnz =: (~:&0 # ]) BN 51x NB. Bernoulli Nos. Non-Zero
drl =: 3 : 0 NB. Display Rational (and Extended) number component Lengths
NB. SYNTAX: drl scalar_rational_no
([ [: dsp@('drl: '&,)@":'#;._2 @(=&'_'@{. }. ]) @,&'r'@":)"0 y.NB. Display N & D lengths
NB. y. NB. Pass y. on as result
)
flg =: (3 : 0)"0 NB. Floor of (base 10) logarithm of rational
NB. SYNTAX: flg rat_no
if. y.>:1 do. _1 + # 0": <. y. NB. Log base 10 of y.
elseif. y.>0 do. (-@# + =&'1'@{. *. */@:=&'0'@}.) 0": >. % y. NB. Log base 10 of y.
elseif. y.=0 do. x: __ return. NB. Case for zero argument
elseif. do. 0j1 [ dsp 'flg: ERROR: negative argument' return.
end.
)
l2r =: 3 : 0 NB. Literal (J format float or integer) value to Rational number
NB. SYNTAX: l2rx literal_value NB. y.= _. is not handled
'A E'=. 2{. <;._2 y.,'e0e' NB. Arg and Exp of equivalent float
'i f'=. 2{. #;._2(>A),'..' NB. Integer and Fraction lengths
". (([ #~ '.'&~:)>A),((0>.".E )#'0'),'r1',(f-0<.".E)#'0' NB. Rational result
)
NB. x.>=0 Return 1 if two rationals in y. differ within specified dp, otherwise 0
rdf =: [: -.@(*/)@:=/ (j.@[ ": ])"0 NB. Rat.Diff; x: dps_acc; y: two rational nos
rr =: (3 : 0)"0 NB. Round Rational
T rr y. NB. Monad definition
:
NB. SYNTAX: [dec_or_sig_places] rr rat_no
if. 0=y. do. y. return. end. NB. Return rat 0 if zero
if. 0<:x. do. l2r (j.x.) ": y. NB. Round to x. decimal places (x.<=0)
else. l2r (j.0>.-x.+>:flg|y.)": y. NB. Round to -x. sig figs
end. NB. %&ln10 60 la is better than 10^.
)
NB. Constants: lnu to lnz, e 1, ln 2, ln 10 and pi to 80 dps (last place(s) are in doubt)
NB. (ea lnu) to (ea lnz) have exact rational representations of 33r32 to 1048577r1048576
lnu =:l2r'0.03077165866675368837102820759677216409169673995889035634986199530109980728337209'
lnx =:l2r'0.00097608597305545889596082490801718667261183433378453623775859827440037212402588'
lny =:l2r'0.00003051711247318637856906951416899468124304989470989364165971885272442917955964'
lnz =:l2r'0.00000095367386165918823390841551496333614360314807097930285811225469343381836078'
e1 =:l2r'2.71828182845904523536028747135266249775724709369995957496696762772407663035354759'
ln2 =:l2r'0.69314718055994530941723212145817656807550013436025525412068000949339362196969472'
ln10=:l2r'2.30258509299404568401799145468436420760110148862877297603332790096757260967735248'
pi =:l2r'3.14159265358979323846264338327950288419716939937510582097494459230781640628620900'
cdp =: -@(10&^.)@(x:^:_1&.%)@+:@|@- NB. Common dps in pair of rational numbers
Cdp =: -@%&ln10@(50&la)@+:@|@- NB. Common dps in pair of rational numbers
NB. Operations with Rational numbers
ca =: (3 : 0)"0
T ca y. NB. Monad definition
:
x.&sa y.+-:pi NB. Cosine; ? accuracy for small |y.
)
ea =: (3 : 0)"0 NB. Rational Exponential Approximation
T ea y. NB. Monad definition
:
NB. SYNTAX: [dps_accuracy] ea exp_rat_arg
y.=.(x:y.)-m=.ctr*<.0.5+y. NB. Scale y. m times by 1
y.=.(x:y.)-lnu* mu=.ctr*<.0.5+y.%lnu NB. Scale y. mu times by 1r32 approx
y.=.(x:y.)-lnx* mx=.ctr*<.0.5+y.%(x:^:_1)lnx NB. Scale y. mx times by c. 1r1024
y.=.(x:y.)-lny* my=.ctr*<.0.5+y.%(x:^:_1)lny NB. Scale y. my times by c. 1r32768
NB. y.=.(x:y.)-lnz* mz=.ctr*<.0.5+y.%(x:^:_1)lnzNB. Scale y. mz times by c. 1r1048576
NB. dsp m, mu, mx, my, mz
PS =. 0r1 NB. Initalise accumulator, Previous Sum
r =. 0 NB. Initial r of series
t =. 1 NB. Initial term of series
while. x. rdf PS, R=.PS + t do. NB. Test Previous Sum against Current Sum
r =. >:r NB. Increment term number
t =. x. rr t*y.%r NB. Update t to next term
PS=. x. rr R NB. Partial result
end.
if. diag do. dsp 'ea terms: ', ":r end. NB. Diagnostic
R =. R * 33r32^mu NB. Reverse scaling 33r32
R =. R * 1025r1024^mx NB. Reverse scaling 1025r1024
R =. R * 32769r32768^my NB. Reverse scaling 32769r32768
NB. R =. R * 1048577r1048576^mz NB. Reverse scaling 1048577r1048576
fi=: x.&rr@*&e1 :.(x.&rr@%&e1) NB. Define obverse function
x. rr fi ^: m R NB. Reverse scaling of result
)
la =: (3 : 0)"0 NB. Rational Natural logarithm Approximation
T la y. NB. Monad definition
:
NB. SYNTAX: [dps_accuracy] la ln_rat_arg
if. diag *. 0=y. do. dsp 'la terms: 0' end.
if. 0=y. do. __ return. end.
y. =. y.* 2x^ - m=.ctr*<.0.5+-/%&ln2^.2 x: y. NB. Scale y. m times by 2
y. =. y.* 33r32^ - mu=.ctr*<.0.5+-/%&(x:^:_1 lnu)^.2 x: y.NB. Scale y. mu times by 1.01
y. =. y.* 1025r1024^ - mx=.ctr*<.0.5+-/%&(x:^:_1 lnx)^.2 x: y.NB. Scale y. mx times by 1.01
y. =. x. rr y.* 32769r32768^ - my=.ctr*<.0.5+-/%&(x:^:_1 lny)^.2 x: y.NB. Scale y. my times b
y 1.01
NB.dsp mu,mx,my,mz
x=.(<: % >:)y. NB. x =. (y.-1)/(y.+1)
R =. 0 NB. Initial accumulator
r =. 1 NB. Initial term no.
i =. x NB.
x2=. x. rr *:x NB.
while. x. rdf R , R1 =. x. rr R + i % r do.
i =. x.rr i*x2 NB.
r =. r+2 NB.
R =. R1 NB.
end.
if. diag do. dsp 'la terms: ', ": -:r+1 end. NB. Diagnostic
NB. R=. x. rr (mh*lnh) + +:R
R=. (mu*lnu) + +:R
R=. (mx*lnx) + R
R=. (my*lny) + R
NB. R=. (mz*lnz) + R
x. rr (m*ln2) + R NB. Reverse scaling
)
ra =: (3 : 0)"0 NB. Rational Approximation to square_root (Newton's Method)
T ra y. NB. Monad definition
:
NB. SYNTAX: [dps_accuracy] ra Root_arg
(10^x:x.) %~ <.@%: y.*10^x:+:x. NB. Root as Roger Hui's pi
)
sa =: (3 : 0)"0 NB. Rational Sine Approximation
T sa y. NB. Monad definition
:
NB. SYNTAX: [dps_accuracy] sa sin_arg
a =. -&(+:pi) ^: (m=.<.0.5+ -:pi%~ y.) y. NB. Scale y. m times by 2pi
ps =. *a NB. Post calculation result multiplier
a =. | a NB. Modulus of a
if. a>-:pi do. a=.pi-a end. NB. Supplementary angle test
if. s=.a>pi%4 do.a=.a--:pi end.NB. Set s for sin: 0 or cos: 1
":m, s, a =. x. rr a NB. Truncate a, diagnostic
t =. s{a, 1x NB. Initial term: sin: angle; cos: 1
i =. - *: a NB. - angle squared
R =. 0 * a NB. Clear accumulator
r=.-.s NB. Initial term no., r, cos: 0; sin: 1
while. x. rdf R , R1 =. x. rr R+t do.
t =. x. rr t*i % (<:r)*r=.r+2 NB. Update r and term
R =. R1 NB. Updte R (? should truncate here also)
if. 2>10|r do. R=.x. rr R [ t=.x. rr t end.
end.
if. diag do. dsp 'sa terms: ', ":<.-:r+1 end. NB. Diagnostic
ps*R NB. Correct result sign if necessary
)
NB. Calculation of pi
pa =: (3 : 0)"0 NB. Rational Pi Approximation
T pa y. NB. Monad definition
:
NB. SYNTAX dps_accuracy pa ''
n %~ <.@o. n=.10^x: x. NB. Based on Roger Hui's Code
)
NB. Calculation of Gamma constant (Approximately 0.577)
ge =: (3 : 0)"0 NB. Knuth's algorithm for Euler's gamma
T ge y. NB. Monad definition
:
NB. SYNTAX: dps_accuracy ga n
n=. x: y. NB. Use n=1000 for upto 80 decimal places
R=.0r1 NB. Accumulator for sum of reciprocals
c =. fc =. 5 NB. Full count for inner loop
'nc lc'=.(>.@%&c ; c&|&.<:) n NB. Number_of_Cycles, Last_Count
for_j. i.nc do. NB. Outer loop
if. j=nc-1 do. c=.lc end. NB. Adjust count for last cycle
R=.x.rr R+x.&rr@+/%(1+j*fc)+i.c NB. Inner loop
end.
k =: +:>:i.15 NB. Bernoulli numbers for positive even k
x. rr R - +/(x. la n), % (2*n), n (-@%&(15{.2 }.bnnz)@] * ^) k
)
NB. Knuth's Constants, Constant Table
KC=: ,: '0) Sqrt 2'; 'ra 2'
KC=:KC, '1) Sqrt 3'; 'ra 3'
KC=:KC, '2) Sqrt 5'; 'ra 5'
KC=:KC, '3) Sqrt 10'; 'ra 10'
KC=:KC, '4) Cube root 2'; 'ea (%&3) la 2'
KC=:KC, '5) Cube root 3'; 'ea (%&3) la 3'
KC=:KC, '6) Fourth root 2'; 'ra ra 2'
KC=:KC, '7) ln 2'; 'la 2'
KC=:KC, '8) ln 3'; 'la 3'
KC=:KC, '9) ln 10'; 'la 10'
KC=:KC, '10) 1/ln 2'; '% la 2'
KC=:KC, '11) 1/ln 10'; '% la 10'
KC=:KC, '12) pi'; 'pa 0'
KC=:KC, '13) pi/180'; 'pi%180 '
KC=:KC, '14) 1/pi'; '% pi '
KC=:KC, '15) pi squared'; '*: pi '
KC=:KC, '16) Sqrt pi, Gamma 1/2'; 'ra pi '
KC=:KC, '17) Gamma 1/3'; 'G1r3=.fa _2r3'
KC=:KC, '18) Gamma 2/3'; '+:pi%G1r3*ra 3'
KC=:KC, '19) e'; 'ea 1'
KC=:KC, '20) 1/e'; 'ea _1'
KC=:KC, '21) e squared'; 'ea 2'
KC=:KC, '22) Euler''s const gamma'; 'EC=.ge 1e2'
KC=:KC, '23) ln pi'; 'la pi '
KC=:KC, '24) phi, Golden Ratio'; '-:>: ra 5'
KC=:KC, '25) e to power gamma'; 'ea EC'
KC=:KC, '26) e to power pi/4'; 'ea (%&4) pa 0'
KC=:KC, '27) sin 1'; 's1=.sa 1'
KC=:KC, '28) cos 1'; 'ra 1-*:s1'
KC=:KC, '29) zeta 3, Sigma 1/n^3'; '''Not calculated'''
KC=:KC, '30) ln phi'; 'la -:>: ra 5'
KC=:KC, '31) 1/ln phi'; '% la -:>: ra 5'
s5 =: 3 : 0
if. '.'e.y. do. (,&'.'@>@{. , ,@(_5&((,&' ')\) )@>@{:) @(<;._2)@,&'.' y.
else. y.
end.
)
RCT =: 3 : 0 NB. Run Constants Table, eg KC
T RCT y. NB. Monad definition
:
NB. SYNTAX: RCT KC
t=. 'Nominal calculation precision, T, ', (": T=:x.), ' decimal places. '
dsp t, 'Results correct to ', (":T-5), ' decimal places', LF
V =: (10#' '), 'D.Places:', 7}. , (10&#@(' '"_) , [ , '0'"_)"0 (<.T%10){.'123456789'
dsp 'LN Constant Time Memory J Code', V
for_i. _2[\,y. do. NB. Allow single item of shape 2
t=. 6!:2 's=. 7!:2 ''r=. ".c=.>{:i''' NB. Time & memory to execute J code
if. 2~:3!:0 r do. v=.s5 (3+j.~T)": r else. c=.'' [ v=.r end.
dsp (23{.>{.i), (6j3":t), (7":s), ' ', (16{.c), v
end.
0 0$''
)
RCTD =: 3 : 0 NB. Guide to no. of decimal places correct in result
(5*>:i.10) RCTD y. NB. Monad definition; 34 sec runtime
:
NB. SYNTAX: T_vals RCTD CT NB. Run Constant Table Diagnostic
dsp 'Guesstimate of number of decimal places correct (has flaws)', LF
R =: (0, #y.=._2[\,y.)$'' NB. Allow single item of shape 2
for_T. x. do.
T =: T NB. Set global T
R1 =. 0$''
for_i. y. do.
t=. 6!:2 's=. 7!:2 ''r=.".c=.>1{i''' NB. Time & memory to execute J code
if. 2~:3!:0 r do. v=.s5 ((5+T)j.T)": r else. c=.'' [ v=.r end.
$ R1 =. R1, <r
end.
R =: R, R1
end.
dsp 'LN Constant T=', 3}. 5":x.
(>{."1 y.),. 5j1": ,. (2& (cdp/\)"1@:> :: (0:/)"1) |:R
)
tc=: {.&bnnz *&}.&}. $&1 _1 % i. * <:@i. NB. Coeffs for Gamma (tc N)
ts=: ($&1@-&2@[ # |@tc@[) % [: */\ ] , $&1 @-&2@[ }.@# *:@] NB. Exp. terms (N ts n)
fa =: (3 : 0)"0 NB. Factorial y. approximation (Knuth p115)
T fa y. NB. Monad definition
:
NT=.51x NB. Arbitrary number of term in series
'y. z' =. (>:y.),(<:# y.)}.45x NB. Initialise z
N=.x:(<:@+:@(+/) % =&0@{:) x. >- 10^. NT ts x:^:_1 z NB. Determine No of terms +2 ???
n=. y+(z-<.y=.>:y.) NB. 25 <= n < 26. Term 41 <= 10^.51
NB. m is ln of correction divisor, +/fract+integers
m=. -(1&| x.&rr@ la@|@(_30&(x.&rr@*/\))@:+ <.@<: }. i.@>:@(z"_))y
k=.+:>:i.<.-:N+2 NB. Values of k for each term
C=. ({&bnnz@>:@-: % [ * <:)k NB. Modify bn coeffs
R=.x.&rr@+/ m, C % */\ n,(<.-:N-1)#*:n NB. Divide coeffs by powers of n and sum
x.&rr(ra +: x.&rr pi*n) * x.&rr ea R + n*la x.&rr n%e1 NB. Root and exponential terms
)