BUG: special: ellip(k,e)inc nan and double expected value

Both ellipkinc and ellipeinc returned nan for a particular value
due to cancelation in a tangent update formula.  Double the
expected value was returned in a small region around the value
returning nan due to using a fragile argument reduction formula.

Fixes gh-3550

This does seem to slow the functions down some.

New
```
In [4]: m = np.abs(np.random.random(2e3))

In [5]: phi = (np.random.random(2e3)+1)*np.pi

In [6]: timeit s.ellipkinc(phi, m[:,None])
1 loops, best of 3: 3.24 s per loop

In [7]: timeit s.ellipeinc(phi, m[:,None])
1 loops, best of 3: 4.5 s per loop
```

Old (using the same `m` and `phi` as above)
```
In [8]: timeit s.ellipkinc(phi, m[:,None])
1 loops, best of 3: 2.34 s per loop

In [9]: timeit s.ellipeinc(phi, m[:,None])
1 loops, best of 3: 3.56 s per loop
```

But correct is better than fast. If anyone has a better replacement for the `mod` variable update in ellik.c and ellie.c, I'll happily change it.

Is it enough to use the "new" `mod` update formula only in the `<= MACHEP` branch?