Selberg integral

Selberg integral is a multi-dimensional integral involving . A nice review can be found in ArXiv article by P. Forrester and S. Warnaar.

An example of Selberg integral was considered by Mehta:

Attempting to verify this beautiful formula using Mathematica for small values of n:

In[1] := Table[Timing[Integrate[
Product[((x[i]-x[j])^2),{i,n},{j,i+1,n}]*
   Exp[-n Sum[x[i]^2,{i,n}]],
Sequence @@ Table[{x[i], -Infinity, Infinity}, {i, n}]]
], {n, 2, 5}]

Out[1] = {{0.86, Pi/4}, {5.703, Pi^(3/2)/(54 Sqrt[3])},
 {20.875, (9 Pi^2)/131072},
 {147.078, (27 Pi^(5/2))/(195312500 Sqrt[5])}}

In[2]:= Table[(2Pi)^(n/2)/(2n)^(n^2/2) Product[k!,{k,n}],{n,2,5}]

Out[2]= {Pi/4, Pi^(3/2)/(54 Sqrt[3]), (9 Pi^2)/131072,
(27 Pi^(5/2))/(195312500 Sqrt[5])}

In order to achieve more efficient evaluation, and check the formula for higher values of n, we note that Integrate[x^(m-1) Exp[-n x^2],{x,-Infinity,Infinity}]==Sin[Pi m/2]^2 Gamma[m/2]/n^(m/2). We can then integrate variable one by one, carefully expanding only when needed with Expand[expr,patt]:

In[1] := Table[Timing[(Pi/n)^(n/2)*
   Fold[Expand[#1, #2] //. {#2^(m_) :> 
    Cos[Pi*m/2]^2*Pochhammer[1/2, m/2]/n^(m/2)} /. {#2->0} &,
   Product[(x[i]-x[j])^2, {i,n}, {j,i+1,n}],
   Table[x[k], {k, n}]]], {n, 2, 7}]

Out[1] = {{0., Pi/4},
  {0., Pi^(3/2)/(54 Sqrt[3])},
  {0.015, (9 Pi^2)/131072},
  {0.125, (27 Pi^(5/2))/(195312500 Sqrt[5])},
  {2.016, (25 Pi^3)/3343537668096},
  {34.594, (273375 Pi^(7/2))/(875799914882589322976 Sqrt[7])}}

In[2]:= %[[All, 2]] ==
 Table[(2 Pi)^(n/2)/(2 n)^(n^2/2) *
  Array[Factorial, n, 1, Times], {n,2, 7}]

Out[2]= True

My 2GB RAM machine runs out of memory for n=8 for both methods, but
the speed-up of the latter method is clear.

Finding curcuit resistance

One of the advanced physics problems which had made an early impression on me, was to find a resistance of an n-steps ladder made out of identical resistors, and to solve this problem for an infinite steps ladder.
Solving of this problem in Mathematica makes it easy to deeper explore involved mathematical concepts, and these include continued fractions, algebraic numbers, recurrence equations to name a few. For instance, the resistance of an n-steps ladder is
Despite appearance of algebraic numbers in the answer, resistance is a rational number for every natural integer n. First few of these read 3 R, 11 R/4, 41 R/15 and so on.

The solution can be downloaded as a Mathematica Player notebook. One would need freely available Mathematica Player to view it and evaluate it. Enjoy!

Divisibility by 7

Some divisibility criteria are being taught in elementary school, some in middle school. The following is the summary of what a high schooler ought to know.

1. A number is divisible by 2, 5, or 10 if its last digit is divisible by, correspondingly, 2, 5, or 10.
2. A number is divisible by 3, or 9 if the total of its digits is divisible by, correspondingly, 3, or 9.
3. A number is divisible by 4 if its last 2 digits are divisible by 4.
4. A number is divisible by 8 if its last 3 digits are divisible by 8.
   

One might wonder when is the number divisible by 7 ? To answer this question it is useful to remember some properties of modular arithmetic. These are

1. (a*b) mod n = (a mod n)*(b mod n) mod n
2. (a+b) mod n = (a mod n) + (b mod n) mod n

Thus for n = d[0] + d[1]*10^1 + d[2]*10^2 + ... + d[k]*10^k,

n mod 7 = sum( d[k] (10^k mod 7), over k )

Quick check in Mathematica shows that 10^k mod 7 is periodic with period 6:

In[30]:= Table[Mod[PowerMod[10, k, 7], 7, -3], {k, 0, 12}]

Out[30]= {1, 3, 2, -1, -3, -2, 1, 3, 2, -1, -3, -2, 1}

This gives the divisibility by 7 test. The number is divisible by 7 if the cyclic sum of its digits (starting with rightmost one) with {1,3,2,-1,-3,-2} is divisible by 7. The following is a Mathematica code to compute this cyclic sum:

DivisibilityBySevenTest[n_Integer] :=
With[{dig = Mod[Reverse[IntegerDigits[n]],7]},
   dig.PadRight[{},Length[dig], {1,3,2,-1,-3,-2}]
]

Consider an easy example. Let n = 28. Since 8 mod 7 = 1, we have (8 mod 7)*1+2*3 = 7, and, of course, 28 = 7*4.

In[47]:= DivisibilityBySevenTest[28]

Out[47]= 7

Now take n=2009. We do (9 mod 7)*1 + 0*3+0*2+2*(-1) = 2-2 = 0.

In[48]:= DivisibilityBySevenTest[2009]

Out[48]= 0

Indeed:

In[46]:= 2009/7

Out[46]= 287

When manually applying the divisibility test for large numbers, it might be easier to use the fact that 1000^k mod 7 is (-1)^(k-1). Thus, we might first form an alternating sum of triples of digits, and test the result afterwards. For example, let n = 123,456,789. It is divisible by 7 if 789 - 456+123 = 333+123 = 456 is divisible by 7. To check this, we compute 6*1+5*3+4*2 = 6+15+8=29 = 1 mod 7. Thus the number n is not divisible by 7. Indeed,

In[50]:= QuotientRemainder[123456789, 7]

Out[50]= {17636684, 1}

Applying the divisibility by 7 directly we get

In[53]:= DivisibilityBySevenTest[123456789]

Out[53]= -13