BUG: truncnorm.cdf returns incorrect values at tail

[fehiepsi]

Describe your issue.

The truncated normal distribution cdf method returns incorrect values when truncated boundaries lie at tail.

Reproducing Code Example

import scipy.stats as ss

ss.truncnorm(13., 15., 0., 1.).cdf(14.)

### Error message

```shell
Expected `1.` but got `0.`

SciPy/NumPy/Python version information

1.7.3 1.21.4 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

[WarrenWeckesser]

The problem is in the master branch, too:

In [23]: import scipy

In [24]: scipy.__version__
Out[24]: '1.8.0.dev0+2101.c255b7d'

In [25]: from scipy.stats import truncnorm

In [26]: truncnorm.cdf(14, 13, 15)
Out[26]: 0.0

I just added an implementation of several methods for the truncated standard normal distribution to mpsci, which can be used to compute the expected value:

In [27]: from mpsci.distributions import truncnorm as truncnorm_mp

In [28]: float(truncnorm_mp.cdf(14, 13, 15))
Out[28]: 0.9999987259565644

I'm working on an overhaul of truncnorm (originally in #13164) that should fix this issue, and fix the slowdowns that have been reported elsewhere.

[tirthasheshpatel]

gh-14753 is a similar report that proposes to use minimax tiling to resolve the inaccurate tails issue of truncnorm.