I don't think this bit is correct.  It's not a problem to have probability values == 0.  e.g. extend my sawtooth example to a double sawtooth. You don't want to remove the point in the middle.  In fact, I don't think you even need to build the probability array at all.  Just build dcdf and then check that all of those values are nonnegative.  

Here's another example that you should check gives reasonable results, a square wave pattern:
```
sq = galsim.DistDeviate(rng=rng,x=[-2,-1,-0.999999,0,0.999999,1,2],p=[1,1,0,0,0,1,1],npoints=5)
```
The mean and variance should be 0 and 7/3.  And it has the advantage that with npoints=5, the linear interpolation is correct, so we should expect to get the correct answers without a lot of points.  I'm currently getting 
```
>>> vals = [sq() for i in range(1000000)]
>>> numpy.mean(vals)
0.49986963964590431
>>> numpy.var(vals)
1.4179576777624152
```
because it's throwing out points with p=0 or dcdf=0.  It's a tricky one to get right.  You'll want to throw out the x=0 point.  But you need to keep both x=-1 and x=1.  And to make sure the LookupTable doesn't complain about repeated elements, you'll need to bump up the cdf value for x=1 slightly.  I recommend waiting until after normalizing cdf to have cdf[-1]=1, and then add 2.23e-16 to the second one.