resolvent

Functions, if you are referring to monads and dyads, can be defined using "3 : '...'" or "3 : 0\n...\n)" using "3" for monads and "4" for dyads.

As for rank information, functions do carry that rank. I believe you are not defining any rank, it which case it becomes infinite rank and depending on usage, it will be applied incorrectly.

For tacit recursion, there is the "$:" operator, which allows a tacit function to call itself anonymously. You may also need the agenda operator "@." to define your base case with a gerund "`".

The "foreigns" table, while useful, seems like a cludgy way to introduce functions that don't or are difficult to introduce cohesively into J's notation.

The LCG is the improved extension which causes you to have to compute values x_0 = seed to x_n in order to calculate x_{n+1}. What I am talking about is just a mapping of index x -> f(x) which is what the OP seemed to have wanted. In that case, `a` only needs to be relatively prime. This cannot account for every permutation of n since the number of values relatively prime to n, the totient, has an upper bound of n, which is the possible values for a. The number of values for b is also n, so at most n^2 possible sequences are generated, which is less than n! for n > 3.

For example, with n = 5: Let a = 3 and b = 2. x = [0, 1, 2, 3, 4], a * x = [0, 3, 6, 9, 12], a * x + b = [2, 5, 8, 11, 14], a * x + b mod n = [2, 0, 3, 1, 4]

So you want a mapping of [1, n] to a permutation of [1, n] without having to generate the list [1, n] and shuffling it. An affine transformation modular n should work. Instead of [1, n], look at [0, n). Find a value `a` in [0, n) such that gcd(a, n) = 1 and pick a random integer b in [0, n). Then the `random` number at each position is `a * x + b mod n`.

This is simple and fast, but is not secure at all. You can solve for a, b by solving the linear congruence. It also does not generate every permutation of `n`. For n = 5, only 20 sequences can be found out of 5! = 120.

The difference lies in the definition of `EllipticK` which uses quantity `k` whereas the paper uses `k^2`. Also, symbolically integrating in Mathematica is much slower. Instead use `NIntegrate` to skip straight to a numerical value.

Print[ NIntegrate[ 1/(1 - Cos[x] Cos[y] Cos[z]), {x, 0, Pi}, {y, 0, Pi}, {z, 0, Pi}, Method -> "LocalAdaptive"] / Pi^3 // Timing]

{0.121543, 1.3932012917627028}

Print[ NIntegrate[ Pi Csc[x/2]^2 EllipticK[-Cos[x] Csc[x/2]^4], {x, 0, Pi}] / Pi^3 // Timing]

{0.021917, 1.3932039076998044}

Print[ Gamma[1/4]^4 / (4 Pi^3) // N // Timing]

{0.000029, 1.3932039296856769}

Print[ 4 / Pi^2 EllipticK[1/2]^2 // N // Timing]

{0.000375, 1.3932039296856764}

For `NIntegrate` the default method is `GlobalAdaptive`. You can try different methods and rules to optimize for the function and bounds.