This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Consider the quotient q=(Subscript[z, 1])/(Subscript[z, 2]), where the (Subscript[z, 1]), (Subscript[z, 2]) are random complex numbers chosen with uniform probability from the square {1+i,2+i,2+2i,1+2i} in the complex plane. Calculate the probability distribution function p(Re(q),Im(q)) analytically.

{qRe, qIm} = 
Simplify[ComplexExpand[#[(x1 + I y1)/(x2 + I y2)], 
                    TargetFunctions -> {Re, Im}]]& /@ {Re, Im}

It is straightforward to calculate exact expressions for the averages of the real and imaginary parts of q. Intuitively one would expect the average value of q to be approximately 1. This is the case to a remarkable precision.

{qrAv, qiAv} = 
(Integrate[#, {y2, 1, 2}, {x2, 1, 2}, {y1, 1, 2}, {x1, 1, 2},
           GenerateConditions -> False] // 
           TrigToExp // FullSimplify)& /@ {qRe, qIm}

N[%]

It is also possible to calculate exactly the variances for the real and imaginary parts. This time, we obtain considerably more complicated expressions containing polylogarithmic functions. From the numerical values of the variances, we expect a half width of about 0.2 for the real and the imaginary part.

(FullSimplify[Integrate[FullSimplify[
 (* split off last integration and simplify before *)
 Integrate[#, {x2, 1, 2}, {y1, 1, 2}, {x1, 1, 2},
           GenerateConditions -> False]], {y2, 1, 2}] /. 
            xA -> qrAv]& /@ (({qRe, qIm} - {qrAv, qiAv})^2)) //
                     (* for a shorter result *) TraditionalForm

N[Sqrt[%], 22]  // N

Using quantifier elimination it is also easily possible to calculate bounds for the range of nonvanishing values for the quotient q. We find 1/2≤Re(q)≤2 and 1-φ≤Im(q)≤φ-1.

Function[reim, Function[lg, 
Exists[{x1, y1, x2, y2}, 
      Element[{x1, y1, x2, y2, q}, Reals] && 
              1 < x1 < 2 && 1 < y1 < 2 &&
              1 < x2 < 2 && 1 < y2 < 2, lg[q, reim]] // 
               Resolve] /@ {Less, Greater}] /@ {qRe, qIm}

Next, we carry out a numerical simulation that uses ten million pairs of random numbers. We use the compiled function qCF for carry simulation.

qCF = Compile[{{n, _Integer}, {pp, _Integer}},
        Module[{m = Table[0, {pp + 1}, {pp + 1}], q, qr, qi, 
                i, j, q0 = {0, -1}, λq = {2, 2}, norm}, 
            Do[(* form the quotient of the two random numbers *)
               q = Random[Complex, {1 + I, 2 + 2 I}]/ 
                   Random[Complex, {1 + I, 2 + 2 I}];
               {qr, qi} = {Re[q], Im[q]}; 
               (* discretize and count *)
               {j, i} = Round[({qr, qi} - q0)/λq pp] + 1;
               m[[i, j]] = m[[i, j]] + 1, {n}]; 
               norm = Plus @@ Flatten[Sign[m]];
               (* return normalized result *)
               m/n norm]];

The following graphic of the modelled probability distribution agrees quantitatively with the above-calculated parameters.

ListPlot3D[qCF[10^7, 60], PlotRange -> All,
           MeshRange -> {{0, 5/2}, {-1, 1}}, Mesh -> False]

Now, we calculate the exact probability distribution function. The bivariate distribution is given as the following four-fold integral:

p((Subscript[q, r]),(Subscript[q, i]))=(∫_1^2)(∫_1^2)(∫_1^2)(∫_1^2)δ((Subscript[q, r])-(((Subscript[x, 1]) (Subscript[x, 2])+(Subscript[y, 1]) (Subscript[y, 2]))/((Subscript[x, 2]^2)+(Subscript[y, 2]^2))))δ((Subscript[q, i])-(((Subscript[x, 2]) (Subscript[y, 1])-(Subscript[x, 1]) (Subscript[y, 2]))/((Subscript[x, 2]^2)+(Subscript[y, 2]^2))))d(Subscript[x, 1])d(Subscript[y, 1])d(Subscript[x, 2])d(Subscript[y, 2]).

Using the well known integration rules for the Dirac delta function, it is straightforward to carry out the integration over (Subscript[x, 1]) and (Subscript[y, 1]).

(* carry out integration over x1 *)
solx1 = x1 /.  Solve[qr - qRe == 0, x1][[1]];
fd1 = -D[qr - qRe, x1];
int1 = 1/fd1 UnitStep[solx1 - 1] UnitStep[2 - solx1];
(* carry out integration over y1; independent factor *)
qiP = (qi - qIm /. x1 -> solx1);
soly1 = y1 /. Solve[qiP == 0, y1][[1]] ;
fd2 = -D[qiP, y1];
int2 = 1/fd2 UnitStep[soly1 - 1] UnitStep[2 - soly1] // Simplify;
(* unite results *)
int3 = (int1 /. y1 -> soly1) int2 // Simplify

This leaves us with the following double integral:

p((Subscript[q, r]),(Subscript[q, i]))=(∫_1^2)(∫_1^2)((Subscript[x, 2]^2)+(Subscript[y, 2]^2)) θ((Subscript[q, r])(Subscript[x, 2])-(Subscript[q, i]) (Subscript[y, 2])-1) θ(2-(Subscript[q, r]) (Subscript[x, 2])+(Subscript[q, i]) (Subscript[y, 2])) θ(2-(Subscript[q, i])(Subscript[x, 2])-(Subscript[q, r])(Subscript[y, 2])) θ((Subscript[q, i])(Subscript[x, 2])+(Subscript[q, r]) (Subscript[y, 2])-1)d(Subscript[x, 2])d(Subscript[y, 2]).

We subdivide the (Subscript[x, 2]),(Subscript[y, 2]),(Subscript[q, r]),(Subscript[q, i])-space into cells using cylindrical algebraic decomposition. Because we are interested only in cells of finite 4D Lebesgue measure, we use the function Experimental`GenericCylindricalAlgebraicDecomposition. This function ignores boundaries having a vanishing 4D Lebesgue measure. Inside each cell, all four arguments of the Heaviside functions will have arguments of definite signs. This will allow to ignore the Heaviside theta function for process of integration.

dC[{x2_, y2_}, {qr_, qi_}] := 
   -1 + qr x2 - qi y2 > 0 &&  2 - qr x2 + qi y2 > 0 &&
    2 - qi x2 - qr y2 > 0 && -1 + qi x2 + qr y2 > 0 &&
	1 < x2 < 2 && 1 < y2 < 2

cad = ((Experimental`GenericCylindricalAlgebraicDecomposition[
                    dC[{x2, y2}, {qr, qi}], {qr, qi, x2, y2}][[1]] /. 
                    (* expand all cells to use explicit description *)
                       Inequality -> inequality) // LogicalExpand) /. 
                                            inequality -> Inequality;
Length[cad]

We obtained 229 cells. Here are the parametrizations for the first two cells.

Take[cad, 2]

Now it is straightforward to carry out the remaining two integrations over (Subscript[x, 2]) and (Subscript[y, 2]). The function doThex2y2Integration integrates (Subscript[x, 2]^2)+(Subscript[y, 2]^2) over each cell and adds Heaviside functions to restrict the (Subscript[q, r]) and (Subscript[q, i]) values appropriately.

doThex2y2Integration[a_And] := 
Module[{ix2, iy2, iqr, iqi, sInt1, dsInt1, sInt2, dsInt2, res},
 (* the variable ranges *)
 {ix2, iy2, iqr, iqi} = Cases[a, HoldPattern[
            Inequality[_, _, #, _, _]]][[1]]& /@ {x2, y2, qr, qi};
 (* carry out the x3=2 and y2 integration *)
 sInt1 = Integrate[x2^2 + y2^2, y2];
 dsInt1 = (sInt1 /. y2 -> iy2[[-1]]) - (sInt1 /. y2 -> iy2[[+1]]);
 sInt2 = Integrate[dsInt1, x2];
 dsInt2 = (sInt2 /. x2 -> ix2[[-1]]) - (sInt2 /. x2 -> ix2[[+1]]);
 (* canonicalize result *)
 res = Factor //@ Together[dsInt2];
 If[(iqr && iqi), Evaluate @ res, 0]]

Here is the result for the first cell.

doThex2y2Integration[cad[[1]]]

It is now a matter to seconds to obtain a closed-form expression for the distribution function.

q𝒫[qr_, qi_] = Plus @@ (doThex2y2Integration /@ cad);

q𝒫[qr_, qi_] = Plus @@
(Function[p, If @@ {p[[1, 1]], Together[Plus @@ (#[[2]]& /@ p)], 0}] /@ 
    Split[List @@ (doThex2y2Integration /@ cad), (#1[[1]] === #2[[1]])&]);

Piecewise the distribution function is a rational function. Overall, we obtain 45 different rational functions.

𝓅s = Union[(List @@ q𝒫[qr, qi]) /. If[_, i_, _] :> i];
𝓅s // Length

Here are three of them.

Take[𝓅s, 3]

Some of the rational function are quite complicated—they are quotients of two total degree 15 polynomials.

{First[#], Length[#]}& /@ Split[Sort[
   Exponent[({Numerator[#], Denominator[#]}& /@ 𝓅s) /. 
                              {qr -> C, qi -> C}, C]]]

The global maximum of p((Subscript[q, r]),(Subscript[q, i])) is p(1,0)=14/3.

To visualize the various rational functions, we number each of the 45 rational functions.

q𝒫C[qr_, qi_] = Evaluate //@ (q𝒫[qr, qi] /. If[cond_, int_, 0] -> 
                   If[cond, {int, Position[𝓅s, int][[1, 1]]}, {0, 0}]);
                   
(* unite If's with identical first condition *)                   
q𝒫C[qr_, qi_] = q𝒫C[qr, qi] //. 
    If[Inequality[qrl_, Less, qr, Less, qru_] && qip1_, r1__] +
    If[Inequality[qrl_, Less, qr, Less, qru_] && qip2_, r2__] :>
    If[Inequality[qrl,  Less, qr, Less, qru], 
       If[qip1, r1] + If[qip2, r2], {0, 0}];

By continuity, we could now continue the distribution functions to the boundaries. For a visualization, it is sufficient to avoid hitting the boundaries. We achieve this by using a slightly asymmetric plot domain.

data = Module[{ε = 1/Pi 10^-6, pp = 241}, N[#, 22]& @ 
       Table[q𝒫C[N[qr, 22], N[qi, 22]], {qi, -1, 1 + ε, (2 + ε)/pp},
                          {qr, 0, 5/2 + ε, (5/2 + ε)/pp}]];

Here is a plot of the resulting probability distribution function. We use a different (randomly chosen) color for each of the 45 rational functions covering the plane.

randomColor[0] := RGBColor[0, 0, 0];
randomColor[x_] := randomColor[x] = ToColor[Hue[Random[]], RGBColor]

Show[Graphics3D[
 ListPlot3D[Map[First, data, {2}], 
            Map[randomColor[Round[#]]&, Map[Last, data, {2}], {2}], 
            PlotRange -> All, MeshRange -> {{0, 5/2}, {-1, 1}},
            Mesh -> False, DisplayFunction -> Identity]] /. 
                                  h_Hue :> SurafeColor[h, h, 2.8],
     DisplayFunction -> $DisplayFunction]