<h3 dir="auto">Describe your issue.</h3>
<p dir="auto">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.</p>
<p dir="auto">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.</p>
<p dir="auto">If one does not have Mathematica installed, one can go to <a href="https://www.wolframalpha.com/" rel="nofollow">https://www.wolframalpha.com/</a>. There one can generate the output for one of the examples below by entering<br>
<code class="notranslate">z =  (5+10*Pi*I); N[Hypergeometric1F1[1/3,1+1/3,z], 10]</code><br>
in the input cell.</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&quot;&quot;&quot;

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.
&quot;&quot;&quot;

# 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)&lt;tol, msg

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


&quot;&quot;&quot;

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:
&quot;&quot;&quot;

# 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))&lt;tol , msg


&quot;&quot;&quot;
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]
&quot;&quot;&quot;

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

assert abs(result1-result1_mathematica)&lt;tol, 'result1 not eqal result1_mathematica'
assert abs(result2-result2_mathematica)&lt;tol, 'result2 not eqal result2_mathematica'"><pre class="notranslate"><span class="pl-s">"""</span>
<span class="pl-s"></span>
<span class="pl-s">EXAMPLE NUMBER 1:</span>
<span class="pl-s"></span>
<span class="pl-s"></span>
<span class="pl-s"></span>
<span class="pl-s">Code Mathematica:</span>
<span class="pl-s"></span>
<span class="pl-s">s =  (-4.559190954155+51.659216953928I);</span>
<span class="pl-s">sconj=Conjugate[s];</span>
<span class="pl-s">Hypergeometric1F1[1/2,3/2,-s]</span>
<span class="pl-s">Hypergeometric1F1[1/2,3/2,-sconj]</span>
<span class="pl-s"></span>
<span class="pl-s"></span>
<span class="pl-s">Output Mathematica:</span>
<span class="pl-s">Out[12]= 1.\[VeryThinSpace]+1.99922*^-11 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">Out[13]= 1.\[VeryThinSpace]-1.99922*^-11 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">For the given value of s we see that Mathematica computes that the function </span>
<span class="pl-s">value equals 1 for all practical purposes.</span>
<span class="pl-s">"""</span>

<span class="pl-c"># Code Python:</span>
<span class="pl-k">import</span> <span class="pl-s1">scipy</span>.<span class="pl-s1">special</span> <span class="pl-k">as</span> <span class="pl-s1">sc</span>
<span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
 
<span class="pl-c"># Parameters, arguments etc:</span>
<span class="pl-s1">p</span><span class="pl-c1">=</span><span class="pl-c1">2</span> 
<span class="pl-s1">s</span> <span class="pl-c1">=</span> <span class="pl-c1">-</span><span class="pl-c1">4.559190954155</span> <span class="pl-c1">+</span> <span class="pl-c1">51.659216953928</span><span class="pl-c1">*</span><span class="pl-c1">1j</span>

<span class="pl-c"># Evaluation:</span>
<span class="pl-s1">result</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-s1">p</span>, <span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-s1">p</span> <span class="pl-c1">+</span> <span class="pl-c1">1</span>, <span class="pl-c1">-</span><span class="pl-s1">s</span>) 
<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">msg</span> <span class="pl-c1">=</span> <span class="pl-s">f'The above  output of <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result</span><span class="pl-kos">}</span></span> from scipy is far from 1, i.e. in large disagreement with Mathematica'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result</span> <span class="pl-c1">-</span> <span class="pl-c1">1</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s1">msg</span>

<span class="pl-c"># print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')</span>


<span class="pl-s">"""</span>
<span class="pl-s"></span>
<span class="pl-s">EXAMPLE NUMBER 2:</span>
<span class="pl-s"></span>
<span class="pl-s">Because the parameters 1/p and 1/p+1 are real, we know that if we replace </span>
<span class="pl-s">s by its conjugate, the output should be the conjugate of the first output. </span>
<span class="pl-s">This turns out not to be the case:</span>
<span class="pl-s">"""</span>

<span class="pl-c"># z complex:</span>
<span class="pl-s1">z</span> <span class="pl-c1">=</span> <span class="pl-c1">5</span><span class="pl-c1">+</span><span class="pl-c1">10</span><span class="pl-c1">*</span><span class="pl-s1">np</span>.<span class="pl-s1">pi</span><span class="pl-c1">*</span><span class="pl-c1">1j</span>
<span class="pl-s1">result1</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-c1">1</span><span class="pl-c1">+</span><span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-s1">z</span>)
<span class="pl-s1">result2</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-c1">1</span><span class="pl-c1">+</span><span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-s1">np</span>.<span class="pl-en">conj</span>(<span class="pl-s1">z</span>))
<span class="pl-en">print</span>(<span class="pl-s">f'sc.hyp1f1(1/3, 1+1/3, z) = <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result1</span><span class="pl-kos">}</span></span>'</span>)
<span class="pl-en">print</span>(<span class="pl-s">f'sc.hyp1f1(1/3, 1+1/3, np.conj(z)) = <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result2</span><span class="pl-kos">}</span></span>'</span>)
<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">msg</span> <span class="pl-c1">=</span> <span class="pl-s">f'np.conj(result2) does not equal result1 as it should'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result1</span><span class="pl-c1">-</span><span class="pl-s1">np</span>.<span class="pl-en">conj</span>(<span class="pl-s1">result2</span>))<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span> , <span class="pl-s1">msg</span>


<span class="pl-s">"""</span>
<span class="pl-s">Code Mathematica:</span>
<span class="pl-s"></span>
<span class="pl-s">z =  (5+10*Pi*I);</span>
<span class="pl-s">zconj=Conjugate[z];</span>
<span class="pl-s">N[Hypergeometric1F1[1/3,1+1/3,z], 10]</span>
<span class="pl-s">N[Hypergeometric1F1[1/3,1+1/3,zconj], 10]</span>
<span class="pl-s"></span>
<span class="pl-s">Output Mathematica:</span>
<span class="pl-s">Out[15]= 0.449118650-1.390784314 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">Out[16]= 0.449118650+1.390784314 \[ImaginaryI]</span>
<span class="pl-s">"""</span>

<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">result1_mathematica</span> <span class="pl-c1">=</span> <span class="pl-c1">0.449118650</span><span class="pl-c1">-</span><span class="pl-c1">1.390784314j</span>
<span class="pl-s1">result2_mathematica</span> <span class="pl-c1">=</span> <span class="pl-c1">0.449118650</span><span class="pl-c1">+</span><span class="pl-c1">-</span><span class="pl-c1">1.390784314j</span>

<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result1</span><span class="pl-c1">-</span><span class="pl-s1">result1_mathematica</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s">'result1 not eqal result1_mathematica'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result2</span><span class="pl-c1">-</span><span class="pl-s1">result2_mathematica</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s">'result2 not eqal result2_mathematica'</span></pre></div>
<h3 dir="auto">Reproducing Code Example</h3>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&quot;&quot;&quot;
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.
&quot;&quot;&quot;

# 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)&lt;tol, msg

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


&quot;&quot;&quot;
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:
&quot;&quot;&quot;

# 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))&lt;tol , msg


&quot;&quot;&quot;
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]
&quot;&quot;&quot;

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

assert abs(result1-result1_mathematica)&lt;tol, 'result1 not eqal result1_mathematica'
assert abs(result2-result2_mathematica)&lt;tol, 'result2 not eqal result2_mathematica'"><pre class="notranslate"><span class="pl-s">"""</span>
<span class="pl-s">EXAMPLE NUMBER 1:</span>
<span class="pl-s"></span>
<span class="pl-s">Code Mathematica:</span>
<span class="pl-s"></span>
<span class="pl-s">s =  (-4.559190954155+51.659216953928I);</span>
<span class="pl-s">sconj=Conjugate[s];</span>
<span class="pl-s">Hypergeometric1F1[1/2,3/2,-s]</span>
<span class="pl-s">Hypergeometric1F1[1/2,3/2,-sconj]</span>
<span class="pl-s"></span>
<span class="pl-s">Output Mathematica:</span>
<span class="pl-s">Out[12]= 1.\[VeryThinSpace]+1.99922*^-11 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">Out[13]= 1.\[VeryThinSpace]-1.99922*^-11 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">For the given value of s we see that Mathematica computes that the function </span>
<span class="pl-s">value equals 1 for all practical purposes.</span>
<span class="pl-s">"""</span>

<span class="pl-c"># Code Python:</span>
<span class="pl-k">import</span> <span class="pl-s1">scipy</span>.<span class="pl-s1">special</span> <span class="pl-k">as</span> <span class="pl-s1">sc</span>
<span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
 
<span class="pl-c"># Parameters, arguments etc:</span>
<span class="pl-s1">p</span><span class="pl-c1">=</span><span class="pl-c1">2</span> 
<span class="pl-s1">s</span> <span class="pl-c1">=</span> <span class="pl-c1">-</span><span class="pl-c1">4.559190954155</span> <span class="pl-c1">+</span> <span class="pl-c1">51.659216953928</span><span class="pl-c1">*</span><span class="pl-c1">1j</span>

<span class="pl-c"># Evaluation:</span>
<span class="pl-s1">result</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-s1">p</span>, <span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-s1">p</span> <span class="pl-c1">+</span> <span class="pl-c1">1</span>, <span class="pl-c1">-</span><span class="pl-s1">s</span>) 
<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">msg</span> <span class="pl-c1">=</span> <span class="pl-s">f'The above  output of <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result</span><span class="pl-kos">}</span></span> from scipy is far from 1, i.e. in large disagreement with Mathematica'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result</span> <span class="pl-c1">-</span> <span class="pl-c1">1</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s1">msg</span>

<span class="pl-c"># print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')</span>


<span class="pl-s">"""</span>
<span class="pl-s">EXAMPLE NUMBER 2:</span>
<span class="pl-s"></span>
<span class="pl-s">Because the parameters 1/p and 1/p+1 are real, we know that if we replace </span>
<span class="pl-s">s by its conjugate, the output should be the conjugate of the first output. </span>
<span class="pl-s">This turns out not to be the case:</span>
<span class="pl-s">"""</span>

<span class="pl-c"># z complex:</span>
<span class="pl-s1">z</span> <span class="pl-c1">=</span> <span class="pl-c1">5</span><span class="pl-c1">+</span><span class="pl-c1">10</span><span class="pl-c1">*</span><span class="pl-s1">np</span>.<span class="pl-s1">pi</span><span class="pl-c1">*</span><span class="pl-c1">1j</span>
<span class="pl-s1">result1</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-c1">1</span><span class="pl-c1">+</span><span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-s1">z</span>)
<span class="pl-s1">result2</span> <span class="pl-c1">=</span> <span class="pl-s1">sc</span>.<span class="pl-en">hyp1f1</span>(<span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-c1">1</span><span class="pl-c1">+</span><span class="pl-c1">1</span><span class="pl-c1">/</span><span class="pl-c1">3</span>, <span class="pl-s1">np</span>.<span class="pl-en">conj</span>(<span class="pl-s1">z</span>))
<span class="pl-en">print</span>(<span class="pl-s">f'sc.hyp1f1(1/3, 1+1/3, z) = <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result1</span><span class="pl-kos">}</span></span>'</span>)
<span class="pl-en">print</span>(<span class="pl-s">f'sc.hyp1f1(1/3, 1+1/3, np.conj(z)) = <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">result2</span><span class="pl-kos">}</span></span>'</span>)
<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">msg</span> <span class="pl-c1">=</span> <span class="pl-s">f'np.conj(result2) does not equal result1 as it should'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result1</span><span class="pl-c1">-</span><span class="pl-s1">np</span>.<span class="pl-en">conj</span>(<span class="pl-s1">result2</span>))<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span> , <span class="pl-s1">msg</span>


<span class="pl-s">"""</span>
<span class="pl-s">Code Mathematica:</span>
<span class="pl-s"></span>
<span class="pl-s">z =  (5+10*Pi*I);</span>
<span class="pl-s">zconj=Conjugate[z];</span>
<span class="pl-s">N[Hypergeometric1F1[1/3,1+1/3,z], 10]</span>
<span class="pl-s">N[Hypergeometric1F1[1/3,1+1/3,zconj], 10]</span>
<span class="pl-s"></span>
<span class="pl-s">Output Mathematica:</span>
<span class="pl-s">Out[15]= 0.449118650-1.390784314 \[ImaginaryI]</span>
<span class="pl-s"></span>
<span class="pl-s">Out[16]= 0.449118650+1.390784314 \[ImaginaryI]</span>
<span class="pl-s">"""</span>

<span class="pl-s1">tol</span> <span class="pl-c1">=</span> <span class="pl-c1">1e-6</span>
<span class="pl-s1">result1_mathematica</span> <span class="pl-c1">=</span> <span class="pl-c1">0.449118650</span><span class="pl-c1">-</span><span class="pl-c1">1.390784314j</span>
<span class="pl-s1">result2_mathematica</span> <span class="pl-c1">=</span> <span class="pl-c1">0.449118650</span><span class="pl-c1">+</span><span class="pl-c1">-</span><span class="pl-c1">1.390784314j</span>

<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result1</span><span class="pl-c1">-</span><span class="pl-s1">result1_mathematica</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s">'result1 not eqal result1_mathematica'</span>
<span class="pl-k">assert</span> <span class="pl-en">abs</span>(<span class="pl-s1">result2</span><span class="pl-c1">-</span><span class="pl-s1">result2_mathematica</span>)<span class="pl-c1">&lt;</span><span class="pl-s1">tol</span>, <span class="pl-s">'result2 not eqal result2_mathematica'</span></pre></div>
<h3 dir="auto">Error message</h3>
<div class="highlight highlight-source-shell notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
&lt;ipython-input-47-f81c740ef128&gt; in &lt;module&gt;
     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'
---&gt; 42 assert abs(result - 1)&lt;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"><pre class="notranslate">---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<span class="pl-k">&lt;</span>ipython-input-47-f81c740ef<span class="pl-k">128&gt;</span> <span class="pl-k">in</span> <span class="pl-k">&lt;</span>module<span class="pl-k">&gt;</span>
     40 tol = 1e-6
     41 msg = f<span class="pl-s"><span class="pl-pds">'</span>The above  output of {result} from scipy is far from 1, i.e. in large disagreement with Mathematica<span class="pl-pds">'</span></span>
---<span class="pl-k">&gt;</span> 42 assert abs(result - 1)<span class="pl-k">&lt;</span>tol, msg
     43 
     44 <span class="pl-c"><span class="pl-c">#</span> print(f'sc.hyp1f1(1/p, 1/p + 1, -s) = {result}')</span>

AssertionError: The above  output of (0.833788272795157+0.18152681828629416j) from scipy is far from 1, i.e. <span class="pl-k">in</span> large disagreement with Mathematica</pre></div>
<h3 dir="auto">SciPy/NumPy/Python version information</h3>
<p dir="auto">import sys, scipy, numpy; print(scipy.<strong>version</strong>, numpy.<strong>version</strong>, sys.version_info)</p>
<p dir="auto">I am using the <code class="notranslate">quad</code> function from <code class="notranslate">scipy.integrate v0.19.1</code> to integrate functions with a square-root like singularity at each end of the integration interval such as for example</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1)"><pre class="notranslate"><code class="notranslate">In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1)
</code></pre></div>
<p dir="auto">(I use the <code class="notranslate">sqrt</code> function from <code class="notranslate">numpy v1.12.0</code>) which immediately yields the correct result pi:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="Out[1]: (3.141592653589591, 6.200897573194197e-10)"><pre class="notranslate"><code class="notranslate">Out[1]: (3.141592653589591, 6.200897573194197e-10)
</code></pre></div>
<p dir="auto">According to the documentation of the <code class="notranslate">quad</code> function the keyword <code class="notranslate">points</code> should be used to indicate the locations of singularities or discontinuities of the integrand, but if I indicate the points <code class="notranslate">[1, -1]</code> where the above integrand is singluar I get a Warning and <code class="notranslate">nan</code> as result:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="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)"><pre class="notranslate"><code class="notranslate">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)
</code></pre></div>
<p dir="auto">Why does <code class="notranslate">quad</code> produce these problems if the singularities of the integrand are specified and just runs fine if these points are not indicated?</p>
<p dir="auto"><strong>EDIT:</strong> 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</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2))"><pre class="notranslate"><code class="notranslate">quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2))
</code></pre></div>
<p dir="auto">and everything is is fine. The above problem with the keyword <code class="notranslate">points</code> is probably due to the fact that a different QUADPACK routine is used if <code class="notranslate">points</code> are specified.</p>