Issue 1:
Describe your issue.

There is a serious bug in the implementation of Kummer's hypergeometric function. If the first two parameters a, b are real, we have that hyp1f1(a, b, z) should be the conjugate of hyp1f1(a, b, np.conj(z)). This is not always the case. Examples are provided below.

Mathematica is one of the industry standards for doing scientific computations. For many arguments hyp1f1(a, b, z) does not agree with Mathematica's output given the same input. Examples are provided below.

If one does not have Mathematica installed, one can go to https://www.wolframalpha.com/. There one can generate the output for one of the examples below by entering

z = (5+10*Pi*I); N[Hypergeometric1F1[1/3,1+1/3,z], 10]

in the input cell.

"""

EXAMPLE NUMBER 1:

Code Mathematica:

s = (-4.559190954155+51.659216953928I);
sconj=Conjugate[s];
Hypergeometric1F1[1/2,3/2,-s]
Hypergeometric1F1[1/2,3/2,-sconj]

Output Mathematica:
Out[12]= 1.\[VeryThinSpace]+1.99922*^-11 \[ImaginaryI]

Out[13]= 1.\[VeryThinSpace]-1.99922*^-11 \[ImaginaryI]

For the given value of s we see that Mathematica computes that the function
value equals 1 for all practical purposes.
"""

# Code Python:
import scipy.special as sc
import numpy as np

# Parameters, arguments etc:
p=2
s = -4.559190954155 + 51.659216953928*1j

# Evaluation:
result = sc.hyp1f1(1/p, 1/p + 1, -s)
tol = 1e-6
msg = f'The above output of {result} from scipy is far from 1, i.e. in large disagreement with Mathematica'
assert abs(result - 1)<tol, msg

# print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')

"""

EXAMPLE NUMBER 2:

Because the parameters 1/p and 1/p+1 are real, we know that if we replace
s by its conjugate, the output should be the conjugate of the first output.
This turns out not to be the case:
"""

# z complex:
z = 5+10*np.pi*1j
result1 = sc.hyp1f1(1/3, 1+1/3, z)
result2 = sc.hyp1f1(1/3, 1+1/3, np.conj(z))
print(f'sc.hyp1f1(1/3, 1+1/3, z) = {result1}')
print(f'sc.hyp1f1(1/3, 1+1/3, np.conj(z)) = {result2}')
tol = 1e-6
msg = f'np.conj(result2) does not equal result1 as it should'
assert abs(result1-np.conj(result2))<tol , msg

"""
Code Mathematica:

z = (5+10*Pi*I);
zconj=Conjugate[z];
N[Hypergeometric1F1[1/3,1+1/3,z], 10]
N[Hypergeometric1F1[1/3,1+1/3,zconj], 10]

Output Mathematica:
Out[15]= 0.449118650-1.390784314 \[ImaginaryI]

Out[16]= 0.449118650+1.390784314 \[ImaginaryI]
"""

tol = 1e-6
result1_mathematica = 0.449118650-1.390784314j
result2_mathematica = 0.449118650+-1.390784314j

assert abs(result1-result1_mathematica)<tol, 'result1 not eqal result1_mathematica'
assert abs(result2-result2_mathematica)<tol, 'result2 not eqal result2_mathematica'

Reproducing Code Example

"""
EXAMPLE NUMBER 1:

Code Mathematica:

s = (-4.559190954155+51.659216953928I);
sconj=Conjugate[s];
Hypergeometric1F1[1/2,3/2,-s]
Hypergeometric1F1[1/2,3/2,-sconj]

Output Mathematica:
Out[12]= 1.\[VeryThinSpace]+1.99922*^-11 \[ImaginaryI]

Out[13]= 1.\[VeryThinSpace]-1.99922*^-11 \[ImaginaryI]

For the given value of s we see that Mathematica computes that the function
value equals 1 for all practical purposes.
"""

# Code Python:
import scipy.special as sc
import numpy as np

# Parameters, arguments etc:
p=2
s = -4.559190954155 + 51.659216953928*1j

# Evaluation:
result = sc.hyp1f1(1/p, 1/p + 1, -s)
tol = 1e-6
msg = f'The above output of {result} from scipy is far from 1, i.e. in large disagreement with Mathematica'
assert abs(result - 1)<tol, msg

# print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')

"""
EXAMPLE NUMBER 2:

Because the parameters 1/p and 1/p+1 are real, we know that if we replace
s by its conjugate, the output should be the conjugate of the first output.
This turns out not to be the case:
"""

# z complex:
z = 5+10*np.pi*1j
result1 = sc.hyp1f1(1/3, 1+1/3, z)
result2 = sc.hyp1f1(1/3, 1+1/3, np.conj(z))
print(f'sc.hyp1f1(1/3, 1+1/3, z) = {result1}')
print(f'sc.hyp1f1(1/3, 1+1/3, np.conj(z)) = {result2}')
tol = 1e-6
msg = f'np.conj(result2) does not equal result1 as it should'
assert abs(result1-np.conj(result2))<tol , msg

"""
Code Mathematica:

z = (5+10*Pi*I);
zconj=Conjugate[z];
N[Hypergeometric1F1[1/3,1+1/3,z], 10]
N[Hypergeometric1F1[1/3,1+1/3,zconj], 10]

Output Mathematica:
Out[15]= 0.449118650-1.390784314 \[ImaginaryI]

Out[16]= 0.449118650+1.390784314 \[ImaginaryI]
"""

tol = 1e-6
result1_mathematica = 0.449118650-1.390784314j
result2_mathematica = 0.449118650+-1.390784314j

assert abs(result1-result1_mathematica)<tol, 'result1 not eqal result1_mathematica'
assert abs(result2-result2_mathematica)<tol, 'result2 not eqal result2_mathematica'

Error message

---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-47-f81c740ef128> in <module>
40 tol = 1e-6
41 msg = f'The above output of {result} from scipy is far from 1, i.e. in large disagreement with Mathematica'
---> 42 assert abs(result - 1)<tol, msg
43
44 # print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')

AssertionError: The above output of (0.833788272795157+0.18152681828629416j) from scipy is far from 1, i.e. in large disagreement with Mathematica

SciPy/NumPy/Python version information

import sys, scipy, numpy; print(scipy.version, numpy.version, sys.version_info)

Issue 2:
I am using the quad function from scipy.integrate v0.19.1 to integrate functions with a square-root like singularity at each end of the integration interval such as for example

In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1)

(I use the sqrt function from numpy v1.12.0) which immediately yields the correct result pi:

Out[1]: (3.141592653589591, 6.200897573194197e-10)

According to the documentation of the quad function the keyword points should be used to indicate the locations of singularities or discontinuities of the integrand, but if I indicate the points [1, -1] where the above integrand is singluar I get a Warning and nan as result:

In [2]: quad(lambda x: 1/sqrt(1-x**2), -1, 1, points=[-1, 1])

RuntimeWarning: divide by zero encountered in double_scalars
IntegrationWarning: Extremely bad integrand behavior occurs at some points of the integration interval.

Out[2]: (nan, nan)

Why does quad produce these problems if the singularities of the integrand are specified and just runs fine if these points are not indicated?

EDIT: I figured out, that the correct way to handle the singularities in the above integral is to use a proper weighting function. For the above integral I can use

quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2))

and everything is is fine. The above problem with the keyword points is probably due to the fact that a different QUADPACK routine is used if points are specified.