code: def sample(N, D):
      
      scale = 31
      result = np.zeros([N, D])
      
      if D > len(directions) + 1:
          raise ValueError("Error in Sobol sequence: not enough dimensions")
  
      L = int(math.ceil(math.log(N) / math.log(2)))
  
      if L > scale:
          raise ValueError("Error in Sobol sequence: not enough bits")
  
      for i in range(D):
          V = np.zeros(L + 1, dtype=long)
  
          if i == 0:
              for j in range(1, L + 1):
                  V[j] = 1 << (scale - j)  # all m's = 1
          else:
              m = np.array(directions[i - 1], dtype=int)
              a = m[0]
              s = len(m) - 1
  
              # The following code discards the first row of the ``m`` array
              # Because it has floating point errors, e.g. values of 2.24e-314
              if L <= s:
                  for j in range(1, L + 1):
                      V[j] = m[j] << (scale - j)
              else:
                  for j in range(1, s + 1):
                      V[j] = m[j] << (scale - j)
                  for j in range(s + 1, L + 1):
                      V[j] = V[j - s] ^ (V[j - s] >> s)
                      for k in range(1, s):
                          V[j] ^= ((a >> (s - 1 - k)) & 1) * V[j - k]
  
          X = long(0)
          for j in range(1, N):
              X ^= V[index_of_least_significant_zero_bit(j - 1)]
              result[j][i] = float(X / math.pow(2, scale))
  
      return result
docs: Generate (N x D) numpy array of Sobol sequence samples