I would like to construct a Finite Element basis by using a generalized Vandermonde matrix. The idea is to compute the values of a suitable modal basis ('prime basis') at a set of points in reference space and express the values of Finite Element basis functions as a linear combination of the members of the prime basis as explained in the references  or .
I have the following problem. The primal basis of choice on triangles/tetrahedra is Dubiner basis and uses Jacobi polynomials. When I evaluate this basis, some entries in the Vandermonde matrix are nearly zero ( ~ 10^-17), but it appears that they should be in fact exactly zero. Later on, when I use the Vandermonde matrix to compute the values of Lagrange shape functions, some of the resulting values are again slightly perturbed. The consequence of this is that the sum of all Lagrange shape functions at certain points is not exactly 1, and the sum of derivatives is not exactly 0.
I am wondering if anyone has ever encountered a problem like this and if there's a remedy. I searched in books and on google, but with no success. As I would like to use higher-order shape functions in a CFD code, I'm afraid that the errors in the shape function evaluation will introduce instabilities.
To give a bit of implementation details: I implemented the computation of Dubiner modes in C++ in two different ways:
a) 'direct' implementation using a recursive formula for Jacobi polynomials as explained in  and . This is known to fail at the vertex (-1,1) of the reference triangle due to a singularity in the collapsed coordinate transformation.
b) a 'workaround' implementation described in . In accordance with the observations made in this paper, the errors obtained by b) are slightly smaller than in case a).
 Jan S. Hesthaven, Tim Warburton: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications
 George Em Karniadakis, Spencer J. Sherwin: Spectral/hp Element Methods for Computational Fluid Dynamics
 Robert C. Kirby: Singularity-free evaluation of collapsed-coordinate orthogonal polynomials