Compute number of matrices with appropriate sums When multiplying monomials in the Milnor basis for the Steenrod algebra, part of the algorithm involves enumerating certain "allowable matrices".

Given two lists of nonnegative integers r1, ... ,rm and  s1, ... ,sn, a matrix of nonnegative integers X 



is allowable if


The sum of the jth column is less than or equal to sj:


The sum of the ith row weighted by powers of 2 is less than or equal to ri:




Task

Write a program which takes a pair of lists r1, ... ,rm and s1, s1, ... ,sn and computes the number of allowable matrices for these lists. Your program may optionally take m and n as additional arguments if need be.


These numbers may be input in any format one likes, for instance grouped into lists or encoded in unary, or anything else.
Output should be a positive integer
Standard loopholes apply.


Scoring

This is code golf: Shortest solution in bytes wins.

Examples:

For [2] and [1], there are two allowable matrices:



For [4] and [1,1] there are three allowable matrices:



For [2,4] and [1,1] there are five allowable matrices:



Test cases:

   Input: [1], [2]
   Output: 1

   Input: [2], [1]
   Output: 2

   Input: [4], [1,1]
   Output: 3

   Input: [2,4], [1,1]   
   Output: 5      

   Input: [3,5,7], [1,2]
   Output: 14

   Input: [7, 10], [1, 1, 1]
   Output: 15       

   Input: [3, 6, 16, 33], [0, 1, 1, 1, 1]
   Output: 38      

   Input: [7, 8], [3, 3, 1]
   Output: 44

   Input: [2, 6, 15, 18], [1, 1, 1, 1, 1]
   Output: 90       

   Input: [2, 6, 7, 16], [1, 3, 2]
   Output: 128

   Input: [2, 7, 16], [3, 3, 1, 1]
   Output: 175
JavaScript (ES7), 163 bytes



f=([R,...x],s)=&gt;1/R?[...Array(R**s.length)].reduce((k,_,n)=&gt;(a=s.map((_,i)=&gt;n/R**i%R|0)).some(c=&gt;(p+=c&lt;&lt;++j)&gt;R,p=j=0)?k:k+f(x,s.map((v,i)=&gt;v-a[i])),0):!/-/.test(s)


Test cases

NB: I've removed the two most time-consuming test cases from this snippet, but they should pass as well.



f=([R,...x],s)=&gt;1/R?[...Array(R**s.length)].reduce((k,_,n)=&gt;(a=s.map((_,i)=&gt;n/R**i%R|0)).some(c=&gt;(p+=c&lt;&lt;++j)&gt;R,p=j=0)?k:k+f(x,s.map((v,i)=&gt;v-a[i])),0):!/-/.test(s)

console.log(f([1],[2]))            // 1
console.log(f([2],[1]))            // 2
console.log(f([4],[1,1]))          // 3
console.log(f([2,4],[1,1]   ))     // 5      
console.log(f([3,5,7],[1,2]))      // 14
console.log(f([7,10],[1,1,1]))     // 15       
console.log(f([7,8],[3,3,1]))      // 44
console.log(f([2,6,7,16],[1,3,2])) // 128
console.log(f([2,7,16],[3,3,1,1])) // 175




Commented

f = (                               // f = recursive function taking:
  [R,                               //   - the input array r[] splitted into:
      ...x],                        //     R = next element / x = remaining elements
  s                                 //   - the input array s[]
) =&gt;                                //
  1 / R ?                           // if R is defined:
    [...Array(R**s.length)]         //   for each n in [0, ..., R**s.length - 1],
    .reduce((k, _, n) =&gt;            //   using k as an accumulator:
      (a =                          //     build the next combination a[] of
        s.map((_, i) =&gt;             //     N elements in [0, ..., R - 1]
          n / R**i % R | 0          //     where N is the length of s[]
        )                           //
      ).some(c =&gt;                   //     for each element c in a[]:
        (p += c &lt;&lt; ++j)             //       increment j; add c * (2**j) to p
        &gt; R,                        //       exit with a truthy value if p &gt; R
        p = j = 0                   //       start with p = j = 0
      ) ?                           //     end of some(); if truthy:
        k                           //       just return k unchanged
      :                             //     else:
        k +                         //       add to k the result of
        f(                          //       a recursive call to f() with:
          x,                        //         the remaining elements of r[]
          s.map((v, i) =&gt; v - a[i]) //         s[] updated by subtracting the values of a[]
        ),                          //       end of recursive call
      0                             //     initial value of the accumulator k
    )                               //   end of reduce()
  :                                 // else:
    !/-/.test(s)                    //   return true if there's no negative value in s[]