MONSTER functions in Wolfram Language

Oleg Marichev, Daniel McDonald
WRI

Join the Conversation #WolframTechConf

Introduction to Monster Functions

What makes a function a monster?

They depend on arbitrarily large sets of parameters {(Subscript[ν, 1]),(Subscript[ν, 2]),…,(Subscript[ν, n])} and variables {(Subscript[z, 1]),(Subscript[z, 2]),…,(Subscript[z, m])}.

They are defined in terms of sums/products/integrals.

They generalize a wide array of more specialized functions.

And/Or Oleg says so.

Example: Generalized Hypergeometric Function

?HypergeometricPFQ

(Subscript[, p])(Subscript[F, q])((Subscript[a, 1]),…,(Subscript[a, p]);(Subscript[b, 1]),…,(Subscript[b, q]);z)==(∑_(k=0)^∞)(((∏_(j=1)^p)(Subscript[((Subscript[a, j])), k]) (z^k))/((∏_(j=1)^q)(Subscript[((Subscript[b, j])), k]) k!))/;q≥p∨q==p-1∧\[LeftBracketingBar]z\[RightBracketingBar]<1∨q==p-1∧\[LeftBracketingBar]z\[RightBracketingBar]==1∧Re((∑_(j=1)^(p-1))(Subscript[b, j])-(∑_(j=1)^p)(Subscript[a, j]))>0

HypergeometricPFQ[{1,1,3},{2,2},x]

HypergeometricPFQ[{-3,a},{b,c},x]

Meijer G-Function

Definition

?MeijerG

Defined through Mellin-Barnes contour integral:

(Subscript[G, p,q]^(m,n))(z\[VerticalSeparator]({{(Subscript[a, 1]),…,(Subscript[a, n]),(Subscript[a, n+1]),…,(Subscript[a, p])}, {(Subscript[b, 1]),…,(Subscript[b, m]),(Subscript[b, m+1]),…,(Subscript[b, q])}}))==(1/(2 π I))(∫_ℒ)((((∏_(k=1)^m)Γ(s+(Subscript[b, k]))) (∏_(k=1)^n)Γ(1-(Subscript[a, k])-s))/(((∏_(k=n+1)^p)Γ(s+(Subscript[a, k])) )(∏_(k=m+1)^q)Γ(1-(Subscript[b, k])-s)))(z^(-s))\[DifferentialD]s/;m∈ℕ∧n∈ℕ∧p∈ℕ∧q∈ℕ∧m≤q∧n≤p

The infinite contour of integration ℒ separates the poles of Γ(1-(Subscript[a, k])-s)at s=1-(Subscript[a, k])+j, j∈ℕ from the poles of Γ((Subscript[b, i])+ s) at s=-(Subscript[b, i])-l, l∈ℕ. Such a contour always exists in the cases  (Subscript[a, k])-(Subscript[b, i])-1∉ℕ.

Generalizing Other Functions

MeijerG[{{},{}},{{v},{-v}},z]

(side note on Bessel functions: http://blog.wolfram.com/2016/05/16/new-derivatives-of-the-bessel-functions-have-been-discovered-with-the-help-of-the-wolfram-language/ )

MeijerG[{{(1/2),(1/2)},{}},{{0},{0}},x]

MeijerG[{{(1/2),(1/2),(1/2)},{}},{{0},{-(3/2),-(3/2)}},-x]

MeijerGReduce[HypergeometricPFQ[{a},{b1,b2},x],x]

Usefulness in Integration

Integrate[KelvinKer[1/2,x],x]

Integrate[x Exp[-((x^2)+(1/x))],{x,0,Infinity}]

Would Oleg like to discuss with you more about MeijerG after our talk?

Yes.

Bell Polynomials

?BellY

Partial Bell Polynomials

(Subscript[Y, n,k])((Subscript[x, 1]),…,(Subscript[x, n-k+1]))= (∑_((Subscript[m, 1])=0)^n)⋯(∑_((Subscript[m, n])=0)^n)((n!)/((Subscript[m, 1])!⋯  (Subscript[m, n])!))Boole[(Subscript[m, 1])+2 (Subscript[m, 2])+⋯+n (Subscript[m, n])==n∧(Subscript[m, 1])+(Subscript[m, 2])+⋯+(Subscript[m, n])==k] (∏_(s=1)^n)((((Subscript[x, s])/s!))^(Subscript[m, s]))

Table[BellY[n,k,Array[x,n-k+1]],{n,1,7},{k,1,n}]//TableForm

Faà di Bruno’s formula

Chain rule gives an expression for the first derivative of a composition of two functions:

D[f[g[z]],z]

One way to generalize the chain rule is to allow higher-order derivatives.

Faà di Bruno’s formula expresses the n(^th) derivative of a composition of two functions:

(Subscript[D, t]^n)f(g(t))==(∑_(k=0)^n)(f^((k)))(g(t)) (Subscript[Y, n,k])(g'(t),g''(t),…,(g^((n-k+1)))(t)).

Faà di Bruno's formula for the third derivative of f(g(z)):

Sum[Derivative[k][f][g[z]] BellY[3,k,{g'[z],g''[z],Derivative[3][g][z]}],{k,1,3}]

D[f[g[z]],{z,3}]

The Wolfram Language now uses Faà di Bruno’s formula for unspecified higher-order derivatives of functional compositions (see http://blog.wolfram.com/2017/05/25/new-in-the-wolfram-language-enhanced-derivatives/):

D[f[g[x]],{x,n}]

Applying Faà di Bruno’s formula twice for compositions of three functions is a bit unwieldy though...

D[f[g[h[x]]],{x,n}]

Generalized Bell Polynomials

Complicated Recursive Definition

(BellY1[({{(Subscript[z, 1,1]), (Subscript[z, 1,2]), …, (Subscript[z, 1,m])}, {(Subscript[z, 2,1]), (Subscript[z, 2,2]), …, (Subscript[z, 2,m])}, {⋮, ⋮, ⋮, ⋮}, {(Subscript[z, n,1]), (Subscript[z, n,2]), …, (Subscript[z, n,m])}})])==(∑_(s=1)^n)(∑_((Subscript[s, 1])=0)^n)(∑_((Subscript[s, 2])=0)^n)… (∑_((Subscript[s, n])=0)^n)(KroneckerDeltaSeq[s,(∑_(j=1)^n)(Subscript[s, j])]) (KroneckerDeltaSeq[n,(∑_(j=1)^n)j (Subscript[s, j])]) ((n!)/((∏_(j=1)^n)(Subscript[s, j])!)) (Subscript[z, 1,s]) (((((BellY1[({{(Subscript[z, 1,2]), …, (Subscript[z, 1,m])}})])/(1!))))^(Subscript[s, 1])) (((((BellY1[({{(Subscript[z, 1,2]), …, (Subscript[z, 1,m])}, {(Subscript[z, 2,2]), …, (Subscript[z, 2,m])}})])/(2!))))^(Subscript[s, 2])) … (((((BellY1[({{(Subscript[z, 1,2]), …, (Subscript[z, 1,m])}, {(Subscript[z, 2,2]), …, (Subscript[z, 2,m])}, {⋮, ⋮, ⋮}, {(Subscript[z, n,2]), …, (Subscript[z, n,m])}})])/(n!))))^(Subscript[s, n]))/; (BellY1[({{(Subscript[z, 1,1])}, {(Subscript[z, 2,1])}, {⋮}, {(Subscript[z, n,1])}})])==(Subscript[z, n,1])

(tabellY=Table[BellY[Array[(Subscript[z, #1,#2])&,{n,m}]],{n,3},{m,3}])//TableForm

Y?

Reverse-engineered to express the n(^th) derivative of a composition of m functions (in a sense, the definition above comes from just applying Faà di Bruno’s formula m-1 times):

(Subscript[D, t]^n)(Subscript[f, 1])((Subscript[f, 2])((Subscript[…f, m])(t)…))==Y(({{(Subscript[f, 1])'((Subscript[f, 2])(… (Subscript[f, m])(t)…)), (Subscript[f, 2])'((Subscript[f, 3])(… (Subscript[f, m])(t)…)), ⋯, (Subscript[f, m])'(t)}, {(Subscript[f, 1])''((Subscript[f, 2])(… (Subscript[f, m])(t)…)), (Subscript[f, 2])''((Subscript[f, 3])(… (Subscript[f, m])(t)…)), ⋯, (Subscript[f, m])''(t)}, {⋮, ⋮, ⋱, ⋮}, {((Subscript[f, 1])^((n)))((Subscript[f, 2])(… (Subscript[f, m])(t)…)), ((Subscript[f, 2])^((n)))((Subscript[f, 3])(… (Subscript[f, m])(t)…)), ⋯, ((Subscript[f, m])^((n)))(t)}}))

Table[D[(Composition@@Array[(Subscript[f, #])&,m])@x,{x,n}],{n,1,3},{m,1,3}]//TableForm

Table[tabellY[[n,m]]/.((Subscript[z, i_,j_]):>(D[(Subscript[f, j])[x],{x,i}]/.x:>(Composition@@Array[(Subscript[f, #])&,m-j,j+1])[x])),{n,3},{m,3}]//Expand//TableForm

In case you were wondering...

Generalized partial Bell polynomials can be defined in terms of generalized Bell polynomials:

(BellYMatrix[n, k, ({{(Subscript[z, 1,1]), (Subscript[z, 1,2]), …, (Subscript[z, 1,m])}, {⋮, ⋮, ⋮, ⋮}, {(Subscript[z, k,1]), (Subscript[z, k,2]), …, (Subscript[z, k,m])}, {⋮, ⋮, ⋮, ⋮}, {(Subscript[z, n,1]), (Subscript[z, n,2]), …, (Subscript[z, n,m])}})])==(BellY1[({{0, (Subscript[z, 1,1]), (Subscript[z, 1,2]), …, (Subscript[z, 1,m])}, {⋮, ⋮, ⋮, ⋮, ⋮}, {0, (Subscript[z, k-1,1]), (Subscript[z, k-1,2]), …, (Subscript[z, k-1,m])}, {1, (Subscript[z, k,1]), (Subscript[z, k,2]), …, (Subscript[z, k,m])}, {0, (Subscript[z, k+1,1]), (Subscript[z, k+1,2]), …, (Subscript[z, k+1,m])}, {⋮, ⋮, ⋮, ⋮, ⋮}, {0, (Subscript[z, n,1]), (Subscript[z, n,2]), …, (Subscript[z, n,m])}})])/;n∈ℤ∧k∈ℤ∧m∈ℤ∧n≥k≥1∧m≥1

Quick Computation of Generalized Bell Polynomials

A Recursive Functional Identity

Recall how we defined Y to compute (Subscript[D, t]^n)(Subscript[f, 1])((Subscript[f, 2])((Subscript[⋯f, m])(t)⋯)):

Y(({{(Subscript[f, 1])'((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])'((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])'(t)}, {(Subscript[f, 1])''((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])''((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])''(t)}, {⋮, ⋮, ⋱, ⋮}, {((Subscript[f, 1])^((n)))((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), ((Subscript[f, 2])^((n)))((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, ((Subscript[f, m])^((n)))(t)}}))=(Subscript[D, t]^n)(Subscript[f, 1])((Subscript[f, 2])((Subscript[⋯f, m])(t)⋯))

Since composition is associative, we can get the following recursive identity for Y by letting 1≤(Subscript[m, 1])<(Subscript[m, 2])<⋯<(Subscript[m, p])<m and noting that (Subscript[f, 1])∘(Subscript[f, 2])∘⋯∘(Subscript[f, m])=((Subscript[f, 1])∘(Subscript[f, 2])∘⋯∘(Subscript[f, (Subscript[m, 1])]))∘((Subscript[f, (Subscript[m, 1])+1])∘(Subscript[f, (Subscript[m, 1])+2])∘⋯∘(Subscript[f, (Subscript[m, 2])]))∘⋯∘((Subscript[f, (Subscript[m, p])+1])∘(Subscript[f, (Subscript[m, p])+2])∘⋯∘(Subscript[f, m])).

Y(({{(Subscript[f, 1]^1)((Subscript[t, 1])), (Subscript[f, 2]^1)((Subscript[t, 2])), ⋯, (Subscript[f, m]^1)((Subscript[t, m]))}, {(Subscript[f, 1]^2)((Subscript[t, 1])), (Subscript[f, 2]^2)((Subscript[t, 2])), ⋯, (Subscript[f, m]^2)((Subscript[t, m]))}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[f, 1]^n)((Subscript[t, 1])), (Subscript[f, 2]^n)((Subscript[t, 2])), ⋯, (Subscript[f, m]^n)((Subscript[t, m]))}}))=Y(({{Y(({{(Subscript[f, 1]^1)((Subscript[t, 1])), (Subscript[f, 2]^1)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^1)((Subscript[t, (Subscript[m, 1])]))}})) , Y(({{(Subscript[f, (Subscript[m, 1])+1]^1)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^1)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^1)((Subscript[t, (Subscript[m, 2])]))}})) , ⋯, Y(({{(Subscript[f, (Subscript[m, p])+1]^1)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^1)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^1)((Subscript[t, m]))}})) }, {Y(({{(Subscript[f, 1]^1)((Subscript[t, 1])), (Subscript[f, 2]^1)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^1)((Subscript[t, (Subscript[m, 1])]))}, {(Subscript[f, 1]^2)((Subscript[t, 1])), (Subscript[f, 2]^2)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^2)((Subscript[t, (Subscript[m, 1])]))}})) , Y(({{(Subscript[f, (Subscript[m, 1])+1]^1)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^1)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^1)((Subscript[t, (Subscript[m, 2])]))}, {(Subscript[f, (Subscript[m, 1])+1]^2)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^2)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^2)((Subscript[t, (Subscript[m, 2])]))}})) , ⋯, Y(({{(Subscript[f, (Subscript[m, p])+1]^1)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^1)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^1)((Subscript[t, m]))}, {(Subscript[f, (Subscript[m, p])+1]^2)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^2)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^2)((Subscript[t, m]))}})) }, {⋮, ⋮, ⋱, ⋮}, {Y(({{(Subscript[f, 1]^1)((Subscript[t, 1])), (Subscript[f, 2]^1)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^1)((Subscript[t, (Subscript[m, 1])]))}, {(Subscript[f, 1]^2)((Subscript[t, 1])), (Subscript[f, 2]^2)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^2)((Subscript[t, (Subscript[m, 1])]))}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[f, 1]^n)((Subscript[t, 1])), (Subscript[f, 2]^n)((Subscript[t, 2])), ⋯, (Subscript[f, (Subscript[m, 1])]^n)((Subscript[t, (Subscript[m, 1])]))}})) , Y(({{(Subscript[f, (Subscript[m, 1])+1]^1)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^1)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^1)((Subscript[t, (Subscript[m, 2])]))}, {(Subscript[f, (Subscript[m, 1])+1]^2)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^2)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^2)((Subscript[t, (Subscript[m, 2])]))}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[f, (Subscript[m, 1])+1]^n)((Subscript[t, (Subscript[m, 1])+1])), (Subscript[f, (Subscript[m, 1])+2]^n)((Subscript[t, (Subscript[m, 1])+2])), ⋯, (Subscript[f, (Subscript[m, 2])]^n)((Subscript[t, (Subscript[m, 2])]))}})) , ⋯, Y(({{(Subscript[f, (Subscript[m, p])+1]^1)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^1)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^1)((Subscript[t, m]))}, {(Subscript[f, (Subscript[m, p])+1]^2)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^2)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^2)((Subscript[t, m]))}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[f, (Subscript[m, p])+1]^n)((Subscript[t, (Subscript[m, p])+1])), (Subscript[f, (Subscript[m, p])+2]^n)((Subscript[t, (Subscript[m, p])+2])), ⋯, (Subscript[f, m]^n)((Subscript[t, m]))}})) }}))

Setting (Subscript[m, 1])=(Subscript[m, p])=m/2 (assuming m is even), we get the following simplified identity.

Y(({{(Subscript[f, 1])'((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])'((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])'(t)}, {(Subscript[f, 1])''((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])''((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])''(t)}, {⋮, ⋮, ⋱, ⋮}, {((Subscript[f, 1])^((n)))((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), ((Subscript[f, 2])^((n)))((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, ((Subscript[f, m])^((n)))(t)}}))=Y(({{Y(({{(Subscript[f, 1])'((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])'((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m/2])'((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}})) , Y(({{(Subscript[f, m/2+1])'((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, m/2+2])'((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])'(t)}}))}, {Y(({{(Subscript[f, 1])'((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])'((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m/2])'((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}, {(Subscript[f, 1])''((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])''((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m/2])''((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}})) , Y(({{(Subscript[f, m/2+1])'((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, m/2+2])'((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])'(t)}, {(Subscript[f, m/2+1])''((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, m/2+2])''((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])''(t)}}))}, {⋮, ⋮}, {Y(({{(Subscript[f, 1])'((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])'((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m/2])'((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}, {(Subscript[f, 1])''((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, 2])''((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m/2])''((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}, {⋮, ⋮, ⋱, ⋮}, {((Subscript[f, 1])^((n)))((Subscript[f, 2])(⋯ (Subscript[f, m])(t)⋯)), ((Subscript[f, 2])^((n)))((Subscript[f, 3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, ((Subscript[f, m/2])^((n)))((Subscript[f, m/2+1])(⋯ (Subscript[f, m])(t)⋯))}})) , Y(({{(Subscript[f, m/2+1])'((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, m/2+2])'((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])'(t)}, {(Subscript[f, m/2+1])''((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), (Subscript[f, m/2+2])''((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, (Subscript[f, m])''(t)}, {⋮, ⋮, ⋱, ⋮}, {((Subscript[f, m/2+1])^((n)))((Subscript[f, m/2+2])(⋯ (Subscript[f, m])(t)⋯)), ((Subscript[f, m/2+2])^((n)))((Subscript[f, m/2+3])(⋯ (Subscript[f, m])(t)⋯)), ⋯, ((Subscript[f, m])^((n)))(t)}}))}}))

A Recursive Program (except this is WL so it’s probably more proper to use bottom-up dynamic programming instead of top-down recursion)

We can use this identity (along with Faà di Bruno’s formula as the base case) to create the following implementation of Y.

(Note that (Subscript[f, 1])((Subscript[f, 2])((Subscript[⋯f, m])(t)⋯))=(Subscript[f, 1])((Subscript[f, 2])((Subscript[⋯f, m])((Subscript[f, m+1])((Subscript[⋯f, M])(t)⋯))⋯)) for (Subscript[f, m+1])(t)=(Subscript[f, m+2])(t)=⋯=(Subscript[f, M])(t)=t, which allows us to add functions to the composition until M is a power of 2.)

bellY[matrix_]:=Module[{rownum,colnum,depth,newcolnum,coldiff,idmat,newmat,fdb,reduce},
{rownum,colnum}=Dimensions[matrix];

(*To simplify our calculation, we pad our matrix with columns consisting of a 1 followed by 0s (corresponding to composition with the identity function (Subscript[f, i])(t)=t) until the number of columns is a power of 2 greater than 1.*)
depth=⌈Log[2,colnum]⌉;
newcolnum=(2^depth);
coldiff=newcolnum-colnum;
idmat=Join[ConstantArray[1,{1,coldiff}],ConstantArray[0,{rownum-1,coldiff}]];
newmat=Join[matrix,idmat,2] ;
(*This is our new, padded matrix.*)

(*Implement Faà di Bruno's formula for a matrix having exactly two columns.*)
fdb[twocolmat_,ord_]:=(∑_(k=1)^ord)twocolmat[[k,1]]BellY[ord,k,twocolmat[[1;;ord-k+1,2]]] ;

(*Apply Faà di Bruno's formula to consecutive columns of the given matrix to halve the number of columns.*)
reduce[mat_,dim_]:=Table[fdb[mat[[1;;i,2j-1;;2j]],i],{i,dim[[1]]},{j,dim[[2]]}] ;

(*Apply the reduce function to iteratively halve the number of columns until one remains, then takes the last element of that column, which is the n(\[Null]^th) derivative.*)
Fold[reduce,newmat,Table[{rownum,newcolnum/(2^k)},{k,depth}]][[rownum,1]]

]

We now compare to the existing BellY function.

matrix=RandomReal[{.5,1.5},{19,33}];
bellY[matrix]
BellY[matrix]

Спасибо за то, что вы слушали!