1use runmat_builtins::{
4 BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
5 BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
6 BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
7 BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
8 BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
9 BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
10};
11use runmat_macros::runtime_builtin;
12use runmat_value::{LogicalArray, Value};
13
14use crate::builtins::common::spec::{
15 BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
16 ReductionNaN, ResidencyPolicy, ShapeRequirements,
17};
18use crate::builtins::common::tensor;
19use crate::builtins::math::optim::common::{call_function, value_to_scalar};
20use crate::builtins::math::optim::type_resolvers::numerical_integral_type;
21use crate::{build_runtime_error, BuiltinResult, RuntimeError};
22
23const NAME: &str = "quad";
24const DEFAULT_TOL: f64 = 1.0e-6;
25const MAX_DEPTH: usize = 30;
26const MAX_FUN_EVALS: usize = 100_000;
27
28const INTEGER_BOUND_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
29 id: "quad-integer-bound",
30 mode: BuiltinExtensionMode::RunMatOnly,
31 description: "quad with native-class integer integration bounds is a RunMat extension",
32 error_identifier: Some("RunMat:compatibility:QuadIntegerBoundExtension"),
33};
34const INTEGER_TOLERANCE_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
35 id: "quad-integer-tolerance",
36 mode: BuiltinExtensionMode::RunMatOnly,
37 description: "quad with a native-class integer tolerance is a RunMat extension",
38 error_identifier: Some("RunMat:compatibility:QuadIntegerToleranceExtension"),
39};
40const INTEGER_CALLBACK_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
41 id: "quad-integer-callback-result",
42 mode: BuiltinExtensionMode::RunMatOnly,
43 description: "quad with a native-class integer integrand result is a RunMat extension",
44 error_identifier: Some("RunMat:compatibility:QuadIntegerCallbackExtension"),
45};
46const LOGICAL_NUMERIC_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
47 id: "quad-logical-numeric",
48 mode: BuiltinExtensionMode::RunMatOnly,
49 description: "quad with logical bounds, tolerance, or integrand values is a RunMat extension",
50 error_identifier: Some("RunMat:compatibility:QuadLogicalNumericExtension"),
51};
52const RESIDENT_INPUT_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
53 id: "quad-resident-input",
54 mode: BuiltinExtensionMode::RunMatOnly,
55 description: "quad host fallback for explicit gpuArray values is a RunMat extension",
56 error_identifier: Some("RunMat:compatibility:QuadResidentInputExtension"),
57};
58pub const EXTENSIONS: [BuiltinExtensionDescriptor; 5] = [
59 INTEGER_BOUND_EXTENSION,
60 INTEGER_TOLERANCE_EXTENSION,
61 INTEGER_CALLBACK_EXTENSION,
62 LOGICAL_NUMERIC_EXTENSION,
63 RESIDENT_INPUT_EXTENSION,
64];
65
66const INTEGER_BOUND_INPUT: [BuiltinIntegerInputCapability; 1] = [BuiltinIntegerInputCapability {
67 name: "a or b",
68 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
69 availability: BuiltinIntegerInputAvailability::RunMatOnly,
70 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
71 notes: "The compatibility target documents single and double finite limits; typed integer bounds are gated and cross binary64 exactly.",
72}];
73const INTEGER_TOLERANCE_INPUT: [BuiltinIntegerInputCapability; 1] =
74 [BuiltinIntegerInputCapability {
75 name: "tol",
76 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
77 availability: BuiltinIntegerInputAvailability::RunMatOnly,
78 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
79 notes: "The compatibility target documents single and double tolerance; typed integer tolerance is a checked RunMat extension.",
80 }];
81const INTEGER_TRACE_INPUT: [BuiltinIntegerInputCapability; 1] = [BuiltinIntegerInputCapability {
82 name: "trace",
83 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
84 availability: BuiltinIntegerInputAvailability::Documented,
85 scalar_double: BuiltinIntegerScalarDoubleRule::Allowed,
86 notes: "The documented nonzero trace toggle reads integer scalars structurally without floating conversion.",
87}];
88const INTEGER_PARAMETER_INPUT: [BuiltinIntegerInputCapability; 1] =
89 [BuiltinIntegerInputCapability {
90 name: "p1, p2, ...",
91 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
92 availability: BuiltinIntegerInputAvailability::Documented,
93 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
94 notes: "Additional parameterized-function arguments are passed to the integrand unchanged and retain exact native integer storage.",
95 }];
96const INTEGER_CALLBACK_INPUT: [BuiltinIntegerInputCapability; 1] =
97 [BuiltinIntegerInputCapability {
98 name: "fun result",
99 classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
100 availability: BuiltinIntegerInputAvailability::RunMatOnly,
101 scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
102 notes: "Integer integrand results are gated and must cross the quadrature's binary64 boundary exactly.",
103 }];
104pub const INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 5] = [
105 BuiltinIntegerCapabilityDescriptor { form: "q = quad(fun, integer_a, integer_b, ___)", inputs: &INTEGER_BOUND_INPUT, computation_domain: BuiltinIntegerComputationDomain::FloatingPoint, output_class: BuiltinIntegerOutputClassRule::Double, overflow: BuiltinIntegerOverflowRule::Error, backend: BuiltinIntegerBackendRule::GatherFallback, overload: BuiltinIntegerOverloadKind::Multiple, notes: "Integer limits are RunMat-only and enter adaptive Simpson arithmetic only after exact checked conversion." },
106 BuiltinIntegerCapabilityDescriptor { form: "q = quad(fun, a, b, integer_tol, ___)", inputs: &INTEGER_TOLERANCE_INPUT, computation_domain: BuiltinIntegerComputationDomain::FloatingPoint, output_class: BuiltinIntegerOutputClassRule::Double, overflow: BuiltinIntegerOverflowRule::Error, backend: BuiltinIntegerBackendRule::GatherFallback, overload: BuiltinIntegerOverloadKind::StructuralParameter, notes: "Integer tolerance is a RunMat-only checked floating control." },
107 BuiltinIntegerCapabilityDescriptor { form: "q = quad(fun, a, b, tol, integer_trace, ___)", inputs: &INTEGER_TRACE_INPUT, computation_domain: BuiltinIntegerComputationDomain::Structural, output_class: BuiltinIntegerOutputClassRule::Double, overflow: BuiltinIntegerOverflowRule::NotApplicable, backend: BuiltinIntegerBackendRule::GpuRestricted, overload: BuiltinIntegerOverloadKind::StructuralParameter, notes: "Trace compares exact integer zero/nonzero state and never materializes the selector as floating point." },
108 BuiltinIntegerCapabilityDescriptor { form: "q = quad(fun, a, b, tol, trace, integer_p1, ___)", inputs: &INTEGER_PARAMETER_INPUT, computation_domain: BuiltinIntegerComputationDomain::Structural, output_class: BuiltinIntegerOutputClassRule::Double, overflow: BuiltinIntegerOverflowRule::NotApplicable, backend: BuiltinIntegerBackendRule::HostAndGpu, overload: BuiltinIntegerOverloadKind::Multiple, notes: "Parameter arguments are callback payload and preserve native integer class/storage exactly." },
109 BuiltinIntegerCapabilityDescriptor { form: "q = quad(integer_returning_fun, a, b, ___)", inputs: &INTEGER_CALLBACK_INPUT, computation_domain: BuiltinIntegerComputationDomain::FloatingPoint, output_class: BuiltinIntegerOutputClassRule::Double, overflow: BuiltinIntegerOverflowRule::Error, backend: BuiltinIntegerBackendRule::GatherFallback, overload: BuiltinIntegerOverloadKind::Multiple, notes: "RunMat-only integer integrand results convert exactly before quadrature arithmetic." },
110];
111
112const QUAD_OUTPUT_Q: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
113 name: "q",
114 ty: BuiltinParamType::NumericScalar,
115 arity: BuiltinParamArity::Required,
116 default: None,
117 description: "Numerical integral estimate.",
118}];
119
120const QUAD_OUTPUT_Q_FCNT: [BuiltinParamDescriptor; 2] = [
121 BuiltinParamDescriptor {
122 name: "q",
123 ty: BuiltinParamType::NumericScalar,
124 arity: BuiltinParamArity::Required,
125 default: None,
126 description: "Numerical integral estimate.",
127 },
128 BuiltinParamDescriptor {
129 name: "fcnt",
130 ty: BuiltinParamType::NumericScalar,
131 arity: BuiltinParamArity::Required,
132 default: None,
133 description: "Number of integrand evaluations.",
134 },
135];
136
137const QUAD_INPUTS_CORE: [BuiltinParamDescriptor; 3] = [
138 BuiltinParamDescriptor {
139 name: "fun",
140 ty: BuiltinParamType::Any,
141 arity: BuiltinParamArity::Required,
142 default: None,
143 description: "Scalar integrand callback.",
144 },
145 BuiltinParamDescriptor {
146 name: "a",
147 ty: BuiltinParamType::Any,
148 arity: BuiltinParamArity::Required,
149 default: None,
150 description: "Lower integration bound.",
151 },
152 BuiltinParamDescriptor {
153 name: "b",
154 ty: BuiltinParamType::Any,
155 arity: BuiltinParamArity::Required,
156 default: None,
157 description: "Upper integration bound.",
158 },
159];
160
161const QUAD_INPUTS_TOL_TRACE_ARGS: [BuiltinParamDescriptor; 6] = [
162 BuiltinParamDescriptor {
163 name: "fun",
164 ty: BuiltinParamType::Any,
165 arity: BuiltinParamArity::Required,
166 default: None,
167 description: "Scalar integrand callback.",
168 },
169 BuiltinParamDescriptor {
170 name: "a",
171 ty: BuiltinParamType::Any,
172 arity: BuiltinParamArity::Required,
173 default: None,
174 description: "Lower integration bound.",
175 },
176 BuiltinParamDescriptor {
177 name: "b",
178 ty: BuiltinParamType::Any,
179 arity: BuiltinParamArity::Required,
180 default: None,
181 description: "Upper integration bound.",
182 },
183 BuiltinParamDescriptor {
184 name: "tol",
185 ty: BuiltinParamType::NumericScalar,
186 arity: BuiltinParamArity::Optional,
187 default: Some("1e-6"),
188 description: "Absolute error tolerance. Empty uses the default.",
189 },
190 BuiltinParamDescriptor {
191 name: "trace",
192 ty: BuiltinParamType::Any,
193 arity: BuiltinParamArity::Optional,
194 default: Some("false"),
195 description: "Nonzero value prints legacy [fcnEvals, a, b-a, Q] trace rows.",
196 },
197 BuiltinParamDescriptor {
198 name: "p",
199 ty: BuiltinParamType::Any,
200 arity: BuiltinParamArity::Variadic,
201 default: None,
202 description: "Additional arguments forwarded to the integrand.",
203 },
204];
205
206const QUAD_SIGNATURES: [BuiltinSignatureDescriptor; 4] = [
207 BuiltinSignatureDescriptor {
208 label: "q = quad(fun, a, b)",
209 inputs: &QUAD_INPUTS_CORE,
210 outputs: &QUAD_OUTPUT_Q,
211 },
212 BuiltinSignatureDescriptor {
213 label: "q = quad(fun, a, b, tol, trace, p1, p2, ...)",
214 inputs: &QUAD_INPUTS_TOL_TRACE_ARGS,
215 outputs: &QUAD_OUTPUT_Q,
216 },
217 BuiltinSignatureDescriptor {
218 label: "[q, fcnt] = quad(fun, a, b)",
219 inputs: &QUAD_INPUTS_CORE,
220 outputs: &QUAD_OUTPUT_Q_FCNT,
221 },
222 BuiltinSignatureDescriptor {
223 label: "[q, fcnt] = quad(fun, a, b, tol, trace, p1, p2, ...)",
224 inputs: &QUAD_INPUTS_TOL_TRACE_ARGS,
225 outputs: &QUAD_OUTPUT_Q_FCNT,
226 },
227];
228
229const QUAD_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
230 code: "RM.QUAD.INVALID_ARGUMENT",
231 identifier: Some("RunMat:quad:InvalidArgument"),
232 when: "Tolerance, trace flag, or argument grammar is invalid.",
233 message: "quad: invalid argument",
234};
235
236const QUAD_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
237 code: "RM.QUAD.INVALID_INPUT",
238 identifier: Some("RunMat:quad:InvalidInput"),
239 when: "Bounds, integrand values, or adaptive solver semantics are invalid.",
240 message: "quad: invalid input",
241};
242
243const QUAD_ERROR_TOO_MANY_OUTPUTS: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
244 code: "RM.QUAD.TOO_MANY_OUTPUTS",
245 identifier: Some("RunMat:quad:TooManyOutputs"),
246 when: "`quad` is called with more than two requested output arguments.",
247 message: "quad: too many output arguments",
248};
249
250const QUAD_ERRORS: [BuiltinErrorDescriptor; 3] = [
251 QUAD_ERROR_INVALID_ARGUMENT,
252 QUAD_ERROR_INVALID_INPUT,
253 QUAD_ERROR_TOO_MANY_OUTPUTS,
254];
255
256pub const QUAD_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
257 signatures: &QUAD_SIGNATURES,
258 output_mode: BuiltinOutputMode::ByRequestedOutputCount,
259 completion_policy: BuiltinCompletionPolicy::Public,
260 errors: &QUAD_ERRORS,
261};
262
263fn quad_error_with_detail(
264 error: &'static BuiltinErrorDescriptor,
265 detail: impl AsRef<str>,
266) -> RuntimeError {
267 let detail = detail.as_ref();
268 let message = if detail.starts_with("quad:") {
269 detail.to_string()
270 } else {
271 format!("{}: {detail}", error.message)
272 };
273 let mut builder = build_runtime_error(message).with_builtin(NAME);
274 if let Some(identifier) = error.identifier {
275 builder = builder.with_identifier(identifier);
276 }
277 builder.build()
278}
279
280fn quad_map_error(err: RuntimeError, fallback: &'static BuiltinErrorDescriptor) -> RuntimeError {
281 if err.identifier().is_some() {
282 err
283 } else {
284 quad_error_with_detail(fallback, err.message())
285 }
286}
287
288fn validate_requested_outputs() -> BuiltinResult<()> {
289 if matches!(crate::output_count::current_output_count(), Some(n) if n > 2) {
290 return Err(quad_error_with_detail(
291 &QUAD_ERROR_TOO_MANY_OUTPUTS,
292 "quad: too many output arguments; maximum is 2",
293 ));
294 }
295 Ok(())
296}
297
298#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::optim::quad")]
299pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
300 name: "quad",
301 op_kind: GpuOpKind::Custom("legacy-adaptive-simpson"),
302 supported_precisions: &[],
303 broadcast: BroadcastSemantics::None,
304 provider_hooks: &[],
305 constant_strategy: ConstantStrategy::InlineLiteral,
306 residency: ResidencyPolicy::GatherImmediately,
307 nan_mode: ReductionNaN::Include,
308 two_pass_threshold: None,
309 workgroup_size: None,
310 accepts_nan_mode: false,
311 notes: "Host adaptive Simpson solver. Callback computations may use GPU-aware builtins, but the adaptive integration loop runs on the CPU.",
312};
313
314#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::optim::quad")]
315pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
316 name: "quad",
317 shape: ShapeRequirements::Any,
318 constant_strategy: ConstantStrategy::InlineLiteral,
319 elementwise: None,
320 reduction: None,
321 emits_nan: false,
322 notes:
323 "Legacy adaptive quadrature repeatedly invokes user code and terminates fusion planning.",
324};
325
326#[runtime_builtin(
327 name = "quad",
328 category = "math/optim",
329 summary = "Approximate finite scalar definite integrals using legacy adaptive Simpson quadrature.",
330 keywords = "quad,numerical integration,adaptive simpson,quadrature,function handle",
331 accel = "sink",
332 type_resolver(numerical_integral_type),
333 descriptor(crate::builtins::math::optim::quad::QUAD_DESCRIPTOR),
334 extensions(crate::builtins::math::optim::quad::EXTENSIONS),
335 integer_capabilities(crate::builtins::math::optim::quad::INTEGER_CAPABILITIES),
336 builtin_path = "crate::builtins::math::optim::quad"
337)]
338async fn quad_builtin(
339 function: Value,
340 a: Value,
341 b: Value,
342 rest: Vec<Value>,
343) -> BuiltinResult<Value> {
344 validate_requested_outputs()?;
345 ensure_quad_extensions(&a, &b, &rest)?;
346 let options = QuadOptions::parse(rest)
347 .await
348 .map_err(|err| quad_map_error(err, &QUAD_ERROR_INVALID_ARGUMENT))?;
349 let a = scalar_real("lower bound", a)
350 .await
351 .map_err(|err| quad_map_error(err, &QUAD_ERROR_INVALID_INPUT))?;
352 let b = scalar_real("upper bound", b)
353 .await
354 .map_err(|err| quad_map_error(err, &QUAD_ERROR_INVALID_INPUT))?;
355
356 let result = if a == b {
357 QuadResult {
358 q: 0.0,
359 func_count: 0,
360 }
361 } else {
362 let sign = if b < a { -1.0 } else { 1.0 };
363 let lo = a.min(b);
364 let hi = a.max(b);
365 let mut result = integrate_quad(&function, lo, hi, &options)
366 .await
367 .map_err(|err| quad_map_error(err, &QUAD_ERROR_INVALID_INPUT))?;
368 result.q *= sign;
369 result
370 };
371
372 finalize(result)
373}
374
375fn ensure_quad_extensions(a: &Value, b: &Value, rest: &[Value]) -> BuiltinResult<()> {
376 for bound in [a, b] {
377 if crate::builtins::common::validation::value_contains_native_integer_class(bound) {
378 crate::compatibility::ensure_builtin_extension_enabled(&INTEGER_BOUND_EXTENSION, NAME)?;
379 }
380 if is_logical_numeric(bound) {
381 crate::compatibility::ensure_builtin_extension_enabled(
382 &LOGICAL_NUMERIC_EXTENSION,
383 NAME,
384 )?;
385 }
386 }
387 if let Some(tol) = rest.first() {
388 if crate::builtins::common::validation::value_contains_native_integer_class(tol) {
389 crate::compatibility::ensure_builtin_extension_enabled(
390 &INTEGER_TOLERANCE_EXTENSION,
391 NAME,
392 )?;
393 }
394 if is_logical_numeric(tol) {
395 crate::compatibility::ensure_builtin_extension_enabled(
396 &LOGICAL_NUMERIC_EXTENSION,
397 NAME,
398 )?;
399 }
400 }
401 if crate::builtins::common::validation::value_contains_explicit_gpu(a)
402 || crate::builtins::common::validation::value_contains_explicit_gpu(b)
403 || rest
404 .iter()
405 .any(|value| crate::builtins::common::validation::value_contains_explicit_gpu(value))
406 {
407 crate::compatibility::ensure_builtin_extension_enabled(&RESIDENT_INPUT_EXTENSION, NAME)?;
408 }
409 Ok(())
410}
411
412fn is_logical_numeric(value: &Value) -> bool {
413 matches!(value, Value::Bool(_) | Value::LogicalArray(_))
414 || matches!(value, Value::GpuTensor(handle) if runmat_accelerate_api::handle_is_logical(handle))
415}
416
417async fn prepare_quad_floating_value(
418 label: &str,
419 value: Value,
420 integer_extension: &'static BuiltinExtensionDescriptor,
421) -> BuiltinResult<Value> {
422 if crate::builtins::common::validation::value_contains_native_integer_class(&value) {
423 crate::compatibility::ensure_builtin_extension_enabled(integer_extension, NAME)?;
424 if !crate::builtins::common::validation::native_integer_value_is_exact_f64_async(&value)
425 .await?
426 {
427 return Err(quad_error_with_detail(
428 &QUAD_ERROR_INVALID_ARGUMENT,
429 format!("integer {label} must be exactly representable as double"),
430 ));
431 }
432 }
433 if is_logical_numeric(&value) {
434 crate::compatibility::ensure_builtin_extension_enabled(&LOGICAL_NUMERIC_EXTENSION, NAME)?;
435 }
436 if crate::builtins::common::validation::value_contains_explicit_gpu(&value) {
437 crate::compatibility::ensure_builtin_extension_enabled(&RESIDENT_INPUT_EXTENSION, NAME)?;
438 }
439 crate::dispatcher::gather_if_needed_async(&value).await
440}
441
442struct QuadOptions {
443 tol: f64,
444 trace: bool,
445 extra_args: Vec<Value>,
446}
447
448impl QuadOptions {
449 async fn parse(rest: Vec<Value>) -> BuiltinResult<Self> {
450 let mut values = rest.into_iter();
451 let tol = match values.next() {
452 Some(value) => parse_optional_tol(value).await?,
453 None => DEFAULT_TOL,
454 };
455 let trace = match values.next() {
456 Some(value) => parse_optional_trace(value).await?,
457 None => false,
458 };
459 Ok(Self {
460 tol,
461 trace,
462 extra_args: values.collect(),
463 })
464 }
465}
466
467async fn parse_optional_tol(value: Value) -> BuiltinResult<f64> {
468 let value =
469 prepare_quad_floating_value("tolerance", value, &INTEGER_TOLERANCE_EXTENSION).await?;
470 if is_empty_value(&value) {
471 return Ok(DEFAULT_TOL);
472 }
473 let tol = scalar_real_sync("tolerance", value, &QUAD_ERROR_INVALID_ARGUMENT)?;
474 if tol > 0.0 {
475 Ok(tol)
476 } else {
477 Err(quad_error_with_detail(
478 &QUAD_ERROR_INVALID_ARGUMENT,
479 "tolerance must be a positive finite scalar",
480 ))
481 }
482}
483
484async fn parse_optional_trace(value: Value) -> BuiltinResult<bool> {
485 let value = crate::dispatcher::gather_if_needed_async(&value).await?;
486 if is_empty_value(&value) {
487 return Ok(false);
488 }
489 scalar_nonzero_sync("trace", value, &QUAD_ERROR_INVALID_ARGUMENT)
490}
491
492fn is_empty_value(value: &Value) -> bool {
493 match value {
494 Value::Tensor(tensor) => tensor::tensor_values_f64(tensor).is_empty(),
495 Value::LogicalArray(LogicalArray { data, .. }) => data.is_empty(),
496 _ => false,
497 }
498}
499
500async fn scalar_real(label: &str, value: Value) -> BuiltinResult<f64> {
501 let value = prepare_quad_floating_value(label, value, &INTEGER_BOUND_EXTENSION).await?;
502 scalar_real_sync(label, value, &QUAD_ERROR_INVALID_INPUT)
503}
504
505fn scalar_nonzero_sync(
506 label: &str,
507 value: Value,
508 error: &'static BuiltinErrorDescriptor,
509) -> BuiltinResult<bool> {
510 if let Some(integer) = tensor::scalar_integer_value(&value) {
511 return Ok(!integer.is_zero());
512 }
513 match value {
514 Value::Bool(flag) => Ok(flag),
515 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => Ok(data[0] != 0),
516 other => Ok(scalar_real_sync(label, other, error)? != 0.0),
517 }
518}
519
520fn scalar_real_sync(
521 label: &str,
522 value: Value,
523 error: &'static BuiltinErrorDescriptor,
524) -> BuiltinResult<f64> {
525 let parsed = match value {
526 Value::Num(n) => n,
527 Value::Int(i) => i.to_f64(),
528 Value::Bool(flag) => {
529 if flag {
530 1.0
531 } else {
532 0.0
533 }
534 }
535 Value::Tensor(tensor) if tensor::is_scalar_tensor(&tensor) => {
536 tensor::tensor_value_f64(&tensor, 0)
537 }
538 Value::LogicalArray(LogicalArray { data, .. }) if data.len() == 1 => {
539 if data[0] != 0 {
540 1.0
541 } else {
542 0.0
543 }
544 }
545 other => {
546 return Err(quad_error_with_detail(
547 error,
548 format!("{label} must be a finite real scalar, got {other:?}"),
549 ))
550 }
551 };
552 if parsed.is_finite() {
553 Ok(parsed)
554 } else {
555 Err(quad_error_with_detail(
556 error,
557 format!("{label} must be finite"),
558 ))
559 }
560}
561
562#[derive(Clone, Copy)]
563struct QuadResult {
564 q: f64,
565 func_count: usize,
566}
567
568async fn integrate_quad(
569 function: &Value,
570 a: f64,
571 b: f64,
572 options: &QuadOptions,
573) -> BuiltinResult<QuadResult> {
574 let fa = call_integrand(function, a, &options.extra_args).await?;
575 let c = midpoint(a, b);
576 let fc = call_integrand(function, c, &options.extra_args).await?;
577 let fb = call_integrand(function, b, &options.extra_args).await?;
578 let whole = simpson(a, b, fa, fc, fb);
579 let mut func_count = 3usize;
580 let mut trace = options.trace.then_some(QuadTrace);
581 let q = adaptive_simpson(
582 function,
583 &options.extra_args,
584 SimpsonState {
585 a,
586 b,
587 fa,
588 fc,
589 fb,
590 whole,
591 tol: options.tol,
592 depth: MAX_DEPTH,
593 },
594 &mut func_count,
595 &mut trace,
596 )
597 .await?;
598 Ok(QuadResult { q, func_count })
599}
600
601#[derive(Clone, Copy)]
602struct SimpsonState {
603 a: f64,
604 b: f64,
605 fa: f64,
606 fc: f64,
607 fb: f64,
608 whole: f64,
609 tol: f64,
610 depth: usize,
611}
612
613#[async_recursion::async_recursion(?Send)]
614async fn adaptive_simpson(
615 function: &Value,
616 extra_args: &[Value],
617 state: SimpsonState,
618 func_count: &mut usize,
619 trace: &mut Option<QuadTrace>,
620) -> BuiltinResult<f64> {
621 if *func_count + 2 > MAX_FUN_EVALS {
622 return Err(quad_error_with_detail(
623 &QUAD_ERROR_INVALID_INPUT,
624 "exceeded maximum function evaluations",
625 ));
626 }
627
628 let c = midpoint(state.a, state.b);
629 let d = midpoint(state.a, c);
630 let e = midpoint(c, state.b);
631 let fd = call_integrand(function, d, extra_args).await?;
632 let fe = call_integrand(function, e, extra_args).await?;
633 *func_count += 2;
634
635 let left = simpson(state.a, c, state.fa, fd, state.fc);
636 let right = simpson(c, state.b, state.fc, fe, state.fb);
637 let refined = left + right;
638 let error = refined - state.whole;
639 if let Some(trace) = trace {
640 trace.record(*func_count, state.a, state.b, refined, error);
641 }
642 if error.abs() <= 15.0 * state.tol {
643 return Ok(refined + error / 15.0);
644 }
645 if state.depth == 0 {
646 return Err(quad_error_with_detail(
647 &QUAD_ERROR_INVALID_INPUT,
648 "adaptive Simpson quadrature did not converge",
649 ));
650 }
651
652 let left_value = adaptive_simpson(
653 function,
654 extra_args,
655 SimpsonState {
656 a: state.a,
657 b: c,
658 fa: state.fa,
659 fc: fd,
660 fb: state.fc,
661 whole: left,
662 tol: state.tol * 0.5,
663 depth: state.depth - 1,
664 },
665 func_count,
666 trace,
667 )
668 .await?;
669 let right_value = adaptive_simpson(
670 function,
671 extra_args,
672 SimpsonState {
673 a: c,
674 b: state.b,
675 fa: state.fc,
676 fc: fe,
677 fb: state.fb,
678 whole: right,
679 tol: state.tol * 0.5,
680 depth: state.depth - 1,
681 },
682 func_count,
683 trace,
684 )
685 .await?;
686 Ok(left_value + right_value)
687}
688
689fn midpoint(a: f64, b: f64) -> f64 {
690 a + (b - a) * 0.5
691}
692
693fn simpson(a: f64, b: f64, fa: f64, fm: f64, fb: f64) -> f64 {
694 (b - a) * (fa + 4.0 * fm + fb) / 6.0
695}
696
697async fn call_integrand(function: &Value, x: f64, extra_args: &[Value]) -> BuiltinResult<f64> {
698 let mut args = Vec::with_capacity(1 + extra_args.len());
699 args.push(Value::Num(x));
700 args.extend(extra_args.iter().cloned());
701 let value = call_function(function, args).await?;
702 let value =
703 prepare_quad_floating_value("integrand result", value, &INTEGER_CALLBACK_EXTENSION).await?;
704 value_to_scalar(NAME, value)
705}
706
707struct QuadTrace;
708
709impl QuadTrace {
710 fn record(&mut self, func_count: usize, a: f64, b: f64, q: f64, _err: f64) {
711 crate::console::record_console_line(
712 crate::console::ConsoleStream::Stdout,
713 format!(
714 " {func_count:>5} {a:13.6e} {width:13.6e} {q:13.6e}",
715 width = b - a,
716 ),
717 );
718 }
719}
720
721fn finalize(result: QuadResult) -> BuiltinResult<Value> {
722 let q = Value::Num(result.q);
723 let fcnt = Value::Num(result.func_count as f64);
724 match crate::output_count::current_output_count() {
725 None => Ok(q),
726 Some(0) => Ok(Value::OutputList(Vec::new())),
727 Some(1) => Ok(crate::output_count::output_list_with_padding(1, vec![q])),
728 Some(2) => Ok(crate::output_count::output_list_with_padding(
729 2,
730 vec![q, fcnt],
731 )),
732 Some(_) => Err(quad_error_with_detail(
733 &QUAD_ERROR_TOO_MANY_OUTPUTS,
734 "quad: too many output arguments; maximum is 2",
735 )),
736 }
737}
738
739#[cfg(test)]
740mod tests {
741 use super::*;
742 use crate::builtins::common::test_support;
743 use futures::executor::block_on;
744 use runmat_accelerate_api::HostTensorView;
745 use runmat_value::{IntegerStorage, Tensor};
746 use std::sync::Arc;
747
748 #[test]
749 fn quad_integrates_sine_with_default_tolerance() {
750 let result = block_on(quad_builtin(
751 Value::FunctionHandle("sin".into()),
752 Value::Num(0.0),
753 Value::Num(std::f64::consts::PI),
754 Vec::new(),
755 ))
756 .expect("quad");
757 match result {
758 Value::Num(value) => assert!((value - 2.0).abs() < 1.0e-6),
759 other => panic!("unexpected value {other:?}"),
760 }
761 }
762
763 #[test]
764 fn quad_respects_tighter_tolerance_on_polynomial() {
765 let _invoker = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
766 |_function, args, requested_outputs| {
767 assert_eq!(requested_outputs, 1);
768 let x = match &args[0] {
769 Value::Num(value) => *value,
770 other => panic!("expected x, got {other:?}"),
771 };
772 Box::pin(async move { Ok(Value::Num(x * x)) })
773 },
774 )));
775
776 let result = block_on(quad_builtin(
777 Value::BoundFunctionHandle {
778 name: "square".to_string(),
779 function: 7,
780 },
781 Value::Num(0.0),
782 Value::Num(1.0),
783 vec![Value::Num(1.0e-10)],
784 ))
785 .expect("quad");
786 match result {
787 Value::Num(value) => assert!((value - (1.0 / 3.0)).abs() < 1.0e-10),
788 other => panic!("unexpected value {other:?}"),
789 }
790 }
791
792 #[test]
793 fn quad_two_outputs_include_function_count() {
794 let _guard = crate::output_count::push_output_count(Some(2));
795 let result = block_on(quad_builtin(
796 Value::FunctionHandle("sin".into()),
797 Value::Num(0.0),
798 Value::Num(std::f64::consts::PI),
799 Vec::new(),
800 ))
801 .expect("quad");
802 match result {
803 Value::OutputList(outputs) => {
804 assert_eq!(outputs.len(), 2);
805 assert!(matches!(&outputs[0], Value::Num(value) if (value - 2.0).abs() < 1.0e-6));
806 assert!(matches!(&outputs[1], Value::Num(fcnt) if *fcnt >= 5.0));
807 }
808 other => panic!("unexpected value {other:?}"),
809 }
810 }
811
812 #[test]
813 fn quad_trace_records_rows() {
814 crate::console::reset_thread_buffer();
815 let result = block_on(quad_builtin(
816 Value::FunctionHandle("sin".into()),
817 Value::Num(0.0),
818 Value::Num(std::f64::consts::PI),
819 vec![Value::Num(1.0e-6), Value::Num(1.0)],
820 ))
821 .expect("quad");
822 assert!(matches!(result, Value::Num(_)));
823
824 let joined = crate::console::take_thread_buffer()
825 .into_iter()
826 .map(|entry| entry.text)
827 .collect::<String>();
828 let first_row: Vec<&str> = joined
829 .lines()
830 .next()
831 .expect("expected at least one trace row")
832 .split_whitespace()
833 .collect();
834 assert_eq!(first_row.len(), 4, "{joined}");
835 assert_eq!(first_row[0], "5", "{joined}");
836 assert!((first_row[1].parse::<f64>().unwrap() - 0.0).abs() < 1.0e-12);
837 assert!(
838 (first_row[2].parse::<f64>().unwrap() - std::f64::consts::PI).abs() < 1.0e-6,
839 "{joined}"
840 );
841 }
842
843 #[test]
844 fn quad_tol_and_trace_read_typed_integer_storage_exactly() {
845 let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
846 crate::console::reset_thread_buffer();
847 let tol = Tensor::new_integer(IntegerStorage::U16(vec![1]), vec![1, 1]).expect("tol");
848 let trace = Tensor::new_integer(IntegerStorage::U16(vec![1]), vec![1, 1]).expect("trace");
849
850 let result = block_on(quad_builtin(
851 Value::FunctionHandle("sin".into()),
852 Value::Num(0.0),
853 Value::Num(1.0),
854 vec![Value::Tensor(tol), Value::Tensor(trace)],
855 ))
856 .expect("quad");
857 assert!(matches!(result, Value::Num(_)));
858 assert!(
859 !crate::console::take_thread_buffer().is_empty(),
860 "typed trace value should enable trace output"
861 );
862 }
863
864 #[test]
865 fn quad_strict_mode_rejects_integer_tolerance() {
866 let _strict = crate::compatibility::push_runmat_extensions_enabled(false);
867 let tolerance = Tensor::new_integer(IntegerStorage::U16(vec![1]), vec![1, 1]).unwrap();
868
869 let error = block_on(quad_builtin(
870 Value::FunctionHandle("sin".into()),
871 Value::Num(0.0),
872 Value::Num(1.0),
873 vec![Value::Tensor(tolerance)],
874 ))
875 .expect_err("integer tolerance is a RunMat-only extension");
876
877 assert_eq!(
878 error.identifier(),
879 INTEGER_TOLERANCE_EXTENSION.error_identifier
880 );
881 }
882
883 #[test]
884 fn quad_rejects_wide_integer_bound_before_float_conversion() {
885 let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
886 let bound =
887 Tensor::new_integer(IntegerStorage::U64(vec![(1_u64 << 53) + 1]), vec![1, 1]).unwrap();
888
889 let error = block_on(quad_builtin(
890 Value::FunctionHandle("sin".into()),
891 Value::Num(0.0),
892 Value::Tensor(bound),
893 Vec::new(),
894 ))
895 .expect_err("wide integer bound cannot cross exactly");
896
897 assert!(error.message().contains("exactly representable"));
898 }
899
900 #[test]
901 fn quad_integer_trace_uses_exact_zero_test_in_strict_mode() {
902 let _strict = crate::compatibility::push_runmat_extensions_enabled(false);
903 crate::console::reset_thread_buffer();
904 let trace = Tensor::new_integer(IntegerStorage::U64(vec![u64::MAX]), vec![1, 1]).unwrap();
905
906 let result = block_on(quad_builtin(
907 Value::FunctionHandle("sin".into()),
908 Value::Num(0.0),
909 Value::Num(1.0),
910 vec![
911 Value::Tensor(Tensor::zeros(vec![0, 0])),
912 Value::Tensor(trace),
913 ],
914 ))
915 .expect("documented integer trace selector");
916
917 assert!(matches!(result, Value::Num(_)));
918 assert!(!crate::console::take_thread_buffer().is_empty());
919 }
920
921 #[test]
922 fn quad_passes_wide_integer_extra_parameter_exactly() {
923 let _strict = crate::compatibility::push_runmat_extensions_enabled(false);
924 let expected = u64::MAX;
925 let _invoker = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
926 move |_function, args, _requested_outputs| {
927 let Value::Tensor(parameter) = &args[1] else {
928 panic!("expected integer parameter")
929 };
930 assert!(matches!(
931 parameter.numeric_value_at(0),
932 Some(runmat_value::NumericScalar::U64(value)) if value == expected
933 ));
934 let Value::Num(x) = args[0] else {
935 panic!("expected scalar quadrature point")
936 };
937 Box::pin(async move { Ok(Value::Num(x)) })
938 },
939 )));
940 let parameter =
941 Tensor::new_integer(IntegerStorage::U64(vec![expected]), vec![1, 1]).unwrap();
942
943 let result = block_on(quad_builtin(
944 Value::BoundFunctionHandle {
945 name: "parameterized".to_string(),
946 function: 902,
947 },
948 Value::Num(0.0),
949 Value::Num(1.0),
950 vec![
951 Value::Tensor(Tensor::zeros(vec![0, 0])),
952 Value::Num(0.0),
953 Value::Tensor(parameter),
954 ],
955 ))
956 .expect("exact callback parameter");
957
958 assert!(matches!(result, Value::Num(value) if (value - 0.5).abs() < 1.0e-6));
959 }
960
961 #[test]
962 fn quad_automatic_resident_bound_gathers_but_explicit_bound_is_gated() {
963 test_support::with_test_provider(|provider| {
964 let values = [1.0];
965 let shape = [1, 1];
966 let automatic = provider
967 .upload(&HostTensorView {
968 data: &values,
969 shape: &shape,
970 })
971 .expect("automatic upload");
972 let automatic =
973 automatic.with_provenance(runmat_accelerate_api::GpuHandleProvenance::Automatic);
974 let result = block_on(quad_builtin(
975 Value::FunctionHandle("sin".into()),
976 Value::Num(0.0),
977 Value::GpuTensor(automatic),
978 Vec::new(),
979 ))
980 .expect("automatic resident bound gathers");
981 assert!(matches!(result, Value::Num(_)));
982
983 let explicit = provider
984 .upload(&HostTensorView {
985 data: &values,
986 shape: &shape,
987 })
988 .expect("explicit upload");
989 let explicit =
990 explicit.with_provenance(runmat_accelerate_api::GpuHandleProvenance::Explicit);
991 let _strict = crate::compatibility::push_runmat_extensions_enabled(false);
992 let error = block_on(quad_builtin(
993 Value::FunctionHandle("sin".into()),
994 Value::Num(0.0),
995 Value::GpuTensor(explicit),
996 Vec::new(),
997 ))
998 .expect_err("explicit resident bound is gated before fallback");
999 assert_eq!(
1000 error.identifier(),
1001 RESIDENT_INPUT_EXTENSION.error_identifier
1002 );
1003 });
1004 }
1005
1006 #[test]
1007 fn quad_forwards_extra_arguments_after_tol_and_trace() {
1008 let _invoker = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
1009 |function, args, requested_outputs| {
1010 assert_eq!(function, 42);
1011 assert_eq!(requested_outputs, 1);
1012 assert_eq!(args.len(), 2);
1013 let x = match &args[0] {
1014 Value::Num(value) => *value,
1015 other => panic!("expected x, got {other:?}"),
1016 };
1017 let scale = match &args[1] {
1018 Value::Num(value) => *value,
1019 other => panic!("expected scale, got {other:?}"),
1020 };
1021 Box::pin(async move { Ok(Value::Num(scale * x)) })
1022 },
1023 )));
1024
1025 let result = block_on(quad_builtin(
1026 Value::BoundFunctionHandle {
1027 name: "scaled_line".to_string(),
1028 function: 42,
1029 },
1030 Value::Num(0.0),
1031 Value::Num(2.0),
1032 vec![
1033 Value::Tensor(Tensor::zeros(vec![0, 0])),
1034 Value::Tensor(Tensor::zeros(vec![0, 0])),
1035 Value::Num(3.0),
1036 ],
1037 ))
1038 .expect("quad");
1039 match result {
1040 Value::Num(value) => assert!((value - 6.0).abs() < 1.0e-8),
1041 other => panic!("unexpected value {other:?}"),
1042 }
1043 }
1044
1045 #[test]
1046 fn quad_handles_oscillatory_integrand() {
1047 let result = block_on(quad_builtin(
1048 Value::FunctionHandle("sin".into()),
1049 Value::Num(0.0),
1050 Value::Num(2.0 * std::f64::consts::PI),
1051 vec![Value::Num(1.0e-8)],
1052 ))
1053 .expect("quad");
1054 match result {
1055 Value::Num(value) => assert!(value.abs() < 1.0e-8),
1056 other => panic!("unexpected value {other:?}"),
1057 }
1058 }
1059
1060 #[test]
1061 fn quad_handles_integrable_endpoint_shape() {
1062 let _invoker = crate::user_functions::install_semantic_function_invoker(Some(Arc::new(
1063 |_function, args, _requested_outputs| {
1064 let x = match &args[0] {
1065 Value::Num(value) => *value,
1066 other => panic!("expected x, got {other:?}"),
1067 };
1068 Box::pin(async move { Ok(Value::Num(x.sqrt())) })
1069 },
1070 )));
1071
1072 let result = block_on(quad_builtin(
1073 Value::BoundFunctionHandle {
1074 name: "sqrt_fn".to_string(),
1075 function: 9,
1076 },
1077 Value::Num(0.0),
1078 Value::Num(1.0),
1079 vec![Value::Num(1.0e-7)],
1080 ))
1081 .expect("quad");
1082 match result {
1083 Value::Num(value) => assert!((value - (2.0 / 3.0)).abs() < 1.0e-6),
1084 other => panic!("unexpected value {other:?}"),
1085 }
1086 }
1087
1088 #[test]
1089 fn quad_reversed_bounds_negate_result() {
1090 let result = block_on(quad_builtin(
1091 Value::FunctionHandle("sin".into()),
1092 Value::Num(std::f64::consts::PI),
1093 Value::Num(0.0),
1094 Vec::new(),
1095 ))
1096 .expect("quad");
1097 match result {
1098 Value::Num(value) => assert!((value + 2.0).abs() < 1.0e-6),
1099 other => panic!("unexpected value {other:?}"),
1100 }
1101 }
1102
1103 #[test]
1104 fn quad_rejects_more_than_two_outputs() {
1105 let _guard = crate::output_count::push_output_count(Some(3));
1106 let err = block_on(quad_builtin(
1107 Value::FunctionHandle("sin".into()),
1108 Value::Num(0.0),
1109 Value::Num(1.0),
1110 Vec::new(),
1111 ))
1112 .expect_err("too many outputs should fail");
1113 assert_eq!(err.identifier(), Some("RunMat:quad:TooManyOutputs"));
1114 }
1115
1116 #[test]
1117 fn quad_descriptor_signatures_cover_legacy_forms() {
1118 let labels: Vec<&str> = QUAD_DESCRIPTOR
1119 .signatures
1120 .iter()
1121 .map(|signature| signature.label)
1122 .collect();
1123 assert_eq!(
1124 labels,
1125 vec![
1126 "q = quad(fun, a, b)",
1127 "q = quad(fun, a, b, tol, trace, p1, p2, ...)",
1128 "[q, fcnt] = quad(fun, a, b)",
1129 "[q, fcnt] = quad(fun, a, b, tol, trace, p1, p2, ...)",
1130 ]
1131 );
1132
1133 let codes: Vec<&str> = QUAD_DESCRIPTOR
1134 .errors
1135 .iter()
1136 .map(|error| error.code)
1137 .collect();
1138 assert_eq!(
1139 codes,
1140 vec![
1141 "RM.QUAD.INVALID_ARGUMENT",
1142 "RM.QUAD.INVALID_INPUT",
1143 "RM.QUAD.TOO_MANY_OUTPUTS",
1144 ]
1145 );
1146 }
1147}