diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a1e2b2..0c24b92 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## 12.1.0 2026-08-17 + +* Rust + * Add allocation-free `hyp2f1_scalar`, `hyp2f1`, and `hyp2f1_par` for fully complex parameters and arguments +* Python + * Add a `complex128` array binding for `hyp2f1` + ## 12.0.0 2026-08-11 * Rust diff --git a/Cargo.lock b/Cargo.lock index cd09805..69b5481 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -128,7 +128,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "cfsem" -version = "12.0.0" +version = "12.1.0" dependencies = [ "criterion", "faer", @@ -136,6 +136,8 @@ dependencies = [ "itertools 0.15.0", "libm", "nalgebra", + "ndarray", + "num-complex", "num-traits", "num_cpus", "numpy", @@ -1020,6 +1022,7 @@ dependencies = [ "portable-atomic", "portable-atomic-util", "rawpointer", + "rayon", ] [[package]] diff --git a/Cargo.toml b/Cargo.toml index 6dd1e19..2a77c1d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "12.0.0" +version = "12.1.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" @@ -20,12 +20,16 @@ unsafe_code = "forbid" [dependencies] pyo3 = { version = "0.29.2", features = ["extension-module", "abi3-py310", "generate-import-lib"], optional = true } numpy = { version="0.29.0", optional=true } # This must match pyo3 version! +# Keep this version aligned with numpy's ndarray dependency. Enabling rayon on +# the shared crate instance provides Zip::par_for_each to the Python bindings. +ndarray = { version="0.17.2", features=["rayon"], optional=true } nalgebra = "^0.34.2" rayon = "^1.12.0" num_cpus = "^1.17.0" libm = "^0.2" num-traits = { version = "0.2.19", features = ["libm"] } +num-complex = "0.4.6" faer = "0.24.4" faer-traits = "0.24.0" @@ -46,7 +50,7 @@ codegen-units = 16 [features] default = [] -python = ["pyo3", "numpy"] +python = ["pyo3", "numpy", "ndarray"] [[bench]] name = "linear_filament" @@ -64,5 +68,9 @@ harness = false name = "boundary_element" harness = false +[[bench]] +name = "hyp2f1" +harness = false + [package.metadata.docs.rs] rustdoc-args = [ "--html-in-header", "katex-header.html" ] diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md new file mode 100644 index 0000000..65717f1 --- /dev/null +++ b/THIRD_PARTY_NOTICES.md @@ -0,0 +1,32 @@ +# Third-Party Notices + +## HypergeometricFunctions.jl + +Parts of the private complex gamma-difference helpers and hypergeometric +continuation recurrences in `src/math/hyp2f1.rs` are adapted from +[HypergeometricFunctions.jl](https://github.com/JuliaMath/HypergeometricFunctions.jl), +version 0.3.30, under the following license: + +> MIT License +> +> Copyright (c) 2018-2023 Richard Mikael Slevinsky and other contributors: +> +> https://github.com/JuliaMath/HypergeometricFunctions.jl/graphs/contributors +> +> Permission is hereby granted, free of charge, to any person obtaining a copy +> of this software and associated documentation files (the "Software"), to deal +> in the Software without restriction, including without limitation the rights +> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +> copies of the Software, and to permit persons to whom the Software is +> furnished to do so, subject to the following conditions: +> +> The above copyright notice and this permission notice shall be included in all +> copies or substantial portions of the Software. +> +> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +> SOFTWARE. diff --git a/benches/hyp2f1.py b/benches/hyp2f1.py new file mode 100644 index 0000000..83c9854 --- /dev/null +++ b/benches/hyp2f1.py @@ -0,0 +1,265 @@ +"""Compare Python-level hyp2f1 throughput for cfsem, mpmath, and SciPy. + +The complex workload exercises several numerical regions and compares cfsem's +array interface with mpmath's scalar complex implementation. The real workload +uses ``z < 1`` so that SciPy and cfsem evaluate the same real branch. + +By default, parallel cfsem processes 13,107,200 values, SciPy processes +6,553,600 values, serial cfsem processes 655,360 values, and mpmath processes +a matching 1,024-value prefix. Run the benchmark with:: + + uv run python benches/hyp2f1.py + +The unequal sizes keep each timed call reasonably short despite the large +throughput difference. The table reports each implementation's sample count. +Use ``--repeats 3`` when more stable timings are worth the wait. +""" + +from __future__ import annotations + +import argparse +import gc +import statistics +import time +from collections.abc import Callable, Sequence +from dataclasses import dataclass +from typing import TypeVar + +import mpmath as mp +import numpy as np +from scipy.special import hyp2f1 as scipy_hyp2f1 + +import cfsem + +DEFAULT_NATIVE_SIZE = 10 * (1 << 16) +DEFAULT_PARALLEL_SIZE = 20 * DEFAULT_NATIVE_SIZE +DEFAULT_SCIPY_SIZE = 10 * DEFAULT_NATIVE_SIZE +DEFAULT_MPMATH_SIZE = 1 << 10 + +# These points cover the direct series, a terminating polynomial, Pfaff's +# transformation, expansions near one and infinity, Taylor continuation, and +# a near-integer parameter difference. +COMPLEX_CASES = np.asarray( + [ + (0.5 + 0.25j, 1.25 - 0.5j, 2.0 + 0.75j, 0.1 + 0.2j), + (-2.0 + 0.0j, 1.2 + 0.4j, 3.5 - 0.2j, 2.0 + 0.5j), + (0.7 + 0.2j, 1.3 - 0.1j, 2.4 + 0.3j, -3.0 + 0.4j), + (0.4 + 0.2j, 1.1 + 0.3j, 2.5 + 0.5j, 0.98 + 0.03j), + (0.4 + 0.2j, 1.1 + 0.3j, 2.7 - 0.2j, 4.0 + 2.0j), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.5 + 0.866_025_403_784_438_6j), + (0.4 + 0.2j, 0.9 - 0.1j, 3.3 + 0.100_000_001j, 0.97 + 0.02j), + ], + dtype=np.complex128, +) + +# SciPy's real-valued interface does not return the complex limiting value for +# real z > 1. Keep this comparison below the cut while retaining direct, +# transformed, near-one, polynomial, and moderate-parameter cases. +REAL_CASES = np.asarray( + [ + (0.5, 1.1, 2.4, 0.2), + (1.2, 0.7, 2.8, -0.5), + (0.4, 1.1, 2.5, 0.95), + (0.7, 1.3, 3.1, -3.0), + (-3.0, 1.4, 2.5, 0.8), + (12.5, 9.25, 17.0, 0.4), + ], + dtype=np.float64, +) + +T = TypeVar("T") + + +@dataclass(frozen=True) +class Timing: + """Median wall time and derived element throughput for one implementation.""" + + name: str + seconds: float + size: int + + @property + def values_per_second(self) -> float: + return self.size / self.seconds + + +def tiled_arguments(cases: np.ndarray, size: int, dtype: np.dtype) -> tuple[np.ndarray, ...]: + """Tile case rows into four contiguous argument arrays of exactly ``size`` elements.""" + + indices = np.arange(size) % len(cases) + return tuple(np.ascontiguousarray(cases[indices, column], dtype=dtype) for column in range(4)) + + +def time_call(name: str, function: Callable[[], T], size: int, repeats: int) -> tuple[Timing, T]: + """Time a no-argument callable while excluding cyclic-GC bookkeeping.""" + + samples = [] + result: T | None = None + gc_was_enabled = gc.isenabled() + gc.disable() + try: + for _ in range(repeats): + start = time.perf_counter() + result = function() + samples.append(time.perf_counter() - start) + finally: + if gc_was_enabled: + gc.enable() + assert result is not None + return Timing(name, statistics.median(samples), size), result + + +def print_timings(title: str, timings: Sequence[Timing]) -> None: + """Print a compact throughput table.""" + + print(f"\n{title}") + print(f"{'implementation':<30} {'values':>10} {'seconds':>10} {'values/s':>15}") + print(f"{'-' * 30} {'-' * 10} {'-' * 10} {'-' * 15}") + for timing in timings: + print( + f"{timing.name:<30} {timing.size:>10,} " + f"{timing.seconds:>10.4f} {timing.values_per_second:>15,.0f}" + ) + + +def mpmath_arguments(arguments: tuple[np.ndarray, ...]) -> tuple[tuple[mp.mpc, ...], ...]: + """Convert NumPy inputs once so mpmath conversion is outside the timed loop.""" + + return tuple(tuple(mp.mpc(value) for value in argument) for argument in arguments) + + +def time_cfsem_parallel( + cases: np.ndarray, size: int, validation_size: int, repeats: int +) -> tuple[Timing, np.ndarray]: + """Time parallel cfsem and retain only the prefix needed for validation.""" + + arguments = tiled_arguments(cases, size, np.dtype(np.complex128)) + out = np.empty(size, dtype=np.complex128) + cfsem.hyp2f1(*(argument[: len(cases)] for argument in arguments), par=True) + timing, values = time_call( + "cfsem (parallel)", + lambda: cfsem.hyp2f1(*arguments, par=True, out=out), + size, + repeats, + ) + return timing, values[:validation_size].copy() + + +def benchmark_complex(parallel_size: int, native_size: int, mpmath_size: int, repeats: int, dps: int) -> None: + """Benchmark complex128 cfsem arrays against scalar mpmath evaluation.""" + + mp.mp.dps = dps + mp_arguments = mpmath_arguments(tiled_arguments(COMPLEX_CASES, mpmath_size, np.dtype(np.complex128))) + + # Initialize mpmath caches without evaluating the full workload twice. + mp.hyp2f1(*(argument[0] for argument in mp_arguments)) + + parallel, validation_values = time_cfsem_parallel(COMPLEX_CASES, parallel_size, mpmath_size, repeats) + + serial_arguments = tiled_arguments(COMPLEX_CASES, native_size, np.dtype(np.complex128)) + serial_out = np.empty(native_size, dtype=np.complex128) + serial, _ = time_call( + "cfsem (serial)", + lambda: cfsem.hyp2f1(*serial_arguments, par=False, out=serial_out), + native_size, + repeats, + ) + mpmath, mpmath_values = time_call( + "mpmath (scalar loop)", + lambda: [mp.hyp2f1(a, b, c, z) for a, b, c, z in zip(*mp_arguments, strict=True)], + mpmath_size, + repeats, + ) + + sample = np.linspace(0, mpmath_size - 1, min(mpmath_size, 32), dtype=np.intp) + reference = np.asarray([complex(mpmath_values[index]) for index in sample]) + np.testing.assert_allclose(validation_values[sample], reference, rtol=3e-10, atol=3e-11) + print_timings("Complex inputs", [parallel, serial, mpmath]) + + +def benchmark_real(parallel_size: int, native_size: int, scipy_size: int, repeats: int) -> None: + """Benchmark cfsem and SciPy for wholly real input arrays below the cut.""" + + parallel, validation_values = time_cfsem_parallel(REAL_CASES, parallel_size, scipy_size, repeats) + + complex_arguments = tiled_arguments(REAL_CASES, native_size, np.dtype(np.complex128)) + serial_out = np.empty(native_size, dtype=np.complex128) + serial, _ = time_call( + "cfsem (serial)", + lambda: cfsem.hyp2f1(*complex_arguments, par=False, out=serial_out), + native_size, + repeats, + ) + + real_arguments = tiled_arguments(REAL_CASES, scipy_size, np.dtype(np.float64)) + scipy_hyp2f1(*(argument[: len(REAL_CASES)] for argument in real_arguments)) + scipy, scipy_values = time_call( + "scipy.special.hyp2f1", lambda: scipy_hyp2f1(*real_arguments), scipy_size, repeats + ) + + np.testing.assert_allclose(validation_values.real, scipy_values, rtol=3e-10, atol=3e-11) + np.testing.assert_allclose(validation_values.imag, 0.0, atol=3e-11) + print_timings("Real inputs (z < 1)", [parallel, serial, scipy]) + + +def positive_integer(value: str) -> int: + """Parse a strictly positive command-line integer.""" + + parsed = int(value) + if parsed <= 0: + raise argparse.ArgumentTypeError("must be a positive integer") + return parsed + + +def main() -> None: + """Parse benchmark settings and run both workloads.""" + + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter + ) + parser.add_argument( + "--native-size", + "--size", + dest="native_size", + type=positive_integer, + default=DEFAULT_NATIVE_SIZE, + help="elements in each serial cfsem workload", + ) + parser.add_argument( + "--mpmath-size", + type=positive_integer, + default=DEFAULT_MPMATH_SIZE, + help="elements in the mpmath workload", + ) + parser.add_argument( + "--parallel-size", + type=positive_integer, + default=DEFAULT_PARALLEL_SIZE, + help="elements in each parallel cfsem workload", + ) + parser.add_argument( + "--scipy-size", + type=positive_integer, + default=DEFAULT_SCIPY_SIZE, + help="elements in the SciPy workload", + ) + parser.add_argument("--repeats", type=positive_integer, default=1, help="timed passes per implementation") + parser.add_argument("--dps", type=positive_integer, default=17, help="mpmath decimal precision") + args = parser.parse_args() + if args.mpmath_size > args.native_size: + parser.error("--mpmath-size must not exceed --native-size") + if args.parallel_size < max(args.native_size, args.scipy_size): + parser.error("--parallel-size must not be smaller than --native-size or --scipy-size") + + print( + f"hyp2f1 throughput: parallel size={args.parallel_size:,}, " + f"serial size={args.native_size:,}, " + f"SciPy size={args.scipy_size:,}, " + f"mpmath size={args.mpmath_size:,}, repeats={args.repeats}, mpmath dps={args.dps}" + ) + benchmark_complex(args.parallel_size, args.native_size, args.mpmath_size, args.repeats, args.dps) + benchmark_real(args.parallel_size, args.native_size, args.scipy_size, args.repeats) + + +if __name__ == "__main__": + main() diff --git a/benches/hyp2f1.rs b/benches/hyp2f1.rs new file mode 100644 index 0000000..f072bcc --- /dev/null +++ b/benches/hyp2f1.rs @@ -0,0 +1,126 @@ +#![allow(clippy::all)] // Criterion requires black-boxed benchmark inputs and outputs. + +use cfsem::math::{hyp2f1, hyp2f1_par, hyp2f1_scalar}; +use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; +use num_complex::Complex64; +use std::hint::black_box; +use std::time::Duration; + +type Case = (&'static str, Complex64, Complex64, Complex64, Complex64); + +fn cases() -> [Case; 7] { + [ + ( + "direct", + Complex64::new(0.5, 0.25), + Complex64::new(1.25, -0.5), + Complex64::new(2.0, 0.75), + Complex64::new(0.1, 0.2), + ), + ( + "polynomial", + Complex64::new(-2.0, 0.0), + Complex64::new(1.2, 0.4), + Complex64::new(3.5, -0.2), + Complex64::new(2.0, 0.5), + ), + ( + "pfaff", + Complex64::new(0.7, 0.2), + Complex64::new(1.3, -0.1), + Complex64::new(2.4, 0.3), + Complex64::new(-0.5, 0.1), + ), + ( + "one", + Complex64::new(0.4, 0.2), + Complex64::new(1.1, 0.3), + Complex64::new(2.5, 0.5), + Complex64::new(0.98, 0.03), + ), + ( + "infinity", + Complex64::new(0.4, 0.2), + Complex64::new(1.1, 0.3), + Complex64::new(2.7, -0.2), + Complex64::new(4.0, 2.0), + ), + ( + "taylor", + Complex64::new(0.7, 0.2), + Complex64::new(1.2, -0.3), + Complex64::new(2.1, 0.1), + Complex64::new(0.5, 0.866_025_403_784_438_6), + ), + ( + "near-integer", + Complex64::new(0.4, 0.2), + Complex64::new(0.9, -0.1), + Complex64::new(3.3, 0.100_000_001), + Complex64::new(0.97, 0.02), + ), + ] +} + +fn bench_scalar_regions(criterion: &mut Criterion) { + let mut group = criterion.benchmark_group("hyp2f1 scalar regions"); + group.sample_size(20); + group.measurement_time(Duration::from_secs(2)); + for (label, a, b, c, z) in cases() { + group.bench_function(label, |bencher| { + bencher.iter(|| { + black_box(hyp2f1_scalar( + black_box(a), + black_box(b), + black_box(c), + black_box(z), + )) + }); + }); + } + group.finish(); +} + +fn bench_vectors(criterion: &mut Criterion) { + let mut group = criterion.benchmark_group("hyp2f1 vectors"); + group.sample_size(10); + group.measurement_time(Duration::from_secs(2)); + let regions = cases(); + for length in [1, 64, 4096, 65_536] { + let mut a = Vec::with_capacity(length); + let mut b = Vec::with_capacity(length); + let mut c = Vec::with_capacity(length); + let mut z = Vec::with_capacity(length); + for index in 0..length { + let (_, ai, bi, ci, zi) = regions[index % regions.len()]; + a.push(ai); + b.push(bi); + c.push(ci); + z.push(zi); + } + let mut out = vec![Complex64::ZERO; length]; + + group.bench_with_input( + BenchmarkId::new("scalar loop", length), + &length, + |bench, _| { + bench.iter(|| { + for index in 0..length { + out[index] = hyp2f1_scalar(a[index], b[index], c[index], z[index]); + } + black_box(&out); + }); + }, + ); + group.bench_with_input(BenchmarkId::new("serial", length), &length, |bench, _| { + bench.iter(|| black_box(hyp2f1(&a, &b, &c, &z, &mut out).unwrap())); + }); + group.bench_with_input(BenchmarkId::new("parallel", length), &length, |bench, _| { + bench.iter(|| black_box(hyp2f1_par(&a, &b, &c, &z, &mut out).unwrap())); + }); + } + group.finish(); +} + +criterion_group!(hyp2f1_benches, bench_scalar_regions, bench_vectors); +criterion_main!(hyp2f1_benches); diff --git a/cfsem/__init__.py b/cfsem/__init__.py index 618037a..f0ab930 100644 --- a/cfsem/__init__.py +++ b/cfsem/__init__.py @@ -58,7 +58,7 @@ ) from cfsem.types import Array3xN -from .cfsem import ellipe, ellipk +from .cfsem import ellipe, ellipk, hyp2f1 from .cfsem import DimensionalityError MU_0 = 4.0 * np.pi * 1e-7 * (1.0 + 5.5e-10) @@ -109,6 +109,7 @@ "self_inductance_distributed_axisymmetric_conductor", "ellipe", "ellipk", + "hyp2f1", "rotate_filaments_about_path", "vector_potential_linear_filament", "vector_potential_linear_filament_hierarchical", diff --git a/cfsem/cfsem.pyi b/cfsem/cfsem.pyi index 92149d4..1c9e658 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -1,9 +1,11 @@ from typing import TypeAlias, TypedDict -from numpy import float32, float64, int64, uint64 +from numpy import complex128, float32, float64, int64, uint64 from numpy.typing import NDArray FloatArray: TypeAlias = NDArray[float64] +ComplexArray: TypeAlias = NDArray[complex128] +ComplexInput: TypeAlias = complex | complex128 | ComplexArray Float32Array: TypeAlias = NDArray[float32] IntArray: TypeAlias = NDArray[int64] UIntArray: TypeAlias = NDArray[uint64] @@ -210,6 +212,15 @@ def body_force_density_linear_filament( ) -> ArrayTriple: ... def ellipe(x: float) -> float: ... def ellipk(x: float) -> float: ... +def hyp2f1( + a: ComplexInput, + b: ComplexInput, + c: ComplexInput, + z: ComplexInput, + par: bool = True, + *, + out: ComplexArray | None = None, +) -> ComplexArray: ... def filament_helix_path( path: ArrayTriple, helix_start_offset: tuple[float, float, float], diff --git a/docs/python/math.md b/docs/python/math.md index bb2a000..5183c23 100644 --- a/docs/python/math.md +++ b/docs/python/math.md @@ -5,3 +5,5 @@ ::: cfsem.ellipe ::: cfsem.ellipk + +::: cfsem.hyp2f1 diff --git a/pyproject.toml b/pyproject.toml index 932cf7e..9930451 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,6 +29,8 @@ dev = [ "pytest-cov == 7.1.0", "ruff == 0.6.2", "ty == 0.0.42", + # Reference generation and benchmarks + "mpmath == 1.3.0", # Docs "mkdocs == 1.6.0", "mkdocs-material == 9.5.25", @@ -45,6 +47,7 @@ features = ["python"] module-name = "cfsem.cfsem" version = "cargo" profile = "release" +include = [{ path = "THIRD_PARTY_NOTICES.md", format = ["sdist", "wheel"] }] [tool.ruff] line-length = 110 diff --git a/src/math.rs b/src/math.rs index 48e80b1..491b542 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,5 +1,9 @@ //! Pure-math functions supporting physics calculations. +mod hyp2f1; + +pub use hyp2f1::{hyp2f1, hyp2f1_par, hyp2f1_scalar}; + use core::ops::{Add, Div, Mul, Sub}; use num_traits::{Float, FromPrimitive}; diff --git a/src/math/hyp2f1.rs b/src/math/hyp2f1.rs new file mode 100644 index 0000000..052bd6f --- /dev/null +++ b/src/math/hyp2f1.rs @@ -0,0 +1,1947 @@ +//! Gauss hypergeometric function with complex parameters and argument. +//! +//! The scalar implementation is an allocation-free polyalgorithm. It combines +//! the defining Gauss series with Euler and Pfaff transformations, stabilized +//! connection expansions about `z = 1` and `z = infinity`, and Taylor +//! continuation through the region in which none of those series converges +//! rapidly. The same scalar kernel backs the serial and Rayon-parallel array +//! interfaces. +//! +//! Complex powers use their principal values. Consequently, values on the +//! conventional branch cut `[1, infinity)` depend on the sign of `z.im`, +//! including signed zero. See [`hyp2f1_scalar`] for the complete public +//! contract and references. + +use num_complex::Complex64; +use rayon::prelude::*; + +use crate::macros::check_length; + +const NAN: Complex64 = Complex64::new(f64::NAN, f64::NAN); +const REL_TOL: f64 = 8.0 * f64::EPSILON; +const DIRECT_RADIUS: f64 = 0.9; +const TAYLOR_INNER_ANCHOR_RADIUS: f64 = 0.875; +const TAYLOR_OUTER_ANCHOR_RADIUS: f64 = 1.1; +const MAX_SERIES_ITERATIONS: usize = 10_000; +const MAX_TAYLOR_ITERATIONS: usize = 512; +const MAX_PARAMETER_ITERATIONS: usize = 10_000; +const LANCZOS_G_MINUS_HALF: f64 = 4.242_187_5; +const LOG_SQRT_TWO_PI: f64 = 0.918_938_533_204_672_7; +const LANCZOS_COEFFICIENTS: [f64; 15] = [ + 0.999_999_999_999_997_1, + 57.156_235_665_862_92, + -59.597_960_355_475_49, + 14.136_097_974_741_746, + -0.491_913_816_097_620_2, + 0.000_033_994_649_984_811_89, + 0.000_046_523_628_927_048_58, + -0.000_098_374_475_304_879_65, + 0.000_158_088_703_224_912_5, + -0.000_210_264_441_724_104_88, + 0.000_217_439_618_115_212_64, + -0.000_164_318_106_536_763_9, + 0.000_084_418_223_983_852_74, + -0.000_026_190_838_401_581_067, + 0.000_003_689_918_265_953_162_5, +]; + +/// Evaluates Gauss's hypergeometric function on its principal branch. +/// +/// This computes +/// +/// `2F1(a, b; c; z) = sum((a)_n (b)_n z^n / ((c)_n n!), n = 0..infinity)` +/// +/// by selecting among the defining series, Euler/Pfaff transformations, +/// stabilized connection expansions about `z = 1` and `z = infinity`, and a +/// Taylor continuation fallback. All four arguments may be complex, and `a` +/// and `b` are interchangeable. +/// +/// Complex powers take their principal values. On the conventional branch cut +/// `z` in `[1, infinity)`, positive and negative zero in `z.im` therefore select +/// the upper and lower limiting values, respectively. +/// +/// A complex NaN is returned for non-finite inputs, mathematical singularities +/// (including a nonpositive-integer `c` unless the series terminates before its +/// pole), or failure of an internal expansion to converge within its limit. +/// +/// # References +/// +/// \[1\] N. Michel and M. V. Stoitsov, “Fast computation of the Gauss +/// hypergeometric function with all its parameters complex with +/// application to the Pöschl–Teller–Ginocchio potential wave functions,” +/// *Computer Physics Communications*, vol. 178, no. 7, pp. 535–551, +/// Apr. 2008, doi: +/// [10.1016/j.cpc.2007.11.007](https://doi.org/10.1016/j.cpc.2007.11.007). +/// +/// \[2\] NIST Digital Library of Mathematical Functions, “§15.2 Definitions +/// and Analytical Properties,” NIST. Accessed: Aug. 17, 2026. \[Online\]. +/// Available: +/// +/// \[3\] JuliaMath, “HypergeometricFunctions.jl,” ver. 0.3.30, GitHub. +/// Accessed: Aug. 17, 2026. \[Online\]. Available: +/// +/// +/// \[4\] SciPy Developers, “scipy.special.hyp2f1,” *SciPy API Reference*. +/// Accessed: Aug. 17, 2026. \[Online\]. Available: +/// +#[inline] +pub fn hyp2f1_scalar(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> Complex64 { + if !finite(a) || !finite(b) || !finite(c) || !finite(z) { + return NAN; + } + if z == Complex64::ZERO || a == Complex64::ZERO || b == Complex64::ZERO { + return Complex64::ONE; + } + + if let Some(degree) = terminating_degree(a, b) { + if let Some(pole_degree) = negative_integer_degree(c) + && degree > pole_degree + { + return NAN; + } + let result = terminating_series(a, b, c, z, degree); + return if result.converged { result.value } else { NAN }; + } + if is_nonpositive_integer(c) { + return NAN; + } + if z == Complex64::ONE { + let balance = c - a - b; + if balance.re <= 0.0 { + return NAN; + } + if terminating_degree(c - a, c - b).is_some() { + return Complex64::ZERO; + } + return gamma_ratio(&[c, balance], &[c - a, c - b]); + } + if c == a { + return (-b * complex_log1p(-z)).exp(); + } + if c == b { + return (-a * complex_log1p(-z)).exp(); + } + let result = general_evaluation(a, b, c, z); + if result.converged { result.value } else { NAN } +} + +/// Evaluates Gauss's hypergeometric function elementwise on equal-length, +/// contiguous slices. +/// +/// Each output element is `hyp2f1_scalar(a[i], b[i], c[i], z[i])`. See +/// [`hyp2f1_scalar`] for the mathematical definition, branch convention, +/// failure policy, implementation notes, and references. +/// +/// Returns `Err("Length mismatch")` without modifying `out` if any input slice +/// has a different length from `out`. +pub fn hyp2f1( + a: &[Complex64], + b: &[Complex64], + c: &[Complex64], + z: &[Complex64], + out: &mut [Complex64], +) -> Result<(), &'static str> { + check_length!(out.len(), a, b, c, z); + for i in 0..out.len() { + out[i] = hyp2f1_scalar(a[i], b[i], c[i], z[i]); + } + Ok(()) +} + +/// Evaluates Gauss's hypergeometric function elementwise in parallel on +/// equal-length, contiguous slices. +/// +/// Rayon dynamically partitions the elementwise work so that expensive +/// continuation cases do not pin an entire worker. See [`hyp2f1_scalar`] for +/// the mathematical definition, branch convention, failure policy, +/// implementation notes, and references. +/// +/// Returns `Err("Length mismatch")` without modifying `out` if any input slice +/// has a different length from `out`. +pub fn hyp2f1_par( + a: &[Complex64], + b: &[Complex64], + c: &[Complex64], + z: &[Complex64], + out: &mut [Complex64], +) -> Result<(), &'static str> { + check_length!(out.len(), a, b, c, z); + out.par_iter_mut() + .zip(a.par_iter()) + .zip(b.par_iter()) + .zip(c.par_iter()) + .zip(z.par_iter()) + .for_each(|((((out, a), b), c), z)| { + *out = hyp2f1_scalar(*a, *b, *c, *z); + }); + Ok(()) +} + +/// A value returned by one candidate expansion together with convergence +/// status. +/// +/// Keeping failure information separate from the value lets the path selector +/// map all internal singularities and iteration-limit failures to the public +/// complex-NaN policy in one place. +#[derive(Clone, Copy, Debug)] +struct EvalOutcome { + value: Complex64, + converged: bool, +} + +impl EvalOutcome { + #[inline] + const fn success(value: Complex64) -> Self { + Self { + value, + converged: true, + } + } + + #[inline] + const fn failure() -> Self { + Self { + value: NAN, + converged: false, + } + } +} + +// Branch-aware complex primitives. These small helpers preserve precision or +// signed-zero information that the straightforward formulas can lose. + +#[inline] +fn finite(z: Complex64) -> bool { + z.re.is_finite() && z.im.is_finite() +} + +#[inline] +fn complex_abs(z: Complex64) -> f64 { + z.re.hypot(z.im) +} + +/// Computes `ln(1 + z)` without first rounding `1 + z` near the origin. +/// +/// The `atan2` expression also preserves which side of the negative real axis +/// was approached, which is required for the branch cut of `hyp2f1`. +#[inline] +fn complex_log1p(z: Complex64) -> Complex64 { + if z == -Complex64::ONE { + return Complex64::new(f64::NEG_INFINITY, z.im); + } + let quadratic = 2.0 * z.re + z.re * z.re + z.im * z.im; + if quadratic > -1.0 { + Complex64::new(0.5 * quadratic.ln_1p(), z.im.atan2(1.0 + z.re)) + } else { + (Complex64::ONE + z).ln() + } +} + +#[inline] +fn complex_expm1(z: Complex64) -> Complex64 { + let (sin_y, cos_y) = z.im.sin_cos(); + Complex64::new( + z.re.exp_m1() * cos_y - 2.0 * (0.5 * z.im).sin().powi(2), + z.re.exp() * sin_y, + ) +} + +#[inline] +fn complex_pow(base: Complex64, exponent: Complex64) -> Complex64 { + (exponent * base.ln()).exp() +} + +/// Computes `1 / z` with scaled complex division to reduce avoidable overflow +/// and underflow when the real and imaginary components have unlike sizes. +#[inline] +fn complex_inverse(z: Complex64) -> Complex64 { + if z.re.abs() >= z.im.abs() { + let ratio = z.im / z.re; + let denominator = z.re + z.im * ratio; + Complex64::new(1.0 / denominator, -ratio / denominator) + } else { + let ratio = z.re / z.im; + let denominator = z.im + z.re * ratio; + Complex64::new(ratio / denominator, -1.0 / denominator) + } +} + +/// Forms `1 - z` while retaining the negated sign of an exactly zero +/// imaginary part. Ordinary complex subtraction rounds `+0.0 - +0.0` back to +/// `+0.0`, which loses the upper/lower-lip distinction on the branch cut. +#[inline] +fn one_minus(z: Complex64) -> Complex64 { + Complex64::new(1.0 - z.re, -z.im) +} + +/// Forms the Pfaff argument `z / (z - 1) = 1 + 1 / (z - 1)` without allowing +/// addition of the real unit to erase the inverse's signed imaginary zero. +#[inline] +fn pfaff_argument(z: Complex64) -> Complex64 { + let inverse = complex_inverse(z - 1.0); + Complex64::new(1.0 + inverse.re, inverse.im) +} + +// Private gamma machinery. The connection formulas need complex gamma, +// digamma, and differences of nearly equal reciprocal-gamma values; keeping +// them local avoids adding a second public special-function API. + +#[inline] +fn is_nonpositive_integer(z: Complex64) -> bool { + z.im == 0.0 && z.re <= 0.0 && z.re.is_finite() && z.re == z.re.trunc() +} + +#[inline] +fn nearest_integer_difference(z: Complex64) -> (i32, Complex64) { + let nearest = z.re.round().clamp(i32::MIN as f64, i32::MAX as f64) as i32; + (nearest, z - nearest as f64) +} + +#[inline] +fn sin_cos_pi_real(x: f64) -> (f64, f64) { + let nearest = x.round(); + let remainder = x - nearest; + let sign = if (nearest % 2.0).abs() == 1.0 { + -1.0 + } else { + 1.0 + }; + let (sine, cosine) = (core::f64::consts::PI * remainder).sin_cos(); + (sign * sine, sign * cosine) +} + +#[inline] +fn sin_pi(z: Complex64) -> Complex64 { + let y = core::f64::consts::PI * z.im; + let (sine, cosine) = sin_cos_pi_real(z.re); + Complex64::new(sine * y.cosh(), cosine * y.sinh()) +} + +#[inline] +fn cos_pi(z: Complex64) -> Complex64 { + let y = core::f64::consts::PI * z.im; + let (sine, cosine) = sin_cos_pi_real(z.re); + Complex64::new(cosine * y.cosh(), -sine * y.sinh()) +} + +fn log_sin_pi(z: Complex64) -> Complex64 { + let y = core::f64::consts::PI * z.im; + if y.abs() < 20.0 { + return sin_pi(z).ln(); + } + // For large imaginary parts, evaluating sinh/cosh before taking the log + // would overflow. Use the leading exponential form directly instead. + let (sine, cosine) = sin_cos_pi_real(z.re); + Complex64::new( + y.abs() - core::f64::consts::LN_2, + (z.im.signum() * cosine).atan2(sine), + ) +} + +#[inline] +fn cot_pi(z: Complex64) -> Complex64 { + let two_x = 2.0 * core::f64::consts::PI * (z.re - z.re.round()); + let two_y = 2.0 * core::f64::consts::PI * z.im; + if two_y.abs() > 350.0 { + return Complex64::new(0.0, -z.im.signum()); + } + let denominator = two_y.cosh() - two_x.cos(); + Complex64::new(two_x.sin() / denominator, -two_y.sinh() / denominator) +} + +#[inline] +fn sinc_pi(z: Complex64) -> Complex64 { + if z == Complex64::ZERO { + Complex64::ONE + } else { + sin_pi(z) / (core::f64::consts::PI * z) + } +} + +fn lanczos_sum(z: Complex64) -> Complex64 { + let mut sum = Complex64::new(LANCZOS_COEFFICIENTS[0], 0.0); + for (index, coefficient) in LANCZOS_COEFFICIENTS.iter().enumerate().skip(1) { + sum += coefficient / (z + (index - 1) as f64); + } + sum +} + +/// Principal complex log-gamma from a 15-term Lanczos approximation, with the +/// reflection formula used to move arguments out of the left half-plane. +fn log_gamma(z: Complex64) -> Complex64 { + if is_nonpositive_integer(z) { + return Complex64::new(f64::INFINITY, f64::NAN); + } + if z.re < 0.5 { + return Complex64::new(core::f64::consts::PI.ln(), 0.0) + - log_sin_pi(z) + - log_gamma(Complex64::ONE - z); + } + let shifted = z + LANCZOS_G_MINUS_HALF; + Complex64::new(LOG_SQRT_TWO_PI, 0.0) + (z - 0.5) * shifted.ln() - shifted + lanczos_sum(z).ln() +} + +#[inline] +fn gamma(z: Complex64) -> Complex64 { + log_gamma(z).exp() +} + +#[inline] +fn reciprocal_gamma(z: Complex64) -> Complex64 { + if is_nonpositive_integer(z) { + Complex64::ZERO + } else { + (-log_gamma(z)).exp() + } +} + +fn gamma_ratio(numerator: &[Complex64], denominator: &[Complex64]) -> Complex64 { + let numerator_log = numerator + .iter() + .copied() + .map(log_gamma) + .fold(Complex64::ZERO, |sum, value| sum + value); + let denominator_log = denominator + .iter() + .copied() + .map(log_gamma) + .fold(Complex64::ZERO, |sum, value| sum + value); + (numerator_log - denominator_log).exp() +} + +fn pochhammer(z: Complex64, order: i32) -> Complex64 { + let mut result = Complex64::ONE; + if order >= 0 { + for n in 0..order { + result *= z + n as f64; + } + } else { + for n in order..0 { + result /= z + n as f64; + } + } + result +} + +fn digamma(mut z: Complex64) -> Complex64 { + if is_nonpositive_integer(z) { + return NAN; + } + // Shifting is especially accurate next to a pole, where the reflection + // formula subtracts a large cotangent term. Reserve reflection for inputs + // so far left that a long recurrence would be needlessly expensive. + if z.re < -64.0 { + return digamma(Complex64::ONE - z) - core::f64::consts::PI * cot_pi(z); + } + let mut result = Complex64::ZERO; + while z.re < 8.0 { + result -= Complex64::ONE / z; + z += 1.0; + } + let inverse = Complex64::ONE / z; + let inverse_squared = inverse * inverse; + let correction = 1.0 / 12.0 + + inverse_squared + * (-1.0 / 120.0 + + inverse_squared + * (1.0 / 252.0 + + inverse_squared * (-1.0 / 240.0 + inverse_squared * (5.0 / 660.0)))); + result + z.ln() - 0.5 * inverse - inverse_squared * correction +} + +fn lanczos_ratio(z: Complex64, epsilon: Complex64) -> Complex64 { + let mut numerator = Complex64::ZERO; + let mut denominator = Complex64::new(LANCZOS_COEFFICIENTS[0], 0.0); + for (index, coefficient) in LANCZOS_COEFFICIENTS.iter().enumerate().skip(1) { + let offset = (index - 1) as f64; + let inverse = Complex64::ONE / (z + offset); + numerator += coefficient * inverse / (z + epsilon + offset); + denominator += coefficient * inverse; + } + numerator / denominator +} + +/// Evaluates `(Gamma(z + epsilon) / Gamma(z) - 1) / epsilon` in a form that +/// remains finite and accurate as `epsilon` tends to zero. +fn log_gamma_difference_over_epsilon(z: Complex64, epsilon: Complex64) -> Complex64 { + let shifted = z + epsilon; + let base = z - 0.5; + let lanczos_argument = base + 4.742_187_5; + if z.re >= 0.5 { + if shifted == z { + return base / lanczos_argument + lanczos_argument.ln() + - 1.0 + - lanczos_ratio(z, epsilon); + } + let difference = base * complex_log1p(epsilon / lanczos_argument) + + epsilon * (lanczos_argument + epsilon).ln() + - epsilon + + complex_log1p(-epsilon * lanczos_ratio(z, epsilon)); + return complex_expm1(difference) / epsilon; + } + + let tangent = sin_pi(z) / cos_pi(z); + if shifted == z { + return log_gamma_difference_over_epsilon(Complex64::ONE - z, epsilon) + - core::f64::consts::PI / tangent; + } + let mut value = (cos_pi(epsilon) + sin_pi(epsilon) / tangent) + * log_gamma_difference_over_epsilon(Complex64::ONE - z, -epsilon) + + 0.5 * epsilon * (core::f64::consts::PI * sinc_pi(0.5 * epsilon)).powi(2) + - core::f64::consts::PI * sinc_pi(epsilon) / tangent; + value /= Complex64::ONE - epsilon * value; + value +} + +/// Evaluates `(1/Gamma(z) - 1/Gamma(z + epsilon)) / epsilon`, including the +/// removable limits at gamma poles and at `epsilon = 0`. +fn gamma_difference_ratio(z: Complex64, epsilon: Complex64) -> Complex64 { + let shifted = z + epsilon; + if complex_abs(epsilon) > 0.1 { + return (reciprocal_gamma(z) - reciprocal_gamma(shifted)) / epsilon; + } + if shifted == z { + if is_nonpositive_integer(z) { + let integer = z.re as i32; + let sign = if (integer + 1).rem_euclid(2) == 0 { + 1.0 + } else { + -1.0 + }; + return sign * gamma(Complex64::new((1 - integer) as f64, 0.0)); + } + return digamma(z) * reciprocal_gamma(z); + } + if is_nonpositive_integer(z) { + return -reciprocal_gamma(shifted) / epsilon; + } + if is_nonpositive_integer(shifted) { + return reciprocal_gamma(z) / epsilon; + } + let (z_integer, _) = nearest_integer_difference(z); + let (shifted_integer, _) = nearest_integer_difference(shifted); + if complex_abs(z + (z_integer.abs() as f64)) + < complex_abs(shifted + (shifted_integer.abs() as f64)) + { + log_gamma_difference_over_epsilon(z, epsilon) * reciprocal_gamma(shifted) + } else { + log_gamma_difference_over_epsilon(shifted, -epsilon) * reciprocal_gamma(z) + } +} + +fn pochhammer_difference_ratio(z: Complex64, epsilon: Complex64, order: i32) -> Complex64 { + debug_assert!(order >= 0); + if order == 0 { + return Complex64::ZERO; + } + let pole_index = -z.re.round() as i32; + // Rounding locates a possible zero factor, but the pole-limit formula is + // valid only when the base itself is exactly a nonpositive integer. + let contains_pole = is_nonpositive_integer(z) && (0..order).contains(&pole_index); + if epsilon == Complex64::ZERO { + if contains_pole { + let mut product = Complex64::ONE; + for index in 0..order { + if index != pole_index { + product *= z + index as f64; + } + } + return product; + } + let mut reciprocal_sum = Complex64::ZERO; + for index in 0..order { + reciprocal_sum += Complex64::ONE / (z + index as f64); + } + return pochhammer(z, order) * reciprocal_sum; + } + if contains_pole { + let mut shifted_product = Complex64::ONE; + let mut log_sum = Complex64::ZERO; + for index in 0..order { + if index != pole_index { + shifted_product *= z + epsilon + index as f64; + log_sum += complex_log1p(epsilon / (z + index as f64)); + } + } + return shifted_product + pochhammer(z, order) * complex_expm1(log_sum) / epsilon; + } + let mut log_sum = Complex64::ZERO; + for index in 0..order { + log_sum += complex_log1p(epsilon / (z + index as f64)); + } + pochhammer(z, order) * complex_expm1(log_sum) / epsilon +} + +#[inline] +fn exponential_difference_ratio(z: Complex64, epsilon: Complex64) -> Complex64 { + if epsilon == Complex64::ZERO { + z + } else { + complex_expm1(epsilon * z) / epsilon + } +} + +#[inline] +fn integer_sign(power: i32) -> f64 { + if power.rem_euclid(2) == 0 { 1.0 } else { -1.0 } +} + +// Stable connection formula about z = 1. Writing c - a - b = m + epsilon +// separates the finite polynomial part from the infinite tail. The paired +// beta/gamma recurrence evaluates their cancellation before epsilon reaches +// machine precision, rather than subtracting two singular connection terms. + +fn one_alpha_zero( + a: Complex64, + b: Complex64, + c: Complex64, + m: i32, + epsilon: Complex64, +) -> Complex64 { + if epsilon == Complex64::ZERO { + integer_sign(m) * gamma(Complex64::new(m as f64, 0.0)) * gamma(c) + / (gamma(a + m as f64) * gamma(b + m as f64)) + } else { + gamma(c) + / (epsilon + * gamma(1.0 - m as f64 - epsilon) + * gamma(a + m as f64 + epsilon) + * gamma(b + m as f64 + epsilon)) + } +} + +fn one_beta_zero( + a: Complex64, + b: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> Complex64 { + let mf = m as f64; + if complex_abs(epsilon) > 0.1 { + return (pochhammer(a, m) * pochhammer(b, m) + / (gamma(1.0 - epsilon) + * gamma(a + mf + epsilon) + * gamma(b + mf + epsilon) + * gamma(Complex64::new(mf + 1.0, 0.0))) + - complex_pow(w, epsilon) / (gamma(a) * gamma(b) * gamma(mf + 1.0 + epsilon))) + * gamma(c) + * complex_pow(w, Complex64::new(mf, 0.0)) + / epsilon; + } + ((gamma_difference_ratio(Complex64::ONE, -epsilon) / gamma(Complex64::new(mf + 1.0, 0.0)) + + gamma_difference_ratio(Complex64::new(mf + 1.0, 0.0), epsilon)) + / (gamma(a + mf + epsilon) * gamma(b + mf + epsilon)) + - (gamma_difference_ratio(a + mf, epsilon) / gamma(b + mf + epsilon) + + gamma_difference_ratio(b + mf, epsilon) / gamma(a + mf)) + / gamma(mf + 1.0 + epsilon) + - exponential_difference_ratio(w.ln(), epsilon) + / (gamma(a + mf) * gamma(b + mf) * gamma(mf + 1.0 + epsilon))) + * gamma(c) + * pochhammer(a, m) + * pochhammer(b, m) + * complex_pow(w, Complex64::new(mf, 0.0)) +} + +fn one_gamma_zero( + a: Complex64, + b: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> Complex64 { + let mf = m as f64; + gamma(c) * pochhammer(a, m) * pochhammer(b, m) * complex_pow(w, Complex64::new(mf, 0.0)) + / (gamma(a + mf + epsilon) + * gamma(b + mf + epsilon) + * gamma(Complex64::new(mf + 1.0, 0.0)) + * gamma(1.0 - epsilon)) +} + +/// Evaluates the finite portion shared by the `z = 1` and `z = infinity` +/// connection formulas after their formula-specific initial term is known. +fn connection_finite_part( + a: Complex64, + second: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, + mut term: Complex64, +) -> EvalOutcome { + debug_assert!(m > 0); + let mut sum = CompensatedSum::new(term); + for n in 0..(m - 1) { + let nf = n as f64; + let denominator = (nf + 1.0) * (1.0 - m as f64 - epsilon + nf); + if denominator == Complex64::ZERO { + return EvalOutcome::failure(); + } + term *= (a + nf) * (second + nf) * w / denominator; + sum.add(term); + if !finite(term) || !finite(sum.value()) { + return EvalOutcome::failure(); + } + } + EvalOutcome::success(sum.value()) +} + +/// Evaluates the infinite tail shared by both connection formulas after their +/// formula-specific beta and gamma seeds are known. +fn connection_infinite_part( + a: Complex64, + second: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, + mut beta: Complex64, + mut gamma_term: Complex64, +) -> EvalOutcome { + let mf = m as f64; + let mut sum = CompensatedSum::new(beta); + let mut small_terms = 0; + for n in 0..MAX_SERIES_ITERATIONS { + let nf = n as f64; + let amn = a + mf + nf; + let second_mn = second + mf + nf; + let shifted_a = amn + epsilon; + let shifted_second = second_mn + epsilon; + let denominator = (mf + nf + 1.0 + epsilon) * (nf + 1.0); + let correction_denominator = (mf + nf + 1.0 + epsilon) * (nf + 1.0 - epsilon); + if denominator == Complex64::ZERO || correction_denominator == Complex64::ZERO { + return EvalOutcome::failure(); + } + beta = shifted_a * shifted_second * w * beta / denominator + + (amn * second_mn / (mf + nf + 1.0) - amn - second_mn - epsilon + + shifted_a * shifted_second / (nf + 1.0)) + * gamma_term + / correction_denominator; + sum.add(beta); + gamma_term *= amn * second_mn * w / ((mf + nf + 1.0) * (nf + 1.0 - epsilon)); + if !finite(beta) || !finite(gamma_term) || !finite(sum.value()) { + return EvalOutcome::failure(); + } + let sum_abs = complex_abs(sum.value()); + if sum_abs > 0.0 && complex_abs(beta) <= REL_TOL * sum_abs { + small_terms += 1; + if small_terms >= 2 { + return EvalOutcome::success(sum.value()); + } + } else { + small_terms = 0; + } + } + EvalOutcome::failure() +} + +fn one_finite_part( + a: Complex64, + b: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> EvalOutcome { + if m <= 0 { + return EvalOutcome::success(Complex64::ZERO); + } + connection_finite_part(a, b, w, m, epsilon, one_alpha_zero(a, b, c, m, epsilon)) +} + +fn one_infinite_part( + a: Complex64, + b: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> EvalOutcome { + connection_infinite_part( + a, + b, + w, + m, + epsilon, + one_beta_zero(a, b, c, w, m, epsilon), + one_gamma_zero(a, b, c, w, m, epsilon) * w, + ) +} + +fn one_expansion(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> EvalOutcome { + let (m, epsilon) = nearest_integer_difference(c - a - b); + if m < 0 || m as usize > MAX_PARAMETER_ITERATIONS { + return EvalOutcome::failure(); + } + let w = one_minus(z); + let finite_part = one_finite_part(a, b, c, w, m, epsilon); + let infinite_part = one_infinite_part(a, b, c, w, m, epsilon); + if !finite_part.converged || !infinite_part.converged { + return EvalOutcome::failure(); + } + let value = integer_sign(m) * (finite_part.value + infinite_part.value) / sinc_pi(epsilon); + if finite(value) { + EvalOutcome::success(value) + } else { + EvalOutcome::failure() + } +} + +// Stable connection formula about z = infinity. This has the same +// m + epsilon construction as the z = 1 expansion, now for b - a, and is +// evaluated in the reciprocal coordinate w = 1 / z. + +fn infinity_alpha_zero(a: Complex64, c: Complex64, m: i32, epsilon: Complex64) -> Complex64 { + if epsilon == Complex64::ZERO { + integer_sign(m) * gamma(Complex64::new(m as f64, 0.0)) * gamma(c) + / (gamma(a + m as f64) * gamma(c - a)) + } else { + gamma(c) + / (epsilon + * gamma(1.0 - m as f64 - epsilon) + * gamma(a + m as f64 + epsilon) + * gamma(c - a)) + } +} + +fn infinity_beta_zero( + a: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> Complex64 { + let mf = m as f64; + let d = 1.0 - c + a; + if complex_abs(epsilon) > 0.1 { + return (pochhammer(a, m) * pochhammer(d, m) + / (gamma(1.0 - epsilon) + * gamma(a + mf + epsilon) + * gamma(c - a) + * gamma(Complex64::new(mf + 1.0, 0.0))) + - complex_pow(-w, epsilon) * pochhammer(d + epsilon, m) + / (gamma(a) * gamma(c - a - epsilon) * gamma(mf + 1.0 + epsilon))) + * gamma(c) + * complex_pow(w, Complex64::new(mf, 0.0)) + / epsilon; + } + ((pochhammer(d + epsilon, m) * gamma_difference_ratio(Complex64::ONE, -epsilon) + - pochhammer_difference_ratio(d, epsilon, m) / gamma(1.0 - epsilon)) + / (gamma(c - a) * gamma(a + mf + epsilon) * gamma(Complex64::new(mf + 1.0, 0.0))) + + pochhammer(d + epsilon, m) + * ((gamma_difference_ratio(Complex64::new(mf + 1.0, 0.0), epsilon) + / gamma(a + mf + epsilon) + - gamma_difference_ratio(a + mf, epsilon) / gamma(mf + 1.0 + epsilon)) + / gamma(c - a) + - (gamma_difference_ratio(c - a, -epsilon) + - exponential_difference_ratio(-(-w).ln(), -epsilon) / gamma(c - a - epsilon)) + / (gamma(mf + 1.0 + epsilon) * gamma(a + mf)))) + * gamma(c) + * pochhammer(a, m) + * complex_pow(w, Complex64::new(mf, 0.0)) +} + +fn infinity_gamma_zero( + a: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> Complex64 { + let mf = m as f64; + gamma(c) + * pochhammer(a, m) + * pochhammer(1.0 - c + a, m) + * complex_pow(w, Complex64::new(mf, 0.0)) + / (gamma(a + mf + epsilon) + * gamma(c - a) + * gamma(Complex64::new(mf + 1.0, 0.0)) + * gamma(1.0 - epsilon)) +} + +fn infinity_finite_part( + a: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> EvalOutcome { + if m <= 0 { + return EvalOutcome::success(Complex64::ZERO); + } + connection_finite_part( + a, + 1.0 - c + a, + w, + m, + epsilon, + infinity_alpha_zero(a, c, m, epsilon), + ) +} + +fn infinity_infinite_part( + a: Complex64, + c: Complex64, + w: Complex64, + m: i32, + epsilon: Complex64, +) -> EvalOutcome { + connection_infinite_part( + a, + 1.0 - c + a, + w, + m, + epsilon, + infinity_beta_zero(a, c, w, m, epsilon), + infinity_gamma_zero(a, c, w, m, epsilon) * w, + ) +} + +/// Evaluates the reciprocal-coordinate connection expansion, swapping `a` +/// and `b` first so that the integer part of `b - a` is nonnegative. +fn infinity_expansion( + mut a: Complex64, + mut b: Complex64, + c: Complex64, + z: Complex64, +) -> EvalOutcome { + if (b - a).re < 0.0 { + core::mem::swap(&mut a, &mut b); + } + let (m, epsilon) = nearest_integer_difference(b - a); + if m < 0 || m as usize > MAX_PARAMETER_ITERATIONS { + return EvalOutcome::failure(); + } + let w = complex_inverse(z); + let finite_part = infinity_finite_part(a, c, w, m, epsilon); + let infinite_part = infinity_infinite_part(a, c, w, m, epsilon); + if !finite_part.converged || !infinite_part.converged { + return EvalOutcome::failure(); + } + let value = integer_sign(m) * complex_pow(-w, a) * (finite_part.value + infinite_part.value) + / sinc_pi(epsilon); + if finite(value) { + EvalOutcome::success(value) + } else { + EvalOutcome::failure() + } +} + +// Direct and terminating series use compensated addition because intermediate +// terms can be much larger than the final answer for complex parameters. + +/// Kahan-style compensated accumulation applied componentwise by complex +/// arithmetic. +#[derive(Clone, Copy, Debug)] +struct CompensatedSum { + sum: Complex64, + correction: Complex64, +} + +impl CompensatedSum { + #[inline] + const fn new(value: Complex64) -> Self { + Self { + sum: value, + correction: Complex64::ZERO, + } + } + + #[inline] + fn add(&mut self, value: Complex64) { + let adjusted = value - self.correction; + let next = self.sum + adjusted; + self.correction = (next - self.sum) - adjusted; + self.sum = next; + } + + #[inline] + fn value(self) -> Complex64 { + self.sum + } +} + +#[inline] +fn negative_integer_degree(z: Complex64) -> Option { + if is_nonpositive_integer(z) && -z.re <= usize::MAX as f64 { + Some((-z.re) as usize) + } else { + None + } +} + +#[inline] +fn terminating_degree(a: Complex64, b: Complex64) -> Option { + match (negative_integer_degree(a), negative_integer_degree(b)) { + (Some(a_degree), Some(b_degree)) => Some(a_degree.min(b_degree)), + (Some(degree), None) | (None, Some(degree)) => Some(degree), + (None, None) => None, + } +} + +fn terminating_series( + a: Complex64, + b: Complex64, + c: Complex64, + z: Complex64, + degree: usize, +) -> EvalOutcome { + if degree > MAX_PARAMETER_ITERATIONS { + return EvalOutcome::failure(); + } + let mut term = Complex64::ONE; + let mut sum = CompensatedSum::new(Complex64::ONE); + for n in 0..degree { + let nf = n as f64; + let denominator = (c + nf) * (nf + 1.0); + if denominator == Complex64::ZERO { + return EvalOutcome::failure(); + } + term *= (a + nf) * (b + nf) * z / denominator; + // Once finite-precision arithmetic underflows a polynomial term to + // exact zero, every later recurrence term remains zero. This avoids + // walking enormous formal degrees for a result already determined. + if term == Complex64::ZERO { + return EvalOutcome::success(sum.value()); + } + sum.add(term); + if !finite(term) || !finite(sum.value()) { + return EvalOutcome::failure(); + } + } + EvalOutcome::success(sum.value()) +} + +/// Sums the defining Gauss series up to the supplied iteration limit. +#[inline] +fn direct_series_with_limit( + a: Complex64, + b: Complex64, + c: Complex64, + z: Complex64, + iteration_limit: usize, +) -> EvalOutcome { + let mut term = Complex64::ONE; + let mut sum = CompensatedSum::new(Complex64::ONE); + let mut small_terms = 0; + + for n in 0..iteration_limit { + let nf = n as f64; + let denominator = (c + nf) * (nf + 1.0); + if denominator == Complex64::ZERO { + return EvalOutcome::failure(); + } + term *= (a + nf) * (b + nf) * z / denominator; + sum.add(term); + let value = sum.value(); + let term_abs = complex_abs(term); + if !finite(term) || !finite(value) { + return EvalOutcome::failure(); + } + + let value_abs = complex_abs(value); + let converged = + term == Complex64::ZERO || (value_abs > 0.0 && term_abs <= REL_TOL * value_abs); + if converged { + small_terms += 1; + if small_terms >= 2 && n >= 1 { + return EvalOutcome::success(value); + } + } else { + small_terms = 0; + } + } + EvalOutcome::failure() +} + +#[inline] +fn direct_series(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> EvalOutcome { + direct_series_with_limit(a, b, c, z, MAX_SERIES_ITERATIONS) +} + +// Region selection and Taylor continuation. Each transformed candidate is +// scored by the modulus of its local series variable; the Taylor path fills +// the compact gap left when every candidate lies outside DIRECT_RADIUS. + +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +enum EvalPath { + Direct, + PfaffDirect, + Infinity, + PfaffInfinity, + One, + PfaffOne, + Taylor, +} + +/// Return every transformed path and the modulus of its local series variable. +fn path_candidates(z: Complex64) -> [(EvalPath, f64); 6] { + let one_minus_z = one_minus(z); + let inverse_z = complex_inverse(z); + let pfaff_z = pfaff_argument(z); + [ + (EvalPath::Direct, complex_abs(z)), + (EvalPath::PfaffDirect, complex_abs(pfaff_z)), + // Prefer the Pfaff form when the two infinity coordinates tie (notably + // at z = 2); its connection terms cancel less severely for large + // parameters in the mpmath reference grid. + (EvalPath::PfaffInfinity, complex_abs(one_minus(inverse_z))), + (EvalPath::Infinity, complex_abs(inverse_z)), + (EvalPath::One, complex_abs(one_minus_z)), + ( + EvalPath::PfaffOne, + complex_abs(complex_inverse(one_minus_z)), + ), + ] +} + +fn select_path(z: Complex64) -> EvalPath { + // The ordering is also the deterministic tie-break policy at region + // boundaries, which prevents small roundoff changes from switching paths. + let mut selected = EvalPath::Taylor; + let mut selected_modulus = f64::INFINITY; + for (path, modulus) in path_candidates(z) { + if modulus <= DIRECT_RADIUS { + let tie = 32.0 * f64::EPSILON * selected_modulus.max(modulus).max(1.0); + if selected == EvalPath::Taylor || modulus + tie < selected_modulus { + selected = path; + selected_modulus = modulus; + } + } + } + selected +} + +fn forced_anchor_evaluation(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> EvalOutcome { + // Anchor construction guarantees a non-Taylor path. Keep that condition + // explicit so a future selector change fails rather than recursing. + debug_assert_ne!(select_path(z), EvalPath::Taylor); + general_evaluation_impl(a, b, c, z, false) +} + +fn taylor_continuation_with_limit( + a: Complex64, + b: Complex64, + c: Complex64, + z: Complex64, + iteration_limit: usize, +) -> EvalOutcome { + let z_abs = complex_abs(z); + if z_abs == 0.0 { + return EvalOutcome::success(Complex64::ONE); + } + let preferred_radius = if z_abs < 1.0 { + TAYLOR_INNER_ANCHOR_RADIUS + } else { + TAYLOR_OUTER_ANCHOR_RADIUS + }; + let preferred_anchor = (preferred_radius / z_abs) * z; + // The outer anchor is close to the target and usually admits a transformed + // series. If it remains in the Taylor gap, move inside the direct disk; + // this guarantees a stable, nonrecursive selector path. + let z0 = if select_path(preferred_anchor) == EvalPath::Taylor { + (TAYLOR_INNER_ANCHOR_RADIUS / z_abs) * z + } else { + preferred_anchor + }; + let q0_outcome = forced_anchor_evaluation(a, b, c, z0); + let shifted_outcome = forced_anchor_evaluation(a + 1.0, b + 1.0, c + 1.0, z0); + if !q0_outcome.converged || !shifted_outcome.converged || c == Complex64::ZERO { + return EvalOutcome::failure(); + } + + // q0 and q1 are the function and its first derivative at z0. Higher + // derivatives follow from the hypergeometric differential equation. + let mut q0 = q0_outcome.value; + let mut q1 = a * b * shifted_outcome.value / c; + let delta = z - z0; + let mut delta_power = delta; + let mut sum = CompensatedSum::new(q0); + sum.add(q1 * delta); + let differential_denominator = z0 * one_minus(z0); + if differential_denominator == Complex64::ZERO || !finite(sum.value()) { + return EvalOutcome::failure(); + } + + let mut small_terms = 0; + for n in 0..iteration_limit { + let nf = n as f64; + let q2 = ((nf * (2.0 * z0 - 1.0) - c + (a + b + 1.0) * z0) * q1 + + (a + nf) * (b + nf) * q0 / (nf + 1.0)) + / (differential_denominator * (nf + 2.0)); + delta_power *= delta; + let term = q2 * delta_power; + sum.add(term); + let value = sum.value(); + if !finite(q2) || !finite(term) || !finite(value) { + return EvalOutcome::failure(); + } + let term_abs = complex_abs(term); + let value_abs = complex_abs(value); + if value_abs > 0.0 && term_abs <= REL_TOL * value_abs { + small_terms += 1; + if small_terms >= 2 { + return EvalOutcome::success(value); + } + } else { + small_terms = 0; + } + q0 = q1; + q1 = q2; + } + EvalOutcome::failure() +} + +#[inline] +fn taylor_continuation(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> EvalOutcome { + taylor_continuation_with_limit(a, b, c, z, MAX_TAYLOR_ITERATIONS) +} + +fn general_evaluation(a: Complex64, b: Complex64, c: Complex64, z: Complex64) -> EvalOutcome { + general_evaluation_impl(a, b, c, z, true) +} + +/// Evaluate one selector path, including its path-specific Pfaff prefactor. +fn evaluate_path( + path: EvalPath, + a: Complex64, + b: Complex64, + c: Complex64, + z: Complex64, + allow_taylor: bool, +) -> EvalOutcome { + let (result, prefactor) = match path { + EvalPath::Direct => (direct_series(a, b, c, z), Complex64::ONE), + EvalPath::PfaffDirect => ( + direct_series(a, c - b, c, pfaff_argument(z)), + (-a * complex_log1p(-z)).exp(), + ), + EvalPath::Infinity => (infinity_expansion(a, b, c, z), Complex64::ONE), + EvalPath::PfaffInfinity => ( + infinity_expansion(a, c - b, c, pfaff_argument(z)), + (-a * complex_log1p(-z)).exp(), + ), + EvalPath::One => (one_expansion(a, b, c, z), Complex64::ONE), + EvalPath::PfaffOne => ( + one_expansion(a, c - b, c, pfaff_argument(z)), + (-a * complex_log1p(-z)).exp(), + ), + EvalPath::Taylor if allow_taylor => (taylor_continuation(a, b, c, z), Complex64::ONE), + EvalPath::Taylor => (EvalOutcome::failure(), Complex64::ONE), + }; + if !result.converged { + return result; + } + let value = prefactor * result.value; + if finite(value) { + EvalOutcome::success(value) + } else { + EvalOutcome::failure() + } +} + +fn general_evaluation_impl( + mut a: Complex64, + mut b: Complex64, + c: Complex64, + z: Complex64, + allow_taylor: bool, +) -> EvalOutcome { + // Euler's transformation makes Re(c-a-b) nonnegative; symmetry in a and b + // then gives the stable parameter ordering assumed by the expansions. + let balance = c - a - b; + let euler_a = c - a; + let euler_b = c - b; + + // Euler's identity converts c-a and c-b into numerator parameters. When + // either is a nonpositive integer, evaluate the resulting polynomial + // before selecting a connection path: those paths contain gamma poles + // even though the complete hypergeometric value has a removable limit. + if let Some(degree) = terminating_degree(euler_a, euler_b) { + let result = terminating_series(euler_a, euler_b, c, z, degree); + if !result.converged { + return result; + } + let value = (balance * complex_log1p(-z)).exp() * result.value; + return if finite(value) { + EvalOutcome::success(value) + } else { + EvalOutcome::failure() + }; + } + + let mut outer_prefactor = Complex64::ONE; + if balance.re < 0.0 { + outer_prefactor = (balance * complex_log1p(-z)).exp(); + a = euler_a; + b = euler_b; + } + if (b - a).re < 0.0 { + core::mem::swap(&mut a, &mut b); + } + + let path = select_path(z); + let mut result = evaluate_path(path, a, b, c, z, allow_taylor); + if !result.converged { + // A small local coordinate predicts convergence, not numerical + // success. If the preferred expansion encounters a removable pole or + // unstable seed, try every other admissible series before reporting + // failure. The normal successful path still pays no retry cost. + for (fallback, modulus) in path_candidates(z) { + if fallback != path && modulus <= DIRECT_RADIUS { + result = evaluate_path(fallback, a, b, c, z, allow_taylor); + if result.converged { + break; + } + } + } + } + if !result.converged { + return result; + } + let value = outer_prefactor * result.value; + if finite(value) { + EvalOutcome::success(value) + } else { + EvalOutcome::failure() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn assert_close(actual: Complex64, expected: Complex64, tolerance: f64) { + let scale = complex_abs(expected).max(1.0); + assert!( + complex_abs(actual - expected) <= tolerance * scale, + "actual={actual:?}, expected={expected:?}, tolerance={tolerance}" + ); + } + + fn parse_gamma_fixture() -> Vec<[f64; 8]> { + let fixture = include_str!("../../test/data/complex_gamma_reference.csv"); + let mut rows = Vec::new(); + for (line_index, line) in fixture.lines().enumerate() { + if line.starts_with('#') || line_index == 1 || line.is_empty() { + continue; + } + let fields: Vec<_> = line.split(',').collect(); + assert_eq!( + fields.len(), + 8, + "complex_gamma_reference.csv row {} has wrong field count", + line_index + 1 + ); + let mut row = [0.0; 8]; + for (field_index, field) in fields.iter().enumerate() { + row[field_index] = field.parse().unwrap_or_else(|error| { + panic!( + "complex_gamma_reference.csv row {}, field {} is not f64: {error}", + line_index + 1, + field_index + 1 + ) + }); + } + rows.push(row); + } + rows + } + + #[derive(Debug)] + struct HypFixtureRow<'a> { + a: Complex64, + b: Complex64, + c: Complex64, + z: Complex64, + expected: Complex64, + tolerance: f64, + label: &'a str, + } + + fn parse_hyp_fixture() -> Vec> { + let fixture = include_str!("../../test/data/hyp2f1_reference.csv"); + let mut rows = Vec::new(); + for (line_index, line) in fixture.lines().enumerate() { + if line.starts_with('#') || line_index == 1 || line.is_empty() { + continue; + } + let fields: Vec<_> = line.split(',').collect(); + assert_eq!( + fields.len(), + 12, + "hyp2f1_reference.csv row {} has wrong field count", + line_index + 1 + ); + let mut values = [0.0; 11]; + for (field_index, field) in fields[..11].iter().enumerate() { + values[field_index] = field.parse().unwrap_or_else(|error| { + panic!( + "hyp2f1_reference.csv row {}, field {} is not f64: {error}", + line_index + 1, + field_index + 1 + ) + }); + } + rows.push(HypFixtureRow { + a: Complex64::new(values[0], values[1]), + b: Complex64::new(values[2], values[3]), + c: Complex64::new(values[4], values[5]), + z: Complex64::new(values[6], values[7]), + expected: Complex64::new(values[8], values[9]), + tolerance: values[10], + label: fields[11], + }); + } + rows + } + + #[test] + fn complex_gamma_and_digamma_match_reference() { + for [ + z_re, + z_im, + gamma_re, + gamma_im, + rgamma_re, + rgamma_im, + digamma_re, + digamma_im, + ] in parse_gamma_fixture() + { + let z = Complex64::new(z_re, z_im); + assert_close(gamma(z), Complex64::new(gamma_re, gamma_im), 8e-13); + assert_close( + reciprocal_gamma(z), + Complex64::new(rgamma_re, rgamma_im), + 8e-13, + ); + assert_close(digamma(z), Complex64::new(digamma_re, digamma_im), 8e-13); + } + } + + #[test] + fn elementary_helpers_preserve_small_complex_increments() { + let z = Complex64::new(1e-12, -2e-12); + assert_close(complex_log1p(z).exp() - 1.0, z, 2e-16); + assert_close(complex_expm1(z), z, 2e-16); + assert_close(sinc_pi(Complex64::ZERO), Complex64::ONE, 0.0); + assert_close( + complex_pow(Complex64::new(-2.0, 0.0), Complex64::new(0.5, 0.0)), + Complex64::new(0.0, 2.0_f64.sqrt()), + 2e-15, + ); + + let upper_cut = Complex64::new(1.3, 0.0); + let lower_cut = Complex64::new(1.3, -0.0); + assert!(one_minus(upper_cut).im.is_sign_negative()); + assert!(one_minus(lower_cut).im.is_sign_positive()); + assert!(pfaff_argument(upper_cut).im.is_sign_negative()); + assert!(pfaff_argument(lower_cut).im.is_sign_positive()); + } + + #[test] + fn gamma_helpers_obey_identities() { + let z = Complex64::new(-0.3, 0.7); + assert_close(gamma(z + 1.0), z * gamma(z), 2e-13); + assert_close( + gamma(z) * gamma(1.0 - z), + core::f64::consts::PI / sin_pi(z), + 3e-13, + ); + assert_close(gamma_ratio(&[z + 1.0], &[z]), z, 3e-13); + assert_eq!(reciprocal_gamma(Complex64::new(-4.0, 0.0)), Complex64::ZERO); + } + + #[test] + fn stabilized_difference_helpers_have_finite_zero_limits() { + let z = Complex64::new(1.2, -0.4); + let tiny = Complex64::new(0.0, 1e-9); + assert_close( + gamma_difference_ratio(Complex64::ONE, -tiny), + Complex64::new(-0.577_215_664_901_532_9, -6.558_780_715_202_539e-10), + 2e-14, + ); + assert_close( + gamma_difference_ratio(Complex64::new(3.0, 0.0), tiny), + Complex64::new(0.461_392_167_549_233_57, -1.141_492_155_637_234e-10), + 2e-14, + ); + assert_close( + gamma_difference_ratio(z, Complex64::ZERO), + digamma(z) * reciprocal_gamma(z), + 2e-14, + ); + assert_close( + pochhammer_difference_ratio(z, Complex64::ZERO, 3), + 3.0 * z * z + 6.0 * z + 2.0, + 2e-14, + ); + assert_close( + pochhammer_difference_ratio(Complex64::new(-0.5, 0.0), Complex64::ZERO, 2), + Complex64::ZERO, + 2e-14, + ); + // Rounding -1.3 locates the nearby factor at index 1, but it is not an + // exact zero. The zero-epsilon limit is the full derivative of (z)_3. + assert_close( + pochhammer_difference_ratio(Complex64::new(-1.3, 0.0), Complex64::ZERO, 3), + Complex64::new(-0.73, 0.0), + 2e-14, + ); + assert_eq!(exponential_difference_ratio(z, Complex64::ZERO), z); + } + + #[test] + fn transformed_terminating_and_exact_integer_infinity_cases_match_mpmath() { + assert_close( + hyp2f1_scalar( + Complex64::new(3.0, 0.0), + Complex64::new(0.5, 0.0), + Complex64::new(2.0, 0.0), + Complex64::new(0.95, 0.02), + ), + Complex64::new(20.007_751_880_108_815, 11.419_384_341_915_018), + 3e-14, + ); + assert_close( + hyp2f1_scalar( + Complex64::new(-0.5, 0.0), + Complex64::new(-0.5, 0.0), + Complex64::ONE, + Complex64::new(1.9, 0.0), + ), + Complex64::new(1.566_746_598_608_614_2, 0.060_785_011_674_705_52), + 3e-14, + ); + } + + #[test] + fn real_parameter_infinity_cases_near_pochhammer_zeros_match_mpmath() { + // These cases exercise the infinity expansion with b-a an integer and + // 1-c+a near, but not equal to, a nonpositive integer. Reference values + // were generated with mpmath 1.3.0 at 50 decimal digits. + let cases = [ + ( + 0.3, + 2.3, + 0.9, + Complex64::new(4.0, 0.5), + Complex64::new(0.247_877_616_659_031_27, 0.297_556_752_826_444_6), + ), + ( + 0.5, + 2.5, + 1.2, + Complex64::new(5.0, 1.0), + Complex64::new(0.024_577_322_512_199_99, 0.222_224_334_815_688_8), + ), + ( + 0.3, + 3.3, + 0.9, + Complex64::new(4.0, 0.5), + Complex64::new(0.213_813_644_541_913_16, 0.277_550_319_435_210_5), + ), + ]; + + for (a, b, c, z, expected) in cases { + assert_close( + hyp2f1_scalar( + Complex64::new(a, 0.0), + Complex64::new(b, 0.0), + Complex64::new(c, 0.0), + z, + ), + expected, + 3e-14, + ); + } + } + + #[test] + fn positive_and_negative_integer_connection_edges_match_mpmath() { + // Both c-a and c-b are negative integers in the reported removable- + // pole case. + assert_close( + hyp2f1_scalar( + Complex64::new(-2.2, 0.3), + Complex64::new(-2.2, 0.3), + Complex64::new(-3.2, 0.3), + Complex64::new(2.0, 0.5), + ), + Complex64::new(-0.189_614_379_742_522_02, 0.129_705_129_001_328_46), + 3e-13, + ); + + let c = Complex64::new(-3.25, 0.25); + let z = Complex64::new(2.0, 0.5); + assert_close( + hyp2f1_scalar( + Complex64::new(-2.25, 0.25), + Complex64::new(-2.25, 0.25), + c, + z, + ), + Complex64::new(-0.201_770_442_240_660_46, 0.169_046_932_792_871_47), + 3e-13, + ); + + // Positive integer differences must continue through the ordinary + // connection machinery rather than being mistaken for gamma poles. + assert_close( + hyp2f1_scalar( + Complex64::new(-4.25, 0.25), + Complex64::new(-4.25, 0.25), + c, + z, + ), + Complex64::new(51.557_704_317_178_434, 269.233_357_645_869_8), + 3e-13, + ); + + assert_eq!( + hyp2f1_scalar( + Complex64::new(-2.25, 0.25), + Complex64::new(-2.25, 0.25), + c, + Complex64::ONE, + ), + Complex64::ZERO + ); + } + + #[test] + fn terminating_series_stops_after_term_underflows_to_zero() { + assert_eq!( + hyp2f1_scalar( + Complex64::new(-(MAX_PARAMETER_ITERATIONS as f64), 0.0), + Complex64::ONE, + Complex64::ONE, + Complex64::new(1e-300, 0.0), + ), + Complex64::ONE + ); + } + + #[test] + fn parameter_magnitude_iteration_caps_report_failure() { + let excessive = MAX_PARAMETER_ITERATIONS + 1; + assert!( + !terminating_series( + Complex64::new(-(excessive as f64), 0.0), + Complex64::ONE, + Complex64::ONE, + Complex64::new(1e-300, 0.0), + excessive, + ) + .converged + ); + + let a = Complex64::new(0.3, 0.0); + let b = Complex64::new(0.4, 0.0); + let excessive = excessive as f64; + assert!(!one_expansion(a, b, a + b + excessive, Complex64::new(0.5, 0.1)).converged); + assert!( + !infinity_expansion( + a, + a + excessive, + Complex64::new(0.9, 0.0), + Complex64::new(4.0, 0.5), + ) + .converged + ); + + let public_value = hyp2f1_scalar( + Complex64::new(-1_000_000_000.0, 0.0), + Complex64::ONE, + Complex64::ONE, + Complex64::new(1e-300, 0.0), + ); + assert!(public_value.re.is_nan() && public_value.im.is_nan()); + } + + #[test] + fn direct_and_polynomial_fixtures_match_reference() { + for row in parse_hyp_fixture() + .into_iter() + .filter(|row| matches!(row.label, "direct" | "polynomial")) + { + assert_close( + hyp2f1_scalar(row.a, row.b, row.c, row.z), + row.expected, + row.tolerance, + ); + } + } + + #[test] + fn transformed_fixtures_match_reference() { + for row in parse_hyp_fixture().into_iter().filter(|row| { + matches!( + row.label, + "pfaff" + | "infinity" + | "infinity-equal-m0" + | "one" + | "one-balanced-m0" + | "one-near-integer" + | "infinity-near-integer" + | "euler" + | "moderate" + | "scipy-1561" + | "upper-cut" + | "lower-cut" + | "pfaff-upper-cut" + | "pfaff-lower-cut" + ) + }) { + assert_close( + hyp2f1_scalar(row.a, row.b, row.c, row.z), + row.expected, + row.tolerance, + ); + } + } + + #[test] + fn taylor_gap_fixtures_match_reference() { + for row in parse_hyp_fixture() + .into_iter() + .filter(|row| row.label.starts_with("taylor")) + { + assert_eq!(select_path(row.z), EvalPath::Taylor); + assert_close( + hyp2f1_scalar(row.a, row.b, row.c, row.z), + row.expected, + row.tolerance, + ); + } + } + + #[test] + fn broad_mpmath_grid_matches_reference() { + let rows: Vec<_> = parse_hyp_fixture() + .into_iter() + .filter(|row| row.label.starts_with("grid-")) + .collect(); + assert_eq!(rows.len(), 406); + let mut failures = Vec::new(); + for row in rows { + let actual = hyp2f1_scalar(row.a, row.b, row.c, row.z); + let scale = complex_abs(row.expected).max(1.0); + let relative_error = complex_abs(actual - row.expected) / scale; + if relative_error > row.tolerance { + failures.push(format!( + "{}: actual={:?}, expected={:?}, error={}, tolerance={}", + row.label, actual, row.expected, relative_error, row.tolerance + )); + } + } + assert!(failures.is_empty(), "{}", failures.join("\n")); + } + + #[test] + fn taylor_cap_and_differential_equation_residual() { + let a = Complex64::new(0.7, 0.2); + let b = Complex64::new(1.2, -0.3); + let c = Complex64::new(2.1, 0.1); + let z = Complex64::new(0.5, 0.866_025_403_784_438_6); + assert!(!taylor_continuation_with_limit(a, b, c, z, 1).converged); + + let step = 1e-4; + let center = hyp2f1_scalar(a, b, c, z); + let plus = hyp2f1_scalar(a, b, c, z + step); + let minus = hyp2f1_scalar(a, b, c, z - step); + let first = (plus - minus) / (2.0 * step); + let second = (plus - 2.0 * center + minus) / (step * step); + let residual = z * (1.0 - z) * second + (c - (a + b + 1.0) * z) * first - a * b * center; + let scale = complex_abs(a * b * center).max(1.0); + assert!( + complex_abs(residual) <= 5e-7 * scale, + "residual={residual:?}" + ); + } + + #[test] + fn every_transformed_selector_path_is_reachable() { + assert_eq!(select_path(Complex64::new(0.2, 0.1)), EvalPath::Direct); + assert_eq!( + select_path(Complex64::new(-0.5, 0.1)), + EvalPath::PfaffDirect + ); + assert_eq!(select_path(Complex64::new(3.0, 4.0)), EvalPath::Infinity); + assert_eq!( + select_path(Complex64::new(1.2, 0.4)), + EvalPath::PfaffInfinity + ); + assert_eq!(select_path(Complex64::new(0.95, 0.05)), EvalPath::One); + assert_eq!(select_path(Complex64::new(-3.0, 0.4)), EvalPath::PfaffOne); + } + + #[test] + fn failed_selected_connection_path_retries_admissible_direct_series() { + let a = Complex64::new(0.2, 0.0); + let b = Complex64::new(0.3, 0.0); + let c = Complex64::new(48.5, 0.0); + let z = Complex64::new(0.6, 0.0); + + assert_eq!(select_path(z), EvalPath::One); + assert!(!one_expansion(a, b, c, z).converged); + assert!(direct_series(a, b, c, z).converged); + assert_close( + hyp2f1_scalar(a, b, c, z), + Complex64::new(1.000_749_430_996_130_3, 0.0), + 3e-14, + ); + } + + #[test] + fn symmetry_euler_and_pfaff_identities_agree() { + let a = Complex64::new(0.4, 0.2); + let b = Complex64::new(0.9, -0.1); + let c = Complex64::new(2.4, 0.3); + let z = Complex64::new(-0.7, 0.2); + let value = hyp2f1_scalar(a, b, c, z); + assert_close(hyp2f1_scalar(b, a, c, z), value, 5e-13); + + let euler = ((c - a - b) * complex_log1p(-z)).exp() * hyp2f1_scalar(c - a, c - b, c, z); + assert_close(euler, value, 2e-12); + + let transformed_z = pfaff_argument(z); + let pfaff = (-a * complex_log1p(-z)).exp() * hyp2f1_scalar(a, c - b, c, transformed_z); + assert_close(pfaff, value, 2e-12); + } + + #[test] + fn exceptional_cases_follow_documented_precedence() { + let a = Complex64::new(1.2, 0.3); + let b = Complex64::new(0.7, -0.2); + assert_eq!(hyp2f1_scalar(a, b, a, Complex64::ZERO), Complex64::ONE); + assert_close( + hyp2f1_scalar(a, b, a, Complex64::new(0.2, -0.1)), + (-b * complex_log1p(Complex64::new(-0.2, 0.1))).exp(), + 2e-14, + ); + + let valid = hyp2f1_scalar( + Complex64::new(-2.0, 0.0), + b, + Complex64::new(-2.0, 0.0), + Complex64::new(2.0, 0.5), + ); + assert!(finite(valid)); + let invalid = hyp2f1_scalar( + Complex64::new(-3.0, 0.0), + b, + Complex64::new(-2.0, 0.0), + Complex64::new(0.2, 0.0), + ); + assert!(invalid.re.is_nan() && invalid.im.is_nan()); + + let non_finite = hyp2f1_scalar(Complex64::new(f64::INFINITY, 0.0), b, a, Complex64::ZERO); + assert!(non_finite.re.is_nan() && non_finite.im.is_nan()); + } + + #[test] + fn argument_unity_uses_gauss_ratio() { + let a = Complex64::new(0.2, 0.1); + let b = Complex64::new(0.3, -0.2); + let c = Complex64::new(2.0, 0.4); + let expected = gamma_ratio(&[c, c - a - b], &[c - a, c - b]); + assert_close(hyp2f1_scalar(a, b, c, Complex64::ONE), expected, 2e-14); + let divergent = hyp2f1_scalar(a, b, a + b, Complex64::ONE); + assert!(divergent.re.is_nan() && divergent.im.is_nan()); + } + + #[test] + fn direct_series_iteration_cap_reports_failure() { + let result = direct_series_with_limit( + Complex64::new(0.7, 0.2), + Complex64::new(1.1, -0.3), + Complex64::new(2.4, 0.1), + Complex64::new(0.8, 0.1), + 1, + ); + assert!(!result.converged); + assert!(result.value.re.is_nan() && result.value.im.is_nan()); + } + + #[test] + fn zero_argument_is_one() { + let value = hyp2f1_scalar( + Complex64::new(1.0, 2.0), + Complex64::new(3.0, -1.0), + Complex64::new(4.0, 0.5), + Complex64::new(0.0, -0.0), + ); + assert_eq!(value, Complex64::ONE); + } + + #[test] + fn length_mismatch_is_transactional() { + let input = [Complex64::ONE]; + let mut out = [Complex64::new(7.0, 8.0); 2]; + assert_eq!( + hyp2f1(&input, &input, &input, &input, &mut out), + Err("Length mismatch") + ); + assert_eq!(out, [Complex64::new(7.0, 8.0); 2]); + + assert_eq!( + hyp2f1_par(&input, &input, &input, &input, &mut out), + Err("Length mismatch") + ); + assert_eq!(out, [Complex64::new(7.0, 8.0); 2]); + } + + #[test] + fn vector_apis_match_scalar_across_regions() { + let rows: Vec<_> = parse_hyp_fixture().into_iter().take(12).collect(); + let a: Vec<_> = rows.iter().map(|row| row.a).collect(); + let b: Vec<_> = rows.iter().map(|row| row.b).collect(); + let c: Vec<_> = rows.iter().map(|row| row.c).collect(); + let z: Vec<_> = rows.iter().map(|row| row.z).collect(); + let expected: Vec<_> = (0..rows.len()) + .map(|index| hyp2f1_scalar(a[index], b[index], c[index], z[index])) + .collect(); + let mut serial = vec![NAN; rows.len()]; + let mut parallel = vec![NAN; rows.len()]; + hyp2f1(&a, &b, &c, &z, &mut serial).unwrap(); + hyp2f1_par(&a, &b, &c, &z, &mut parallel).unwrap(); + assert_eq!(serial, expected); + assert_eq!(parallel, expected); + } + + #[test] + fn vector_apis_cover_empty_singleton_subslice_and_nan() { + let empty: [Complex64; 0] = []; + let mut empty_out = []; + assert_eq!( + hyp2f1(&empty, &empty, &empty, &empty, &mut empty_out), + Ok(()) + ); + assert_eq!( + hyp2f1_par(&empty, &empty, &empty, &empty, &mut empty_out), + Ok(()) + ); + + let a = [Complex64::new(0.5, 0.2); 3]; + let b = [Complex64::new(1.1, -0.3); 3]; + let c = [Complex64::new(2.4, 0.1); 3]; + let z = [ + Complex64::new(99.0, 0.0), + Complex64::new(0.3, 0.2), + Complex64::new(f64::NAN, 0.0), + ]; + let mut out = [Complex64::new(7.0, 8.0); 3]; + hyp2f1(&a[1..], &b[1..], &c[1..], &z[1..], &mut out[1..]).unwrap(); + assert_eq!(out[0], Complex64::new(7.0, 8.0)); + assert_eq!(out[1], hyp2f1_scalar(a[1], b[1], c[1], z[1])); + assert!(out[2].re.is_nan() && out[2].im.is_nan()); + + let mut singleton = [NAN]; + hyp2f1_par(&a[1..2], &b[1..2], &c[1..2], &z[1..2], &mut singleton).unwrap(); + assert_eq!(singleton[0], out[1]); + } +} diff --git a/src/python.rs b/src/python.rs index 179880e..282b5d2 100644 --- a/src/python.rs +++ b/src/python.rs @@ -1,10 +1,14 @@ use numpy::Element as NumpyElement; -use numpy::PyArray1; -use numpy::borrow::{PyReadonlyArray1, PyReadonlyArray2, PyReadonlyArray3, PyReadwriteArray1}; +use numpy::borrow::{ + PyReadonlyArray1, PyReadonlyArray2, PyReadonlyArray3, PyReadonlyArrayDyn, PyReadwriteArray1, + PyReadwriteArrayDyn, +}; +use numpy::ndarray::{ArrayView0, ArrayViewD, IxDyn, Zip}; +use numpy::{Complex64, PyArray1, PyArrayDyn, PyArrayMethods, PyUntypedArray}; use pyo3::create_exception; use pyo3::exceptions; use pyo3::prelude::*; -use pyo3::types::{PyDict, PyTuple}; +use pyo3::types::{PyComplex, PyComplexMethods, PyDict, PyTuple}; use std::ffi::CString; use std::fmt::Debug; @@ -3315,6 +3319,190 @@ fn ellipk(x: f64) -> f64 { math::ellipk(x) } +/// A borrowed NumPy `complex128` array or a scalar complex value. +/// +/// Scalars become zero-dimensional ndarray views so that `Zip` can broadcast +/// them without allocating expanded parameter arrays. +enum Hyp2f1Input<'py> { + Array(PyReadonlyArrayDyn<'py, Complex64>), + Scalar(Complex64), +} + +impl<'py> Hyp2f1Input<'py> { + /// Extract a strictly complex scalar or `complex128` NumPy array. + fn extract(value: &Bound<'py, PyAny>, name: &str) -> PyResult { + if let Ok(array) = value.cast::>() { + return Ok(Self::Array(array.try_readonly()?)); + } + if value.cast::().is_ok() { + return Err(exceptions::PyTypeError::new_err(format!( + "{name} must have dtype complex128" + ))); + } + if let Ok(value) = value.cast::() { + return Ok(Self::Scalar(Complex64::new(value.real(), value.imag()))); + } + Err(exceptions::PyTypeError::new_err(format!( + "{name} must be a complex scalar or complex128 ndarray" + ))) + } + + /// Return the shape that this input requires, ignoring scalar inputs. + fn nonscalar_shape(&self) -> Option<&[usize]> { + use numpy::PyUntypedArrayMethods as _; + + match self { + Self::Array(array) if array.ndim() != 0 => Some(array.shape()), + Self::Array(_) | Self::Scalar(_) => None, + } + } + + /// Return an ndarray view suitable for `Zip::and_broadcast`. + fn view(&self) -> ArrayViewD<'_, Complex64> { + match self { + Self::Array(array) => array.as_array(), + Self::Scalar(value) => ArrayView0::from_shape((), std::slice::from_ref(value)) + .expect("one value always forms a scalar view") + .into_dyn(), + } + } +} + +/// Determine and validate the exact elementwise output shape. +/// +/// Non-scalar inputs must have identical shapes. If every input is scalar, a +/// supplied output defines the shape; otherwise the result is zero-dimensional. +fn hyp2f1_output_shape( + inputs: [&Hyp2f1Input<'_>; 4], + out_shape: Option<&[usize]>, +) -> PyResult> { + let input_shape = inputs.iter().find_map(|input| input.nonscalar_shape()); + let shape = input_shape.or(out_shape).unwrap_or_default(); + + if inputs + .iter() + .filter_map(|input| input.nonscalar_shape()) + .any(|input_shape| input_shape != shape) + || out_shape.is_some_and(|out_shape| out_shape != shape) + { + return Err(exceptions::PyValueError::new_err( + "all non-scalar inputs and out must have the same shape", + )); + } + Ok(shape.to_vec()) +} + +/// Conservatively detect whether an array's logical elements may overlap. +/// +/// Axes are considered from smallest to largest absolute byte stride. Each +/// successive axis must step beyond the full byte span reachable through the +/// preceding axes, including the width of one element. Strides which are not a +/// whole number of elements are also rejected because ndarray cannot represent +/// them exactly. A `false` result therefore guarantees non-overlap; unusual +/// injective layouts may conservatively return `true`. +fn array_may_self_overlap(shape: &[usize], strides: &[isize]) -> bool { + if shape.contains(&0) { + return false; + } + + let mut axes: Vec<_> = (0..shape.len()).filter(|&axis| shape[axis] > 1).collect(); + axes.sort_unstable_by_key(|&axis| strides[axis].unsigned_abs()); + + let element_size = std::mem::size_of::(); + let mut preceding_span = element_size - 1; + for axis in axes { + let stride = strides[axis].unsigned_abs(); + if stride % element_size != 0 || stride <= preceding_span { + return true; + } + preceding_span = preceding_span.saturating_add((shape[axis] - 1).saturating_mul(stride)); + } + false +} + +/// Fill an output array using ndarray's shape-aware serial or parallel `Zip`. +fn hyp2f1_fill( + out: &mut PyReadwriteArrayDyn<'_, Complex64>, + inputs: [&Hyp2f1Input<'_>; 4], + par: bool, +) { + let zip = Zip::from(out.as_array_mut()) + .and_broadcast(inputs[0].view()) + .and_broadcast(inputs[1].view()) + .and_broadcast(inputs[2].view()) + .and_broadcast(inputs[3].view()); + let evaluate = + |out: &mut Complex64, &a: &Complex64, &b: &Complex64, &c: &Complex64, &z: &Complex64| { + *out = math::hyp2f1_scalar(a, b, c, z); + }; + if par { + zip.par_for_each(evaluate); + } else { + zip.for_each(evaluate); + } +} + +/// Evaluate Gauss's hypergeometric function elementwise on its principal branch. +/// +/// Each input may be a complex scalar or an arbitrarily strided NumPy +/// `complex128` array of any dimensionality. All non-scalar arrays must have the +/// same shape; scalars are broadcast across that shape without being expanded. +/// Values on the cut `[1, +inf)` distinguish the sign of zero in `z.imag`; +/// mathematical singularities and unsupported numerical failures produce +/// complex NaN values. Accuracy is not guaranteed uniformly for unbounded +/// parameter magnitudes. +/// +/// Pass a writable `complex128` array with the result shape as `out` to reuse +/// its storage. The output may be strided but must not contain overlapping +/// elements, and the same array object is returned. If `out` is `None`, a new +/// C-contiguous output is allocated. If every input is scalar, `out` may define +/// any result shape; without `out`, the result is a zero-dimensional array. +/// `out` must not alias an input array. +#[pyfunction(signature = (a, b, c, z, par = true, *, out = None))] +fn hyp2f1( + py: Python<'_>, + a: &Bound<'_, PyAny>, + b: &Bound<'_, PyAny>, + c: &Bound<'_, PyAny>, + z: &Bound<'_, PyAny>, + par: bool, + out: Option>>, +) -> PyResult>> { + use numpy::PyUntypedArrayMethods as _; + + let a = Hyp2f1Input::extract(a, "a")?; + let b = Hyp2f1Input::extract(b, "b")?; + let c = Hyp2f1Input::extract(c, "c")?; + let z = Hyp2f1Input::extract(z, "z")?; + let inputs = [&a, &b, &c, &z]; + let shape = hyp2f1_output_shape(inputs, out.as_ref().map(|out| out.shape()))?; + + if out + .as_ref() + .is_some_and(|out| array_may_self_overlap(out.shape(), out.strides())) + { + return Err(exceptions::PyValueError::new_err( + "out must have a non-overlapping memory layout", + )); + } + + match out { + Some(out) => { + let returned = out.clone().unbind(); + let mut out = out.try_into_readwrite()?; + hyp2f1_fill(&mut out, inputs, par); + Ok(returned) + } + None => { + let out = PyArrayDyn::::zeros(py, IxDyn(&shape), false); + let returned = out.clone().unbind(); + let mut out = out.try_into_readwrite()?; + hyp2f1_fill(&mut out, inputs, par); + Ok(returned) + } + } +} + /// Python bindings for cfsemrs::physics::flux_density_circular_filament_cartesian #[pyfunction] fn flux_density_circular_filament_cartesian( @@ -4420,6 +4608,7 @@ fn _cfsem<'py>(_py: Python, m: Bound<'py, PyModule>) -> PyResult<()> { // Pure math m.add_function(wrap_pyfunction!(ellipe, m.clone())?)?; m.add_function(wrap_pyfunction!(ellipk, m.clone())?)?; + m.add_function(wrap_pyfunction!(hyp2f1, m.clone())?)?; // Filamentization and meshing m.add_function(wrap_pyfunction!(filament_helix_path, m.clone())?)?; diff --git a/test/data/complex_gamma_reference.csv b/test/data/complex_gamma_reference.csv new file mode 100644 index 0000000..f86b401 --- /dev/null +++ b/test/data/complex_gamma_reference.csv @@ -0,0 +1,11 @@ +# mpmath=1.3.0, dps=100 +z_re,z_im,gamma_re,gamma_im,rgamma_re,rgamma_im,digamma_re,digamma_im +0.2,0.3,1.1707421186241772670427731995774877133177611439867350443841974174626524384821839,-2.1041380778637433694988449377477141897331956013146054753917902281485455711335231,0.20192052797748195097997861960073333432897257868365447879048477159850271449250645,0.36290542969365843686363859945096373031118527496107017357749980548549908515182102,-1.7642419236812976318623888848418710621474451952380796330088773637603439639539229,2.6740928489238811571516103952228291468186468933642911671012956611784263187535521 +0.2,-0.3,1.1707421186241772670427731995774877133177611439867350443841974174626524384821839,2.1041380778637433694988449377477141897331956013146054753917902281485455711335231,0.20192052797748195097997861960073333432897257868365447879048477159850271449250645,-0.36290542969365843686363859945096373031118527496107017357749980548549908515182102,-1.7642419236812976318623888848418710621474451952380796330088773637603439639539229,-2.6740928489238811571516103952228291468186468933642911671012956611784263187535521 +1.2,-2.5,0.065026022252581368736450223419630120216777327145520688506102029368201225794076034,-0.068180793591090160444601318790546650859071106567457486131251354112525226500403309,7.3252215386613878758222749343008917952282654080722927118375271107572566252259719,7.6806084769648273516285361071398752118988209313008859683388836839576457132898973,0.94867561056762661571185627826033808315242932467325841790867317211624938841090312,-1.2944207545567560014379368281306796768911683993526690757981205807340999500235618 +1.2,2.5,0.065026022252581368736450223419630120216777327145520688506102029368201225794076034,0.068180793591090160444601318790546650859071106567457486131251354112525226500403309,7.3252215386613878758222749343008917952282654080722927118375271107572566252259719,-7.6806084769648273516285361071398752118988209313008859683388836839576457132898973,0.94867561056762661571185627826033808315242932467325841790867317211624938841090312,1.2944207545567560014379368281306796768911683993526690757981205807340999500235618 +8.5,1.25,-10946.305869219624039696521568063495573593652870609009793453260666457261747725063,6508.8137621077734802217753670160243325862411712176815132403177542185652596402547,-0.000067492186311940786183195859931359425954058732387897824208670687785976803277829434,-0.000040131718988153937871727767075219816980734597344519368598218973498603933137501047,2.0921058044345760428446831077683063307620943891866678983926521980792223032122571,0.15480384809663781399572731012379175661570702969101351588539273604418785816779724 +-0.3,0.7,-0.84835962739534080496307373449374610532359972025060614302160842421807228525719798,-0.53024136947899736437507009091719396857573180834664906888926309352725537091640482,-0.84762222377677940260859564189077705448592348540953884248681151632022919316515035,0.52978047778644330923963549922574640018109292819705399267194033162446936899241702,0.14296332294659075628593866446651901313813773183826771642301000206990410768689807,2.432038869506476422108099917118336097043512759886743844644777528304176079710315 +-0.3,-0.7,-0.84835962739534080496307373449374610532359972025060614302160842421807228525719798,0.53024136947899736437507009091719396857573180834664906888926309352725537091640482,-0.84762222377677940260859564189077705448592348540953884248681151632022919316515035,-0.52978047778644330923963549922574640018109292819705399267194033162446936899241702,0.14296332294659075628593866446651901313813773183826771642301000206990410768689807,-2.432038869506476422108099917118336097043512759886743844644777528304176079710315 +-4.0,1e-08,0.062754902851325007729976231144098276953343923106123602792586672703973347088204584,-4166666.666666665467627744640760575350710571990200225599146542920518643709379445,0.0000000000000036146824042363217056850820197417928893186913921509858351243371919885017072828878,0.00000024000000000000001462327156079356805785570139605928072251712264826968439281294107,1.5061176684318004751663078555066558347822374640299066456851500276535229216222525,100000000.00000002859319569658052673302102327271779326613916666956459375052838231 +-4.0,-1e-08,0.062754902851325007729976231144098276953343923106123602792586672703973347088204584,4166666.666666665467627744640760575350710571990200225599146542920518643709379445,0.0000000000000036146824042363217056850820197417928893186913921509858351243371919885017072828878,-0.00000024000000000000001462327156079356805785570139605928072251712264826968439281294107,1.5061176684318004751663078555066558347822374640299066456851500276535229216222525,-100000000.00000002859319569658052673302102327271779326613916666956459375052838231 diff --git a/test/data/hyp2f1_reference.csv b/test/data/hyp2f1_reference.csv new file mode 100644 index 0000000..fcf3ee3 --- /dev/null +++ b/test/data/hyp2f1_reference.csv @@ -0,0 +1,437 @@ +# mpmath=1.3.0, dps=100 +a_re,a_im,b_re,b_im,c_re,c_im,z_re,z_im,expected_re,expected_im,rtol,label +0.5,0.25,1.25,-0.5,2.0,0.75,0.1,0.2,1.0509220049050646553612806861072754262290340591982140778855947347796157425062318,0.068955456364755334566033316887706214207138306957600871268042232735954814914304379,2e-13,direct +-2.0,0.0,1.2,0.4,3.5,-0.2,2.0,0.5,0.14250768254131209440624423874138796435807241164393253361141930101659146986541673,-0.1929973079525080270244276242938194881417865229225925849720661141562562975742573,2e-13,polynomial +0.7,0.2,1.3,-0.1,2.4,0.3,-3.0,0.4,0.52896612874970450501316518116761876890855725724974446070296338946233705322327926,-0.0091005805790748396511769301397638051641618571587402970476535616628315955436338592,2e-12,pfaff +0.4,0.2,1.1,0.3,2.7,-0.2,4.0,2.0,0.57114639244136977039487188169017374402468916317568142623628825199492539655835878,0.27728473367458146275845130549395895909589484626276923954587376346175059655954809,3e-12,infinity +0.4,0.2,0.4,0.2,2.7,-0.2,4.0,2.0,0.80303215149342885035938991158352113195945607487103134290541562955041793635824429,0.16109751153992063974089492858924516780027513842266852645317952919403574261117684,8e-11,infinity-equal-m0 +0.4,0.2,1.1,0.3,2.5,0.5,0.98,0.03,1.3040555290794469928916781678631748209447157201219914174836260484842610330030864,0.29759669389184639748068722444002491479522542041789824834072889033744244116134576,3e-12,one +0.4,0.2,0.9,-0.1,1.3,0.1,0.98,0.03,2.0211813513970894291446097180081758783051242394475819855609775570030662626235877,0.7890884263440798291178495809280001221697615939347065504679586855092850849060725,8e-11,one-balanced-m0 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.97,0.02,1.1771526374790107314814426471553640700163718341332435790238436712737844183861888,0.074731605972286677140398904851747292382971409738413322344959337335637184288212665,8e-11,one-near-integer +0.4,0.2,0.9,-0.1,3.300000001,0.1,0.97,0.02,1.1771526373574450371258955582200874556769846395240118508188413229322924291747276,0.074731606017464640300282791932095231616332718600915812134227960930900760864158839,8e-11,one-near-integer +0.4,0.2,0.9,-0.1,3.299999999,0.1,0.97,0.02,1.1771526375241886677770349508472556602429688953229229428580039347549690004189877,0.074731606093852350886510166387651783090430085607621098210046101052895691286000249,8e-11,one-near-integer +0.4,0.2,0.9,-0.1,3.3,0.099999999,0.97,0.02,1.1771526374026230102942893927015874717322743037809074386405568142326989384481106,0.074731606139030330931868799459615331469994855952558034696153025105473120591622486,8e-11,one-near-integer +0.4,0.2,2.400000001,0.200000001,3.1,-0.3,5.0,1.0,0.28606215220276783100053448319888182832766779064759733746743384492964555953088835,0.2870645844673169053729605592847638386874989523789309701072137467391602188426313,8e-11,infinity-near-integer +0.4,0.2,2.400000001,0.2,3.1,-0.3,5.0,1.0,0.28606215210637735129952507526469135197376349222564699275048937075907729009543922,0.28706458456373466861797794392239884651403802442930316948475176124647852708542804,8e-11,infinity-near-integer +0.4,0.2,2.399999999,0.2,3.1,-0.3,5.0,1.0,0.28606215229921289406228096099159244843208703950076968871709892231272159038241531,0.28706458475651564426450399193398608334359107299815648055280250010390021176623087,8e-11,infinity-near-integer +0.4,0.2,2.4,0.199999999,3.1,-0.3,5.0,1.0,0.28606215210640464275893101346015868597051965203294592020764637679793712138308226,0.28706458475654291976050784539294715774441269918997624448039826013313642413588495,8e-11,infinity-near-integer +1.4,0.2,1.2,-0.3,1.1,0.4,0.6,0.2,4.009772502592252205190772759373880352384959753261987214522870279732168849316121,0.60710334319476934493896703466117122663804363951969560291458055588614076879416992,8e-11,euler +12.5,2.0,9.25,-1.5,17.0,0.5,0.4,0.2,-10.790148307965848817444415790924155025893247081893043819574451059866649013835927,22.040961318697734036535840855536559214046965853709552311927419823314041122145826,2e-10,moderate +0.7,0.2,1.2,-0.3,2.1,0.1,0.5,0.8660254037844386,0.94735539412071878930998051596180600424566070808031964704602016196997640623822576,0.470770670933339595282119509293296750858057853601463592736711037994176894020621,5e-12,taylor +0.7,0.2,1.2,-0.3,2.1,0.1,0.5,-0.8660254037844386,0.94443553063972950692923206808357633783682715528204684658277126211718251395748733,-0.44025720419080339207882218115932935043834659107391924040509856677058364735211908,5e-12,taylor +0.7,0.2,1.2,-0.3,2.1,0.1,0.495,0.8573651497465942,0.95072714180554300757277375763211028716197780373771279664346906118600323734003234,0.46718168860770713493465827977160160687679325537650307570914105082124462519015108,5e-12,taylor +0.7,0.2,1.2,-0.3,2.1,0.1,0.495,-0.8573651497465942,0.94737707397929214777052343002093275915356596276363133719448231737596553097121868,-0.43709086250305418536110254684257060614320160341617783498524304910190298009656854,5e-12,taylor +0.7,0.2,1.2,-0.3,2.1,0.1,0.505,0.874685657822283,0.94395666848703384435243425387537497923132448006124637254522471204451297580862889,0.47430335801718034884938507340831877043742151236619345424172358053382763203079079,5e-12,taylor +0.7,0.2,1.2,-0.3,2.1,0.1,0.505,-0.874685657822283,0.94147607205375161897270302152141397791840981174288249428833431800442037628953208,-0.44337546596439116519118947174952318961996923607619662477830528686492639184366551,5e-12,taylor +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,-0.81,-0.016940742038142853771161286512004869167832316293606694655428546922034005755587541,-0.041396732070281265120153563365866517646350590788745833152206586906247807246211583,2e-08,taylor-large-parameter-regression +1.0,0.0,1.0,0.0,4.0,0.0,3.0,4.0,0.4923438400096354501167611672435443424531929971283547504015486399318285080909635,0.60513406166123967627028771297240769292953259802082496424023011849172190033034179,3e-12,scipy-1561 +1.2,0.3,0.7,-0.1,2.8,0.2,2.0,0.0,0.96293041819066552538451760258188905240041471289178403910708613248459998937132174,1.3491213367180762528869290092877330046660667136377606106959059925204975445615717,5e-12,upper-cut +1.2,0.3,0.7,-0.1,2.8,0.2,2.0,-0.0,1.0730901960696174711428494488274671025704823832386993255055043066955419159342041,-1.4909004718780877674806972084977087991970442350541644735594287704538481714093445,5e-12,lower-cut +1.2,0.3,0.7,-0.1,2.8,0.2,1.3,0.0,1.7672425327283729565380221292123008736461547182459402217866376870482660862526959,1.047984138525224070751989794028849508480410844584535981352260350125211271133826,5e-11,pfaff-upper-cut +1.2,0.3,0.7,-0.1,2.8,0.2,1.3,-0.0,1.9568677494350155614530067065043291790434756123598276456142914973685672941178965,-1.0144734953540480341175598785308460892121537763699964461547744197782732774063804,5e-11,pfaff-lower-cut +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,-2.0,0.53847301877698610675077823444468940163891144544451935112334862345894814317793047,-0.079271176529689843918143744389564582287689445662453324751565765022536685863175633,2e-09,grid-generic-000 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,-0.81,0.56239329693339348732216241869980989192147304818291110196689026119416315223469536,-0.013327404337652756113746192747393864610600379064214755628842860605238988993726389,2e-09,grid-generic-001 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,-0.2,0.5658736582123395592322784570080864081948898157917500544983288082029078245795438,0.026113650593518999796075300106392474561022229962242313594263044393412475570135044,2e-09,grid-generic-002 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,0.2,0.56392698160547996564292936135787368730904439032831848086695400546278042983415353,0.052842246040415435724135640188918747739855818366761598499256857003141209733867186,2e-09,grid-generic-003 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,0.81,0.55416582768966255944934726142171254513848716362728570574129365797153591714065636,0.092969120199192882760419676198291464214829854211126261077311483740517490237663695,2e-09,grid-generic-004 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,2.0,0.51480949345493926147656933742708217278043853601496369414934539163364440807797772,0.16039078144611744741761512968591393279679367617807090537748224582097232073785365,2e-09,grid-generic-005 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,-2.0,0.671334158770113238947216428014031323356866168799747861827829141366877650171415,-0.21545775801892020774110194830424779287109165571713538482751377740596560048713525,2e-09,grid-generic-006 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,-0.81,0.7979439938188421856144758911745880599514326076945742044221676228911228053920514,-0.11672541654660591501496911249107731977804096434113198829022546086123068319552146,2e-09,grid-generic-007 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,-0.2,0.85360087991071392542701842863841232269252166164012779618148672088521930710307434,-0.013303933409337363669704258453014224411401765562623608092279356723444613079895806,2e-09,grid-generic-008 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,0.2,0.86666095508735699208961904293898984275680549314133287386368594335569972892115513,0.074738658330897807147288467519865387068117691000035679579881912221311586373562114,2e-09,grid-generic-009 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,0.81,0.83279506714487134127139281353504513655673517668021830591758097534859016325057746,0.21007616847490874396722271622013308614955127465489005437745963415877151858206673,2e-09,grid-generic-010 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,2.0,0.66439402596194413568139285957575907615149698699273285571825940918940804709191106,0.36158983441307875582025254870201138350926312074187445190917655490265262067235963,2e-09,grid-generic-011 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,-2.0,0.69624591744940966779012430809893541461706091873224301734099827643099261775043341,-0.29789006147934398232329964851407087468261458601378553839726186722379076298505758,2e-09,grid-generic-012 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,-0.81,0.89779299592006329631552787683280488050710620704152019840829162210145000952234364,-0.23908673955293880326972385554050502582946930509922198909932880278255002480827298,2e-09,grid-generic-013 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,-0.2,1.0534099046423838591270028451734583682958350591386614610976500046033461592384739,-0.11144035987243722505751147475484358167102332855789863938197152273448148912617286,2e-09,grid-generic-014 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,0.2,1.119743290238804614557335386599311836567768673566063392516833007174731575124244,0.063593294123067165867561180115297667133282186901196578974099521301931477556064534,2e-09,grid-generic-015 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,0.81,1.0282560607926373221579928446954964424440106667241162790017059055147544698651764,0.33905000059385286681654252515793495581996256988018235368654776981739356406819933,2e-09,grid-generic-016 +0.5,0.25,1.25,-0.5,2.0,0.75,0.25,2.0,0.69754891696556241499593265436778020847906368638167194684383583051288310389566819,0.49970216015222079866156021774087893236489606617886112565156443688109380490807957,2e-09,grid-generic-017 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,-2.0,0.69817790343522322037215112872923233475888569217914472544551036325293176585456128,-0.34096047438200276051390187432587613990592023409874159250788875216807244082294418,2e-09,grid-generic-018 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,-0.81,0.92679818819575994790951333913006759016351107551864753723921421651824384513688775,-0.32867032568066534717446770008225990676081609803220581289062752927181099668126597,2e-09,grid-generic-019 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,-0.2,1.1768938109771924108753062481147900799839272094063337292606572239857660491772216,-0.23930769957680687164245242157890017874336114276737283122656557615964517358278243,2e-09,grid-generic-020 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,0.2,1.3545249084778554265743311988034074543950092163945526244483741840134251550984252,0.03650272331622641660014805001565611875220323358604929721444007639314020967836744,2e-09,grid-generic-021 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,0.81,1.1504498383451207526787259176708113947777822509241218904283096143932627317742129,0.46980967934271055471837189341698805277125654357137423588683675238149282359407441,2e-09,grid-generic-022 +0.5,0.25,1.25,-0.5,2.0,0.75,0.59,2.0,0.69550202003845580449370572816954289104349665581478412884871905173245512293776617,0.58071338915761139606704926209147318884846527537741077907314478564985747200607308,2e-09,grid-generic-023 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,-2.0,0.69332602103693834293987184360579424306687628749200166379319286828170743773451566,-0.38124620692142827173813507584756856890624396964716337100813118004114341337834252,2e-09,grid-generic-024 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,-0.81,0.92653791618516824507789699443545595443283821132705092300470725395505882921242995,-0.42355517291478455243461718776944594952628041065958841277881813437455554054902997,2e-09,grid-generic-025 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,-0.2,1.231545903066439308781688700517347489126151816909205037853523296106110532824761,-0.47633370631997037602140544794173093606523214459332359652697994001098956189220025,2e-09,grid-generic-026 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,0.2,1.865919527001623908622487679272459771912857260921969979127121372489548281826045,0.0054082506164883297896490238900619485074757921809877985738004148178888798488734544,2e-09,grid-generic-027 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,0.81,1.2518497759990914967449355063327162373730654004005022461080409408020265605407337,0.67863795325445218961362877857242848501540946151468992852229152924963691864242493,2e-09,grid-generic-028 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,2.0,0.67708068382700146097947708751995301811776081617307313041466118350030871764564376,0.66170536330823802488846782119609511104816491291813450620958007463827441978731288,2e-09,grid-generic-029 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,-2.0,0.67756352332395385707864584061291111494846359234794089796662957730350346827982803,-0.43161594181511335970531389080580382807147913686510665466436661094994579934080708,2e-09,grid-generic-030 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,-0.81,0.88230996711430920664221719196363337008551412008912695509179947401464013441696629,-0.53685563184841030879611449694333156602268471383377948789716586726196098947207956,2e-09,grid-generic-031 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,-0.2,1.0606026901581528651528360987532022129886588120616462718148469631895477559696171,-0.70373887693607514918373555410105010077257001882653818592843237426003095781394507,2e-09,grid-generic-032 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,0.2,2.6259827399322874490504002975813068954104128390086573553901357667708817663426765,1.5870592694800135083953388866626188185522306230995397609180265635871766425798575,2e-09,grid-generic-033 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,0.81,1.2308940661642823677311176191756951489903996058689672645198724325603978727627801,1.0771535288057801579551800210717831625625670744203469368389676812478861069516434,2e-09,grid-generic-034 +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,2.0,0.62388289766869456361710900586619138286443669190728728403001827638552255051624224,0.76692048254513188311829731077980460620158335066592940986358104866061297264297159,2e-09,grid-generic-035 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,-2.0,0.54084331871690494014017224292957347813433505512912359693293171338231315840412392,-0.5725566462195350931891249833688258214511898116607156081599001581140312589557206,2e-09,grid-generic-036 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,-0.81,0.58276744409431626762645547201110458045077109334520355786934450399892465200733043,-0.71064174718113167213757020925029838368092438001365374948468529417411337701539369,2e-09,grid-generic-037 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,-0.2,0.58203765136626831363513369227628771247516369214863346215732600250413872566748681,-0.80695270658876191329461211319612219868050954878154153053487265305268356927496762,2e-09,grid-generic-038 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,0.2,-0.31025756490463737414186363803062122720907680673770801633583456529364884091616574,1.6082111626451416011807502482085296485919280816408892866823477625575443087711224,2e-09,grid-generic-039 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,0.81,-0.015688317380076755810621784373968958091878140226709179837311678724776637507310925,1.3558007347194900177301559682527956814385403441458048454909667605388347316323084,2e-09,grid-generic-040 +0.5,0.25,1.25,-0.5,2.0,0.75,3.0,2.0,0.16842989779571448156523505599301421923795852764945985401403020352594602973545927,0.94290733341778546447085755456792187984853393775150609591550634370725915496447276,2e-09,grid-generic-041 +0.5,0.25,1.25,-0.5,2.0,0.75,-3.0,0.0,0.5653400419622808140161391483540940444779606811108331331817131157651038525278143,0.039452816447344625703529186489944346514614464110405966638345739371583784210769033,2e-09,grid-generic-042 +0.5,0.25,1.25,-0.5,2.0,0.75,-0.5,0.0,0.86337871675501459103260002893739539382125482534053972341576684549809763000916657,0.029384243395073748307793692165686754679908674142502175684503374949579578591224259,2e-09,grid-generic-043 +0.5,0.25,1.25,-0.5,2.0,0.75,0.0,0.0,1.0,0.0,2e-09,grid-generic-044 +0.5,0.25,1.25,-0.5,2.0,0.75,0.5,0.0,1.2276603822573563063143901572166638799223633666321448650728971713220946172684006,-0.096604849349378184060908021456879274621401744772532668138635837298358166782691469,2e-09,grid-generic-045 +0.5,0.25,1.25,-0.5,2.0,0.75,0.9,0.0,1.5405590534882572498219909263081866465386607137010661167708638176406067775542393,-0.49995776420793301676094554387884207708853919772749798102249450686323753001485939,2e-09,grid-generic-046 +0.5,0.25,1.25,-0.5,2.0,0.75,0.99,0.0,1.3409921638447517565204018230602489935083259231250490992990573469953443635053554,-0.8315710234346652916056869297843697195227018778723072419774580615668301492253996,2e-09,grid-generic-047 +0.5,0.25,1.25,-0.5,2.0,0.75,1.01,0.0,-0.31831466956016010691704567633070278574866599054884131245642605412510932820686745,-2.6407166303472442862125375786711849573643799238520234426039566171204721109098255,2e-09,grid-generic-048-upper-cut +0.5,0.25,1.25,-0.5,2.0,0.75,1.01,-0.0,1.3180089884147052620774849600187191765542563266440813144520892307858567388352474,-0.73290185534565503476189363370143476269491741353065830386396622861272977200398396,2e-09,grid-generic-049-lower-cut +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,0.0,3.8976998256716173009701189764245577734367683600470433303401953312142418468342624,2.629081735917990519405600070187473975506647271079872726767549282636407985504821,2e-09,grid-generic-050-upper-cut +0.5,0.25,1.25,-0.5,2.0,0.75,1.3,-0.0,1.1124974071276094386038789153942975804629956042000789761669659487870552456957902,-0.80886930992211921784317224446755126538441177268697273764604979123448346149001608,2e-09,grid-generic-051-lower-cut +0.5,0.25,1.25,-0.5,2.0,0.75,1.9,0.0,0.63636185888828587689175678784747811896748896718142470044958434859176876991479681,2.905746219354396353736549011145595278805763497776759626823246114407441431250197,2e-09,grid-generic-052-upper-cut +0.5,0.25,1.25,-0.5,2.0,0.75,1.9,-0.0,0.843458255631391051105693719512885503032562041001212805213295762626067791542761,-0.86353210343576842882847588474701728418571452578871011284411648018171745683262651,2e-09,grid-generic-053-lower-cut +0.5,0.25,1.25,-0.5,2.0,0.75,2.0,0.0,0.41042196470412212031622177331879555594531877119274973234238739345118438648830776,2.7666173585899250404522351033903992302513175755193220090694055645883554506084607,2e-09,grid-generic-054-upper-cut +0.5,0.25,1.25,-0.5,2.0,0.75,2.0,-0.0,0.81036465421419778364697313112623798518168350869769463206360379428358620399607679,-0.86495670136977319439272103237217692185190522007565375718831115801874594591579057,2e-09,grid-generic-055-lower-cut +0.5,0.25,1.25,-0.5,2.0,0.75,4.0,0.0,-0.5597441235345862346916107585332278964535280461362461185395934301858322148432189,1.1062268486703621347105862828990871356389560400794109346503801769989043040383781,2e-09,grid-generic-056-upper-cut +0.5,0.25,1.25,-0.5,2.0,0.75,4.0,-0.0,0.43795853779675071545971921572929773436449987934494899086518899931065878634681264,-0.79781253471977773569911727456429768011927607758610103908587822576569297000978133,2e-09,grid-generic-057-lower-cut +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,-2.0,0.32013003887954502033324911738772414387009957542667352765340129347381308840745434,-0.19031045415352647376494853651832014291121225041238333157531686748470169334574768,5e-09,grid-integer-c1-000 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,-0.81,0.35256842311482749974774719763729414271279772994676508753193711135660823684952238,-0.12586744416968247331952878685075239788922853506624766929128656145550459942306993,5e-09,grid-integer-c1-001 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,-0.2,0.35683926473962499843648500206839156922178026484064791339413335409290465610135263,-0.085996897406055274603158805872978454956059057919962947284927990994328808827801841,5e-09,grid-integer-c1-002 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,0.2,0.35389821880755685197643098712613395003840997573496074724559492161380782747132592,-0.059644485814732530846810251345749808526597186553114993001921159807752717174074163,5e-09,grid-integer-c1-003 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,0.81,0.34111573160957666367741784098561973714166348731784149863631495328054339043571727,-0.022308256961333870891730760992582507722422192734762051556538096595379220256538703,5e-09,grid-integer-c1-004 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,2.0,0.29708892301845484567129308750508817998291405958153034314411266170733776441370404,0.030578279116469795013133071145101341280116815587849466999397726973119505113994782,5e-09,grid-integer-c1-005 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,-2.0,0.43422827032331109878199509181056499260629960472692406162286136126673465213277118,-0.39187890374033528138805460541676811610840205799092387522009295069441820463503355,5e-09,grid-integer-c1-006 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,-0.81,0.65576019820567603531422448673994827694495561150423403428693336474238922611953188,-0.27913933507782087974254089786328785679045646569191492307148598913982679767936383,5e-09,grid-integer-c1-007 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,-0.2,0.73803766854623436510895592647138775456872615080105678675467943190576593375343348,-0.099150292017699808479807601226485316011803260405411347582422791315162109816862324,5e-09,grid-integer-c1-008 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,0.2,0.72139412851380294753852783722190867148318379718371448055242987391482294314369547,0.045832670235127999983393226133409004663559922143412976658447807274545915263271789,5e-09,grid-integer-c1-009 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,0.81,0.59433607840175291288753562168760221495867110791358439006199854200421723858720636,0.20482025724379065703819582390560106026686505177689189153813977095171919394689055,5e-09,grid-integer-c1-010 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,2.0,0.34272466888176269821316483893659682050471574632897255884117384811694464738423145,0.24664645032320150352889403683981225293778303902538009053245960743808678310333143,5e-09,grid-integer-c1-011 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,-2.0,0.42013160139731350085799793559910636313802044193085987545806553064573175584185262,-0.51984306739730245492690031237464689082778592970942671708579881694603361230420427,5e-09,grid-integer-c1-012 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,-0.81,0.81803734179638664054771478385548895701318104376891514486682700724465719772354634,-0.54306858461651864342366822578905670452372470897679638146395884128021363992281519,5e-09,grid-integer-c1-013 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,-0.2,1.2002138035765922053379011537997939401721375498393247152898363911797615927075492,-0.21839292293230183932602410083111239012085647105672764012139357905316747784320944,5e-09,grid-integer-c1-014 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,0.2,1.1722040679741543473554288418454301289171371171920518229352754031987111700533289,0.25213126511905844320295652873473936183947435307498477082081663557537623326688841,5e-09,grid-integer-c1-015 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,0.81,0.69159414530440123031367116684310703411110610153787841818474210151018722020255055,0.53032875175312768846530433906810394668239364727973547478398300597383711037972592,5e-09,grid-integer-c1-016 +0.5,0.25,1.25,-0.5,1.0,0.0,0.25,2.0,0.27933837624440644617617083446940280760691327280340756410215026600956622981330947,0.36237302941389986851747052369523355309476302701952819963397246363772347454771983,5e-09,grid-integer-c1-017 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,-2.0,0.38904858353284680675520898566069540911109552171273860770128209403151993754374145,-0.58070530430678086256478723017045441915426372365194832571392960215643305466259666,5e-09,grid-integer-c1-018 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,-0.81,0.82127455076607961799369921453642128048138768220821574011208525139935720959008534,-0.77580193022592690093564904248859256984227451568041058906110388541031578966577065,5e-09,grid-integer-c1-019 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,-0.2,1.6968003778010666941399421712162214320424965007016259206963088571230708746133059,-0.55462047042258763385614760754714874004383781316973622077080975122294664896559811,5e-09,grid-integer-c1-020 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,0.2,1.7034869931262064746358124742796811324399165635469077883478436274107087156941117,0.65204854610647198898122777827322305184024944233570108435116952054856029417737349,5e-09,grid-integer-c1-021 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,0.81,0.62478754655366461323772520526585253082765041200441638118176070993962437171995413,0.80455615249306243960940784191531093888494477132520342408222690768369273976291113,5e-09,grid-integer-c1-022 +0.5,0.25,1.25,-0.5,1.0,0.0,0.59,2.0,0.22127858544233468346898007360267148694850860016135456445814927445921693998101843,0.40679354153401332210199924402783249148209383857075192031721236526326073869417189,5e-09,grid-integer-c1-023 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,-2.0,0.34599062699363778769777615108614168790107195657339988715183538564058852893393946,-0.62970573909084345883598375091382464917201961823811604677275430482301683265451177,5e-09,grid-integer-c1-024 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,-0.81,0.69094601485479158474048091857098431908312644685763292616856796439226954936991708,-1.0044334264983985712330774499634220627617165458797922377232141329269089193410281,5e-09,grid-integer-c1-025 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,-0.2,1.9439167953785046879956425130043564970721919052102385433466182611947060498411257,-2.0047519129428408341146002607852909961591243769608365432682288721007107433035129,5e-09,grid-integer-c1-026 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,0.2,2.2655720536224902895361739489498507007333344897030930504789895871550085844613918,2.5662665634251816160020722270572203788328203737841901906722354784556790890952966,5e-09,grid-integer-c1-027 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,0.81,0.37109886278040222912024109804589519623158118017164724988999829680818057031636407,1.0522803131095414230801117112259681389537028083645658909044559933315731371858988,5e-09,grid-integer-c1-028 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,2.0,0.15421643655300752471297971440787083015791041759622092636751772103652616603100432,0.43185933279972784633292265904479641213351675013140794733256886552816322481650332,5e-09,grid-integer-c1-029 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,-2.0,0.27518684015291965120381041519383856217404975271928010101208238736728945016318291,-0.67524663910011290468778455201083554794023548782833263412143750959338347056760296,5e-09,grid-integer-c1-030 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,-0.81,0.37931915989786420083276872382680062059567108631913432344025162321509640337461136,-1.1294028032468538246137461180080985584303361727707397612525743036311705654171054,5e-09,grid-integer-c1-031 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,-0.2,-0.06715621181688166685737316485200418643900450720316793056058949347718483685755863,-1.8591042388459326790479311207965664588070007999646672717746194633240633423413314,5e-09,grid-integer-c1-032 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,0.2,-1.6663001401520661905582912702586322109943931972708609551499533254842871717576974,2.2106976426386505089129458969805120327293627047371695275258782127353561683872636,5e-09,grid-integer-c1-033 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,0.81,-0.14072633689277111750138549638957277181269928013101080431093962084727282544048122,1.054075952977872587741267994750174540198629601712798911207727369149184010302868,5e-09,grid-integer-c1-034 +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,2.0,0.059721099600544685452345866816612244827840973183285179879041998591788831406744179,0.43338546223530020239134627194737711570567073416203288387668772442009729798209911,5e-09,grid-integer-c1-035 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,-2.0,0.0038398829103907636385641146848522337672545053097116210200348423797898905258937401,-0.65410148837295748461120709114127625722076913176864355465177325118006906304498544,5e-09,grid-integer-c1-036 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,-0.81,-0.12578138752557664547468470345535472660263712469557787600029901930613873781422211,-0.77371360810745487086251580329051759568782197932513443477201875812414277489308327,5e-09,grid-integer-c1-037 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,-0.2,-0.23556561193566234054189168335115801664031143634677707045873656372164776683787861,-0.80136385919240992473528052323937352159540368913639412689802648949258336694832515,5e-09,grid-integer-c1-038 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,0.2,-0.42397117694740794412801682620417888383575197352908652949968101426192147760597126,-0.00082562456847999968053926737954255377100754827040231058146999951285485064864188679,5e-09,grid-integer-c1-039 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,0.81,-0.33612428170761611846424014797985053990876111111596014140386662442459704740023533,0.14269805617880479295584832428326463440092278630338635847480791259931722039657407,5e-09,grid-integer-c1-040 +0.5,0.25,1.25,-0.5,1.0,0.0,3.0,2.0,-0.13706488665459489087278127077782096673872892071022400305893120238515286406252211,0.21394015022407325229876857169600680878861848809125803655724643196996516939662879,5e-09,grid-integer-c1-041 +0.5,0.25,1.25,-0.5,1.0,0.0,-3.0,0.0,0.35594190891338755580066746753469347307490110163827134949217851163538201794030277,-0.072735879859274792585568103428003987918276150236357365114346375066461209614321885,5e-09,grid-integer-c1-042 +0.5,0.25,1.25,-0.5,1.0,0.0,-0.5,0.0,0.73820296941441454716475138746764105967279767700436516979697186596076763611367968,-0.025947033014552058316247911307246257604003595957681273688019042234613323839249914,5e-09,grid-integer-c1-043 +0.5,0.25,1.25,-0.5,1.0,0.0,0.0,0.0,1.0,0.0,5e-09,grid-integer-c1-044 +0.5,0.25,1.25,-0.5,1.0,0.0,0.5,0.0,1.6831000677775836920964725475801253034727798280651130975071094786255330417230549,0.027441187461574909381715268022554819185802237018585657057812837991946799455282002,5e-09,grid-integer-c1-045 +0.5,0.25,1.25,-0.5,1.0,0.0,0.9,0.0,5.5429140821036274162999657256865989007445168141548000770295642112830575098380833,-0.65749886581217485536875769396033472345205727584121035348238680479776585158652516,5e-09,grid-integer-c1-046 +0.5,0.25,1.25,-0.5,1.0,0.0,0.99,0.0,25.816046086507279510671989783379990587250765791562613919942922024901604823511162,-16.1202079476002120143706904834809924080027489824309791761243115435563994732313,5e-09,grid-integer-c1-047 +0.5,0.25,1.25,-0.5,1.0,0.0,1.01,0.0,-16.107768012703603611606135590981167881275334675386741196373011702689919710411687,61.44205760610252624975194508516695880217485666447428836710314624690125593546811,5e-09,grid-integer-c1-048-upper-cut +0.5,0.25,1.25,-0.5,1.0,0.0,1.01,-0.0,-13.090562437444035290795321965976968107375229225916173642448771275771145095644025,-4.5497086176066415842874701077523650073623514647531404332168292913701959700918804,5e-09,grid-integer-c1-049-lower-cut +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,0.0,-3.5063850787105349956340293919456613295519873375072586461212499529673969730641155,1.513553517765340917873670542133813572933596004893950763853219686101768472639917,5e-09,grid-integer-c1-050-upper-cut +0.5,0.25,1.25,-0.5,1.0,0.0,1.3,-0.0,-0.62755399885319999890458587248470683864585413454661482032936866209744872790701292,-1.8010759687459307899351593879448094325589370925755855076970382351654207366641291,5e-09,grid-integer-c1-051-lower-cut +0.5,0.25,1.25,-0.5,1.0,0.0,1.9,0.0,-1.2163187610881811949571807927854154164025819299437701002309875261268891036035769,0.10749276996431602261735837727966679864042411456985605648983703930385143852142123,5e-09,grid-integer-c1-052-upper-cut +0.5,0.25,1.25,-0.5,1.0,0.0,1.9,-0.0,-0.32521467679511506916943427145152240983253576602663607445843909091704589223364103,-1.1483613001769208512153380377468304210719424090863820078805229030791456446640831,5e-09,grid-integer-c1-053-lower-cut +0.5,0.25,1.25,-0.5,1.0,0.0,2.0,0.0,-1.0780734471000054268419045217754071887436486772100359097148619395517811390875312,0.062000109992741197662723663027235230614533540079278746456992858662124204912132499,5e-09,grid-integer-c1-054-upper-cut +0.5,0.25,1.25,-0.5,1.0,0.0,2.0,-0.0,-0.31418681738971561836281606306276346763887149967657632605458533500078017033421807,-1.0973459022695113290409225846277474769172261915134691805464392864501897121620047,5e-09,grid-integer-c1-055-lower-cut +0.5,0.25,1.25,-0.5,1.0,0.0,4.0,0.0,-0.22327064521432855667353056113340781625626594327656776280700175693930927203031693,-0.061920232152644847649079292277022637317918096520080475037928543141459652084201263,5e-09,grid-integer-c1-056-upper-cut +0.5,0.25,1.25,-0.5,1.0,0.0,4.0,-0.0,-0.26270201537787414147176419785452939776945039615093753997984454039440207470884068,-0.65512407927478315372198052401827695341746261850657908597453825865098698876413517,5e-09,grid-integer-c1-057-lower-cut +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,-2.0,0.39261058259614961398931150989042097192578633928404831664507392663864906849765012,-0.018574104268091517634590867798289968120517829769636653359747335735514157420372514,5e-09,grid-integer-c2-000 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,-0.81,0.43240762305413525056964441135746368008158791963629246309916862320145895968678777,0.044535550197997991163146687746056527780210700808699421398117682159176648040371243,5e-09,grid-integer-c2-001 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,-0.2,0.44316399615855458206538386159244468142158716803011109026725871279926184897435907,0.085900227767330503880823039238948261311413713436758485139178654613345493191197062,5e-09,grid-integer-c2-002 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,0.2,0.44509603385653570541345246328117936787143909256394235721619107985100345826219116,0.11479909550774695174673115902025568078502362277641670979172923942436425407866173,5e-09,grid-integer-c2-003 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,0.81,0.43961465069241602285647852528068012340040311872629451811944135637729762882345331,0.15892288182693020020877283990785050683464007933955360938804224124217755855127898,5e-09,grid-integer-c2-004 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,2.0,0.40374958252198517575606504989035581166087985976850580840217887206109125336854098,0.23368258565072220631752181720476347005971973840953837477826841091917131587276282,5e-09,grid-integer-c2-005 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,-2.0,0.48654016328053133178207434714899175782610255963950605047418554373993435578655202,-0.2141046194189482329488258159986005090732946480745040317282545621309547162704613,5e-09,grid-integer-c2-006 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,-0.81,0.68660664059029216613895236737843926286310409475763272597141903004450267684527956,-0.1398715555633667774724073844649557985122257639805533310394729008426370783474515,5e-09,grid-integer-c2-007 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,-0.2,0.78584226306667455156566438913680357136511666733496605441287837875890189102879878,-0.0099583527563674387897891072082385768112433961148733127107105937949886522053581104,5e-09,grid-integer-c2-008 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,0.2,0.81089030771207566796067322283222907152424571300504198336682389229780742150527523,0.10836230048343037448638306617276631048222382218826852042482064767339305168825368,5e-09,grid-integer-c2-009 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,0.81,0.76694322954282181805926103009659335443262830950325951934017725935297941590199628,0.28078546153514888846692125626821304072739852073976747439099736618284770329506667,5e-09,grid-integer-c2-010 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,2.0,0.56764200745421960503852814562536516704469162432818879972198377595779616292610631,0.44791109946490696205712667645068790558384612346531086020975588764881636518839948,5e-09,grid-integer-c2-011 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,-2.0,0.46948778950297702467005188850916461448044295272014432818563429816013082294826435,-0.32445207494020555154516712912406079633726933728449514380449414220475657814928734,5e-09,grid-integer-c2-012 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,-0.81,0.78656549669547930800291555920394333710052097062635577663896251825634581062726335,-0.36065483288860745803074766593306695997014190998003577009105701940927311454439309,5e-09,grid-integer-c2-013 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,-0.2,1.0867949356752492642397824396638151157370059066080279237214574480003464647789112,-0.18240510925951226183879749344572643460369222021211485464727810158295587152153755,5e-09,grid-integer-c2-014 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,0.2,1.1759761392359283716493873607692499651554830938554616181607221431488870389758491,0.11243844169496131488606995028017038575910285302369302528738254311862164902813423,5e-09,grid-integer-c2-015 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,0.81,0.98958312740639585347536652488446916526253389932582465503448772282084794458018592,0.46516672037631809904966479490207994504315282155554784059009409876471435067228565,5e-09,grid-integer-c2-016 +1.4,0.2,0.7,-0.3,2.0,0.0,0.25,2.0,0.59386364697835768458655285178521232633057039893071688407367611774102957282875304,0.59739099573006771604903737652192952823697057403602596011696872035829020554182669,5e-09,grid-integer-c2-017 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,-2.0,0.4421188538597769841857689843621221797715852468315869785213329843376596755275,-0.37615352876877241395294982478984095273814240928965035687363356657682960330041413,5e-09,grid-integer-c2-018 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,-0.81,0.78242211428339090009856264711296594592390166527823064497442632451890893119756819,-0.53030467085389771074855836148638696298013857904028765035903288555769915676941033,5e-09,grid-integer-c2-019 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,-0.2,1.3340871687388738884998729951106774360257107344541347505075828079633828214849256,-0.43445818017441965596321265156404300259835179304941287472293561176161875902897862,5e-09,grid-integer-c2-020 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,0.2,1.5568199960116478997317087456259373319947827263738478034993315857723120175104671,0.16560254113612193573105171915772023512985993088554510095049419181010976819705753,5e-09,grid-integer-c2-021 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,0.81,1.0967015330167956706676590161767173999467105088775667757130960505165801118314783,0.6530810412209859502446464616188364084482354041813643639464690259932139508980036,5e-09,grid-integer-c2-022 +1.4,0.2,0.7,-0.3,2.0,0.0,0.59,2.0,0.58385594957194956774214769075947955078042554195718870367403601216556351628740966,0.68076830022638838696840332254981266878563514260739911247627201415052404640861591,5e-09,grid-integer-c2-023 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,-2.0,0.40534305082315784237960881699916920466518151590894366778783616356837294959682531,-0.41835506843174698983624112888691050542423240545766670808883611700026856173773903,5e-09,grid-integer-c2-024 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,-0.81,0.69762726175620006130276783537063977228678287214431795732035837162822150542573802,-0.7011693393174566831702090140484038615741527315574184209102581716465205247880296,5e-09,grid-integer-c2-025 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,-0.2,1.4849251777146799035233280385169698098090598099381932869607008326807981709142657,-1.1293295402829668384817429494193557815470902304426555428296630535243283078128003,5e-09,grid-integer-c2-026 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,0.2,2.2771499428247862581246494325893654536003210483172129254845939884749202874392805,0.61157440861863292797792161262786480558834348383298793505074863972377090458750842,5e-09,grid-integer-c2-027 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,0.81,1.1305345719930790313143026324145744775640473481838623325239657103179661709773532,0.91148656664858984994378576787793674618920091068542795270102322301665113378381632,5e-09,grid-integer-c2-028 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,2.0,0.55670095006262022062021895442807754093660701728939557130983654487214126205689087,0.75983775308829246597017313454521445725879665850339189849102908293178267787397021,5e-09,grid-integer-c2-029 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,-2.0,0.34470027492392501969602996970956039820163987060775996026027405803930211131740628,-0.45912789191207065416501325299536246868972580129375472781591483099735032168991481,5e-09,grid-integer-c2-030 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,-0.81,0.4746909412580954595890985616506313660208413715793289185115746970657157085017667,-0.83393517614864255256989694042595030392372795747682433334471021919632206372018389,5e-09,grid-integer-c2-031 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,-0.2,0.39175112837842770456147452744602660282333127981699895644551177571820450142276889,-1.5355822613246265289107977212658362405605707708544666059836210608697640675438334,5e-09,grid-integer-c2-032 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,0.2,1.6529640984699015220558372326717915265448321948787313473030052002342091462185011,2.2132176008187683333123615701055010040618350287104666556545336338959112994111137,5e-09,grid-integer-c2-033 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,0.81,0.97267088902159194852931479428427741871949040844731003468366077613214685062031552,1.267478484320554108640394647800471181465014346851919087263061976159585497834246,5e-09,grid-integer-c2-034 +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,2.0,0.49412958334128815098544436079566605345731715396586151619300573232637222424557754,0.85551292649656886312889072817750065391606462019788626581668358550133344017181385,5e-09,grid-integer-c2-035 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,-2.0,0.098622894672072623749087688029796370187805732061866451593441620077544221529358096,-0.44249972737250324991732790574629193267018923431472404524246226006523069199607211,5e-09,grid-integer-c2-036 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,-0.81,-0.032785972182369158224736452466275929858074651083448612245663778082146203420372539,-0.5572999914595092060518744040868754236652016686040235183794102771304276075158466,5e-09,grid-integer-c2-037 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,-0.2,-0.15515129360084174392959818480645629105707967620709341905757472690239503287312832,-0.58312729757212669658242823007859786222841909051556776711874005550572510834381174,5e-09,grid-integer-c2-038 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,0.2,-0.25911057219143473909373008940136133147856368253386211294201878832243999199187745,1.4995172137645131327547465657709901844454892140578974364479406346483298722677243,5e-09,grid-integer-c2-039 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,0.81,-0.060952437078351574098949037644388430330443179225193351663356197851181548301082999,1.3184449135380110842868048614028296878216871549861835815137830993165545150526388,5e-09,grid-integer-c2-040 +1.4,0.2,0.7,-0.3,2.0,0.0,3.0,2.0,0.079261687350263436455562274700888696310423438134125859555539790475715907588103418,0.99262386565636989835768359030072452035264043794280233248734438261952388474559102,5e-09,grid-integer-c2-041 +1.4,0.2,0.7,-0.3,2.0,0.0,-3.0,0.0,0.44467182956348244562358655869436281397091742205683861298591547626153911318154951,0.10025602564569671245874922884811369756984722311859742107998220367221895209130092,5e-09,grid-integer-c2-042 +1.4,0.2,0.7,-0.3,2.0,0.0,-0.5,0.0,0.80384100570224827150625587315785997520316869446449701691468509630999218761525413,0.047486913582232518049066970469197083747184314455440845915167311734576758271646306,5e-09,grid-integer-c2-043 +1.4,0.2,0.7,-0.3,2.0,0.0,0.0,0.0,1.0,0.0,5e-09,grid-integer-c2-044 +1.4,0.2,0.7,-0.3,2.0,0.0,0.5,0.0,1.4036709917465354897341997083063369690951743896508855181119483389414568875738903,-0.12807164696593890424781560043795659827644678849100095175197383667714212330341511,5e-09,grid-integer-c2-045 +1.4,0.2,0.7,-0.3,2.0,0.0,0.9,0.0,2.6514864844005103197797315947912643302336458278393982888004597381579913173759899,-0.71734957312137717022932140238163395409628861009683102291928987942038401063551054,5e-09,grid-integer-c2-046 +1.4,0.2,0.7,-0.3,2.0,0.0,0.99,0.0,4.9473888586094187941853204742181500493616417536191386544622194547316166144941823,-2.4346698275247962507407961903513910928122146876997239320811788535414913499153461,5e-09,grid-integer-c2-047 +1.4,0.2,0.7,-0.3,2.0,0.0,1.01,0.0,7.8559732068802529336615077986826889154350376337077687635372663842852725215822005,2.1896865826199168819226283853983765942154272037011564421954237326406224223826694,5e-09,grid-integer-c2-048-upper-cut +1.4,0.2,0.7,-0.3,2.0,0.0,1.01,-0.0,1.7692814020250577316074745977689637971744115634984677755040943658868794814005012,-4.8584018949764971341015300662204068667898890875041312843035915813569858550520101,5e-09,grid-integer-c2-049-lower-cut +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,0.0,1.6796330679750063627149639210270034470672816329395932881424858055427284641841338,3.127880089776983169086766070606915874414188345294802799065746826696706279939421,5e-09,grid-integer-c2-050-upper-cut +1.4,0.2,0.7,-0.3,2.0,0.0,1.3,-0.0,0.01222782470638034825886030276201502109626942280600193202789497466418941199923132,-1.8787780252916480690896674969736947286722238157867731597211211691712201185274289,5e-09,grid-integer-c2-051-lower-cut +1.4,0.2,0.7,-0.3,2.0,0.0,1.9,0.0,0.21265336004777517466430118566784370723377077611324518915383748281028814024825055,2.3258118855416943434071096803756409994152206998252610192358548737023500445612326,5e-09,grid-integer-c2-052-upper-cut +1.4,0.2,0.7,-0.3,2.0,0.0,1.9,-0.0,-0.21365301396437501619890417254562164091486843708358278488760421112466637632358492,-1.0351789921534204001163232025372764192434998138333236651317745527737404750840835,5e-09,grid-integer-c2-053-lower-cut +1.4,0.2,0.7,-0.3,2.0,0.0,2.0,0.0,0.11182087032188755354176708508471280429459578546426297123301777854100005056128615,2.2254534919935751363239833128043381080411409019501816326924292288290209661802234,5e-09,grid-integer-c2-054-upper-cut +1.4,0.2,0.7,-0.3,2.0,0.0,2.0,-0.0,-0.21888878910377086217982194860897647148671622246799264610480059866510924628867681,-0.9652707620488584935791568968696285137757775358636150902561899103781751356517025,5e-09,grid-integer-c2-055-lower-cut +1.4,0.2,0.7,-0.3,2.0,0.0,4.0,0.0,-0.46966551510418656659257628897438359026895281109498488619935220958637053891137723,1.1797685751316721689152805631717245425770037920342443518547998407850336536191933,5e-09,grid-integer-c2-056-upper-cut +1.4,0.2,0.7,-0.3,2.0,0.0,4.0,-0.0,-0.16460458387981631919245097721489463143175166954886571080324039419650119991048538,-0.4162095053172912658390254258000742686717583739375233207615566882936340176980051,5e-09,grid-integer-c2-057-lower-cut +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,-2.0,0.11029141516576803657029685185740732257039530403580333536173486582393238680172366,0.050740060643967738427099149302860205410920180874343031284081909119581138965295423,5e-09,grid-euler-000 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,-0.81,0.12393093283476998337567837497768197316795631326640182489186461552167412208551493,0.096485678230848201558273548000384770909779269972065908368689861931779838199062182,5e-09,grid-euler-001 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,-0.2,0.11971891347332078006795101437146733874502262024181480382573139082641293547518369,0.12782435429014867170255800763027724521850389930545575259213572224135219468529984,5e-09,grid-euler-002 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,0.2,0.11063385967575243613553900875863047876976342054652013306258303280736358139832334,0.14885129368895716490085283191167643015718874234981680851665504379641941479704535,5e-09,grid-euler-003 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,0.81,0.086983831758396572193969066646381009169772896404496953782845034861488751264781558,0.17715671031199071772857187113182261455859069921298102180226189579699818758338861,5e-09,grid-euler-004 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,2.0,0.020424512720872094431510141007168681074302165147418648951988847589515665127341375,0.20324088185597653323591237238542296908228762128281237160008071810804028520223571,5e-09,grid-euler-005 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,-2.0,0.14391112858802754399477302951288409946231856252214238827113017783711302218221903,-0.072678367922641928172861413539928982945031001043466810019187973626960642082949208,5e-09,grid-euler-006 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,-0.81,0.34120729096074860773992189469636277793264466393846192085347354232287345840435009,-0.094212525493713517546235556738034369549045149258287505980914671751435604866421455,5e-09,grid-euler-007 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,-0.2,0.52912479855177396326834724434335831325515305709430275537761923073604136864607297,0.053979221693425279391740747265789743441106754243191520761854837879248670580803215,5e-09,grid-euler-008 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,0.2,0.57189158307887632301630760590704965831853564166795960927667696353433626966770796,0.26957652498010691682314947763770775252248243244088144000293479237749027919981369,5e-09,grid-euler-009 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,0.81,0.37857326721656392348136924565766805407338868203483067702742825007765296443133172,0.55971511378258750082914646134609334572853197622551678181666552870111270302583261,5e-09,grid-euler-010 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,2.0,-0.084736917027655569892777470522518833388015139219570097199910234526161055494193055,0.5202689614156749366294552062065033150646009550840378847972665852495917190664362,5e-09,grid-euler-011 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,-2.0,0.090195838076920861874290243072346186698481915426986602452400235765142944592401517,-0.11643229804851741136691817453169378449406904608344005058453198393975145859535233,5e-09,grid-euler-012 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,-0.81,0.24162810986287703368318361442621217378008053868912270872643579394919480908902405,-0.41298282494202369100790147279416238500355759011044176000088211943808683985746616,5e-09,grid-euler-013 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,-0.2,1.0278469200520959476965185489276545636922874728784932448842139665862709816252963,-0.59651896884885379590730614507676265624259933232145202976213817694889856153709116,5e-09,grid-euler-014 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,0.2,1.6322590456918869664090758081243728504912613071206791648512320381349999098272371,0.32095100610252881141557478341640949082124527838871015600018068290376046712074709,5e-09,grid-euler-015 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,0.81,0.57362675062615475100317989200030737720717408749417083111276758190425564693090599,1.3160625168105963681569781266288263975410344367930192739600745799147775448337687,5e-09,grid-euler-016 +1.4,0.2,1.2,-0.3,1.1,0.4,0.25,2.0,-0.31889704984875939303257811419366887440273357576837961488751511721658583551272012,0.6718238465536092995048818615977405811482370612329458419927239247035204355776928,5e-09,grid-euler-017 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,-2.0,0.056796104015790889691777530235219571156663502992366760136463867507703097669518974,-0.11819491988831822643116424397023853117763844502589638822500712664607815061325469,5e-09,grid-euler-018 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,-0.81,-0.017545831579289031498669015804105935988775476506944815767872578181398611717939755,-0.49659327527055793189894329001718974542428192533412600278726931079447806076280718,5e-09,grid-euler-019 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,-0.2,0.8943452243589219004560133773537811860024384440520465710621420993485908671107876,-2.0557618432961449501580630499364718769101826842703411845593236948158943849442867,5e-09,grid-euler-020 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,0.2,3.8778743725002068118616776943073583990723207845348066974144123578095776941727439,0.57745221699962298809409782328455436450060014291816057233886238381947520068692779,5e-09,grid-euler-021 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,0.81,0.30868924403041517985522087692728995913357603431334066806735181554878510731428095,2.148114536244717356232364255758042688941001259660839287686228000575967697695552,5e-09,grid-euler-022 +1.4,0.2,1.2,-0.3,1.1,0.4,0.59,2.0,-0.49640095763663913143110294696379381784461865410567921422526817411391209373241077,0.69170315247127100368770287916447922857454345804401092560394953470492748838113811,5e-09,grid-euler-023 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,-2.0,0.029957751851397124910266995103223816736958813904102904394810313230393388593337593,-0.10787718145507678152801563794894669946844085443415055166910074374454027123578954,5e-09,grid-euler-024 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,-0.81,-0.26614438783662274725877251589651720991681354132859877034196472827224025652874017,-0.33410282645085792439268906460913373368305850171636100996469049903285109271592175,5e-09,grid-euler-025 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,-0.2,-3.8264216609740378233281548750437502111410428502640185860782497759895872789982356,-2.264698813995735483628731066744627969635336002084368657460766324302585433331124,5e-09,grid-euler-026 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,0.2,11.988319240741247617352676324930722763689760426670039840689951391615979427527006,9.9528581542816729408178494589541296326199534593758542714273224581253563588526457,5e-09,grid-euler-027 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,0.81,-0.80751551004513647789805010108344716510321974806236456011904185889813339370606317,2.8876564361987128568469024387182485712170735051902998779103056759394921516829678,5e-09,grid-euler-028 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,2.0,-0.68058879424143642129592199577291701765863744806778689973526291471919172665967814,0.64738795223477946424028480587953139653328628278197863249288262046875692081315299,5e-09,grid-euler-029 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,-2.0,0.008458167220550258317270231329272809228553011230554837679170479909465908838625515,-0.084496686424048481342620320017514545554248186655688624366505071912906411539402514,5e-09,grid-euler-030 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,-0.81,-0.23788487044407977522773374104077928226218681331713652218799761481224195442211219,-0.019470724828654942963007685487031839791360730556415893890771681213344856297866042,5e-09,grid-euler-031 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,-0.2,0.083770611656774960897059078148550775455289662372538534171264313646983326119224403,0.95501987921183083415284920347295524247595780316641972349787557233030009996873676,5e-09,grid-euler-032 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,0.2,-15.807253424576158187836416886366599367996048436835462960761630710216973379945232,2.1516529079751873302363607110228047274815221669026559736924831329553134946782784,5e-09,grid-euler-033 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,0.81,-2.8540851840496675351669867770546031846261886145703760737285246362223421810966603,2.0235896251189163626266460185317897005818163668169738165790219971174613610446546,5e-09,grid-euler-034 +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,2.0,-0.89459782837591147140628449447100820791395919316528705552882125931608518966834574,0.48513627412900367642821972744424665620543404030871970963167481510587435937264608,5e-09,grid-euler-035 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,-2.0,0.017782888513565311453985236648848784510294396848396290664548389894645387040473168,-0.031039317116173840489600405868216432380236198779893228929694727148956943891933445,5e-09,grid-euler-036 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,-0.81,0.028015837753142916587486064421391524646881216753663541202889458596888590462753095,-0.016455600260430782008673060743795525934045947980285460720679662572088662502823668,5e-09,grid-euler-037 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,-0.2,0.046538268001978451190236131026520918639098575412319238944729129185118763059447388,-0.02202561202404147668128862794679184549634942964143308062486022459450175493029719,5e-09,grid-euler-038 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,0.2,-0.55266812631334861648277320961295174796633160254435281660607058600412417430392432,-1.6827583861773390391312992805337679866254103018702823128861139066448345463337818,5e-09,grid-euler-039 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,0.81,-0.90154123153672036707108831329164754866275693524104435970122094713250843413694025,-1.1321405667409591772732904366112068134605053432860986674468093478386780336107474,5e-09,grid-euler-040 +1.4,0.2,1.2,-0.3,1.1,0.4,3.0,2.0,-0.76836506222723843536430824684527883395545636585732993888449230040888799397967916,-0.39759419180771467148796780677510521164268639848611671212650597907305089478485212,5e-09,grid-euler-041 +1.4,0.2,1.2,-0.3,1.1,0.4,-3.0,0.0,0.11583941932008774783213541627642056129429809859166993326573979031163490881044482,0.13843059813736838355629196875291170855021213875935338619008760586009777785966158,5e-09,grid-euler-042 +1.4,0.2,1.2,-0.3,1.1,0.4,-0.5,0.0,0.56680012656104505882126404748444030266758188239171787134356148365766866222137331,0.15388108612789431333905688433883760769977863567846953617892285262736308115735312,5e-09,grid-euler-043 +1.4,0.2,1.2,-0.3,1.1,0.4,0.0,0.0,1.0,0.0,5e-09,grid-euler-044 +1.4,0.2,1.2,-0.3,1.1,0.4,0.5,0.0,2.3689570525742544356229771029785419599380806910211600036009084204810007617704532,-1.1240281129017230048164460407455168464094521880481899371667196432467798302219429,5e-09,grid-euler-045 +1.4,0.2,1.2,-0.3,1.1,0.4,0.9,0.0,5.6735443415410716251071996070613456640803012208527723687881005508183892829501769,-27.577640082806410211358379995081252704333157524160467685028138769162718490530219,5e-09,grid-euler-046 +1.4,0.2,1.2,-0.3,1.1,0.4,0.99,0.0,-742.41637737709104357388141531251691616309449576394843157927184392759575479516211,-494.81072971286418917439114864896597261663974995315844989311124772328574515973795,5e-09,grid-euler-047 +1.4,0.2,1.2,-0.3,1.1,0.4,1.01,0.0,-2360.389012530573798518697367400267879037978656100007811845734438118201239551429,3596.3143748264234408722451924581877435759267940613240137825225733908800663378184,5e-09,grid-euler-048-upper-cut +1.4,0.2,1.2,-0.3,1.1,0.4,1.01,-0.0,102.09693700278319558337483904880770067051767022576511228375554866249205755088762,-155.64830594913405929759358789190529245661426289567120544340734658486966179860707,5e-09,grid-euler-049-lower-cut +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,0.0,-21.76783656383528424726924415029105809510430584787760429285293074237145116946548,-16.072990676318679493874496539046499696302283496416721344012290116604663959441585,5e-09,grid-euler-050-upper-cut +1.4,0.2,1.2,-0.3,1.1,0.4,1.3,-0.0,1.0175163513596397952024547750547378115440024155248769893502873443636180264891592,0.49439523977199633843618177212179458248493056233030644526401297830239556065228159,5e-09,grid-euler-051-lower-cut +1.4,0.2,1.2,-0.3,1.1,0.4,1.9,0.0,-2.5763786928403553357139067932088881852081560621096157005319483792303130543462306,-4.9287484271722742302184211653074225668834086248514411389197547804306386381384676,5e-09,grid-euler-052-upper-cut +1.4,0.2,1.2,-0.3,1.1,0.4,1.9,-0.0,0.16661670202930388408436300397488051837821295505017656676026946859112697024627586,0.060766621726430688379690363299250913086816999391087476719772702913935884645015488,5e-09,grid-euler-053-lower-cut +1.4,0.2,1.2,-0.3,1.1,0.4,2.0,0.0,-2.0476301451364429281766443487031670833001208373245871266691171669447045789530397,-4.3372029561491993013300591374768766904622159541522930597067270902109667546023499,5e-09,grid-euler-054-upper-cut +1.4,0.2,1.2,-0.3,1.1,0.4,2.0,-0.0,0.14133711187836531191965565560518946362365041602747288221523759142573629214707452,0.040962217121098022043311965904371168779691889505187664462766423996841933528476903,5e-09,grid-euler-055-lower-cut +1.4,0.2,1.2,-0.3,1.1,0.4,4.0,0.0,-0.044951327624327221276542116592141838125751053529513694583168911707512675573317904,-1.0756830948980185676328056855188572622107332449513807852123972743229898201741205,5e-09,grid-euler-056-upper-cut +1.4,0.2,1.2,-0.3,1.1,0.4,4.0,-0.0,0.030538722989833397412597327710201073150425591667145754894244848744072308316345005,-0.037509564931920996071759126219801486098870046591679158810593161026621633279145152,5e-09,grid-euler-057-lower-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,-2.0,0.80308056073514666308182844339457465112653079406475948128026168078037019316890655,-0.14337581620664273045876128678183015672728389564331165578153040072969783104887535,8e-09,grid-near-one-integer-000 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,-0.81,0.8056527413972034554076055729418591312300939315588912436353935294762653415912515,-0.093502473030975278241345143852436167295372881685709936948829535497688806164870011,8e-09,grid-near-one-integer-001 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,-0.2,0.80199438121843628282957123529380589277873679602175496750951036972362666278916473,-0.067494294162405793539582543914292860829370688414029193869819242574737052787354203,8e-09,grid-near-one-integer-002 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,0.2,0.79772195580785801080757020734241737508691784908827443718738330092418570506152279,-0.05083624416081585027968671816521775892893435408111698381489500103584409748849393,8e-09,grid-near-one-integer-003 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,0.81,0.78857758528572626130171868687977549996261885718123282907591740869329849314915324,-0.026651936244021449778427983395645594529515774801895137687553871557829492607913958,8e-09,grid-near-one-integer-004 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,2.0,0.76366178978339007886665315740234952796769653518375951577535036572292482316469939,0.014220136800743425749133474819760090594596751522418720177509681983947628204285978,8e-09,grid-near-one-integer-005 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,-2.0,0.9327516996531043111798014845050414897582339410347294487192940816125031472377189,-0.183660066763575427823600747473468525554033891974624894879369052326442833518867,8e-09,grid-near-one-integer-006 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,-0.81,0.95812567402239983716874331196080420741809227820070087490974125191474794428608076,-0.090641077510022556531365499190174355942788099312806927803256642574981848007012301,8e-09,grid-near-one-integer-007 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,-0.2,0.95396534772011931282823483262921921095272244023351875544151621258685869941447964,-0.034641234659958526883581394185526014705112282657635753628112709308798305908581875,8e-09,grid-near-one-integer-008 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,0.2,0.94275625472228053630432694255465934515090566835684645276674517512831233768178042,0.0010262906447898502663279868428935991114621916053067902360082816744803046312184193,8e-09,grid-near-one-integer-009 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,0.81,0.91518572128345026310398737189095207120655792705357084583433375830816445239156048,0.047998449868879084681170257038603776738798392918663495744828937333041694975318984,8e-09,grid-near-one-integer-010 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,2.0,0.84785667571864815619791237149978610801182852156526052795646897216187792539725036,0.10690724216089190267795576302298567183832730383777477671775411603471903223145097,8e-09,grid-near-one-integer-011 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,-2.0,0.98433165029541623939565085693737252792831695764174913413189109961487084783221235,-0.21492410706219997103247719148467349189468592598372362499420784092086060945825566,8e-09,grid-near-one-integer-012 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,-0.81,1.0374475621998179167223579341064584562421981600883108541389146185565121692469892,-0.10074804147326732092506500453143442444025113364300903390330273508966920934639022,8e-09,grid-near-one-integer-013 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,-0.2,1.0390158219604501772090264690634467763795585830884859512939765145330180244363422,-0.017390621729869434778860109894552184652851144784525076536184713479599662432237674,8e-09,grid-near-one-integer-014 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,0.2,1.0201164019203417473580615129754509323003636831767851590617628700326059401097413,0.03720963596179017739671949496097382715061032946103124945059573928059105200941481,8e-09,grid-near-one-integer-015 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,0.81,0.97029747763294284061004771366014095607852804462413142965843317557308976334392555,0.099418304585104924531056622934428781791382750788299792539728839920410373230727193,8e-09,grid-near-one-integer-016 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.25,2.0,0.87085886395446068136410534254850284377864225749299516734191642113701777653037216,0.1555681091285443492164763380682620916537896510874760117675683300318893922771664,8e-09,grid-near-one-integer-017 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,-2.0,1.0085488620976020916654990848094497644688414544915024160160623512091219163824908,-0.23489492980155650337856838935608669458237610685832797518136266824582341179190404,8e-09,grid-near-one-integer-018 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,-0.81,1.0832388949028905883403863209966951868018419294222486833360882274317144681437483,-0.11322535400671112733063832425310186568929410584889538490219229772880964726544713,8e-09,grid-near-one-integer-019 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,-0.2,1.0952260152412086313003899623371388395488394081190965587346136634786005791849884,-0.0081867146001052473156637198289069109075676739926528708194277241574422151721698408,8e-09,grid-near-one-integer-020 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,0.2,1.0683021198492990229689000303613509636445852277340208281233041522090998768098089,0.064991051936944403369357459589134332893454865496263670936075267073557069232385935,8e-09,grid-near-one-integer-021 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,0.81,0.99672537701719597479982605452727956171080944189575773223738465718470728097945334,0.13444773913910378089364660018163759382390850308696313185602136016532860925575061,8e-09,grid-near-one-integer-022 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.59,2.0,0.87841149890520607381837384527750709755753264079063853966737307883544735193403897,0.18198652625382952562856190983617154654622429776347902914240738133507254190404148,8e-09,grid-near-one-integer-023 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,-2.0,1.0300993031358607413778920535324866225180917466900431929326444089900757144686224,-0.25690838302404719734118230092681304192334494040507716893744559563811842529622196,8e-09,grid-near-one-integer-024 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,-0.81,1.13038567495743840677483817281005762752859968893800741457420393801786428625776,-0.13380417490822043974027548704529343695287281442062103377240870848464192439788481,8e-09,grid-near-one-integer-025 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,-0.2,1.1679948201198444518958778478985044262127862288645265949144637501897540944903482,-0.0027846376926365891891347082736199357307050595254052289152667209846640923730210075,8e-09,grid-near-one-integer-026 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,0.2,1.1240082958473123862794123966368502945126466951081733244536043777513939750278575,0.10774635202821580024767637728995245866809591143954878884723431027573785800770595,8e-09,grid-near-one-integer-027 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,0.81,1.0182251594686900151847581720543227472119985714997018730302392085584524538755059,0.17583256987508580366837603469386458321890150418311238971681371181360637331121388,8e-09,grid-near-one-integer-028 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,2.0,0.88270508113239074572575552777620280690046413866864366966948842186634851915103029,0.20822547428330639908165530749173668260089879153691690470478417443113653292701767,8e-09,grid-near-one-integer-029 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,-2.0,1.0556222431685503588750493454609338487420806118155623264937775570232455182747846,-0.29088002660064849507313414807378495000376858375810068368558566503634616622729439,8e-09,grid-near-one-integer-030 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,-0.81,1.1928692706001086794095945127535731400964171148800135054996437566393285506836056,-0.18013418216892235739431206936520984934113289215815307649886713832139770638206004,8e-09,grid-near-one-integer-031 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,-0.2,1.303519922192672447025516447306293994141148735607166064356462292786101440779108,-0.046376409204778911543104849497988046014626279970657703652952808865043067488638739,8e-09,grid-near-one-integer-032 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,0.2,1.1845005244352441594479297300484335157909387177936753828528374375980292150882582,0.22267783824183016031167885114322148909971106922539085048340666034981330300374999,8e-09,grid-near-one-integer-033 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,0.81,1.0340672583567554810639164454876961562320605461281855909018372831197397028481856,0.242441210585327528662398980809435727938296215112518505315622397414359109183235,8e-09,grid-near-one-integer-034 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,2.0,0.88352725645331641346214936837845978308363896913326236439993837416663573823785723,0.24426517301351951667286230697707444531217527505405003931491965242737054032861295,8e-09,grid-near-one-integer-035 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,-2.0,1.098476116581965760602877331518753081699544145772773372770472842595880762672249,-0.47761785950186298324334709111353975108235979879544717367688712121003367679657968,8e-09,grid-near-one-integer-036 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,-0.81,1.2650618901022916619365367848756303522777474980801001688940472764535681280012683,-0.50701577282423674061748834993762669269077953434170145331608585260428849244206312,8e-09,grid-near-one-integer-037 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,-0.2,1.381917055725305836855616869449309644119474425569934412109095009573540885070713,-0.55037111392004318448707051456129461839546970700185862224150669349119831898986252,8e-09,grid-near-one-integer-038 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,0.2,0.95459291558237993794958744500760734838604431773189540904467649480197400087063672,0.56408282758624916638602872180658894169393025567905174855957531649308639035049233,8e-09,grid-near-one-integer-039 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,0.81,0.90938496589585876766533282512873048509506105491570754559078287664734771246638276,0.48212832949052970798628469163617587946216350218640688781618704657685673025527945,8e-09,grid-near-one-integer-040 +0.4,0.2,0.9,-0.1,3.3,0.100000001,3.0,2.0,0.82309182440980761787897353041379129537598622145200470942829169861580043629488413,0.38585021181206307615454582321646438032295629633213492088960730307754553971151518,8e-09,grid-near-one-integer-041 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-3.0,0.0,0.80003879452732808287639091179661388613219755048239738339702641213813013877623741,-0.059105057559660194729560551624641996965198329682864665216614926725527008575850973,8e-09,grid-near-one-integer-042 +0.4,0.2,0.9,-0.1,3.3,0.100000001,-0.5,0.0,0.94917629711020990209916902907995768507211062925148807944335602561494949587127578,-0.016489686779704333585807353840880710683371519774187705410699648480946426758493721,8e-09,grid-near-one-integer-043 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.0,0.0,1.0,0.0,8e-09,grid-near-one-integer-044 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.5,0.0,1.0696907727462898875965092721563148010408368611333692833625618694918361542079067,0.024282994649894001177548281381111945695476869644797455912086399346439502399740292,8e-09,grid-near-one-integer-045 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.9,0.0,1.1574210375559813040592643945148740773360945152740005240523756520654033695538416,0.057916608334608091422501314414292368480938715909168489835506557560002050692494783,8e-09,grid-near-one-integer-046 +0.4,0.2,0.9,-0.1,3.3,0.100000001,0.99,0.0,1.1883169059256471097431136820616360735054208649954084749290388913011246839334648,0.070822323078149814667102233183267742628133849153633778903122942876810791539746013,8e-09,grid-near-one-integer-047 +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.01,0.0,1.197070289602854823561363716903748938919161781731096299245116452609355692223778,0.074882602504335497074438735212351523415835637961927804932617135644526704351166845,8e-09,grid-near-one-integer-048-upper-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.01,-0.0,1.197262597643259992760025641895424496698353522778557466470004332661366225347255,0.074532277685031795325974514804918733400580778785295037420449338456972413424153746,8e-09,grid-near-one-integer-049-lower-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,0.0,1.2700022255786815054852348331294845243763946383438044624208937809455259561153805,0.21125926510833606553735641980958665833035472822882933855706335711793058019961052,8e-09,grid-near-one-integer-050-upper-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.3,-0.0,1.3610499793092289428162108922229911223420327290132524046673559467294799105891806,0.033395986360889673487625104629407959563184246586577631199145062440667064582422844,8e-09,grid-near-one-integer-051-lower-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.9,0.0,1.1923819042931760904249360675990601255573091528642831143021511858585880670284262,0.44631726952609124981686165310417195226541269506054286366757087378320098538243406,8e-09,grid-near-one-integer-052-upper-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,1.9,-0.0,1.5005452160116712229207520562417108222983480917586126105796084617451255854211607,-0.22828305048963118665584639372827256753350196433800095911170397793495762852482657,8e-09,grid-near-one-integer-053-lower-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,2.0,0.0,1.170898096468203089243705319924751365800726344560770736613294398862179152002696,0.47123682849471226227283140079567060581689447294346097396807266555361257148693308,8e-09,grid-near-one-integer-054-upper-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,2.0,-0.0,1.5038098425109692956591809830810801284050669250923593669245301668124683989624669,-0.26955329281234982113664997111773554696727876252339435904844868960081612473713909,8e-09,grid-near-one-integer-055-lower-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,4.0,0.0,0.8176886631302220493185518032783840891697959022293106037519592465372629433658382,0.62320130438881068558698542642945261167974224652782618930902536240134784799775757,8e-09,grid-near-one-integer-056-upper-cut +0.4,0.2,0.9,-0.1,3.3,0.100000001,4.0,-0.0,1.2926052288321706229057676689050378906787042857468592459161768348810826628348966,-0.73534816689883877545304039382480752009760457627896975220365909324001941157222292,8e-09,grid-near-one-integer-057-lower-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,-2.0,0.61090133439759568132005949054102235445958371820776246533369296607979298728848082,-0.31247818535655920374296584962301059307935205375588554717319365159356980568233893,8e-09,grid-near-infinity-integer-000 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,-0.81,0.62617200004578891047711502543093549400287033039574159914299685531663834396422538,-0.23591207378902322471883009701733523752523557136339413633830484609075426493227329,8e-09,grid-near-infinity-integer-001 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,-0.2,0.62269626608611301173605528537245794374110535127887405002243974999696477829011869,-0.19552697564305164287291918142030198603722690278581824347545805994624719384014158,8e-09,grid-near-infinity-integer-002 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,0.2,0.61631512553109209675011131130290746743826738270229465321707644906347802170489473,-0.1702999696973375777112533142627195926214734734257928922897133146417135856649365,8e-09,grid-near-infinity-integer-003 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,0.81,0.60136668137481766442248135946979660961381880603564544747902931941145732300459215,-0.13534401673193843605714844839336760439489018511016459132914722231352157172735006,8e-09,grid-near-infinity-integer-004 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,2.0,0.56118488874406166612920243800908473838907245448943186304072462996481441472430939,-0.083018453087866116373916440686140938554377325119768363705723035057859915576009466,8e-09,grid-near-infinity-integer-005 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,-2.0,0.82251353404606830700962492495294173869739041283271350966815153809713483591654109,-0.44619118367595210148470280719906743046362801666836639433619239417314193223232019,8e-09,grid-near-infinity-integer-006 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,-0.81,0.92112187133138389324236269690493543984253051625962538140975907956770394820970201,-0.24875521842456415484007011857062107833972910937965093296000264628788714743940951,8e-09,grid-near-infinity-integer-007 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,-0.2,0.90495713863955457036807367656559862431870149368055340826493189729965763502906622,-0.11456244756952373518462357222932516424999253485551755669494977598561084516413446,8e-09,grid-near-infinity-integer-008 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,0.2,0.86248318366133474755848203165880800049850995743930317015973474417904271621968808,-0.039665488427799862573018966362310832809971692026008017214454241920001016618892497,8e-09,grid-near-infinity-integer-009 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,0.81,0.77636326925300983869108943989297765011663222713714890858532104652289182694888926,0.032768559051759803785261896381697162929893965939573233224372728986192976860994215,8e-09,grid-near-infinity-integer-010 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,2.0,0.63350792007226659009696991754972100852547865036181644929571761789100375363643051,0.073450287271304150963761132531926079919486332624337186883025584982815056980403956,8e-09,grid-near-infinity-integer-011 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,-2.0,0.90422524980911097863563502433129053283118293345838867932705619768348968289891836,-0.55769685507040071611438604559260201391618400807403631928023335804051601651316096,8e-09,grid-near-infinity-integer-012 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,-0.81,1.1409694456172904469959151730545189644832312238862590917405178480299210316629458,-0.29776955179172162813170164654421349871258250904052293545954173758385855912038315,8e-09,grid-near-infinity-integer-013 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,-0.2,1.1322280214719739864595563826319268430872354082588970263016038978689265453812205,-0.017240583779913119122578782208647450694783762535440310457655421815143863859400363,8e-09,grid-near-infinity-integer-014 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,0.2,1.0090995782892402236584992603598212943709120330351522717552754592105151976525729,0.1220674891543732844806526076494438477185039373985134887443780714046036647413271,8e-09,grid-near-infinity-integer-015 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,0.81,0.81798906311651628280920598583159966954784606877733509524586571153766535029218094,0.17797768786698946992156563597256877950190485455429016078439289815770578381224863,8e-09,grid-near-infinity-integer-016 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.25,2.0,0.62515963523465453281047025801648085599476155404550933354961178335851675885050722,0.14956310747469070818171373528899870651186759681376516298990305041532559885988338,8e-09,grid-near-infinity-integer-017 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,-2.0,0.93336967763927667581842667578112428184081661716068197478258106803621404339923716,-0.62971557800283140217353038930674151123552927000428579585713468570032630969279006,8e-09,grid-near-infinity-integer-018 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,-0.81,1.2919642594728628395897579307052002070783221634983691959068269710503703583733333,-0.37348630122961955916376523164524071388438682856372205796637089547762547823439142,8e-09,grid-near-infinity-integer-019 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,-0.2,1.3435633488634566724386572946111692383800752257272258711602998968053912230950311,0.083166494620728562321378171114056811030746482901821447971311515728388086162930459,8e-09,grid-near-infinity-integer-020 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,0.2,1.0784737700435960783944217109967967029970999217547998254034616102267910591504889,0.29110961341537357233002218011134402500731389796693549527699266531667634742921513,8e-09,grid-near-infinity-integer-021 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,0.81,0.80586438279421075514844436072264159569226217568558733013143918732636559727705186,0.26912106760200233022961025401235647624019813812763061294203386491493314504301822,8e-09,grid-near-infinity-integer-022 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.59,2.0,0.61075989616432837818828025410886209901424720540199112370700572938693965336430454,0.18455982212392606816843951504600840394807054776975610843607422776946721668865292,8e-09,grid-near-infinity-integer-023 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,-2.0,0.94877977886724916458047560915535190582147546129552282204507415915652450399990855,-0.70697281493687403036100002466492772217398902390477875050919113260591995327464309,8e-09,grid-near-infinity-integer-024 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,-0.81,1.4451847143467423002885814672729612560865106413172687557784242633536632181233787,-0.52188853918472539331154450251015212156982026111318359099096385227336295229630832,8e-09,grid-near-infinity-integer-025 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,-0.2,1.8199818872954431240213054406809925472150461677130170937389306019808990548590594,0.23494230616132751297292762048868526697445467391798735471841916616186099511242759,8e-09,grid-near-infinity-integer-026 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,0.2,1.0318434990591937787805400055619845357965509321511619617260162038465423945996327,0.56216983751049273728754962469205556140829111556591776955872165810295983385231833,8e-09,grid-near-infinity-integer-027 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,0.81,0.76140448236959076876039526482311571840334790887960703568300046901693462221370456,0.35166856026971910132817148675996416136777094844881717397931031790330955372692285,8e-09,grid-near-infinity-integer-028 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,2.0,0.59157354500710545478245450003721540029236335406689037014575247712308977530215556,0.21402868961926467085348439345127033645929008266895490260714120134225814950218321,8e-09,grid-near-infinity-integer-029 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,-2.0,0.94518615945101180723825921331986938140037940256877103535725178634351857457994937,-0.81738901509883393039044237028596698705545047188422769447597916224297597464222544,8e-09,grid-near-infinity-integer-030 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,-0.81,1.5519729883436235322365345430769303352764266252949906234503209887228705341637356,-0.8571932298421839674832053956011970163931578245759465753148225357533778331161315,8e-09,grid-near-infinity-integer-031 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,-0.2,2.754263053931352558669048012252187784804250161517689470358562114944827000241178,-0.89779325715073279341861622616034432968594676719043635644393974366032762435253438,8e-09,grid-near-infinity-integer-032 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,0.2,0.69951394690175403993292163971981206818296224963769859104966240229570014512578448,0.67825642123604439908735658989070381506516271402430638793330140467345924592983888,8e-09,grid-near-infinity-integer-033 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,0.81,0.66569454659822156742588487105323914672375971027681159509282438876338162383690869,0.42246108319418276791207331956825134429482191864201845426814143350588362958088485,8e-09,grid-near-infinity-integer-034 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,2.0,0.55986345404016731905328996489465919472499000033033350070158046172384790100364998,0.24575239273557359601265723198395031549207567246444270502593372007570013512161302,8e-09,grid-near-infinity-integer-035 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,-2.0,0.65399818627326626124950808264097840222481087649012781128023680200596004896618415,-1.1610922447741196271130598637996204770266286281988521062329297000688763644211932,8e-09,grid-near-infinity-integer-036 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,-0.81,0.70968807088352989841234488123754443477444617868542324645177084883866492161822071,-1.5455854244986688237435351133496471310250489214756438421452767106359823048143305,8e-09,grid-near-infinity-integer-037 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,-0.2,0.62167035543548427479889003530882771815286972733494184323101008392906331486015699,-1.8408260522073051171619956184544125009595695610139681302826613542512696921585233,8e-09,grid-near-infinity-integer-038 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,0.2,0.3329589429242914499507650363699228604402660805782121639242614025678396566950144,0.43231714379271366894243516915165622770415960267486568205864933086825652924833877,8e-09,grid-near-infinity-integer-039 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,0.81,0.37915650363871531897894144124034791844364226600936290832681084821665192638475563,0.3838828861672841501499142043908549671324988624033773305550132278542647298436368,8e-09,grid-near-infinity-integer-040 +0.4,0.2,2.400000001,0.2,3.1,-0.3,3.0,2.0,0.4128760633562692819414723315503020418600511025949120089854988006511040553677047,0.28734584443263235847028681290381603264679188887015675951252726161322853172970858,8e-09,grid-near-infinity-integer-041 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-3.0,0.0,0.61988589454784599909093830073631201191181614547067916335992535458407757287913338,-0.18273097276036494015945074473536655816839163743424199251977496218455205926309947,8e-09,grid-near-infinity-integer-042 +0.4,0.2,2.400000001,0.2,3.1,-0.3,-0.5,0.0,0.88622849475587422173910148868570451830351770863744359413853872184160299514246205,-0.074713161725428089506731842475064131521622021998825466553729356064473079606288698,8e-09,grid-near-infinity-integer-043 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.0,0.0,1.0,0.0,8e-09,grid-near-infinity-integer-044 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.5,0.0,1.1797283288285503436838092329590649025438032215101456047078111768755441032073626,0.17027002756287132936424178414708600721583628185456798288769642135376470985120184,8e-09,grid-near-infinity-integer-045 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.9,0.0,1.3950591839060606978893213546012653746323019180789079904580008351597227119324687,0.69520804137183212852631523654300382445611980314889200577377137430734959598402666,8e-09,grid-near-infinity-integer-046 +0.4,0.2,2.400000001,0.2,3.1,-0.3,0.99,0.0,1.1352772066689683898862759420812835217935023886215641012630911034450889845121193,1.1955243400250503384902986983958099620956928118396328749364739807257746813605499,8e-09,grid-near-infinity-integer-047 +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.01,0.0,0.95618818320203661238014153376285349407104638441738446188109001526167807807174283,1.025545449640732976792882133631053060037637426285857701361603822180302552625316,8e-09,grid-near-infinity-integer-048-upper-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.01,-0.0,0.78743021565761643312358536665244703149696108423470412851218118309878226619492162,3.0438278604504763052627661640890890610092235338144355565605120620908843959626457,8e-09,grid-near-infinity-integer-049-lower-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,0.0,0.64471329909295459970276935742487664670388340368435015185308187862931874427609509,0.81035488662822565836450310604747785000692022078188480333128294891389163707133326,8e-09,grid-near-infinity-integer-050-upper-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.3,-0.0,3.8448130013840030588099831165128753980993941058637516436387323785959305624938535,-1.3775183314447100009582016417145793947225613532073490164147471086655956353666824,8e-09,grid-near-infinity-integer-051-lower-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.9,0.0,0.43030217784822924701503702862620807683694915380645061293782718582721007681307284,0.61694084888085705199315955135490561543065880127919017668627252241088079575989431,8e-09,grid-near-infinity-integer-052-upper-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,1.9,-0.0,1.7054698218598970247852593041145519675782846357009565012532094271703158205195028,-2.358567462698995976049033788621689556275922232900296135243887807544359804317473,8e-09,grid-near-infinity-integer-053-lower-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,2.0,0.0,0.41204907329723787335183140262189517994198577096298685295725641206752990925621091,0.59455355122011671165533856679914534083500757156345658502529478933386209538410172,8e-09,grid-near-infinity-integer-054-upper-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,2.0,-0.0,1.5138181945867451736983969562278018845465537831469773616106874493264810698612376,-2.3378555076306686544705989078066338873821816118118080655759832732049275225655771,8e-09,grid-near-infinity-integer-055-lower-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,4.0,0.0,0.27171001478671909740589966174235672369438980745570040535030342003926408665741479,0.36396714671462320423132306127239677238010784143413957251595513599343798537427887,8e-09,grid-near-infinity-integer-056-upper-cut +0.4,0.2,2.400000001,0.2,3.1,-0.3,4.0,-0.0,0.23537033720512011947243942738226563827804111533501175993571383333885069256037003,-1.6376516302642631624747585958148989700767006926537610122058334310344732996496843,8e-09,grid-near-infinity-integer-057-lower-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,-2.0,-0.000028257330002639784230226373791628416110536267461057568390660382328308576369153052,0.000021902629257346833715357786612254327626043951066280086251912884966020731431122343,2e-07,grid-large-taylor-000 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,-0.81,-0.00037563970048227093927551152814086573565810842330016343591767408447501523049715372,-0.0002832486001752997654193678436665506022530392191327602673847522294180133972520265,2e-07,grid-large-taylor-001 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,-0.2,0.00019891580420125382451763738288158864402535357567142874035969385055244697450132743,-0.0017622530712133205634714814114820176286713285561799905476586986031259460508713499,2e-07,grid-large-taylor-002 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,0.2,0.0032467476386838852137826034434530916616997318003719404519961337871480432865646693,-0.0023543933363941061322624187297439277901438165857473693958292326323589875171409577,2e-07,grid-large-taylor-003 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,0.81,0.0078657139982737401296899412829951132135396599156857231715181004714734534610655429,0.0089144552534504429124950715781585331493583327390114902352160581721698771602911401,2e-07,grid-large-taylor-004 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,2.0,-0.023677941544061426939177306238190132008028101252471173479915711519445600093965357,-0.038465407476593912717164248172071968820989530140068064179703871456918392934489524,2e-07,grid-large-taylor-005 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,-2.0,-0.000032298546000524265046343930568195028925995044432666805120017578985289997052647116,0.000040403853279042369576607607414880262467964926988711410379134511282391603120347626,2e-07,grid-large-taylor-006 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,-0.81,0.0024350456080652526384264077399239629497353305750864585963030919642787677367077109,0.0043070525943331019097852585992447544670487943051665418193789843154018130096314942,2e-07,grid-large-taylor-007 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,-0.2,0.22786428795416952998863028557689576761030165536508933623265705310962510842694814,0.20189793810523525163398310049177621517731369679451522744951902129921957594270827,2e-07,grid-large-taylor-008 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,0.2,-2.9328944159127596492911833578124358145101060778476157193984561626537531687489709,13.264136543009975982078991636379954910693682043298395372754659649768476200446583,2e-07,grid-large-taylor-009 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,0.81,1070.9654957496939966966774730401704277728610475808955863586344990874095644347731,435.2545578277850381220627313268895082044100967374247017263019489503363646983426,2e-07,grid-large-taylor-010 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,2.0,-1900.6508923102778804669580963042780638774323368564854229483389313986597106432557,-1562.2120624733070002172527958819461486408939447242909273008299473479611538258988,2e-07,grid-large-taylor-011 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,-2.0,0.000019853670149965765025065450485658444576201719369915147679388679884865596507796776,0.0000039951696205729582446121541625719715726568557864033029260640969871726023365541925,2e-07,grid-large-taylor-012 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,-0.81,-0.00055014634684320314145012989687701723910488002887725796130246952830398122535956489,-0.00061473135636085871617148332387940965571665474928841541530718774718954450580388848,2e-07,grid-large-taylor-013 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,-0.2,15.107382667760265921336456225301369288872556227028351872933087194263513885057026,-10.632074979364895505332037644582122665213648055627565038571948279472402399766961,2e-07,grid-large-taylor-014 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,0.2,622217.51001494653880435888978974223686896622624158322601616412578505779923366031,-76611.2922378910799192930406240255347914288556193257082616892075059063651629386,2e-07,grid-large-taylor-015 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,0.81,-10910670.043939134702976444797755034112364493138897740585970884418642997647514427,-51357772.846473418220064727705055704064447654973316025545137369688461075587220398,2e-07,grid-large-taylor-016 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.25,2.0,382161.41357222831142203334243808529013229091827480579918974500493101096328391325,-122853.02702423836013972749820074669567005697516962231471022925744640652442728361,2e-07,grid-large-taylor-017 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,-2.0,0.0000093677081567662075326561644882801725083753922186131748184999681776416202471662551,-0.0000064559500852447467804556196012370736362123378200201004158216434196810371247655205,2e-07,grid-large-taylor-018 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,-0.81,-0.016940742038142853771161286512004869167832316293606694655428546922034005755587541,-0.041396732070281265120153563365866517646350590788745833152206586906247807246211583,2e-07,grid-large-taylor-019 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,-0.2,74377674.480632852262428928287760501529911561322563107883895265823806869525768447,213975346.91788725899165216088771230310993906563731612304660050229204173275097671,2e-07,grid-large-taylor-020 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,0.2,48078025085.921453563346452089645106490791794834448305889165268467461914306158722,12773808044787.99475580266437971926707334807175581112660556026613321628507500109,2e-07,grid-large-taylor-021 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,0.81,76077808167.047705103285639978396778300118196498335860239113208355609207139447479,16656868885.868003080249556477153959739869748345227070063299501491085266607458476,2e-07,grid-large-taylor-022 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.59,2.0,320231.807541368737196932747333389619169402129951652155447141806374238587516023,3838077.1680195127976887832457084484228548675982582022286947196894417224938017524,2e-07,grid-large-taylor-023 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,-2.0,0.0000020317001264831641065552034724671782892385638069744694218626280230355810631601081,-0.0000060640325844604899190752223032092286305856926141034703936538812521559565722311204,2e-07,grid-large-taylor-024 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,-0.81,-0.071872873496538260967725564861108427222595456405503795493701706227981724387278756,0.31893026328214611479826614543197992036454128571654678965687026509868841912762215,2e-07,grid-large-taylor-025 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,-0.2,-295738121129582.25968434698363100252384181752202113796222884170801344086178278135,-82572622797910.371157151219184460374298049181931244092929360845970889841588247728,2e-07,grid-large-taylor-026 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,0.2,-332504652726107384503581.14346879550370361577516790667312657014029649366931284696,592368107005281541282694.85679232367829195298328388329829169338575754929550634584,2e-07,grid-large-taylor-027 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,0.81,20551061367797.904867433042242342499787458592095288343412201844697028949965123647,-33838653721490.664131197138454703038292918740632745285184249703128981526646283328,2e-07,grid-large-taylor-028 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,2.0,-20980950.141390203958066846479064813092487727498838651652252281992213925983546771,-9920366.4821391089546758089208806321751474575061132676394327841007388601143856661,2e-07,grid-large-taylor-029 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,-2.0,-0.0000011028711961437623775857416998216984149448420988322428956450549985102639776172221,-0.0000027200565811848157370690257211141173619668935255868872208015176501664701573901682,2e-07,grid-large-taylor-030 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,-0.81,-0.017651973061814197818545127787695203717879590613840809776310762087883124792502542,-0.0063938362979585818370108722449145690206651863959905184768204528379350304151271013,2e-07,grid-large-taylor-031 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,-0.2,-172188.27384233571345909182286898160425163423214691749455154733382088926786086337,-82313.230338933373649860182691456137012005081303134768456651228013856669378745782,2e-07,grid-large-taylor-032 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,0.2,20988106226361275070911869.71446235649542701660082976904850996398029607615468255,17006708900740163526312368.542477617798677738992723525552535347024673999496337246,2e-07,grid-large-taylor-033 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,0.81,-1297160150298378.0340947582013133705081725611864951724067846910098184940900178985,1673638628579866.7872753957932218120344578300918387710559828454694094673403177801,2e-07,grid-large-taylor-034 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,2.0,114287719.21016855895974313044370775988908577330547820830347291049411165738954814,58439945.537292390036492229549075119420456366275150268646377833809546989026879132,2e-07,grid-large-taylor-035 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,-2.0,-0.000000097386107599340751197678889708586045774909937362248045228983892034920122452381503,0.000000073659964767593694706890568787503525305466007263549406912885569634633835832448519,2e-07,grid-large-taylor-036 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,-0.81,0.00000021512340045397684550084716681167275618740662813887863524711596691397523121466333,0.000000055333875628584386417345416085368651849822367606703538265843869910891872471644646,2e-07,grid-large-taylor-037 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,-0.2,0.0000001236713732353412687695573802691355708480250350882428368417220698995064604306693,-0.00000020386136654292873352694679407081238486723339394689557253273523443383026697687751,2e-07,grid-large-taylor-038 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,0.2,-1718327181582.5314157077527086909697072125013170067755525315737691145962559905598,-1507882660437.0964845160665623047461730641598706976326023427518852357684804060417,2e-07,grid-large-taylor-039 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,0.81,90529010692.827750773486909363924078227775794279281717821839093918557191329485665,-78489385445.757330322519487336228323316082638410972017025424013173460317790087321,2e-07,grid-large-taylor-040 +12.0,-3.4,10.0,-0.8,-1.6,4.6,3.0,2.0,74313735.91736455656483643126281810470485194814765758142563415308526284004382823,59261756.187867548054078207178689463457058853720622340825720521720316690659404918,2e-07,grid-large-taylor-041 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-3.0,0.0,0.0013197435039269698194927821300093256331139537165138832286389817694525864446548527,-0.0023416601888382540289634301698534000903018390470669048089765264756905016615229925,2e-07,grid-large-taylor-042 +12.0,-3.4,10.0,-0.8,-1.6,4.6,-0.5,0.0,0.99734664731533640727178351969989277810704919787351722699502993006329128886354253,1.6309876057556031551639096365000210983295549842831704315793181574623590578429169,2e-07,grid-large-taylor-043 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.0,0.0,1.0,0.0,2e-07,grid-large-taylor-044 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.5,0.0,-1775141464.4697780236347532920081762521819777435452892770096646225641324643293152,-66325420.009780359724126419673624626576389621659194464647054404410123249174184972,2e-07,grid-large-taylor-045 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.9,0.0,1044726428786281716977929428.0185567489439552676195169466992181462516518466801889,-1993639604704747919001441681.8643100928572883615733755248148826811980451940940461,2e-07,grid-large-taylor-046 +12.0,-3.4,10.0,-0.8,-1.6,4.6,0.99,0.0,-1680564017317898778197732062651713681207477360111003.2656487428747302622503538238,-334977033741570759184490560239036557790398553807477.22416230708256124585018934703,2e-07,grid-large-taylor-047 +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.01,0.0,-809402795148005507409235877234563731539885965296282804940735948.13227748872060159,1825041118795093213882717064840824216666612937206720016118229044.4140457650136188,2e-07,grid-large-taylor-048-upper-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.01,-0.0,-405565158006859075501394440157822811368.70045066423711673500068214680433421328585,-1894568609566551618165567789731557177852.6759591050334323342617890168058982529194,2e-07,grid-large-taylor-049-lower-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,0.0,127424139462222162182155326977.60137765665769848537774225736292353727886679308942,-104894168055467843040305986201.03249833188777795139806001254206902698648670764556,2e-07,grid-large-taylor-050-upper-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.3,-0.0,-40208.981176942351930519936763127293485951469435766339141629325910389226255857114,155039.38413386049624851268064118158577427481866945356758728428132813870728428522,2e-07,grid-large-taylor-051-lower-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.9,0.0,10391902571254055256.004304270966069781435292934854794943340788270262207412202817,12126605386570609430.515680200508924844225039457700844073263565943305474202363001,2e-07,grid-large-taylor-052-upper-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,1.9,-0.0,-0.000015288695887806048807471440434408574427925741091588640439145191850338733076391882,-0.0000073229444407973666836758854882746488530547647172843156674011167955169795226124481,2e-07,grid-large-taylor-053-lower-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,2.0,0.0,65772627271960583.227483950465761496387109265521595400155243611113640013678354518,1988745517124138883.2593595737993739693149619989432463834212980076966347356719009,2e-07,grid-large-taylor-054-upper-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,2.0,-0.0,-0.0000012998049575958622794259383749239976122248231308364677068424146449278262294715907,-0.0000042992343531786909933751995552482894384421736465314706252647534448372590619138513,2e-07,grid-large-taylor-055-lower-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,4.0,0.0,501356220.36441169220228920791704649176300590879045582721346258770725596954965599,3647674561.625566439550081183489545990095970530791181403199357181705632581759307,2e-07,grid-large-taylor-056-upper-cut +12.0,-3.4,10.0,-0.8,-1.6,4.6,4.0,-0.0,0.0000000093270918757223025601411493782478058542265613178776696453791188146000707368590064,-0.000000033874692448282996714357299442805318572641130096720040890204862910017603220346151,2e-07,grid-large-taylor-057-lower-cut diff --git a/test/test_math.py b/test/test_math.py index 67bffdf..7056453 100644 --- a/test/test_math.py +++ b/test/test_math.py @@ -1,10 +1,16 @@ """Test rusteq.math module against scipy/numpy""" +import csv +from pathlib import Path + import numpy as np -from scipy.special import ellipe, ellipk +import pytest +from scipy.special import ellipe, ellipk, hyp2f1 as scipy_hyp2f1 import cfsem +DATA = Path(__file__).parent / "data" + def test_ellipe(): # 64-bit version @@ -16,3 +22,235 @@ def test_ellipk(): # 64-bit version xs = np.linspace(0.0, 1.0 - 1e-7, 100) assert np.allclose(ellipk(xs), np.array([cfsem.ellipk(x) for x in xs])) + + +def test_hyp2f1_complex128_binding(): + a = np.array([0.5 + 0.2j, 1.2 - 0.3j], dtype=np.complex128) + b = np.array([1.1 - 0.1j, 0.7 + 0.4j], dtype=np.complex128) + c = np.array([2.4 + 0.3j, 2.8 - 0.2j], dtype=np.complex128) + z = np.array([0.0 - 0.0j, 2.0 + 0.5j], dtype=np.complex128) + serial = cfsem.hyp2f1(a, b, c, z, False) + parallel = cfsem.hyp2f1(a, b, c, z) + assert serial.dtype == np.complex128 + assert serial.shape == (2,) + assert serial[0] == 1.0 + 0.0j + np.testing.assert_array_equal(parallel, serial) + + +@pytest.mark.parametrize("par", [False, True]) +def test_hyp2f1_writes_optional_output_and_returns_same_array(par): + a = np.array([0.5 + 0.2j, 1.2 - 0.3j], dtype=np.complex128) + b = np.array([1.1 - 0.1j, 0.7 + 0.4j], dtype=np.complex128) + c = np.array([2.4 + 0.3j, 2.8 - 0.2j], dtype=np.complex128) + z = np.array([0.2 + 0.1j, 2.0 + 0.5j], dtype=np.complex128) + expected = cfsem.hyp2f1(a, b, c, z, par=par) + out = np.full(a.shape, np.nan + 1j * np.nan, dtype=np.complex128) + + returned = cfsem.hyp2f1(a, b, c, z, par=par, out=out) + + assert returned is out + np.testing.assert_array_equal(out, expected) + + +def test_hyp2f1_output_shape_mismatch_is_transactional(): + inputs = [np.ones((2, 3), dtype=np.complex128) for _ in range(4)] + out = np.full((3, 2), 7.0 + 8.0j, dtype=np.complex128) + + with pytest.raises(ValueError, match="same shape"): + cfsem.hyp2f1(inputs[0], inputs[1], inputs[2], inputs[3], out=out) + + np.testing.assert_array_equal(out, np.full((3, 2), 7.0 + 8.0j, dtype=np.complex128)) + + +def test_hyp2f1_rejects_non_complex128_output(): + valid = np.ones(3, dtype=np.complex128) + with pytest.raises(TypeError): + cfsem.hyp2f1(valid, valid, valid, valid, out=np.ones(3, dtype=np.float64)) + + +@pytest.mark.parametrize( + "bad", + [ + np.ones(3, dtype=np.float64), + np.ones(3, dtype=np.complex64), + ], +) +def test_hyp2f1_rejects_non_complex128_array_inputs(bad): + valid = np.ones(3, dtype=np.complex128) + with pytest.raises(TypeError): + cfsem.hyp2f1(bad, valid, valid, valid) + + +@pytest.mark.parametrize("bad", [0.5, 1, np.float64(0.5)]) +def test_hyp2f1_requires_complex_scalar_inputs(bad): + arguments = [0.5 + 0.0j, 1.1 + 0.0j, 2.4 + 0.0j, 0.2 + 0.0j] + names = ["a", "b", "c", "z"] + + for index, name in enumerate(names): + invalid = arguments.copy() + invalid[index] = bad + with pytest.raises(TypeError, match=rf"^{name} must be a complex scalar"): + cfsem.hyp2f1(invalid[0], invalid[1], invalid[2], invalid[3]) + + +def test_hyp2f1_rejects_unequal_shapes(): + matrix = np.ones((2, 3), dtype=np.complex128) + vector = np.ones(6, dtype=np.complex128) + with pytest.raises(ValueError, match="same shape"): + cfsem.hyp2f1(matrix, vector, matrix, matrix) + + +@pytest.mark.parametrize("par", [False, True]) +def test_hyp2f1_rejects_self_overlapping_output(par): + backing = np.full(4, 7.0 + 8.0j, dtype=np.complex128) + outputs = [ + np.lib.stride_tricks.as_strided(backing[:1], shape=(4,), strides=(0,), writeable=True), + np.lib.stride_tricks.as_strided( + backing, shape=(4,), strides=(backing.itemsize // 2,), writeable=True + ), + np.lib.stride_tricks.as_strided( + backing[:3], + shape=(2, 2), + strides=(backing.itemsize, backing.itemsize), + writeable=True, + ), + ] + + for out in outputs: + before = backing.copy() + with pytest.raises(ValueError, match="non-overlapping"): + cfsem.hyp2f1(0.5 + 0.2j, 1.1 - 0.1j, 2.4 + 0.3j, 0.2 + 0.1j, par=par, out=out) + np.testing.assert_array_equal(backing, before) + + +@pytest.mark.parametrize("par", [False, True]) +def test_hyp2f1_multidimensional_strided_inputs_and_output(par): + a = np.array([0.5 + 0.2j, 0.7 - 0.1j, 1.2 + 0.3j, 0.8 - 0.2j, 1.1 + 0.1j, 0.6 - 0.4j]).reshape(2, 3) + b = np.asfortranarray( + np.array([1.1 - 0.1j, 0.9 + 0.2j, 0.7 + 0.4j, 1.3 - 0.2j, 0.6 + 0.1j, 1.0 - 0.3j]).reshape(2, 3) + ) + c_storage = np.full((2, 6), np.nan + 1j * np.nan, dtype=np.complex128) + c = c_storage[:, ::2] + c[...] = np.array([2.4 + 0.3j, 2.8 - 0.2j, 3.1 + 0.1j, 2.2 + 0.4j, 2.7 - 0.3j, 3.3 + 0.2j]).reshape(2, 3) + z = np.array([0.1 + 0.1j, 0.2 - 0.2j, 0.4 + 0.1j, 0.3 - 0.1j, 0.5 + 0.2j, 0.6 - 0.1j]).reshape(2, 3)[ + :, ::-1 + ] + out_storage = np.full((2, 6), np.nan + 1j * np.nan, dtype=np.complex128) + out = out_storage[:, ::2] + expected = np.empty(a.shape, dtype=np.complex128) + for index in np.ndindex(a.shape): + expected[index] = cfsem.hyp2f1( + complex(a[index]), complex(b[index]), complex(c[index]), complex(z[index]) + ).item() + + returned = cfsem.hyp2f1(a, b, c, z, par=par, out=out) + + assert returned is out + np.testing.assert_array_equal(out, expected) + + +@pytest.mark.parametrize("par", [False, True]) +def test_hyp2f1_broadcasts_complex_scalars(par): + a = 0.5 + 0.2j + b = np.complex128(1.1 - 0.1j) + c = np.array(2.4 + 0.3j, dtype=np.complex128) + z = np.array([[0.1 + 0.1j, 0.2 - 0.2j], [0.4 + 0.1j, 0.3 - 0.1j]]) + expected = np.empty(z.shape, dtype=np.complex128) + for index in np.ndindex(z.shape): + expected[index] = cfsem.hyp2f1(a, b, c, complex(z[index]), par=False).item() + + actual = cfsem.hyp2f1(a, b, c, z, par=par) + + assert actual.shape == z.shape + np.testing.assert_array_equal(actual, expected) + + +def test_hyp2f1_all_scalar_output_shape(): + arguments = (0.5 + 0.2j, 1.1 - 0.1j, 2.4 + 0.3j, 0.2 + 0.1j) + scalar_result = cfsem.hyp2f1(*arguments) + out = np.empty((2, 3), dtype=np.complex128) + + returned = cfsem.hyp2f1(*arguments, out=out) + + assert scalar_result.shape == () + assert returned is out + np.testing.assert_array_equal(out, np.full(out.shape, scalar_result.item())) + + +def _hyp2f1_reference_rows(): + with (DATA / "hyp2f1_reference.csv").open(newline="") as stream: + records = (line for line in stream if not line.startswith("#")) + yield from csv.DictReader(records) + + +def test_hyp2f1_full_complex_mpmath_reference(): + rows = list(_hyp2f1_reference_rows()) + + def values(prefix): + return np.array( + [complex(float(row[f"{prefix}_re"]), float(row[f"{prefix}_im"])) for row in rows], + dtype=np.complex128, + ) + + actual = cfsem.hyp2f1(values("a"), values("b"), values("c"), values("z")) + expected = values("expected") + for index, row in enumerate(rows): + tolerance = float(row["rtol"]) + np.testing.assert_allclose( + actual[index], + expected[index], + rtol=tolerance, + atol=tolerance, + err_msg=row["label"], + ) + + +def test_hyp2f1_agrees_with_scipy_real_parameter_subset(): + a = np.array([0.5, 1.2, 1.0, -3.0, 0.7], dtype=np.complex128) + b = np.array([1.1, 0.7, 1.0, 1.4, 1.3], dtype=np.complex128) + c = np.array([2.4, 2.8, 4.0, 2.5, 3.1], dtype=np.complex128) + z = np.array([0.2 + 0.1j, 2.0 + 0.4j, 3.0 + 4.0j, 2.0 + 0.5j, -3.0 + 0.2j]) + actual = cfsem.hyp2f1(a, b, c, z) + expected = scipy_hyp2f1(a.real, b.real, c.real, z) + np.testing.assert_allclose(actual, expected, rtol=2e-11, atol=2e-12) + + +def test_hyp2f1_scipy_1561_regression(): + one = np.array([1.0 + 0.0j], dtype=np.complex128) + c = np.array([4.0 + 0.0j], dtype=np.complex128) + z = np.array([3.0 + 4.0j], dtype=np.complex128) + actual = cfsem.hyp2f1(one, one, c, z)[0] + expected = scipy_hyp2f1(1.0, 1.0, 4.0, 3.0 + 4.0j) + np.testing.assert_allclose(actual, expected, rtol=3e-12, atol=3e-12) + + +def test_hyp2f1_conjugation_and_branch_approach(): + a = np.array([1.2 + 0.3j], dtype=np.complex128) + b = np.array([0.7 - 0.1j], dtype=np.complex128) + c = np.array([2.8 + 0.2j], dtype=np.complex128) + z = np.array([0.4 + 0.6j], dtype=np.complex128) + value = cfsem.hyp2f1(a, b, c, z)[0] + conjugate = cfsem.hyp2f1(a.conj(), b.conj(), c.conj(), z.conj())[0] + np.testing.assert_allclose(conjugate, value.conjugate(), rtol=2e-12, atol=2e-12) + + upper_cut = cfsem.hyp2f1(a, b, c, np.array([complex(2.0, 0.0)]))[0] + lower_cut = cfsem.hyp2f1(a, b, c, np.array([complex(2.0, -0.0)]))[0] + assert upper_cut != lower_cut + upper_near = cfsem.hyp2f1(a, b, c, np.array([2.0 + 1e-10j]))[0] + lower_near = cfsem.hyp2f1(a, b, c, np.array([2.0 - 1e-10j]))[0] + np.testing.assert_allclose(upper_near, upper_cut, rtol=2e-9, atol=2e-9) + np.testing.assert_allclose(lower_near, lower_cut, rtol=2e-9, atol=2e-9) + + +def test_hyp2f1_derivative_identity(): + a = np.array([0.5 + 0.2j], dtype=np.complex128) + b = np.array([1.1 - 0.1j], dtype=np.complex128) + c = np.array([2.4 + 0.3j], dtype=np.complex128) + z = np.array([0.3 + 0.2j], dtype=np.complex128) + step = 1e-5 + plus = cfsem.hyp2f1(a, b, c, z + step)[0] + minus = cfsem.hyp2f1(a, b, c, z - step)[0] + numerical = (plus - minus) / (2.0 * step) + shifted = cfsem.hyp2f1(a + 1.0, b + 1.0, c + 1.0, z)[0] + expected = (a[0] * b[0] / c[0]) * shifted + np.testing.assert_allclose(numerical, expected, rtol=2e-9, atol=2e-10) diff --git a/tools/generate_hyp2f1_reference.py b/tools/generate_hyp2f1_reference.py new file mode 100644 index 0000000..f2e21ad --- /dev/null +++ b/tools/generate_hyp2f1_reference.py @@ -0,0 +1,201 @@ +"""Regenerate complex hyp2f1 and gamma reference fixtures with mpmath 1.3.0.""" + +from __future__ import annotations + +import csv +import math +from pathlib import Path + +import mpmath as mp + +MPMATH_VERSION = "1.3.0" +PRECISION = 100 +ROOT = Path(__file__).resolve().parents[1] +DATA = ROOT / "test" / "data" +HypRow = tuple[complex, complex, complex, complex, str, float] + + +def parts(value: mp.mpc) -> tuple[str, str]: + return mp.nstr(value.real, 80), mp.nstr(value.imag, 80) + + +def broad_hyp_rows() -> list[HypRow]: + """Build a broad cross-product of parameter families and argument regions.""" + + parameter_families = [ + ("generic", 0.5 + 0.25j, 1.25 - 0.5j, 2.0 + 0.75j, 2e-9), + ("integer-c1", 0.5 + 0.25j, 1.25 - 0.5j, 1.0 + 0.0j, 5e-9), + ("integer-c2", 1.4 + 0.2j, 0.7 - 0.3j, 2.0 + 0.0j, 5e-9), + ("euler", 1.4 + 0.2j, 1.2 - 0.3j, 1.1 + 0.4j, 5e-9), + ("near-one-integer", 0.4 + 0.2j, 0.9 - 0.1j, 3.3 + 0.100000001j, 8e-9), + ("near-infinity-integer", 0.4 + 0.2j, 2.400000001 + 0.2j, 3.1 - 0.3j, 8e-9), + ("large-taylor", 12.0 - 3.4j, 10.0 - 0.8j, -1.6 + 4.6j, 2e-7), + ] + arguments = [ + complex(real, imag) + for real in (-3.0, -0.5, 0.25, 0.59, 0.9, 1.3, 3.0) + for imag in (-2.0, -0.81, -0.2, 0.2, 0.81, 2.0) + ] + arguments.extend(complex(real, 0.0) for real in (-3.0, -0.5, 0.0, 0.5, 0.9, 0.99)) + for real in (1.01, 1.3, 1.9, 2.0, 4.0): + arguments.extend((complex(real, 0.0), complex(real, -0.0))) + + rows = [] + for family, a, b, c, tolerance in parameter_families: + for index, z in enumerate(arguments): + lip = "" + if z.real > 1.0 and z.imag == 0.0: + lip = "-lower-cut" if math.copysign(1.0, z.imag) < 0.0 else "-upper-cut" + rows.append((a, b, c, z, f"grid-{family}-{index:03d}{lip}", tolerance)) + return rows + + +def hyp_rows() -> list[HypRow]: + points = [ + (0.5 + 0.25j, 1.25 - 0.5j, 2.0 + 0.75j, 0.1 + 0.2j, "direct", 2e-13), + (-2.0 + 0j, 1.2 + 0.4j, 3.5 - 0.2j, 2.0 + 0.5j, "polynomial", 2e-13), + (0.7 + 0.2j, 1.3 - 0.1j, 2.4 + 0.3j, -3.0 + 0.4j, "pfaff", 2e-12), + (0.4 + 0.2j, 1.1 + 0.3j, 2.7 - 0.2j, 4.0 + 2.0j, "infinity", 3e-12), + (0.4 + 0.2j, 0.4 + 0.2j, 2.7 - 0.2j, 4.0 + 2.0j, "infinity-equal-m0", 8e-11), + (0.4 + 0.2j, 1.1 + 0.3j, 2.5 + 0.5j, 0.98 + 0.03j, "one", 3e-12), + (0.4 + 0.2j, 0.9 - 0.1j, 1.3 + 0.1j, 0.98 + 0.03j, "one-balanced-m0", 8e-11), + ( + 0.4 + 0.2j, + 0.9 - 0.1j, + 3.3 + 0.100000001j, + 0.97 + 0.02j, + "one-near-integer", + 8e-11, + ), + (0.4 + 0.2j, 0.9 - 0.1j, 3.300000001 + 0.1j, 0.97 + 0.02j, "one-near-integer", 8e-11), + (0.4 + 0.2j, 0.9 - 0.1j, 3.299999999 + 0.1j, 0.97 + 0.02j, "one-near-integer", 8e-11), + (0.4 + 0.2j, 0.9 - 0.1j, 3.3 + 0.099999999j, 0.97 + 0.02j, "one-near-integer", 8e-11), + ( + 0.4 + 0.2j, + 2.400000001 + 0.200000001j, + 3.1 - 0.3j, + 5.0 + 1.0j, + "infinity-near-integer", + 8e-11, + ), + (0.4 + 0.2j, 2.400000001 + 0.2j, 3.1 - 0.3j, 5.0 + 1.0j, "infinity-near-integer", 8e-11), + (0.4 + 0.2j, 2.399999999 + 0.2j, 3.1 - 0.3j, 5.0 + 1.0j, "infinity-near-integer", 8e-11), + (0.4 + 0.2j, 2.4 + 0.199999999j, 3.1 - 0.3j, 5.0 + 1.0j, "infinity-near-integer", 8e-11), + (1.4 + 0.2j, 1.2 - 0.3j, 1.1 + 0.4j, 0.6 + 0.2j, "euler", 8e-11), + (12.5 + 2.0j, 9.25 - 1.5j, 17.0 + 0.5j, 0.4 + 0.2j, "moderate", 2e-10), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.5 + 0.8660254037844386j, "taylor", 5e-12), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.5 - 0.8660254037844386j, "taylor", 5e-12), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.495 + 0.8573651497465942j, "taylor", 5e-12), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.495 - 0.8573651497465942j, "taylor", 5e-12), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.505 + 0.874685657822283j, "taylor", 5e-12), + (0.7 + 0.2j, 1.2 - 0.3j, 2.1 + 0.1j, 0.505 - 0.874685657822283j, "taylor", 5e-12), + ( + 12.0 - 3.4j, + 10.0 - 0.8j, + -1.6 + 4.6j, + 0.59 - 0.81j, + "taylor-large-parameter-regression", + 2e-8, + ), + (1.0 + 0j, 1.0 + 0j, 4.0 + 0j, 3.0 + 4.0j, "scipy-1561", 3e-12), + (1.2 + 0.3j, 0.7 - 0.1j, 2.8 + 0.2j, 2.0 + 0.0j, "upper-cut", 5e-12), + (1.2 + 0.3j, 0.7 - 0.1j, 2.8 + 0.2j, 2.0 - 0.0j, "lower-cut", 5e-12), + (1.2 + 0.3j, 0.7 - 0.1j, 2.8 + 0.2j, 1.3 + 0.0j, "pfaff-upper-cut", 5e-11), + (1.2 + 0.3j, 0.7 - 0.1j, 2.8 + 0.2j, 1.3 - 0.0j, "pfaff-lower-cut", 5e-11), + ] + return points + broad_hyp_rows() + + +def generate_hyp2f1() -> None: + path = DATA / "hyp2f1_reference.csv" + with path.open("w", newline="") as stream: + stream.write(f"# mpmath={MPMATH_VERSION}, dps={PRECISION}\n") + writer = csv.writer(stream, lineterminator="\n") + writer.writerow( + [ + "a_re", + "a_im", + "b_re", + "b_im", + "c_re", + "c_im", + "z_re", + "z_im", + "expected_re", + "expected_im", + "rtol", + "label", + ] + ) + for a, b, c, z, label, rtol in hyp_rows(): + mz = mp.mpc(z.real, z.imag) + if label.endswith("upper-cut"): + mz = mp.mpc(z.real, mp.mpf("1e-80")) + elif label.endswith("lower-cut"): + mz = mp.mpc(z.real, -mp.mpf("1e-80")) + expected = mp.hyp2f1(mp.mpc(a), mp.mpc(b), mp.mpc(c), mz) + writer.writerow( + [ + a.real, + a.imag, + b.real, + b.imag, + c.real, + c.imag, + z.real, + "-0.0" if label.endswith("lower-cut") else z.imag, + *parts(expected), + rtol, + label, + ] + ) + + +def generate_gamma() -> None: + path = DATA / "complex_gamma_reference.csv" + points = [ + 0.2 + 0.3j, + 0.2 - 0.3j, + 1.2 - 2.5j, + 1.2 + 2.5j, + 8.5 + 1.25j, + -0.3 + 0.7j, + -0.3 - 0.7j, + -4.0 + 1e-8j, + -4.0 - 1e-8j, + ] + with path.open("w", newline="") as stream: + stream.write(f"# mpmath={MPMATH_VERSION}, dps={PRECISION}\n") + writer = csv.writer(stream, lineterminator="\n") + writer.writerow( + [ + "z_re", + "z_im", + "gamma_re", + "gamma_im", + "rgamma_re", + "rgamma_im", + "digamma_re", + "digamma_im", + ] + ) + for z in points: + value = mp.mpc(z) + writer.writerow( + [ + z.real, + z.imag, + *parts(mp.gamma(value)), + *parts(mp.rgamma(value)), + *parts(mp.digamma(value)), + ] + ) + + +if __name__ == "__main__": + if mp.__version__ != MPMATH_VERSION: + raise RuntimeError(f"expected mpmath {MPMATH_VERSION}, found {mp.__version__}") + mp.mp.dps = PRECISION + DATA.mkdir(parents=True, exist_ok=True) + generate_hyp2f1() + generate_gamma() diff --git a/uv.lock b/uv.lock index c4d076a..66c82ba 100644 --- a/uv.lock +++ b/uv.lock @@ -64,6 +64,7 @@ dev = [ { name = "mkdocs" }, { name = "mkdocs-material" }, { name = "mkdocstrings", extra = ["python"] }, + { name = "mpmath" }, { name = "plotly" }, { name = "pytest" }, { name = "pytest-cov" }, @@ -89,6 +90,7 @@ dev = [ { name = "mkdocs", specifier = "==1.6.0" }, { name = "mkdocs-material", specifier = "==9.5.25" }, { name = "mkdocstrings", extras = ["python"], specifier = "==1.0.3" }, + { name = "mpmath", specifier = "==1.3.0" }, { name = "plotly", specifier = ">=6.6.0,<7" }, { name = "pytest", specifier = "==9.0.2" }, { name = "pytest-cov", specifier = "==7.1.0" },