GSoC 2018 Application Sidhant Nagpal: Transforms, Convolution & Linear Recurrence Evaluation

Title

Transforms, Convolution & Linear Recurrence Evaluation

Table of Contents

• About Me

  • Personal Information
  • Personal Background
  • Programming Skills

• Contributions to SymPy

  • Pull Requests
  • Issues

• The Project

  • Problem & Motivation
  • Theory
  • Implementation

• Timeline

• References

About Me

Personal Information

Heading	Details
Name	Sidhant Nagpal
University	Netaji Subhas Institute of Technology, India
Email	sidhantnagpal97@gmail.com
GitHub/Gitter	sidhantnagpal
LinkedIn	sidhantnagpal
Timezone	IST (UTC+5:30)

Personal Background

I am Sidhant Nagpal, currently enrolled in my sixth semester (third year) for a four year bachelors degree (B.E.) in Computer Engineering (COE) at Netaji Subhas Institute of Technology (NSIT) Delhi, India. I have a strong Mathematics and Competitive Programming background. Since higher primary school, I have had a keen interest in Mathematics and Algebra, having participated in numerous Mathematical Olympiads - IMO (International Mathematics Olympiad) , IMMO (International Master Mathematics Olympiad), IOM (International Olympiad of Mathematics), RMO (Regional Mathematics Olympiad). I have also participated in Science Olympiads - KVPY (Kishore Vaigyanik Protsahan Yojana) Scholar, NTSE (National Talent Search Examination) Stage 2. Since the start of college, I have been avid for Algorithms & Data Structures, having solved over a 1000 algorithmic problems on various Online Judges - Codeforces, SPOJ, CodeChef, TopCoder. I have been presented with Ramanujan Award by NSIT for excellent academic performance in the Mathematics Courses which were a part of the curriculum. I have participated and achieved good ranks in Programming Contests - Google Code Jam (reached Round 2), Google Kickstart (Round B rank 60, Round G rank 105), ACM ICPC India Regionals (Amritapuri Onsite rank 41, Gwalior Onsite rank 41).

Programming Skills

• I work on Ubuntu 16.04 Xenial Xerus (LTS). I prefer Sublime Text as my primary text editor, because it is fast, responsive, and minimalistic. I have been using Sublime Text and GNU/Linux (Ubuntu) for about 3 years.
• I have been programming for about 5 years. I am proficient in C/C++ and Python. My GitHub Page has details about all the projects I've taken up.
• I am pretty comfortable in Python having used it for over 2 years. The features I find particularly useful in Python are generator expressions, list/dict comprehension, and argument unpacking.
• The functionality in SymPy that I particularly find very useful is

  • subs function Repetitive evaluation of complicated symbolic expressions and formulae can be very cumbersome, but subs function does it simply by merely substituting the symbols for values.

    >>> (x*z*cos(x*y*z)).subs({x:2, y:pi, z:S(1)/3})
    -1/3

  • @cacheit decorator Seamlessly reduces the complexity of dynamic programming problems. For example, the complexity of recursive evaluation of Pascal's Identity is reduced from exponential to quadratic, by just adding the @cacheit decorator.

    @cacheit
    def ncr(n, r):
      if r>n or r<0: return 0
      if r==0 or r==n: return 1
      return ncr(n-1, r) + ncr(n-1, r-1)

  • diff function Eases the computation of derivatives, especially for multivariate functions, of higher order.

    >>> diff(x*y**2*z**3, x, y, y, z, 3)
    12

• I am familiar with the concept of Version Control Systems, but have only used Git and GitHub, which I find very intuitive and simple to use.

Contributions to SymPy

I started contributing to SymPy about a month back, and in this duration, I learnt how impactful symbolic expressions can be in Computational Mathematics. I introduced myself on the SymPy mailing list, where the developers are active, and went through the Contribution Guidelines on SymPy wiki. Then, I directly hopped on to the SymPy development workflow, familiarizing myself with the codebase, which I found to be well-documented and readable. In case, I had some queries, I asked them on Gitter and got them resolved by active discussions. I started multiple issues which I found while using SymPy. I opened multiple PRs aimed at optimizing existing methods, adding new methods, and fixing open issues. I got timely feedback from the community for the same which helped me streamline the commits for the final merge.

Pull Requests

• (Merged) PR #14321: Inverse trigonometric evaluation for compatible arguments: Added support for evaluation of inverse trigonometric non-symbolic expressions with real arguments, taking into account the range of inverse trigonometric functions.
• (Merged) PR #14544: Fix evaluation of Reciprocal Trigonometric Functions: Updated eval in ReciprocalTrigonometricFunction to fix evaluation of sec and csc functions.
• (Merged) PR #14333: Resolve sign of modular inverse and handle negative modulo: Updated mod_inverse for negative modulo. Updated mod_inverse result to have same sign as the modulo.
• (Merged) PR #14278: Abs of base raised to exponent for complex non-symbol bases: Added support for computation of absolute value of complex non-symbolic powers.
• (Merged) PR #14249: Optimizing _eval_Mod in cores.power for very large powers: Integer modular exponentiation is optimized for very large powers, using Euler Totient Theorem and CRT (Chinese Remainder Theorem).
• (Merged) PR #14314: Simple optimization to speed up LCM: Precede multiplication (of operand) with division (of gcd) to speed up lcm.
• (Open) PR #14347: Matrix inverse in prime fields, gcdex for inv_mod: For matrix inv_mod, changed modular inverse calculation of determinant from Euler Totient Method to Extended Euclidean Algorithm. Added prime field support for Gaussian Elimination.
• (Open) PR #14365: Support gcd and lcm for compatible irrationals

Issues

• (Fixed) Issue #14320: Improper results with inverse trigonometric functions
• (Fixed) Issue #14543: Improper evaluation of Reciprocal Trigonometric Functions
• (Fixed) Issue #14332: Modular inverse for negative modulo and sign resolve
• (Fixed) Issue #14274: Abs((4 + 5*I)**(6 + 7*I))
• (Fixed) Issue #14260: abs(I**n)
• (Partially Fixed) Issue #14262: log(I**I)
• (Open) Issue #14311: Integral 1/(sin(x)**4 + cos(x)**4) can not be computed
• (Open) Issue #14549: Definite Integral for Beta Function hangs
• (Open) Issue #14355: Bidirectional limit of floor(sin(x)/x) at 0, is not correct
• (Open) Issue #14412: Sum(binomial(n, m)*sin(m*x)*cos((n-m)*x), (m, 0, n)).doit()
• (Open) Issue #14478: Incorrect evaluation of limit(x/2*floor(3/x),x,0,'+')
• (Open) Issue #14344: inv_mod in matrices.py uses totient calculation
• (Open) Issue #14345: Inverse of a matrix for prime fields is not optimal
• (Open) Issue #14364: gcd and lcm of similar type of irrationals
• (Open) Issue #14258: Rationalizing instead of taking LCM of the Denominator

The Project

Problem & Motivation

The concept and applications of Transforms and Convolution is very important in Mathematics and Computer Algebra. I have suggested the implementation of modules on the SymPy Group. One important application of this module will be, that it will provide a transform (Number Theoretic Transform) that will be used for the optimal Evaluation of Linear Homogeneous Recurrences with Constant Coefficients (LHRCC) in special prime fields, that I have proposed on the SymPy Group. For generic computation, a useful Divide and Conquer method that uses Cayley Hamilton Theorem will be used to comprehend the evaluation at large values, which will solve this problem more efficiently than direct matrix exponentiation. Another application of this method, would be to reduce large powers of square matrices (not necessarily diagonalizable), if their characteristic polynomial is known. This has an open issue on GitHub. Since it is difficult to obtain and represent closed form expressions of high order recurrences (of specified type), this approach will enable evaluation of these recurrences (even with high order), at very large values (possibly modulo some number, as the values grow exponentially). Reduction of powers of square matrices, subject to their characteristic polynomial will help in representing them as minimal self-polynomial without having to require polynomial division or polynomial modulo algorithms (which have a large constant factor).

The following major modules are to be implemented:

• Transforms

  • Fast Fourier Transform (FFT)
  • Number Theoretic Transform (NTT)
  • Fast Walsh Hadamard Transform (FWHT)
  • Fast Zeta Transform
  • Fast Mobius Transform

• Convolution

  • Linear Convolution
  • Arithmetic Convolution (Bitwise AND, OR, XOR)
  • Subset Convolution (Bitmasks)

• Numerical Evaluation of LHRCC

  • using Cayley Hamilton Theorem (generic method)
  • using Cayley Hamilton Theorem & Number Theoretic Transform (for special prime fields)

Theory

• Transforms

  

  • Fast Fourier Transform

    

  • Number Theoretic Transform

    

  • Fast Walsh Hadamard Transform

    

  • Zeta Transform for sets S'' ⊆ S' ⊆ S or equivalently bitmasks j ⊆ i ⊆ n, denoted by capped f.

    

  • Mobius Transform for sets S'' ⊆ S' ⊆ S or equivalently bitmasks j ⊆ i ⊆ n, denoted by checked f.

    

• Convolution of two list of numbers a and b, denoted by ∗

  • Linear Convolution

    

  • Arithmetic Convolution for op = ^ (XOR), & (AND), | (OR)

    

  • Subset Convolution

    

• Evaluation of LHRCC

  • Cayley Hamilton Theorem (CHT) & Generating Function (GF)

    

Implementation

Functions written in C++ are used to illustrate the main idea for the implementation.

• Transforms for N-length

  • FFT, NTT, FWHT will be implemented using Divide and Conquer Algorithmic Paradigm - O(N lg(N)).

  These transformations will transform a polynomial in coefficient form to point-value form. The following functions illustrate the iterative implementation for the same.

  void fft(vector< complex<double> >& a, bool op){
      /* op = true for inverse fft */
  
      int n = a.size();
      // n <= N, both being powers of 2, N is required for precomputation of roots
      for(int i = 0; i < n; ++i){
          int j = 0, x = i, y = n-1;
          while(y > 0){
              j = j << 1 | (x & 1);
              x >>= 1;
              y >>= 1;
          } // reverse bitmask
          if(i < j) swap(a[i], a[j]);
      }
      for(int h = 2; h <= n; h <<= 1){
          int ut = op? N-N/h: N/h;
          int hf = h >> 1; cx u, v;
          for(int i = 0; i < n; i += h){
              for(int j = i, k = 0; j < i+hf; ++j){
                  u = a[j]; v = w[k] * a[j+hf]; // w stores roots of unity in complex field
                  a[j] = u + v; a[j+hf] = u - v; // butterfly operation in complex field
                  k += ut; k -= (k >= N)*N;
              }
          }
      }
  }
  
  void ntt(vector<int>& a, int md, bool op);
  /* op = true for inverse ntt, md is the prime mod required for ntt */
  /* implementation of ntt is similar to fft (all operations in prime field) */
  /* precalculation of roots involves finding the generator for the prime field */
  
  void fwht(vector<int>& a, bool op){
      /* op = true for inverse fwht */
  
      int n = a.size();
      for(int h = 2; h <= n; h <<= 1){
          int hf = h >> 1; int u, v;
          for(int i = 0; i < n; i += h)
              for(int j = i; j < i+hf; ++j){
                  u = a[j], v = a[j+hf];
  
                  /* butterfly operation */
              }
      }
  }

  • Fast Zeta & Mobius Transform will be implemented using Dynamic Programming Paradigm - O(N lg(N)).
    This is implemented using a technique called Yate's DP (Dynamic Programming).

  void yate(int dp[], int N, int sgn){ // sgn = +1 for zeta, -1 for mobius
    for(int i = 0; i < lg(N); ++i){
      for(int mask = 0; mask < N; ++mask){
        if(mask >> i & 1) continue;
        smask = mask ^ (1<<i); // consider only supermasks of mask
        dp[mask] += sgn*dp[smask]; // sum over supersets
      }
    }
  }
  // for sum over subsets the if condition is negated,
  // hence only submasks of mask are considered

• Convolution for N-length

  • Linear Convolution for complex/prime field will be implemented using FFT/NTT respectively - O(N lg(N)).

  conv(a, b) = ifft(fft(a) fft(b)) (juxtaposition denotes pointwise-multiplication)
  

  • XOR Convolution will be implemented using FWHT - O(N lg(N)).

  conv_xor(a, b) = ifwht(fwht(a) fwht(b)) (juxtaposition denotes pointwise-multiplication)
  

  • AND/OR Convolution will be implemented using Fast Zeta Transform and Fast Mobius Transform - O(N lg(N)).

  fmt(fzt(a) fzt(b)) (juxtaposition denotes pointwise-multiplication)
  will give and/or convolution respectively for superset/subset variant of the zeta, mobius transforms
  

  • Subset Convolution will be implemented using Submask Iteration - O(N ** lg(3)).

  void conv_subset(int a[], int b[], int c[], int N){
    for(int mask = 0; mask < N; ++mask){
      // traverse submasks of mask in lexicographic order
      for(int smask = mask; ; smask = (smask-1) & mask){
        c[mask] += a[smask] * b[mask^smask];
        if(smask == 0) break;
      }
    }
  }
  // Complexity of the subroutine (for N = 2**k) is
  // sum {0 <= i <= k} C(k, i) * 2**i = (1 + 2)**k = 3**k
  //                                  = 3**lg(N) = N**lg(3)

• Evaluation of LHRCC of order k at n

  

  int k; // order of the recurrence
  int md; // for modular reduction
  vector<int> b; // base cases for recurrence
  vector<int> c; // coefficients of recurrence
  vector<int> rec(long long n){
      /* returns the vector of coefficients for representing A^n */
  
      if(n<k){
        /* handle base case */
      }else{
        vector<int> h = rec(n >> 1);
        vector<int> s(k << 1 | (n & 1), 0); // size 2k if n is even, else 2k + 1
        int *u = &s[n&1]; // multiply by A if n is odd
  
        /* minimal polynomial for (A**[n/2])**2 */
  
        /* use CHT for reduction from 2k to k */
  
        s.resize(k);
        return s;
      }
  }

The definition of utility methods in Python is given below,

def conv(a, b, type='linear', mod=None):
    '''
    a, b  :  list of coefficients
    type  :  convolution type in set('linear', 'xor', 'and', 'or', 'subset')
    mod   :  for modular reduction of coefficients of convolution

    for type == 'linear'
      let n be the convolution size (power of 2),
      then if mod is prime and n | (mod - 1) use ntt, else use fft
    '''

def transform(a, type='fourier', mod=None):
    '''
    a     :  list of coefficients
    type  :  transform type in set('fourier', 'walsh_hadamard', 'mobius_subset',
            'mobius_superset', 'zeta_subset', 'zeta_superset')
    mod   :  for modular reduction of coefficients of transform

    for type == 'fourier'
      let n be the transform size (power of 2),
      then if mod is prime and n | (mod - 1) use ntt, else use fft
    '''

def eval_lhrcc(b, c, n, mod=None):
    '''
    lhrcc :  linear homogeneous recurrence with constant coefficients (say h)
    b     :  base cases for recurrence
    c     :  list of coefficients of recurrence
    n     :  evaluation of recurrence at n
    mod   :  modular reduction of value of recurrence at n

    let k = len(b) = len(c) = order of h, then
    h(n) = b[n]                                 ; for 0 <= n < k
    h(n) = sum{0 <= i < k} c[i] * h(n - 1 - i)  ; for n >= k
    '''

The API and usage for SymPy is given below,

• conv

In []: a = [1, 2, 3] # 1x^0 + 2x^1 + 3x^2
In []: b = [4, 5, 6] # 4x^0 + 5x^1 + 6x^2
In []: conv(a, b, type='linear')
Out[]: [4, 13, 28, 27, 18]

The same can be extended for symbolic functions,

In []: conv(1 + 2*x + 3*x**2, 4 + 5*x + 6*x**2, type='linear')
Out[]: 4 + 13*x + 28*x**2 + 27*x**3 + 18*x**4

• transform

In []: a = [1, 2, 3, 4] # indices are 00, 01, 10, 11, in binary
In []: transform(a, type='zeta_subset') # sum over subsets
Out[]: [1, 3, 4, 10] # a[00], a[01]+a[00], a[10]+a[00], a[11]+a[10]+a[01]+a[00]

• eval_lhrcc

In []: b = [1, 1] # f(0) = 1, f(1) = 1
In []: c = [1, 1] # f(n) = 1 * f(n-1) + 1 * f(n-2)
In []: eval_lhrcc(b, c, n=20) # 20th fibonacci number
Out[]: 10946

Timeline

I would be able to work about 40 hours for SymPy per week on an average.

I will devote around 4-5 hours on weekdays and 7-8 hours on weekends. My college reopens in the last week of July, but that would not limit the number of hours I can work, as the beginning of the semester would not involve a hectic schedule. The timeline might seem a bit non-uniform. This is because I have allotted myself sparse work for the period, where I might have other commitments.

Community Bonding Period (April 23 - May 15)

I would spend this time studying the existing modules thoroughly and comprehending the modules required to be implemented. I would also take this time for interacting with mentors and get concrete details of expected work to be completed. I will be having my end semester examination for the first two weeks of May, which will affect my pace slightly during the latter period.

Week 1, 2, 3 (May 16 - June 4)
Implement the modules:

• FFT (Fast Fourier Transform)
• NTT (Number Theoretic Transform)
• FWHT (Fast Walsh Hadamard Transform)

and write basic tests for the same.

Week 4, 5 (June 5 - June 15)
Implement the modules:

• Linear Convolution
• Arithmetic Convolution (XOR)

by extending the transforms (FFT, NTT, FWHT), and write basic tests for the same.

End of Week 5 - Evaluation I

Week 6, 7 (June 16 - June 29)
Write the module for

• Evaluation of Linear Homogeneous Recurrences with Constant Coefficients (General-Purpose Method)

and test evaluation for complex recurrences.

Week 8, 9 (June 30 - July 13)
Implement the modules:

• Fast Zeta Transform
• Fast Mobius Transform

and write basic tests for them.

End of Week 9 - Evaluation II

Week 10, 11 (July 14 - July 30)
Implement the modules for

• Arithmetic Convolution (AND, OR) using Fast Zeta/Mobius Transforms, with appropriate tests.
• Subset Convolution using Submask Traversal, with proper tests.

Week 12 (August 1 - August 6)

• Write proper documentation and special tests for the modules, apart from including basic tests during coding.
• Complete any unfinished work.
• Buffer period.

Final Evaluation Period (August 7 - August 14)

Post GSoC

SymPy is pretty close to my interests and academics, I plan to actively maintain and review my code, and keep contributing to SymPy.
I plan to continue discussing and working on the following areas post GSoC:

• Ability to find coefficients of Rational Generating Functions by reduction to LHRCC.
• Optimizing the evaluation of LHRCC (for specialised prime fields with NTT), as the generic method will fulfill the needs for practical usage without having any restrictions on the modulo.
• Subset Convolution can be further sped upto O(Nlg^2N) using Fast Ranked Mobius/Zeta Transform, for better performance.

References

[1] SymPy Mailing List: Convolution and Transformation Modules

[2] SymPy Mailing List: Numerical Evaluation of Linear Recurrences

[3] SymPy Mailing List: Discussion on coefficient enumeration for Rational Generating Functions

[4] Issue #14266: Issue about Cayley Hamilton Theorem

[5] FFT - engineering.purdue.edu: Fast Fourier Transform & Convolution

[6] Iterative FFT Explanation - staff.ustc.edu.cn: Polynomials and the FFT

[7] NTT - apfloat.org: Number Theoretic Transform

[8] FWHT - wikiwand.com: Fast Walsh Hadamard Transform

[9] Yate's DP - home.iitk.ac.in: Fast Mobius & Zeta Transform, Subset Convolution