Skip to content

Support for hypergeometric functions #1383

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 32 commits into from
Apr 22, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
32 commits
Select commit Hold shift + click to select a range
53c91aa
Hypergeometric1F1
Feb 12, 2025
bbecd08
Black and isort: hypergeom.py
Feb 12, 2025
43d30fa
Apply suggestions from code review
aravindh-krishnamoorthy Feb 12, 2025
2d369c7
Hypergeometric functions based on sympy and mmpmath.
Feb 13, 2025
6eb4e20
Black and isort on hypergeom.py
Feb 13, 2025
52d43fd
Apply suggestions from code review
aravindh-krishnamoorthy Feb 14, 2025
235cb2f
Symbolic by default
Feb 14, 2025
31116cb
Updates
Feb 14, 2025
96405f1
Minor updates
Feb 14, 2025
bbb4c2f
Updates
Feb 14, 2025
7cb353a
Black and iSort
Feb 14, 2025
f89b80c
Remove eval_N_prec for later.
Feb 14, 2025
8c66f2d
Use MeijerG numerical results for HypergeometricU
Feb 14, 2025
1780e70
Remove redundant mpmath_name for HypergeometricU
Feb 14, 2025
e55d74f
Updates to N-evaluation.
Feb 14, 2025
bece7d9
Updates to documentation of HypergeometricU
Feb 14, 2025
18fe6fd
Minor updates to documentation
Feb 14, 2025
ea22016
Minor updates to documentation/2
Feb 14, 2025
2e90531
Merge branch 'master' into hypergeom
aravindh-krishnamoorthy Feb 15, 2025
4a91ca1
Use tracing for sympy and mpmath calls.
Feb 15, 2025
2c563eb
Merge branch 'master' into hypergeom
aravindh-krishnamoorthy Feb 22, 2025
3eb4d18
Numerical evaluation with machine precision `z` as in MMA
Feb 24, 2025
4c646db
Align attributes
aravindh-krishnamoorthy Mar 16, 2025
d1f3775
Fix division by zero returns for numerical evaluations
aravindh-krishnamoorthy Mar 16, 2025
7e417a4
Initial test_hypergeom.py
aravindh-krishnamoorthy Mar 16, 2025
33e26d6
Limit MeijerG expansion to valid cases.
aravindh-krishnamoorthy Mar 16, 2025
c180f4a
Move scalar special case to 1F1
aravindh-krishnamoorthy Mar 16, 2025
313bead
Test cases v2
aravindh-krishnamoorthy Mar 16, 2025
360a029
Add special case for PFQ[b,b,z] + test.
aravindh-krishnamoorthy Mar 16, 2025
a5489f5
Minor update to HypergeometricU doc
aravindh-krishnamoorthy Mar 16, 2025
6bd51d1
Housekeeping: isort and black
aravindh-krishnamoorthy Mar 16, 2025
ea70742
Update mathics/builtin/specialfns/hypergeom.py
aravindh-krishnamoorthy Apr 22, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 0 additions & 33 deletions mathics/builtin/specialfns/bessel.py
Original file line number Diff line number Diff line change
Expand Up @@ -566,39 +566,6 @@ class HankelH2(_Bessel):
sympy_name = "hankel2"


class HypergeometricU(MPMathFunction):
"""
<url>
:Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
:mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
:WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
<dl>
<dt>'HypergeometricU'[$a$, $b$, $z$]
<dd>returns $U(a, b, z)$.
</dl>

>> HypergeometricU[3, 2, 1.]
= 0.105479

Plot 'U'[3, 2, x] from 0 to 2 in steps of 0.5:
>> Plot[HypergeometricU[3, 2, x], {x, 0.5, 2}]
= -Graphics-

We handle this special case:
>> HypergeometricU[0, c, z]
= 1
"""

attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
mpmath_name = "hyperu"
nargs = {3}
rules = {
"HypergeometricU[0, c_, z_]": "1",
}
summary_text = "Tricomi confluent hypergeometric function"
sympy_name = ""


# Kelvin Functions


Expand Down
231 changes: 231 additions & 0 deletions mathics/builtin/specialfns/hypergeom.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,231 @@
"""
Hypergeometric functions

See also <url>
:Chapter 15 Hypergeometric Functions in the Digital Library of Mathematical Functions:
https://dlmf.nist.gov/15</url>.
"""

import mpmath
import sympy

import mathics.eval.tracing as tracing
from mathics.core.attributes import (
A_LISTABLE,
A_N_HOLD_FIRST,
A_NUMERIC_FUNCTION,
A_PROTECTED,
A_READ_PROTECTED,
)
from mathics.core.builtin import MPMathFunction
from mathics.core.convert.mpmath import from_mpmath
from mathics.core.convert.sympy import from_sympy
from mathics.core.evaluation import Evaluation
from mathics.core.number import FP_MANTISA_BINARY_DIGITS
from mathics.core.systemsymbols import SymbolComplexInfinity, SymbolMachinePrecision


class HypergeometricPFQ(MPMathFunction):
"""
<url>
:Generalized hypergeometric function: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function</url> (<url>
:mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
:sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.hyper</url>, <url>
:WMA: https://reference.wolfram.com/language/ref/HypergeometricPFQ.html</url>)
<dl>
<dt>'HypergeometricPFQ'[${a_1, ..., a_p}, {b_1, ..., b_q}, z$]
<dd>returns ${}_p F_q({a_1, ..., a_p}; {b_1, ..., b_q}; z)$.
</dl>
>> HypergeometricPFQ[{2}, {2}, 1]
= E

Result is symbollicaly simplified by default:
>> HypergeometricPFQ[{3}, {2}, 1]
= HypergeometricPFQ[{3}, {2}, 1]
unless a numerical evaluation is explicitly requested:
>> HypergeometricPFQ[{3}, {2}, 1] // N
= 4.07742
>> HypergeometricPFQ[{3}, {2}, 1.]
= 4.07742

The following special cases are handled:
>> HypergeometricPFQ[{}, {}, z]
= 1
>> HypergeometricPFQ[{0}, {b}, z]
= 1
>> Hypergeometric1F1[b, b, z]
= E ^ z
"""

attributes = A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
mpmath_name = "hyper"
nargs = {3}
rules = {
"HypergeometricPFQ[{}, {}, z_]": "1",
"HypergeometricPFQ[{0}, b_, z_]": "1",
"HypergeometricPFQ[b_, b_, z_]": "Exp[z]",
}
summary_text = "compute the generalized hypergeometric function"
sympy_name = "hyper"

def eval(self, a, b, z, evaluation: Evaluation):
"HypergeometricPFQ[a_, b_, z_]"
try:
a_sympy = [e.to_sympy() for e in a]
b_sympy = [e.to_sympy() for e in b]
result_sympy = tracing.run_sympy(
sympy.hyper, a_sympy, b_sympy, z.to_sympy()
)
return from_sympy(result_sympy)
except Exception:
pass

def eval_N(self, a, b, z, prec, evaluation: Evaluation):
"N[HypergeometricPFQ[a_, b_, z_], prec_]"
try:
result_mpmath = tracing.run_mpmath(
mpmath.hyper, a.to_python(), b.to_python(), z.to_python()
)
return from_mpmath(result_mpmath)
except ZeroDivisionError:
return SymbolComplexInfinity
except Exception as ex:
pass

def eval_numeric(self, a, b, z, evaluation: Evaluation):
"HypergeometricPFQ[a:{__?NumericQ}, b:{__?NumericQ}, z_?MachineNumberQ]"
return self.eval_N(a, b, z, SymbolMachinePrecision, evaluation)


class Hypergeometric1F1(MPMathFunction):
"""
<url>
:Kummer confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
:mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#hyper</url>, <url>
:WMA: https://reference.wolfram.com/language/ref/Hypergeometric1F1.html</url>)
<dl>
<dt>'Hypergeometric1F1'[$a$, $b$, $z$]
<dd>returns ${}_1 F_1(a; b; z)$.
</dl>

Result is symbollicaly simplified by default:
>> Hypergeometric1F1[3, 2, 1]
= HypergeometricPFQ[{3}, {2}, 1]
unless a numerical evaluation is explicitly requested:
>> Hypergeometric1F1[3, 2, 1] // N
= 4.07742
>> Hypergeometric1F1[3, 2, 1.]
= 4.07742

Plot 'M'[3, 2, x] from 0 to 2 in steps of 0.5:
>> Plot[Hypergeometric1F1[3, 2, x], {x, 0.5, 2}]
= -Graphics-
Here, plot explicitly requests a numerical evaluation.
"""

attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED
mpmath_name = ""
nargs = {3}
rules = {
"Hypergeometric1F1[a_, b_, z_]": "HypergeometricPFQ[{a},{b},z]",
}
summary_text = "compute Kummer confluent hypergeometric function"
sympy_name = ""


class MeijerG(MPMathFunction):
"""
<url>
:Meijer G-function: https://en.wikipedia.org/wiki/Meijer_G-function</url> (<url>
:mpmath: https://mpmath.org/doc/current/functions/hypergeometric.html#meijerg</url>, <url>
:sympy: https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.hyper.meijerg</url>, <url>
:WMA: https://reference.wolfram.com/language/ref/MeijerG.html</url>)
<dl>
<dt>'MeijerG'[${{a_1, ..., a_n}, {a_{n+1}, ..., a_p}}, {{b_1, ..., b_m}, {b_{m+1}, ..., a_q}}, z$]
<dd>returns $G^{m,n}_{p,q}(z | {a_1, ..., a_p}; {b_1, ..., b_q})$.
</dl>
Result is symbollicaly simplified by default:
>> MeijerG[{{1, 2}, {}}, {{3}, {}}, 1]
= MeijerG[{{1, 2}, {}}, {{3}, {}}, 1]
unless a numerical evaluation is explicitly requested:
>> MeijerG[{{1, 2},{}}, {{3},{}}, 1] // N
= 0.210958
>> MeijerG[{{1, 2},{}}, {{3},{}}, 1.]
= 0.210958
"""

attributes = A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
mpmath_name = "meijerg"
nargs = {3}
rules = {}
summary_text = "compute the Meijer G-function"
sympy_name = "meijerg"

def eval(self, a, b, z, evaluation: Evaluation):
"MeijerG[a_, b_, z_]"
try:
a_sympy = [[e2.to_sympy() for e2 in e1] for e1 in a]
b_sympy = [[e2.to_sympy() for e2 in e1] for e1 in b]
result_sympy = tracing.run_sympy(
sympy.meijerg, a_sympy, b_sympy, z.to_sympy()
)
return from_sympy(result_sympy)
except Exception:
pass

def eval_N(self, a, b, z, prec, evaluation: Evaluation):
"N[MeijerG[a_, b_, z_], prec_]"
try:
result_mpmath = tracing.run_mpmath(
mpmath.meijerg, a.to_python(), b.to_python(), z.to_python()
)
return from_mpmath(result_mpmath)
except ZeroDivisionError:
return SymbolComplexInfinity
except Exception:
pass

def eval_numeric(self, a, b, z, evaluation: Evaluation):
"MeijerG[a:{___List?(AllTrue[#, NumericQ, Infinity]&)}, b:{___List?(AllTrue[#, NumericQ, Infinity]&)}, z_?MachineNumberQ]"
return self.eval_N(a, b, z, SymbolMachinePrecision, evaluation)


class HypergeometricU(MPMathFunction):
"""
<url>
:Confluent hypergeometric function: https://en.wikipedia.org/wiki/Confluent_hypergeometric_function</url> (<url>
:mpmath: https://mpmath.org/doc/current/functions/bessel.html#mpmath.hyperu</url>, <url>
:WMA: https://reference.wolfram.com/language/ref/HypergeometricU.html</url>)
<dl>
<dt>'HypergeometricU'[$a$, $b$, $z$]
<dd>returns $U(a, b, z)$.
</dl>
Result is symbollicaly simplified, where possible:
>> HypergeometricU[3, 2, 1]
= MeijerG[{{-2}, {}}, {{0, -1}, {}}, 1] / 2
>> HypergeometricU[1,4,8]
= HypergeometricU[1, 4, 8]
unless a numerical evaluation is explicitly requested:
>> HypergeometricU[3, 2, 1] // N
= 0.105479
>> HypergeometricU[3, 2, 1.]
= 0.105479

Plot 'U'[3, 2, x] from 0 to 10 in steps of 0.5:
>> Plot[HypergeometricU[3, 2, x], {x, 0.5, 10}]
= -Graphics-

We handle this special case:
>> HypergeometricU[0, b, z]
= 1
"""

attributes = A_LISTABLE | A_NUMERIC_FUNCTION | A_PROTECTED | A_READ_PROTECTED
mpmath_name = "hyperu"
nargs = {3}
rules = {
"HypergeometricU[0, c_, z_]": "1",
"HypergeometricU[a_, b_, z_] /; (a > 0) && (a-b+1 > 0)": "MeijerG[{{1-a},{}},{{0,1-b},{}},z]/Gamma[a]/Gamma[a-b+1]",
}
summary_text = "compute the Tricomi confluent hypergeometric function"
sympy_name = ""
59 changes: 59 additions & 0 deletions test/builtin/specialfns/test_hypergeom.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,59 @@
# -*- coding: utf-8 -*-
"""
Unit tests for mathics.builtins.specialfns.hypergeom
"""
from test.helper import check_evaluation

import pytest


@pytest.mark.parametrize(
("str_expr", "msgs", "str_expected", "fail_msg"),
[
(
"N[HypergeometricPFQ[{4},{},1]]",
None,
"ComplexInfinity",
None,
),
(
"MeijerG[{{},{}},{{0,0},{0,0}},100^2]",
None,
"MeijerG[{{}, {}}, {{0, 0}, {0, 0}}, 10000]",
None,
),
("N[MeijerG[{{},{}},{{0,0},{0,0}},100^2]]", None, "0.000893912", None),
(
"HypergeometricU[{3,1},{2,4},{7,8}]",
None,
"{MeijerG[{{-2}, {}}, {{0, -1}, {}}, 7] / 2, HypergeometricU[1, 4, 8]}",
None,
),
("N[HypergeometricU[{3,1},{2,4},{7,8}]]", None, "{0.00154364, 0.160156}", None),
("HypergeometricU[0,c,z]", None, "1", None),
(
"Hypergeometric1F1[{3,5},{7,1},{4,2}]",
None,
"{HypergeometricPFQ[{3}, {7}, 4], HypergeometricPFQ[{5}, {1}, 2]}",
None,
),
("N[Hypergeometric1F1[{3,5},{7,1},{4,2}]]", None, "{7.12435, 199.505}", None),
("Hypergeometric1F1[b,b,z]", None, "E ^ z", None),
("HypergeometricPFQ[{6},{1},2]", None, "HypergeometricPFQ[{6}, {1}, 2]", None),
("N[HypergeometricPFQ[{6},{1},2]]", None, "354.182", None),
("HypergeometricPFQ[{},{},z]", None, "1", None),
("HypergeometricPFQ[{0},{c1,c2},z]", None, "1", None),
("HypergeometricPFQ[{c1,c2},{c1,c2},z]", None, "E ^ z", None),
],
)
def test_private_hypergeom(str_expr, msgs, str_expected, fail_msg):
""" """
check_evaluation(
str_expr,
str_expected,
to_string_expr=True,
to_string_expected=True,
hold_expected=True,
failure_message=fail_msg,
expected_messages=msgs,
)
Loading