# Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

Compute the GCD and successive quotients for a pair of Laurent polynomials

Contributed by:
Daniel Lichtblau

ResourceFunction["LaurentExtendedGCD"][ computes the extended GCD of the Laurent polynomials |

The polynomials must be univariate and have numeric coefficients.

The degree of a Laurent polynomial is the difference between the highest and lowest powers of the variable.

In contrast to the ring of ordinary polynomials, units are comprised of all monomials (since highest and lowest powers of the variable are the same).

A GCD for a pair of Laurent polynomials is unique up to unit factors.

Again in contrast to the case of ordinary polynomials, the list of quotients is not unique up to unit multiples. This is because we have the option of removing terms from either end in the process of taking a quotient.

ResourceFunction["LaurentExtendedGCD"] uses a method that attempts to keep coefficients as close to unity as possible in absolute value. The standard of proximity is relative; that is, it tries to minimize the maximum of the log of the coefficient absolute values. This is done via a greedy approach: nontrivial reductions typically allow either two or three possible quotients (reduce two coefficients at the top, two at the bottom, or one at each end, if degrees differ by one). The best one at each step is chosen by the criterion above.

ResourceFunction["LaurentExtendedGCD"] takes a Tolerance option. With the setting Tolerance→*tol*, an inexact numeric coefficient *c* is considered to be zero if Chop[*c*,*tol*] is zero.

The tolerance must be a non-negative real value.

When the coefficients are approximate numbers, the default tolerance is 10^{-8}.

The tolerance can also be specified as a pair {*tol*_{1},*tol*_{2}}. In this case *tol*_{1} is used as above while *tol*_{2} is a relative tolerance. In particular, if the coefficient list of a polynomial is *clist*, then Chop[*clist*,Max[Abs[*clist*]]*tol*] is also performed.

An important application of Laurent extended GCD is in wavelet computations, where the quotient sequence is used to create fast lifting filters.

In the wavelet setting one can also obtain a set of valid quotients and remainder using WaveletFilterCoefficients; it need not (and in general will not) be the same as the set given by ResourceFunction["LaurentExtendedGCD"].

Compute a quotient sequence and GCD for two polynomials:

In[1]:= |

Out[1]= |

Create a pair of random Laurent polynomials of degree 10:

In[2]:= |

Out[7]= |

Compute a quotient sequence and GCD:

In[8]:= |

Out[8]= |

Create two exact Laurent polynomials:

In[9]:= |

Out[10]= |

Compute a quotient sequence and GCD:

In[11]:= |

Out[11]= |

LaurentExtendedGCD can return a nontrivial GCD:

In[12]:= |

Out[13]= |

Compute a quotient sequence and GCD:

In[14]:= |

Out[14]= |

LaurentExtendedGCD works with approximate numeric coefficients and gives an equivalent result:

In[15]:= |

Out[15]= |

Compare the difference numerically:

In[16]:= |

Out[16]= |

Form the Laurent pair corresponding to the Daubechies 6 (or db3) wavelet filter coefficients:

In[17]:= |

Out[18]= |

Compute the GCD and quotients:

In[19]:= |

Out[19]= |

A set of quotients and remainder can also be obtained using WaveletFilterCoefficients with "LiftingFilter" and "LiftingLaurentForm":

In[20]:= |

Out[21]= |

Obtain the Daubechies 4 (D4, also called db2) wavelet coefficients:

In[22]:= |

Out[22]= |

Create two corresponding Laurent polynomials:

In[23]:= |

Out[23]= |

Compute the Laurent GCD and list of quotients:

In[24]:= |

Out[24]= |

Create a list of 2×2 matrices that represent the operation of reducing a pair of Laurent polynomials by their quotient:

In[25]:= |

Out[26]= |

Take their product:

In[27]:= |

Out[27]= |

Applying this product to the vector comprised of the GCD and 0 recovers the original pair of Laurent polynomials:

In[28]:= |

Out[28]= |

Invert the reducing matrices:

In[29]:= |

Out[30]= |

Form the reverse-product of these inverses:

In[31]:= |

Out[31]= |

Check that this is in fact the inverse:

In[32]:= |

Out[32]= |

This inverse product matrix provides the Bezout multipliers (top row) and syzygy (bottom row), as can be seen by its action on the input pair:

In[33]:= |

Out[33]= |

Create the quadrature mirror filter from the D4 coefficients:

In[34]:= |

Out[34]= |

A certain upper triangular matrix is needed to recover these Laurent polynomials:

In[35]:= |

Out[35]= |

Form the augmented set of 2×2 matrices required to create the wavelet polyphase matrix:

In[36]:= |

Out[37]= |

Compute the polyphase matrix:

In[38]:= |

Out[38]= |

Show that the second column is indeed the quadrature mirror filter:

In[39]:= |

Out[39]= |

The same matrix can be obtained from WaveletFilterCoefficients:

In[40]:= |

Out[41]= |

Create a function to recover the discrete wavelet transform (DWT) from these 2×2 matrix factors. For this, powers of the Laurent variable *z* become shifts with positive exponents being leftward, and indices are modulo the length of the sequence (that is, using periodic wrapping):

In[42]:= |

Work with {1,2,…,16} split into odds and evens (the so-called "Lazy Wavelet Transform"):

In[43]:= |

Compute the DWT:

In[44]:= |

Out[44]= |

Compare to the forward lift from the LiftingFilterData for the D2 wavelet:

In[45]:= |

Out[45]= |

Get the same result directly from LiftingWaveletTransform:

In[46]:= |

Out[46]= |

Invert the process to recover the original values, using Riffle at the end to undo the Lazy Wavelet Transform:

In[47]:= |

Out[48]= |

This is equivalent to using InverseWaveletTransform on the DiscreteWaveletData object:

In[49]:= |

Out[49]= |

The conventional way to use the quotient matrices and their inverses to form a lifting filter actually makes them alternating upper and lower triangular. For simplicity of exposition, this step was omitted.

Create the even and odd Laurent polynomials that arise from the Daubechies10 wavelet coefficients:

In[50]:= |

Out[52]= |

Compute an extended GCD with quotients coefficients that are reasonably close in magnitude to unity:

In[53]:= |

Out[53]= |

Another good set of quotients for the D10 Laurent pair is given in the reference by Maslen and Abbott; the latter happens to be the one provided by WaveletFilterCoefficients:

In[54]:= |

Out[55]= |

Create the quotient matrices and the scaling value (a factor of * z^{2}* is required for this particular factorization):

In[56]:= |

Out[64]= |

Set up functions that transform from Laurent polynomials to vector operations:

In[65]:= |

Pad a square image to have a power of 2 dimension, removing the average pixel value from each row:

In[66]:= |

Out[72]= |

Compute the DWT on rows and then compute the DWT on the columns of that result:

In[73]:= |

Create a sparse array, chopping values larger than 0.1 to reduce size:

In[74]:= |

Out[75]= |

Create the inverse operations:

In[76]:= |

Perform one inversion and compare to the first level of the previous discrete transform:

In[77]:= |

Out[78]= |

Perform a second inversion and compare to the input:

In[79]:= |

Out[81]= |

Recover a reasonable approximation to the original image:

In[82]:= |

Out[82]= |

- 1.0.0 – 08 February 2022

This work is licensed under a Creative Commons Attribution 4.0 International License