Merge 3.2: issue #14937.

1 files changed, 256 insertions, 0 deletions

diff --git a/Modules/_decimal/libmpdec/literature/matrix-transform.txt b/Modules/_decimal/libmpdec/literature/matrix-transform.txt
new file mode 100644
index 0000000000..ff6219810e
--- /dev/null
+++ b/Modules/_decimal/libmpdec/literature/matrix-transform.txt
@@ -0,0 +1,256 @@
+
+
+(* Copyright (c) 2011 Stefan Krah. All rights reserved. *)
+
+
+The Matrix Fourier Transform:
+=============================
+
+In libmpdec, the Matrix Fourier Transform [1] is called four-step transform
+after a variant that appears in [2]. The algorithm requires that the input
+array can be viewed as an R*C matrix.
+
+All operations are done modulo p. For readability, the proofs drop all
+instances of (mod p).
+
+
+Algorithm four-step (forward transform):
+----------------------------------------
+
+  a := input array
+  d := len(a) = R * C
+  p := prime
+  w := primitive root of unity of the prime field
+  r := w**((p-1)/d)
+  A := output array
+
+  1) Apply a length R FNT to each column.
+
+  2) Multiply each matrix element (addressed by j*C+m) by r**(j*m).
+
+  3) Apply a length C FNT to each row.
+
+  4) Transpose the matrix.
+
+
+Proof (forward transform):
+--------------------------
+
+  The algorithm can be derived starting from the regular definition of
+  the finite-field transform of length d:
+
+            d-1
+           ,----
+           \
+    A[k] =  |  a[l]  * r**(k * l)
+           /
+           `----
+           l = 0
+
+
+  The sum can be rearranged into the sum of the sums of columns:
+
+            C-1     R-1
+           ,----   ,----
+           \       \
+         =  |       |  a[i * C + j] * r**(k * (i * C + j))
+           /       /
+           `----   `----
+           j = 0   i = 0
+
+
+  Extracting a constant from the inner sum:
+
+            C-1           R-1
+           ,----         ,----
+           \             \
+         =  |  r**k*j  *  |  a[i * C + j] * r**(k * i * C)
+           /             /
+           `----         `----
+           j = 0         i = 0
+
+
+  Without any loss of generality, let k = n * R + m,
+  where n < C and m < R:
+
+                C-1                          R-1
+               ,----                        ,----
+               \                            \
+    A[n*R+m] =  |  r**(R*n*j) * r**(m*j)  *  |  a[i*C+j] * r**(R*C*n*i) * r**(C*m*i)
+               /                            /
+               `----                        `----
+               j = 0                        i = 0
+
+
+  Since r = w ** ((p-1) / (R*C)):
+
+     a) r**(R*C*n*i) = w**((p-1)*n*i) = 1
+
+     b) r**(C*m*i) = w**((p-1) / R) ** (m*i) = r_R ** (m*i)
+
+     c) r**(R*n*j) = w**((p-1) / C) ** (n*j) = r_C ** (n*j)
+
+     r_R := root of the subfield of length R.
+     r_C := root of the subfield of length C.
+
+
+                C-1                             R-1
+               ,----                           ,----
+               \                               \
+    A[n*R+m] =  |  r_C**(n*j) * [ r**(m*j)  *   |  a[i*C+j] * r_R**(m*i) ]
+               /                     ^         /
+               `----                 |         `----    1) transform the columns
+               j = 0                 |         i = 0
+                 ^                   |
+                 |                   `-- 2) multiply
+                 |
+                 `-- 3) transform the rows
+
+
+   Note that the entire RHS is a function of n and m and that the results
+   for each pair (n, m) are stored in Fortran order.
+
+   Let the term in square brackets be f(m, j). Step 1) and 2) precalculate
+   the term for all (m, j). After that, the original matrix is now a lookup
+   table with the mth element in the jth column at location m * C + j.
+
+   Let the complete RHS be g(m, n). Step 3) does an in-place transform of
+   length n on all rows. After that, the original matrix is now a lookup
+   table with the mth element in the nth column at location m * C + n.
+
+   But each (m, n) pair should be written to location n * R + m. Therefore,
+   step 4) transposes the result of step 3).
+
+
+
+Algorithm four-step (inverse transform):
+----------------------------------------
+
+  A  := input array
+  d  := len(A) = R * C
+  p  := prime
+  d' := d**(p-2)             # inverse of d
+  w  := primitive root of unity of the prime field
+  r  := w**((p-1)/d)         # root of the subfield
+  r' := w**((p-1) - (p-1)/d) # inverse of r
+  a  := output array
+
+  0) View the matrix as a C*R matrix.
+
+  1) Transpose the matrix, producing an R*C matrix.
+
+  2) Apply a length C FNT to each row.
+
+  3) Multiply each matrix element (addressed by i*C+n) by r**(i*n).
+
+  4) Apply a length R FNT to each column.
+
+
+Proof (inverse transform):
+--------------------------
+
+  The algorithm can be derived starting from the regular definition of
+  the finite-field inverse transform of length d:
+
+                  d-1
+                 ,----
+                 \
+    a[k] =  d' *  |  A[l]  * r' ** (k * l)
+                 /
+                 `----
+                 l = 0
+
+
+  The sum can be rearranged into the sum of the sums of columns. Note
+  that at this stage we still have a C*R matrix, so C denotes the number
+  of rows:
+
+                  R-1     C-1
+                 ,----   ,----
+                 \       \
+         =  d' *  |       |  a[j * R + i] * r' ** (k * (j * R + i))
+                 /       /
+                 `----   `----
+                 i = 0   j = 0
+
+
+  Extracting a constant from the inner sum:
+
+                  R-1                C-1
+                 ,----              ,----
+                 \                  \
+         =  d' *  |  r' ** (k*i)  *  |  a[j * R + i] * r' ** (k * j * R)
+                 /                  /
+                 `----              `----
+                 i = 0              j = 0
+
+
+  Without any loss of generality, let k = m * C + n,
+  where m < R and n < C:
+
+                     R-1                                  C-1
+                    ,----                                ,----
+                    \                                    \
+    A[m*C+n] = d' *  |  r' ** (C*m*i) *  r' ** (n*i)   *  |  a[j*R+i] * r' ** (R*C*m*j) * r' ** (R*n*j)
+                    /                                    /
+                    `----                                `----
+                    i = 0                                j = 0
+
+
+  Since r' = w**((p-1) - (p-1)/d) and d = R*C:
+
+     a) r' ** (R*C*m*j) = w**((p-1)*R*C*m*j - (p-1)*m*j) = 1
+
+     b) r' ** (C*m*i) = w**((p-1)*C - (p-1)/R) ** (m*i) = r_R' ** (m*i)
+
+     c) r' ** (R*n*j) = r_C' ** (n*j)
+
+     d) d' = d**(p-2) = (R*C) ** (p-2) = R**(p-2) * C**(p-2) = R' * C'
+
+     r_R' := inverse of the root of the subfield of length R.
+     r_C' := inverse of the root of the subfield of length C.
+     R'   := inverse of R
+     C'   := inverse of C
+
+
+                     R-1                                      C-1
+                    ,----                                    ,----  2) transform the rows of a^T
+                    \                                        \
+    A[m*C+n] = R' *  |  r_R' ** (m*i) * [ r' ** (n*i) * C' *  |  a[j*R+i] * r_C' ** (n*j) ]
+                    /                           ^            /       ^
+                    `----                       |            `----   |
+                    i = 0                       |            j = 0   |
+                      ^                         |                    `-- 1) Transpose input matrix
+                      |                         `-- 3) multiply             to address elements by
+                      |                                                     i * C + j
+                      `-- 3) transform the columns
+
+
+
+   Note that the entire RHS is a function of m and n and that the results
+   for each pair (m, n) are stored in C order.
+
+   Let the term in square brackets be f(n, i). Without step 1), the sum
+   would perform a length C transform on the columns of the input matrix.
+   This is a) inefficient and b) the results are needed in C order, so
+   step 1) exchanges rows and columns.
+
+   Step 2) and 3) precalculate f(n, i) for all (n, i). After that, the
+   original matrix is now a lookup table with the ith element in the nth
+   column at location i * C + n.
+
+   Let the complete RHS be g(m, n). Step 4) does an in-place transform of
+   length m on all columns. After that, the original matrix is now a lookup
+   table with the mth element in the nth column at location m * C + n,
+   which means that all A[k] = A[m * C + n] are in the correct order.
+
+
+-- 
+
+  [1] Joerg Arndt: "Matters Computational"
+      http://www.jjj.de/fxt/
+  [2] David H. Bailey: FFTs in External or Hierarchical Memory
+      http://crd.lbl.gov/~dhbailey/dhbpapers/
+
+
+