ENH: Added max_gap keyword for series.interpolate
- [x] closes #12187 and #26796 
- [x] tests added / passed
- [x] passes `git diff upstream/master -u -- "*.py" | flake8 --diff`
- [x] whatsnew entry

This PR introduces the new keyword `max_gap` for `interpolate`. For all NaN-gaps which are wider than `max_gap` no interpolation is carried out and all NaNs of the gap are preserved. This is in contrast to using the `limit` kwarg which does not prevent interpolating values in a longer NaN gap.

I added numpy-based implementation that searches for NaN-gaps wider than `max_gap`. In line with the current implementations for NaN handling in `series.interpolate`, a set of NaN-indices that has
to be preserved is generated. This is used in the end, after a full interpolation of all NaN is done, to restore the NaNs gaps that shall not be interpolated.

Test and documentation were also added.

It will need some small PEP8-cleanup and maybe tests using other interpolation methods then `linear` (edit: Done). But before I continue, I would like to get feedback if my approach is in general okay.

This PR might also be extended to close #16457 which is on interpolation directly after resampling.

Example usage:

```python
In [1]: import numpy as np                                                                                                                                                                                                         

In [2]: import pandas as pd                                                                                                                                                                                                        

In [3]: s = pd.Series([np.nan, 1., np.nan, 2., np.nan, np.nan,  
   ...:                5., np.nan, np.nan, np.nan, -1., np.nan, np.nan])                                                                                                                                                           

# Using the new `max_gap` kwarg
In [4]: s.interpolate(max_gap=2, limit_area='inside')                                                                                                                                                                               
Out[4]: 
0     NaN
1     1.0
2     1.5
3     2.0
4     3.0
5     4.0
6     5.0
7     NaN
8     NaN
9     NaN
10   -1.0
11    NaN
12    NaN
dtype: float64

# Compare to the result when using the existing `limit` kwarg
In [5]: s.interpolate(limit=2, limit_area='inside')                                                                                                                                                                                
Out[5]: 
0     NaN
1     1.0
2     1.5
3     2.0
4     3.0
5     4.0
6     5.0
7     3.5
8     2.0
9     NaN
10   -1.0
11    NaN
12    NaN
dtype: float64
```

Timing:
```python
In [6]: %timeit s.interpolate(max_gap=2, limit_area='inside')                                                                                                                                                                       
708 µs ± 14.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [7]: %timeit s.interpolate(limit=2, limit_area='inside')                                                                                                                                                                        
631 µs ± 5.85 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
```
 
The relative speed difference is similar for larger Series.