Skip to content
BranchesTags

Repository files navigation

HyperPrecisionAsir

Overview

HyperPrecisionAsir is an independent Risa/Asir implementation for complete Horn-type multivariate hypergeometric series. It follows the method of Banik--Bera, HyperPrecision (arXiv:2605.30216v2). The algebraic frontend constructs annihilating partial differential equations from neighboring coefficient ratios. Exact repeated differentiation and Gaussian elimination give a full Pfaffian connection

$$ \frac{\partial Y}{\partial x_i}=\Omega_i(x)Y. $$

The numerical layer first dispatches a predefined function to a term recurrence, a total-degree convolution, a neighbor-ratio shell recurrence, or Pfaffian continuation. It restricts a full Pfaffian connection to a piecewise-linear path and transports a factorized full fundamental matrix by local Taylor series. Based closed paths give numerical monodromy matrices. The implementation uses Risa/Asir for symbolic computation and MPFR midpoint arithmetic. It does not require Mathematica, FiniteFlow, or AMFlow.

For Appell's function $F_2$, the algebraic frontend obtains the derivative basis

$$ [F,\theta_xF,\theta_yF,\theta_x\theta_yF], $$

and holonomic rank $4$.

Features

The current implementation provides the following functions.

  1. It represents a general complete Horn series with arbitrary integral weights, including negative Pochhammer indices.
  2. It constructs neighboring coefficient ratios and Euler-type annihilating partial differential equations.
  3. It prolongs the partial differential equations by Euler derivatives and uses exact Gaussian elimination over a rational-function field.
  4. It extracts a finite derivative basis, the holonomic rank, and the full Pfaffian connection.
  5. It verifies the flatness equations of the full Pfaffian connection exactly.
  6. It supports directly supplied full Pfaffian connections.
  7. It provides constructors for generalized ${}_pF_q$, Appell F1–F4, Horn G1–G3 and H1–H7, and Lauricella $F_A$–$F_D$.
  8. It evaluates ${}_pF_q$ by an $O((p+q+1)D)$ term recurrence and routes ${}_2F_1(1/2,1/2;1;z)$ to an arithmetic-geometric mean recurrence on its principal interior branch.
  9. It evaluates Lauricella $F_A$, $F_B$, and $F_C$ by $O(nD^2)$ total-degree convolutions. Appell F2–F4 use the corresponding two-variable kernels, and Appell F1 uses the $F_D$ kernel.
  10. It evaluates Horn G1–G3 and H1–H7 by bivariate neighbor-ratio shells and uses exact Gauss reductions on coordinate axes.
  11. It evaluates Lauricella $F_D$ by a grouped total-degree recurrence or by a closed rank-$(n+1)$ Pfaffian connection, without generic multi-index enumeration or generic rational-function elimination.
  12. It selects the predefined method by convergence and operation gates and admits explicit series, pfaffian, agm, or auto requests where the method applies.
  13. It selects the Lauricella $F_D$ method by measured cost gates and admits an explicit series, pfaffian, or auto method.
  14. It rejects a generic exact prolongation before RREF when a rectangular monomial estimate exceeds two million matrix cells.
  15. It extracts singular factors from the denominators of the full Pfaffian connection and computes the restricted singularities of a path segment.
  16. It constructs canonical paths, conservative safe_opt paths, sampled fast_opt paths, and user paths.
  17. It transports all columns of a full fundamental matrix and stores the local matrices as a FactorizedFundamentalTransport object.
  18. It compares Taylor orders $N$ and $N+4$, evaluates a differential residual, performs an independent reverse-path check, and records precision escalation.
  19. It constructs one-variable based loops, coordinate loops, and multivariate meridian loops.
  20. It returns a numerical monodromy representation with its basepoint, derivative basis, generators, matrices, and computation history.
  21. It evaluates a distinguished Horn-series solution from a series boundary vector and reconstructs epsilon Laurent coefficients by a Cauchy grid.
  22. It provides an arithmetic-geometric mean evaluator for an independent Gauss ${}_2F_1$ regression.

Requirements and installation

The source files require Risa/Asir with MPFR support. If asir is on PATH, run

asir -quiet -f test/runtests.rr

The test runners use the following pinned public container image by default.

ghcr.io/nakanoryunosuke/risa-asir-container@sha256:a572b603f710167f708e976900ba883b44c03ecf51378e87a796b7b127bf0aed

For example, run the regular tests in Docker by

docker run --rm -v "$PWD:/workspace" -w /workspace \
  ghcr.io/nakanoryunosuke/risa-asir-container@sha256:a572b603f710167f708e976900ba883b44c03ecf51378e87a796b7b127bf0aed \
  asir -quiet -f test/runtests.rr

The environment variable RISA_ASIR_IMAGE replaces the default image used by the test and benchmark runners.

The file src/hyperprecision.rr includes its predefined numerical layer by a quoted source-relative include. Thus an absolute load of src/hyperprecision.rr also works when the current directory is not the repository root.

Quick start

The following computation constructs the full Pfaffian connection of Appell's function $F_2$ and evaluates its distinguished series solution.

load("src/hyperprecision.rr")$
setprec(55)$

H = hp_appell_f2(2,3/2,5/4,4,7/3,m,n,x,y)$
PF = hp_pfaffian(H,[tx,ty],2)$
print(hp_pfaffian_summary(PF))$

R = hp_transport(H,PF,[1/10,1/5],s,0,1/64,
                 26,28,50,24,16)$
print(R[0])$

The following computation evaluates a seven-variable Lauricella function. The report records the selected method, the truncation degree, and the number of generic terms avoided.

A = 1/4$
B = [1/4,1/5,1/6,1/7,1/8,1/9,1/10]$
C = 5/4$
X = [1/2,2/5,1/3,1/4,1/5,1/6,1/7]$

R = hp_lauricella_fd_eval_report(A,B,C,X,"auto",18)$
print(R)$

The result of hp_transport has the form

[function_value, final_basis_vector, ode_error_estimate,
 boundary_shell_estimate, accepted_frobenius_leaf_count,
 boundary_series_seconds, transport_seconds]

Thus R[5] measures the boundary-series construction and R[6] measures the subsequent numerical transport. Public predefined evaluation reports have the separate shape [value, method, metadata]; the flat diagnostics are the key–value records described below.

Derivative reports have the shape [value, derivatives, method, metadata], where derivatives contains all ordinary first derivatives in target-coordinate order. Their flat diagnostics contain the scalar fields together with derivatives, derivative_error_estimates, error_status, path_provenance, working_precision, and resource_admission. A derivative diagnostic calls one derivative report; it does not issue a second derivative evaluation.

The following computation transports a full fundamental matrix and computes a Gauss monodromy matrix about $z=1$.

G = hp_gauss_2f1(1/3,1/4,7/6,m,z)$
PG = hp_pfaffian(G,[tz],2)$

Path = hp_plan_path(PG,[1/5],[1/2],"principal","safe_opt",50)$
T = hp_transport_fundamental(PG,Path,s,18,50,12,500,"fast")$
U = hp_factorized_materialize(T)$
print(hp_factorized_diagnostics(T))$

Loop = hp_univariate_loop(PG,[1/5],1,1/8,10,"around_1",50)$
Rho = hp_monodromy(PG,[Loop],s,18,50,12,1000,"fast")$
M1 = hp_monodromy_matrix(Rho,["L",1])$

General predefined dispatch

The unified dispatcher accepts a family name, a parameter list, a target list, a method, a decimal precision, and a polynomial-detour parameter. The following calls use the generalized hypergeometric recurrence, the Appell $F_2$ convolution alias, and a named Horn shell recurrence, respectively.

P = hp_predefined_diagnostics(
    "HypergeometricPFQ",[[1,1],[2]],[1/2],"auto",30,0
)$

A = hp_predefined_diagnostics(
    "AppellF2",[1/3,2/5,3/7,5/6,7/8],
    [1/10,1/20],"auto",30,0
)$

H = hp_predefined_diagnostics(
    "HornH3",[1/3,2/5,5/6],[1/50,3/200],"auto",30,0
)$

Each flat diagnostic record contains value, family, method_used, degree, error_estimate, elapsed_seconds, compressed_dimension, and metadata. The nested metadata records the recurrence, convergence test, operation count, working precision, and branch detour used by the selected method.

For generalized hypergeometric functions, we set $u_0=1$ and use

$$ u_{k+1}=u_k z \frac{\prod_{i=1}^{p}(a_i+k)} {(k+1)\prod_{j=1}^{q}(b_j+k)}. $$

Thus a truncation through degree $D$ requires $O((p+q+1)D)$ arithmetic operations. Equal upper and lower parameters are canceled exactly before a lower-parameter pole or a terminating upper parameter is tested. A terminating polynomial is evaluated at every finite target. A nonterminating ${}_pF_q$ series is rejected when $p>q+1$, and a nonterminating ${}pF{p-1}$ series requires $|z|<1$. The stopping test waits until the index has passed the absolute sizes of the lower parameters and then uses a geometric ratio bound. The exact raw real or complex parameters determine additional guard digits from their distance to a nonpositive integer. For an exact complex parameter, the Euclidean distance is obtained from the squared rational components before the parameter is converted to MPFR. A finite sum or a guarded near-pole sum is repeated at a higher MPFR precision. The evaluator returns a value only after the two working precisions agree at the requested scale. The stability scale is at least one, so an exact value near zero is tested with an absolute tolerance.

The public constructor hp_bf(X,Prec) retains bounded provenance when an exact nonreal complex input is converted to MPFR. Hence a later evaluator can recover both the exact distance from a nonpositive integer and the source precision. The registry contains at most 64 rounded values. If distinct exact sources round to the same value, it retains the closest near-pole witness and the largest source precision, so a collision cannot decrease either guard. The real recurrence hot path does not query this registry. A bigfloat created outside hp_bf has no recoverable per-object precision in the Risa/Asir language; its metadata reports source_precision_digits=0.

For Lauricella $F_A$, $F_B$, and $F_C$, the implementation first forms one univariate coefficient array for each variable and convolves the arrays by total degree. It then multiplies the degree-$k$ shell by $(a)_k$, $1/(c)_k$, or $(a)_k(b)_k$, respectively. The defining-series gates are

$$ \sum_i|x_i|<1 \quad(F_A),\qquad \max_i|x_i|<1 \quad(F_B),\qquad \sum_i\sqrt{|x_i|}<1 \quad(F_C). $$

The nonterminating kernels compare degrees $D$ and $2D$. The reported uncertainty is a heuristic doubled-degree midpoint diagnostic; it is not a ball enclosure. Finite sums and lower parameters near a nonpositive integer also use a second working precision. Appell F1, F2, F3, and F4 use $F_D^{(2)}$, $F_A^{(2)}$, $F_B^{(2)}$, and $F_C^{(2)}$, respectively. Their derivative reports accumulate the value and all ordinary first derivatives in the same convolution state pass when the specialized series applies. A Pfaffian report extracts all derivatives from the same transported state. Shifted-parameter reports are used for a coordinate removed by an exact radial compression. Before guard selection, exact $F_A$ coordinate factors with $b_i=c_i$ are canceled and zero coordinates are removed on a radial path. The same normalization applies to Appell $F_2$. A derivative in a removed coordinate is not assumed to vanish: its shifted-parameter identity is evaluated when it is defined and rejected when its own lower factor is singular. A nonzero Delta activates every coordinate on the polynomial detour, so zero-coordinate compression is deliberately disabled for that path.

Horn G1–G3 and H1–H7 use uncanceled neighboring coefficient ratios, except that an identical upper/lower parameter pair with an identical weight row is canceled exactly before dynamic guard selection. The implementation sums complete total-degree shells and compares degrees $D$ and $2D$. Automatic series dispatch uses a conservative interior gate. On a coordinate axis, each named Horn function is reduced exactly to a Gauss function by the negative-index and duplication identities. Nonnegative terminating upper-weight rows give an exact finite-support test. When these rows bound both indices, the polynomial is evaluated at every finite target without applying the interior gate.

Each public predefined evaluator rounds its value to the requested decimal precision with two output guard digits. The metadata contains output_digits, output_guard_digits, output_rounding_error, precision_stability_error, and an error_estimate_kind label. The total error_estimate includes the observed output-rounding change.

A nonzero Delta always selects Pfaffian continuation. The dispatcher does not replace such a request by a principal defining series, an arithmetic-geometric mean recurrence, or a terminating-polynomial shortcut. The sign and magnitude of Delta are retained in the report. The automatic dispatcher does not enable Euler, Laplace, Bessel, or Mellin–Barnes integrals.

Automatic production dispatch uses predicted arithmetic work before it starts a series or an exact Pfaffian construction. For ${}pF{p-1}$, the model compares the predicted $O(D)$ recurrence work with a precision-dependent scalar-Pfaffian crossover. For Appell functions and the two-variable Lauricella aliases, it compares the $O(D^2)$ convolution work with a measured rank-four crossover. The named Horn model uses a separate two-million-operation admission limit because a neighbor-ratio cell costs more than a convolution cell in the pinned runtime. Exact termination and exact closed reductions precede these comparisons. The metadata records dispatch_reason, projected_series_operations, pfaffian_crossover_operations, and resource_admission.

The exact-system cache stores at most eight Pfaffian systems. It does not store function values, boundary vectors, transports, or error checks. Thus a second target with the same exact parameters reuses Macaulay/RREF construction and still performs the complete numerical computation. The function hp_fast_clear_system_cache() clears this cache. Pfaffian reports separate pfaffian_build_seconds, boundary_series_seconds, and transport_seconds, and they state whether the exact system was found in the cache.

Lauricella FD dispatch

We group the Lauricella $F_D$ series by total degree. We set

$$ Q(t)=\prod_{i=1}^n(1-x_it)^{-b_i}=\sum_{k=0}^{\infty}q_kt^k, . $$

By logarithmic differentiation, we have

$$ kq_k=\sum_{j=1}^k\left(\sum_{i=1}^n b_ix_i^j\right)q_{k-j}, . $$

Thus the implementation computes a degree-$D$ truncation in $O(D^2+nD)$ arithmetic operations. It does not generate the $\binom{D+n}{n}$ multi-indices. The same recurrence evaluates the first derivatives used by the Pfaffian boundary vector. For the adaptive scalar evaluation, we put

$$ q=\max_i|x_i|,\qquad B=\sum_i|b_i|, \qquad M_k=q^k\left|\frac{(a)_k}{(c)_k}\right|\frac{(B)_k}{k!}. $$

The quantity $M_k$ bounds the absolute value of the degree-$k$ shell. For $j&gt;|c|$, the remaining ratios admit the upper bound

$$ \frac{M_{j+1}}{M_j} \leq q\frac{j+|a|}{j-|c|} \max\left{1,\frac{j+B}{j+1}\right}. $$

The adaptive stopping test requires this upper bound to be less than one and uses the resulting geometric bound for the complete omitted tail. Thus a few small shells before a near-pole amplification cannot trigger termination.

The function hp_lauricella_fd_pfaffian(A,B,C,Variables) gives the connection with respect to the basis

$$ \left[F,\frac{\partial F}{\partial x_1},\ldots, \frac{\partial F}{\partial x_n}\right]. $$

The constructor uses the closed Lauricella differential equations and does not call hp_prolong_relations or hp_rref. Its holonomic rank is $n+1$.

The functions hp_lauricella_fd_eval and hp_lauricella_fd_eval_report accept "series", "pfaffian", and "auto". A nonterminating series call requires $\max_i|x_i|&lt;1$; when $a$ is a nonpositive integer, the resulting polynomial is valid at every finite target. The pfaffian convenience method uses a radial path in the interior and a default polynomial detour outside the defining polydisc. The ordinary-derivative chart requires pairwise distinct target coordinates. The lower-level path API admits other nonsingular paths. The auto method uses the grouped series on inexpensive interior points and uses the closed Pfaffian connection near the boundary or outside the defining polydisc. If $a$ is a nonpositive integer, auto evaluates the terminating series through degree $-a$. If $c=a$ and $\max_i|x_i|&lt;1$, the principal auto or series route first cancels the common parameter and uses

$$ F_D^{(n)}(a;b_1,\ldots,b_n;a;x_1,\ldots,x_n) =\prod_{i=1}^n(1-x_i)^{-b_i}. $$

The cancellation also applies when the common parameter is a nonpositive integer: it is performed before the lower-pole test. The product is restricted to the principal germ. A branch-specific request or a nonterminating target outside the defining polydisc retains Pfaffian continuation. For such an outside target, the convenience method adds an upper-half-plane polynomial detour with Delta=1/20; the report records this choice. The function hp_lauricella_fd_eval_report_with_delta preserves a user-specified sign and magnitude of Delta. Equal coordinates inside the polydisc are combined by adding their $b_i$ parameters. Factors with $x_i=0$ or total exponent zero are removed. This principal-germ compression also covers the diagonal reduction to Gauss ${}_2F_1$. The reduced problem uses the $O(D)$ generalized hypergeometric recurrence or its scalar Pfaffian continuation; it does not enter the $O(D^2)$ grouped $F_D$ recurrence. On the calibrated $|z|\leq 0.9$ interior route, an exact diagonal auto request and a forced series request enter the same reduced pFq kernel without rebuilding an $F_D$ selector or report. An exact rational coordinate is compared with $0.9$ exactly. For a pre-existing BF coordinate, the comparison uses the working precision and retains a $p+5$ decimal uncertainty band; a coordinate in that band is returned to the ordinary pFq selector. In particular, rounding cannot move a point outside the calibrated region into the direct series route. For distinct but close coordinates with $\max_i|x_i|\leq0.9$, auto uses the grouped series because the ordinary-derivative connection contains the factors $(x_i-x_j)^{-1}$. Exact coordinate compression is performed before parameter guard selection, so large opposite exponents that sum to zero do not cause a spurious guard-cap failure. This reduction is not applied to a nonzero-Delta path. When $a=c$, the principal auto or series evaluation removes the common parameter from the dynamic magnitude guard and uses the endpoint product; an explicit series request retains series as its reported method. An explicit Pfaffian or nonzero-Delta request retains both the full continuation guard and the requested path in its report.

The flat report returned by hp_lauricella_fd_diagnostics contains value, method_used, degree, error_estimate, elapsed_seconds, and compressed_dimension. For the grouped series, error_estimate is the scalar tail majorant or the doubled-degree difference reported by the corresponding convergence test. For Pfaffian continuation, it is the maximum of the accepted local ODE term estimate and the boundary-shell estimate. The nested metadata entry retains the method-specific data.

The evaluator compares raw coordinates before conversion to MPFR values. Therefore, exact equal coordinates are combined and distinct exact coordinates are retained. The working precision includes guard digits determined by the decimal orders of the parameters and the nonzero coordinate separations. If the absolute-parameter tail majorant is too large, auto compares two fixed grouped sums after the degree has passed the possible amplification near $-c$. A successful comparison is reported as doubled_degree. A rejected comparison falls back to the closed Pfaffian connection, and the metadata records fallback_from and fallback_reason.

The generic Horn evaluator rejects a truncation with more than $10^6$ terms. The named Horn shell evaluator admits at most $2\times10^6$ predicted operations, while the grouped Lauricella evaluator has a separate hard limit of $2\times10^7$ recurrence operations. These gates reject an oversized request before enumeration begins. If an automatically selected series exceeds its operation limit or convergence cap, auto changes to the closed Pfaffian connection when the ordinary-derivative chart is nonsingular. The pinned Risa/Asir runtime does not expose a native adaptive quadrature callback for an Asir function. Therefore, auto selects between the measured grouped-series and closed-Pfaffian paths; it does not run an interpreted quadrature loop.

Full Pfaffian connection

The function hp_horn represents the series

$$ F(x)=\sum_{m\in\mathbb N_0^n} \frac{\prod_r(a_r)_{\mu_r\mathbin{\cdot}m}} {\prod_s(b_s)_{\nu_s\mathbin{\cdot}m}} \frac{x^m}{m!}. $$

The corresponding Risa/Asir input is

H = hp_horn(
    [a,b1,b2],
    [[1,1],[1,0],[0,1]],
    [c1,c2],
    [[1,0],[0,1]],
    [m,n],
    [x,y]
)$

The functions hp_neighbor_ratios(H) and hp_pde_generator(H,[tx,ty]) return the neighboring coefficient ratios and the annihilating partial differential equations. For each coordinate $i$, hp_pde_generator returns [g_i(m),h_i(m),L_i(theta)], where

$$ L_i=h_i(\theta-e_i)-x_i g_i(\theta). $$

The factors $g_i$ and $h_i$ remain uncancelled during this construction. This convention retains boundary factors at resonant parameter values.

The main algebraic APIs are as follows.

API Operation
hp_horn(A,Mu,B,Nu,Indices,Variables) Construct a Horn series
hp_hypergeometric_pfq(Upper,Lower,m,z) Construct a generalized hypergeometric series
hp_gauss_2f1(...) Construct Gauss ${}_2F_1$
hp_appell_f1(...)–hp_appell_f4(...) Construct Appell F1–F4
hp_horn_g1(...)–hp_horn_g3(...), hp_horn_h1(...)–hp_horn_h7(...) Construct Horn G1–G3 and H1–H7
hp_lauricella_fa(...)–hp_lauricella_fc(...) Construct Lauricella $F_A$–$F_C$
hp_lauricella_fd(...) Construct Lauricella $F_D$
hp_pfq_eval(...) Evaluate ${}_pF_q$ by recurrence, AGM, or Pfaffian continuation
hp_appell_f1_eval(...)–hp_appell_f4_eval(...) Evaluate an Appell function through a specialized Lauricella alias
hp_horn_g1_eval(...)–hp_horn_g3_eval(...), hp_horn_h1_eval(...)–hp_horn_h7_eval(...) Evaluate a named Horn function by shells or continuation
hp_lauricella_fa_eval(...)–hp_lauricella_fc_eval(...) Evaluate Lauricella $F_A$–$F_C$ by convolution or continuation
hp_predefined_eval(...) Evaluate a named predefined family through one dispatcher
hp_predefined_diagnostics(...) Return the common flat diagnostic record
hp_predefined_derivative_eval_report(...) Return a value and all first derivatives through one dispatcher
hp_predefined_derivative_diagnostics(...) Return the common flat derivative diagnostic record
hp_pfq_derivative_eval_report(...) Return ${}_pF_q$ and its first derivative
hp_pfq_derivative_diagnostics(...) Return the flat ${}_pF_q$ derivative diagnostic record
hp_appell_f1_derivative_eval_report(...)–hp_appell_f4_derivative_eval_report(...) Return an Appell value and both first derivatives
hp_appell_f1_derivative_diagnostics(...)–hp_appell_f4_derivative_diagnostics(...) Return a flat Appell derivative diagnostic record
hp_lauricella_fa_derivative_eval_report(...)–hp_lauricella_fd_derivative_eval_report(...) Return a Lauricella value and all first derivatives
hp_lauricella_fa_derivative_diagnostics(...)–hp_lauricella_fd_derivative_diagnostics(...) Return a flat Lauricella derivative diagnostic record
hp_horn_derivative_eval_report(...) Return a named Horn value and both first derivatives
hp_horn_derivative_diagnostics(...) Return a flat named-Horn derivative diagnostic record
hp_lauricella_fd_series_value(...) Evaluate $F_D$ by the grouped recurrence
hp_lauricella_fd_pfaffian(...) Construct the closed rank-$(n+1)$ connection
hp_lauricella_fd_eval(...) Evaluate $F_D$ by auto, series, or pfaffian
hp_lauricella_fd_eval_report(...) Evaluate $F_D$ and return method diagnostics
hp_lauricella_fd_diagnostics(...) Return flat value, method, error, timing, and dimension diagnostics
hp_lauricella_fd_eval_report_with_limit(...) Evaluate $F_D$ with an explicit degree cap
hp_lauricella_fd_eval_report_with_delta(...) Continue $F_D$ with a specified polynomial detour
hp_neighbor_ratio(H,i) Compute a neighboring coefficient ratio
hp_pde_generator(H,Theta) Construct Euler-type annihilators
hp_prolong_relations(H,Theta,Depth) Prolong the annihilators
hp_pfaffian_resource_estimate(H,Theta,Depth) Estimate generic exact matrix cells before prolongation
hp_pfaffian(H,Theta,Depth) Construct the full Pfaffian connection
hp_user_pfaffian(Basis,Variables,Omegas) Construct a directly supplied connection
hp_check_integrability(PF) Verify flatness exactly
hp_holonomic_rank(PF) Return the rank
hp_singular_factors(PF) Extract singular factors
hp_restricted_singularities(PF,Paths,t,Prec) Compute restricted singularities

The unified constructor, scalar evaluator, and derivative evaluator use one canonical schema validator. HypergeometricPFQ takes [upper,lower], where both entries are scalar parameter lists, and Gauss2F1 takes [a,b,c]; their targets have dimension one. Appell F1 and F4 take four scalar parameters, whereas F2 and F3 take five; their targets have dimension two. Lauricella $F_A$–$F_D$ take [a,b-list,c-list], [a-list,b-list,c], [a,b,c-list], and [a,b-list,c], respectively, with every list matching the target dimension. Horn G1–G3 take 3, 4, and 2 scalar parameters. Horn H1–H7 take 4, 5, 3, 4, 3, 3, and 4 scalar parameters. A wrong arity, a misplaced scalar/list group, or a target-dimension mismatch raises an explicit schema error before an evaluator indexes the supplied parameters.

Risa/Asir matrices use zero-based indices. In a Pfaffian object, PF[1] is the derivative basis, PF[2] is the list of connection matrices, and PF[7] is the list of variables.

Path restriction uses simultaneous substitution through fresh placeholders. A system variable that has the same name as a path parameter is not captured by a later substitution. Restricted polynomial roots are computed after normalization by the leading coefficient. The solver checks iteration convergence and the residual of every root against the original nonmonic polynomial. A singular factor that vanishes identically on a path is rejected.

The following constructor defines a full Pfaffian connection supplied by a user.

Omega = matrix(1,1)$
Omega[0][0] = 1/(3*(1-z))$
P = hp_user_pfaffian([[0]],[z],[Omega])$

Exact rational inputs such as 1/10 retain the requested working precision. A machine double such as 0.1 does not recover its lost digits when Prec is increased.

Path planning

A path object has the form

["PiecewiseLinearPath", vertices, path_class, planner, metadata]

The planner API is

Path = hp_plan_path(PF,Start,Target,PathClass,Planner,Prec)$

The planners have the following meanings.

  • canonical returns a deterministic coordinate-by-coordinate reference path.
  • safe_opt returns the canonical path without modification. Its metadata records homotopy_status="not_certified_no_change".
  • fast_opt samples the strip between the direct path and the canonical path. It can accept the direct path when the sample avoids the singular factors and the predicted restricted-root step count does not increase. This finite sampling does not certify the path class.
  • hp_user_path(Vertices,PathClass) stores a user path and its path class without modification.

For example, the following two paths have the same complex endpoint.

P1 = hp_plan_path(PF,[1/5],[1/2+@i/10],
                  "principal","safe_opt",50)$
P2 = hp_user_path([[1/5],[2/5+@i/20],[1/2+@i/10]],
                  "principal")$

The metadata of a shortcut accepted by fast_opt contains homotopy_status="sampled_only_not_certified" and certified=0. The automatic meridian connector uses the canonical path. It does not select a fast_opt shortcut implicitly.

The earlier distinguished-solution API remains available. The call

hp_transport(H,PF,Target,t,Delta,T0,Cutoff,
             Order,Prec,GoalDigits,MaxSplit)$

uses a radial path when Delta=0. A nonzero value of Delta adds the polynomial detour

$$ x_i(t)=\mathrm{Target}_i t+ \sqrt{-1},\mathrm{Delta}(i+1)t(1-t). $$

The sign of Delta selects the side of a real singular locus. A negative value corresponds to the default IDelta -> -I convention in the reference paper.

Fundamental transport

The full fundamental transport API is

T = hp_transport_fundamental(PF,Path,t,Order,Prec,
                             GoalDigits,MaxSteps,"fast")$

At an ordinary center, the algorithm expands

$$ A(t)=\sum_{k\geq 0}A_kz^k, \qquad U(t)=\sum_{n\geq 0}U_nz^n, $$

and computes

$$ U_0=I, \qquad U_{n+1}=\frac{1}{n+1}\sum_{k=0}^{n}A_kU_{n-k}. $$

The local step is at most two thirds of the distance to the closest restricted singularity. The algorithm compares orders $N$ and $N+4$ and evaluates the differential residual at the endpoint of each local step.

The result has the form

["FactorizedFundamentalTransport",
 start, end, rank, local_factors, path, history, diagnostics, mode]

The local factors are stored in traversal order. The function hp_factorized_apply(T,V) applies the local factors successively to a vector. The function hp_factorized_materialize(T) returns their matrix product. The function hp_factorized_inverse(T) reverses and inverts the local factors. This algebraic inverse records reverse_error=-1 and reverse_check="not_checked"; it is not an independently transported reverse path.

The computation history of each local factor contains the path segment, local parameter interval, Taylor order, working precision, order-comparison error, differential residual, restricted radius, and number of step halvings. The transport diagnostics contain the maximum local error, accepted step count, reverse-path error, and precision-escalation count.

The reverse-path error is computed from an independent transport along the reversed geometric path. It measures

$$ \left|T_{\gamma^{-1}}T_\gamma-I\right|_{\max}. $$

If this error exceeds the requested threshold, the fast mode repeats the forward and reverse transports once with 20 additional working digits and four additional Taylor orders.

Monodromy

The standard monodromy computation transports a full fundamental matrix around a based closed path. It does not replace a resonant closed-path transport by the exponential of a residue matrix.

For a one-variable singularity $s$, the function hp_univariate_loop connects the basepoint to a polygon centered at $s$, traverses the polygon counterclockwise, and returns to the basepoint.

Loop = hp_univariate_loop(PF,Base,s,Radius,Sides,PathClass,Prec)$
Rho = hp_monodromy(PF,[Loop],t,Order,Prec,
                   GoalDigits,MaxSteps,"fast")$
M = hp_monodromy_matrix(Rho,["L",1])$

For a multivariate singular factor $Q_j(x)$, the function hp_meridian_generator restricts $Q_j$ to coordinate lines through the basepoint. It selects an isolated smooth root $p_j$ and a coordinate direction $v_j$ for which

$$ dQ_j(p_j)[v_j]\ne 0. $$

The function tests every polygon edge by restricting every singular factor to the edge. It also tests radial spokes for enclosed foreign components and additional germs of the target component. The function reduces the meridian radius until these tests pass. The following code constructs a meridian of the divisor $x=1$ for Appell's function $F_3$.

Factors = hp_singular_factors(PF3)$
IX = hp_find(Factors,x-1)$
G = hp_meridian_generator(PF3,[1/5,1/7],IX,1/8,8,50)$
Rho = hp_monodromy(PF3,[G],s,18,50,12,1500,"fast")$
Mx1 = hp_monodromy_matrix(Rho,["D",IX+1])$

The function hp_meridian_generators constructs several meridians from a list of component indices or from the string "all". A numerical monodromy representation has the form

["NumericalMonodromyRepresentation",
 basepoint, basis, generators, matrices, verified_relations,
 generator_set_complete, history, mode]

The value of generator_set_complete is "unknown". The automatic meridians are not asserted to give a presentation of the full fundamental group. Before transport, hp_monodromy requires exact flatness, at least one loop, closed paths of the correct dimension, a common numerical basepoint, and unique generator labels.

The functions hp_matrix_trace, hp_matrix_determinant, and hp_projective_distance provide numerical matrix invariants.

Tests and benchmarks

The PowerShell test commands are

powershell -ExecutionPolicy Bypass -File .\run-tests.ps1
powershell -ExecutionPolicy Bypass -File .\run-tests.ps1 -Monodromy
powershell -ExecutionPolicy Bypass -File .\run-tests.ps1 -Extended
powershell -ExecutionPolicy Bypass -File .\run-tests.ps1 -All

The corresponding POSIX commands are

sh ./run-tests.sh
sh ./run-tests.sh --monodromy
sh ./run-tests.sh --extended
sh ./run-tests.sh --all

The regular test suite verifies Horn ratios, annihilators, full Pfaffian connections for Appell functions F1–F4, exact flatness, distinguished-solution transport, epsilon reconstruction, and the arithmetic-geometric mean regression. It also verifies the seven-variable Lauricella diagonal reduction, an ordinary derivative identity, the closed rank-eight constructor, and agreement between the grouped series and the closed Pfaffian transport at a shared distinct point. The default runner also executes test/general_hypergeometric_fast.rr, test/production_dispatch.rr, and test/phase3_derivatives.rr. The production regression verifies the measured one- and two-variable crossovers, exact-system cache behavior, separated Pfaffian timings, the $F_D$ diagonal-to-Gauss reduction, complete diagnostics, and preservation of an explicit polynomial detour. Four negative scripts verify that expensive F2 and Horn shells and an oversized closed $F_D$ connection are rejected before allocation, and that a canonical predefined family with the wrong parameter schema fails explicitly. The general hypergeometric suite uses independent decimal values for Lauricella $F_A$, $F_B$, and $F_C$ and for Horn G1–G3 and H1–H7. It tests an exact terminating ${}_pF_q$ polynomial after parameter cancellation, a lower-parameter amplification probe, an AGM value after exact cancellation, an explicit nonzero detour, exact coordinate-axis reductions, and ordinary-versus-Euler derivative normalization. It tests large-cancellation ${}_1F_0$, ${}_0F_0$, and Lauricella $F_A$ values, four independent decimal oracles with a lower parameter $-100+10^{-100}$, complex approaches to the same poles, the 4096-digit guard cap, removable zero-parameter derivatives, derivative-domain errors, and an exact named-Horn polynomial outside the interior gate. Additional adversarial cases cover 200-, 300-, and 600-degree terminating cancellation, Euclidean complex guard boundaries, 5000-digit exact FA/F2 and Horn cancellations, radial zero-coordinate compression with every ordinary derivative, the principal $a=c$ and equal-coordinate $F_D/F_1$ reductions, and preservation of nonzero-detour Pfaffian paths. Matching uncancelled and singular-derivative scripts verify fail-closed behavior. The absolute-load regression executes from /tmp. The suite also checks that a generic seven-variable $F_D$ prolongation is rejected before its estimated $31{,}171{,}875$ matrix cells are formed.

The Phase 3 regression checks the negative-integer $a=c$ cancellation against a fixed principal-product oracle and confirms that a nonzero detour remains a Pfaffian path. It preserves a 210-digit hp_bf parameter at distance $10^{-100}$ from a lower pole and compares a cancellation-sensitive degree-101 sum with an independent exact finite sum. At 15, 50, and 100 digits, independent finite-polynomial sums check every component of the pFq, Appell, Lauricella, and Horn derivative reports, including their component error estimates. The canonical schema matrix covers pFq, Gauss, Appell F1–F4, Lauricella $F_A$–$F_D$, and Horn G1–G3 and H1–H7. At the origin, the G3, H6, and H7 reports obtain both first derivatives from exact first coefficients rather than finite differences. The same precisions retain the fixed $q=0.999$ oracle.

The monodromy test suite verifies the following computations.

  • Simultaneous substitution preserves a path parameter, and the verified root solver obtains the roots $1/2$ and $1$ of $2t^2-3t+1$.
  • Noncommuting rank-two local factors are multiplied in traversal order.
  • A directly supplied rank-one connection reaches a complex target along a canonical path and a user path.
  • The Appell $F_2$ safe_opt path equals the canonical path, while an explicitly requested fast_opt path uses fewer measured local Taylor steps.
  • The Gauss ${}_2F_1$ monodromy about $z=1$ has the expected trace and determinant.
  • The resonant Appell function $F_3(1,1,1,1;2;x,y)$ has nonidentity unipotent monodromy about $x=1$, including the numerical relation $(M-I)^2=0$.
  • Open paths, mixed basepoints, duplicate labels, nonflat connections, invalid loop geometry, and identically singular restrictions are rejected.
  • The meridian regression for the factors $x$ and $16x-y-1$ forces radius reduction when the foreign divisor crosses a radial spoke.
  • The Lauricella function $F_D^{(3)}$ gives a flat full Pfaffian connection of rank $4$.

Run the wall-time benchmark by

powershell -ExecutionPolicy Bypass -File .\run-benchmarks.ps1
powershell -ExecutionPolicy Bypass -File .\run-production-benchmarks.ps1

The first command runs the regular benchmark set. The second command also runs the 15-, 50-, and 100-digit production portfolio. It reports the median of five container-plus-empty-Asir startup samples separately from source loading and warm numerical kernels. The first paired index contains every admitted forced method and auto, and the first method rotates with the sample index. Thus host-load drift is shared across the methods. After that first index, a forced candidate that exceeds five seconds and is not the fastest forced method keeps one complete warm sample. The fastest forced method and auto retain at least five paired samples. A nonfastest candidate in the 4.5--5-second boundary band retains one full five-sample round; this prevents a small host-speed change from expanding a three-round measurement abruptly. The benchmark prints all medians, sampling policies, sample counts, and paired samples before it evaluates the timing gate. At 50 digits, competitive candidates use three rounds of five paired samples and the median of the three round medians. The $q=0.95$ pFq control row uses the same three-round rule at 50 digits because its subsecond measurements are sensitive to a bimodal host load. For this row, the output records the automatic dispatch overhead and the method, degree, value, error estimate, and path signature in every round. The forced and automatic signatures must agree, and failures of the timing inequality in at least two of the three rounds are a hard failure. Workloads outside these 50-digit controls use one round. The output also contains an adaptive candidate matrix. Every mathematically applicable series, Pfaffian, AGM, and automatic candidate is labelled admitted or resource_rejected; a hard resource refusal is not reported as inapplicable. An admitted long-running candidate retains one complete warm sample. Scalar and derivative workloads are measured separately with five paired, alternating warm samples at each requested precision.

The benchmark runner includes benchmark/general_hypergeometric_fast.rr. Every dispatch gate uses warmed calls and requires

$$ t_{\mathrm{auto}}\leq 1.25\min_j t_j+0.005\ \text{seconds}, $$

where the index $j$ ranges over the forced methods included at the same target. The production portfolio uses a 0.005-second timer floor. Its fastest forced candidate and automatic route use at least five warm samples; the long-candidate rule above retains the required one-sample evidence for a slower forced route. It compares identical inputs, requested digits, paths, and error checks. A forced method receives no numerical sample only after its operation admission rejects the request before allocation.

The production portfolio retains $q=0.999$ as a fixed-oracle functional and resource-admission regression at 15, 50, and 100 digits. A legacy run that repeated every 100-digit $q=0.999$ forced candidate five times exceeded 13 minutes before the portfolio completed. The adaptive policy now measures the $q=0.999$ workload itself: auto and the fastest forced candidate use five paired samples, while a slower forced candidate above five seconds keeps one complete sample. The timing gate is applied directly to this $q=0.999$ row. The $q=0.995$ row remains a secondary crossover probe and is not a substitute for the primary $q=0.999$ gate. A single host-dependent pinned-container probe at 100 digits required 1.36 seconds for the $q=0.995$ recurrence and 13.09 seconds for scalar Pfaffian continuation. The precision-dependent crossover selects the recurrence in this regime.

All wall times in this section depend on the host. On 2026-08-20, the pinned image gave a 0.447-second median for container plus empty-Asir startup and 0.0679 seconds for source loading. The production portfolio gave the following warm medians. Each row passed the dispatch gate.

Workload 15 digits 50 digits 100 digits
Generic ${}_2F_1$, $q=0.95$ 0.0188 s, series 0.0309 s, series 0.0603 s, series
Generic ${}_2F_1$, near unit 0.477 s, pfaffian, $q=0.999$ 1.62 s, pfaffian, $q=0.999$ 0.598 s, series, $q=0.995$
Appell $F_2$, score $0.95$ 0.421 s, pfaffian 3.98 s, pfaffian 23.9 s, pfaffian
Appell $F_4$, score $0.90$ 0.292 s, pfaffian 2.55 s, pfaffian 20.0 s, pfaffian
Horn $H_3$, score $0.80$ 0.191 s, pfaffian 1.89 s, pfaffian 17.1 s, pfaffian
Seven-variable diagonal $F_D$, $q=0.90$ 0.00943 s, series 0.0215 s, series 0.0384 s, series

At 15 digits, the main checkout selected series for the $F_2$ and $H_3$ rows. We compare the main checkout and this branch at five paired sample indices on the same host. Each sample uses one warmed Asir runtime per checkout, and the checkout order alternates with the sample index. The internal timer excludes container startup and source loading. The following checkout medians are host-dependent.

15-digit workload main checkout Current branch main / current
Appell $F_2$, score $0.95$ 7.40 s, series 0.821 s, pfaffian 9.02
Horn $H_3$, score $0.80$ 2.26 s, series 0.336 s, pfaffian 6.72

At 100 digits, the F2, F4, and H3 series estimates exceed their family-specific admission limits. Their rows compare forced Pfaffian and automatic calls at each of the same five sample indices. On 2026-08-20, one pinned-container run gave the following warm medians.

Workload Forced method auto Selected method
Gauss ${}_2F_1(1/2,1/2;1;1/2)$ 0.0000861 s (AGM) 0.0000920 s agm
Appell $F_2$ convolution 0.0181 s 0.0172 s series
Horn $H_3$ shells 0.152 s 0.151 s series
Three-variable Lauricella $F_C$ 0.0194 s 0.0231 s series

At the Appell $F_2$ benchmark point, the generic degree-60 multi-index sum required 0.308 s and the specialized convolution required 0.0181 s. At the Horn $H_3$ point, the generic degree-60 sum required 0.335 s and the shell recurrence required 0.152 s. Wall times depend on the host.

On 2026-08-19, the pinned container completed the unchanged regular regression suite in 6.018 seconds. The monodromy regression suite, including the explicit differential-residual and meridian-edge checks, completed in 33.267 seconds. The measured Appell $F_2$ step counts were 11 for the canonical path and 9 for the fast_opt path. The fast_opt reverse-path error was less than $8.55\times10^{-12}$. Wall times depend on the host.

The Lauricella benchmark uses the same parameters and target as the quick-start example at 18 requested digits. On 2026-08-19, 20 closed rank-eight constructions required 9.46 seconds, or 0.473 seconds per construction. Twenty grouped-series evaluations required 0.203 seconds, or 0.0101 seconds per evaluation. The warm medians were 0.01108 seconds for forced series and 0.01462 seconds for auto; auto selected series. One closed-Pfaffian transport required 4.61 seconds. The two values differed by $2.21\times10^{-23}$. The quarter-parameter seven-variable oracle uses a five-second dispatch gate in the same benchmark. The terminating-series and $c=a$ product cases required less than 0.00019 seconds each. The corresponding forced Pfaffian and forced series evaluations required 0.0723 seconds and 0.00776 seconds, respectively. The $10^{40}$ signed-parameter cancellation case required 0.301 seconds and used the grouped series with dynamic guard precision. At degree 40, generic seven-variable enumeration would generate 62,891,499 terms. The benchmark is implemented in benchmark/lauricella_fd7.rr and is included in run-benchmarks.ps1.

The extended test reproduces the epsilon Laurent coefficients of the paper's Appell $F_2$ example,

$$ F_2\left(2,\frac32,1+\epsilon;4,-1-\epsilon;3,\frac{11}{3}\right). $$

The pinned container gives

epsilon^-1   0.5149686368228051 - 6.97e-11 i
epsilon^0    0.5286628168944232 - 4.194390017098441 i
epsilon^1  -10.97813788622874   - 4.834942153461091 i

The paper gives 0.5149686376, 0.528662817 - 4.194390019 i, and -10.978138236 - 4.834942296 i, respectively.

The independent arithmetic-geometric mean test uses

$$ {}_2F_1\left(\frac12,\frac12;1;z\right) =\frac{1}{\mathrm{AGM}(1,\sqrt{1-z})}. $$

The regular test suite includes $z=1-10^{-4}$ and $z=1-10^{-8}$. At the latter point, the difference between the two algorithms is approximately $1.19\times10^{-11}$ with 65 MPFR working digits.

Numerical modes and guarantees

The fast mode uses arbitrary-precision MPFR midpoint arithmetic. It records the following diagnostics.

  • It compares local Taylor orders $N$ and $N+4$.
  • It evaluates the last Taylor term and a differential residual.
  • It limits each local step by the restricted singularities.
  • It transports the reversed path independently.
  • It increases the working precision and Taylor order once when the reverse-path threshold fails.

The function hp_numerical_capabilities() returns the corresponding capability flags.

The fast mode does not provide a rigorous enclosure. Its errors are midpoint diagnostics, not certified bounds.

The predefined series layer follows the same convention. The ${}_pF_q$ recurrence reports an analytic ratio-tail estimate after the index has passed the lower-parameter magnitude guard. The Lauricella convolution and named Horn shell kernels label their tail values as heuristic doubled-degree estimates. Exact terminating sums report zero truncation error together with the observed change between two working precisions. The reported total also includes the output-rounding change. None of these quantities is a complex ball enclosure.

The certified mode is not implemented. Risa/Asir does not provide the complex ball arithmetic, restricted-root ball isolation, coefficient balls, and tail bounds required by this mode. A request for "certified" raises an explicit error. The function hp_certified_mode_available() returns 0.

The safe_opt planner leaves the canonical path unchanged because midpoint sampling cannot certify a path class. The fast_opt planner can use a midpoint sampling test, and its metadata records certified=0.

Limitations

  • The algebraic frontend assumes a complete Horn system of finite holonomic rank. It does not prove completeness or the convergence domain.
  • Automatic named-Horn series dispatch uses a conservative interior measure. A rejected point can still belong to the mathematical convergence domain.
  • The exact named-Horn polynomial test uses terminating upper rows with nonnegative weights. A terminating mechanism that requires a cancellation among signed-weight rows is not inferred automatically.
  • Generic Pfaffian construction uses a rectangular monomial estimate with a two-million-cell limit. The estimate prevents large exact RREF jobs; it is not an estimate of the holonomic rank.
  • Automatic dispatch does not use Euler, Laplace, Bessel, or Mellin–Barnes integral representations. These methods are not enabled without a precision-controlled backend and branch-path diagnostics.
  • Generic Pfaffian continuation for generalized hypergeometric functions requires $p=q+1$. The recurrence also evaluates convergent cases with $p\leq q$ and every terminating polynomial.
  • Precision stability is tested up to 4096 working digits. The evaluator raises an error when a finite or near-pole sum does not stabilize before this cap.
  • The prolongation depth Depth is a user parameter. A larger depth can cause expression swell in exact rational-function elimination.
  • The implementation uses Risa/Asir exact Gaussian elimination. It does not use FiniteFlow-style modular reconstruction.
  • The local Frobenius method treats an ordinary point with exponent zero. It does not construct generalized Frobenius or Levelt bases at a singular point.
  • The safe_opt planner does not shorten the canonical path without an interval or exact homotopy certificate.
  • The fast_opt planner uses a finite strip sample and compares the direct path with the canonical path. It does not prove a path-class identity, construct a PRM/A* candidate graph, or move general waypoints.
  • The automatic meridian search requires an isolated smooth point visible on a coordinate line through the basepoint. It does not search arbitrary transverse directions.
  • The automatic meridians do not certify a complete generator set. Braid generators and a presentation of the fundamental group are not constructed.
  • Invariant Hermitian forms, invariant bilinear forms, relation-word certification, and algebraic recognition are not implemented.
  • The polynomial Delta detour does not extrapolate an i0 limit.
  • The scalar grouped-series report uses an analytic tail majorant evaluated in MPFR midpoint arithmetic. The Pfaffian boundary-shell error and the epsilon-reconstruction error are empirical estimates.
  • Epsilon reconstruction requires a user-supplied Laurent range, radius, and sample count. The pole order is not inferred, and the samples run serially.
  • Singular parameter specializations can change the holonomic rank.
  • Complex ball certified mode is not implemented.

License and references

This repository is distributed under GNU GPL version 3 only (GPL-3.0-only). See LICENSE, NOTICE, and THIRD_PARTY_NOTICES.md.

This repository is an independent Risa/Asir port and reimplementation. It does not include the reference Mathematica source. The algorithms and formulas are based on the paper by Banik--Bera and its GPL-3.0 reference implementation. Risa/Asir, OpenXM, and the Docker image are separate distributions.

OpenAI Codex contributed to the initial Risa/Asir implementation, tests, and documentation. NAKANO Ryuosuke specified the mathematical scope, the Risa/Asir container, the independent arithmetic-geometric mean boundary test, the license, and the release requirements.

About

Risa/Asir implementation of HyperPrecision for high-precision Horn-type hypergeometric functions

Resources

Readme
GPL-3.0 license
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases

Packages

Contributors

Languages