diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 444f8dc..ce6dbe0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -65,8 +65,8 @@ jobs: steps: - uses: actions/checkout@v4 - - name: Check feature configuration - run: cargo check ${{ matrix.args }} + - name: Test feature configuration + run: cargo test ${{ matrix.args }} msrv: name: MSRV / Rust 1.85 diff --git a/Cargo.toml b/Cargo.toml index ac615c7..bd0309b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "i_float" -version = "4.1.0" +version = "5.0.0" authors = ["Nail Sharipov "] edition = "2024" rust-version = "1.85" @@ -20,3 +20,11 @@ serde = ["dep:serde", "core"] serde = { version = "^1.0", default-features = false, features = ["derive"], optional = true } glam = { optional = true, version = ">=0.27" } libm = "^0.2" + +[[example]] +name = "cordic_arc_bench" +required-features = ["core"] + +[[example]] +name = "cordic_precision" +required-features = ["core"] diff --git a/README.md b/README.md index 7149d00..0b5d552 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ The crate is `no_std` and supports `i16`, `i32`, and `i64` coordinate types. ```toml [dependencies] -i_float = "4.0" +i_float = "5.0" ``` The default `core` feature exposes the complete numeric and geometry API. @@ -68,21 +68,73 @@ fit. Conversely, an `IntVector` constructed directly from arbitrary wide values is not covered by the point-coordinate bound. Floating-point input should normally be mapped with `FloatPointAdapter`. For an -explicit general-purpose safety margin, use `with_coordinate_bits` with at most -`I::BITS - 3`; algorithms with stronger range analysis may select a larger bit +explicit general-purpose safety margin, use the conservative constructors. Their +`CONSERVATIVE_COORDINATE_BITS = I::BITS - 3` budget reserves an extra bit for +rounding within the arithmetic range; algorithms with stronger range analysis may select a larger bit budget. +### Fast normalization + +`IntVector::fast_normalize()` returns an approximate `UnitIntVector`, or +`None` for a zero vector. It favors speed over precision: about 6, 14, or 30 bits +of direction precision for `i16`, `i32`, or `i64`, respectively. +The length is at most one. Normalization uses `sqr_length()` under the same +arithmetic-range contract as other vector operations on point differences. +It shifts the squared length to retain fractional precision in the reciprocal, +then applies that reciprocal to the original components. + +```rust +use i_float::int::vector::IntVector; + +let direction = IntVector::::new(3, 4).fast_normalize().unwrap(); +let offset = direction * 10; +assert_eq!(offset, IntVector::::new(6, 8)); +assert!(IntVector::::new(0, 0).fast_normalize().is_none()); +``` + +`UnitIntVector::x()` and `y()` return stored integers with scale `2^14`, `2^30`, +or `2^62` for `i16`, `i32`, or `i64`, respectively. Multiplication by a scalar +of type `T` (or `.scale(scalar)`) returns an `IntVector`, rounding to the +nearest integer with midpoint values away from zero. Approximation error in +the direction grows with the magnitude of the scalar. + +## Floating-point coordinate range + +Floating-point geometry supports finite input coordinates and rectangle bounds +within these inclusive limits: + +| Scalar | Maximum absolute coordinate | +| --- | --- | +| `f32` | `2^60` (approximately `1.15e18`) | +| `f64` | `2^500` (approximately `3.27e150`) | + +These limits leave headroom for point differences, their dot and cross products, +squared lengths, and midpoints without overflow. They do not guarantee exact +arithmetic: ordinary floating-point rounding, cancellation, and underflow still +apply. Arbitrary scaling and repeated operations require their own range analysis. + +`FloatPoint::normalize` and `FloatPointMath::normalize` additionally require a +positive, finite, normal squared length: at least `f32::MIN_POSITIVE` or +`f64::MIN_POSITIVE`. A nonzero vector alone is insufficient because squaring tiny +components can underflow. Rescale such vectors before normalizing them. + ## Floating-point adapter `FloatPointAdapter` maps a bounded floating-point coordinate space onto an integer grid. The same adapter converts results back into the original space. +Its input bounds must satisfy the floating-point coordinate range above and +have `min <= max` on each axis. All adapter constructors validate bounds, +including rectangles assembled through public fields. Fallible constructors +return `FloatPointAdapterScaleError::InvalidRect`; infallible constructors panic. +Checked point conversions validate membership in the original rectangle or +the enclosing integer grid. ```rust use i_float::adapter::FloatPointAdapter; use i_float::float::rect::FloatRect; use i_float::int::point::IntPoint; -let bounds = FloatRect::new(-10.0_f64, 10.0, -5.0, 5.0); +let bounds = FloatRect::new(-10.0_f64, 10.0, -5.0, 5.0).unwrap(); let adapter = FloatPointAdapter::<[f64; 2], i32>::new(bounds); let source = [2.5, -1.25]; @@ -94,12 +146,58 @@ assert!((restored[0] - source[0]).abs() <= tolerance); assert!((restored[1] - source[1]).abs() <= tolerance); ``` +Use `new_conservative(rect)` or `with_iter_conservative(iter)` for automatic +scaling with the conservative coordinate range. Use +`try_with_scale_conservative(rect, scale)` or +`try_with_iter_and_scale_conservative(iter, scale)` to validate an explicit scale +against the same budget. The associated `CONSERVATIVE_COORDINATE_BITS` constant +is `I::BITS - 3`: converted coordinates stay within the inclusive range +`[-2^(I::BITS - 3), 2^(I::BITS - 3)]`, leaving an extra bit for rounding inside +the strict point arithmetic range. + Use `with_coordinate_bits` when an algorithm has an explicit coordinate-bit budget. The value controls only the converted coordinate magnitude; it does not prove that every later arithmetic expression is safe. Use `try_with_scale` or `try_with_scale_and_coordinate_bits` when a caller supplies the scale and invalid or unsafe scales must be rejected. +For input iterators, use `with_iter_and_coordinate_bits` or +`try_with_iter_and_scale_and_coordinate_bits`. Both accept any `IntNumber` +implementation, so downstream algorithms can share these constructors while +choosing their own bit budget: + +```rust +use i_float::adapter::FloatPointAdapter; + +let points = [[-3.0_f64, -1.0], [3.0, 1.0]]; +let adapter = FloatPointAdapter::<[f64; 2], i32>::with_iter_and_coordinate_bits( + points.iter(), i32::BITS - 3, +); +let fixed = FloatPointAdapter::<[f64; 2], i32>::try_with_iter_and_scale_and_coordinate_bits( + points.iter(), 100.0, i32::BITS - 3, +)?; +assert_eq!(fixed.dir_scale(), 100.0); +# Ok::<(), i_float::adapter::FloatPointAdapterScaleError>(()) +``` + +`new` and `with_iter` return the adapter directly and panic for invalid bounds +or input points. An empty iterator uses zero bounds and scale one. + +Automatic scales are capped at the largest finite power of two of the scalar +type (`2^127` for `f32`, `2^1023` for `f64`), trading precision for finite scales +on very small bounds. Within the supported coordinate range, automatic scales +have finite reciprocals. Checked constructors preserve explicit scales and return +`ScaleTooSmall` if the reciprocal is non-finite. `with_scale` panics for invalid +scales. `with_coordinate_bits` panics for invalid bounds or a bit budget greater +than `I::BITS - 2`. + +The adapter retains the original floating-point bounds for input validation. +Integer-to-float conversion checks an internal `IntRect` enclosing those bounds +on the selected grid (minimum rounded down, maximum rounded up). Snapping a valid +input can therefore return a grid point just outside the original float bounds. +`rect()` continues to return the original bounds; `try_snap_to_grid` rejects +inputs outside them. + ## Fixed-scale ratios `UnitRatio` represents a value in the inclusive range `0..=1`. Its stored @@ -126,6 +224,20 @@ stored value between zero and `DENOMINATOR`, `from_float` expects a finite value between zero and one, and `from_int` expects `0 <= numerator <= denominator`. These preconditions are checked by debug assertions. +## Integer rotations + +`int::angle::{Angle, Rotation}` provides integer CORDIC angle measurement and +reusable rotation matrices for `UnitIntVector`. Arc subdivision and storage stay +in the consumer. See the [API, accuracy results, and reproducible arc benchmark](docs/cordic.md). +`Rotation::::with_precision(angle, 4)` selects fewer iterations with an angle/16 +error budget (3 gives angle/8). `Rotation` stores coefficients in `I`, sharing +the vector scale (Q14/Q30/Q62); i64 retains Q30 kernel precision. `rotation.angle()` +reports the achieved step for calculating counts and remainders. +For runtime iteration/precision controls, run +`cargo run --release --example cordic_precision -- --sweep`; +or use `--relative 16 --step 20` to select rotation iterations from a relative +angular tolerance. See the [precision experiment](docs/cordic-precision.md). + ## Features | Feature | Default | Description | diff --git a/docs/cordic-precision.md b/docs/cordic-precision.md new file mode 100644 index 0000000..e006c22 --- /dev/null +++ b/docs/cordic-precision.md @@ -0,0 +1,231 @@ +# Runtime precision experiment + +## Relative error of the rotation step + +```sh +cargo run --release --example cordic_precision -- --relative 16 --step 20 +cargo run --release --example cordic_precision -- --relative 8 --step 20 +``` + +`--relative 16` targets angular error no greater than the requested rotation +angle divided by 16 (6.25%); `--relative 8` allows 12.5%. Iterations are chosen +for **each step angle**, rather than setting a single count for every arc. +This mode keeps vectoring/atan2 at the public 29-iteration precision and matrix +coefficients at Q30. Only rotation CORDIC is shortened. It cannot be combined +with explicit `--iterations` or `--q`. + +The selector uses the CORDIC residual bound `2^-(n-1)` radians, converted to +binary angle units, and reserves two binary units for coefficient/table +rounding. It finds the smallest sufficient iteration count with shifts, +leading-bit inspection, and one integer comparison; no runtime logarithm, +float arithmetic, or division is needed. Clockwise angles use their unsigned +magnitude. For extremely tiny angles whose budget is below the fixed-precision +floor, it falls back to 29 iterations without claiming the relative bound. +That floor is outside this experiment's arc-step range. + +| Requested step | Iterations, step/8 | Iterations, step/16 | +|---|---:|---:| +| 45 degrees | 5 | 6 | +| 20 degrees | 6 | 7 | +| 5 degrees | 8 | 9 | +| 1 degree | 10 | 11 | +| 360/1024 degrees | 12 | 13 | + +This is a tolerance on **one matrix's angle**, not a bound on the phase error +of all generated arc points. Repeating a matrix repeats its angular error with +the same sign. With fixed reference counts, after `n` steps the phase error +is approximately `n * step_error`, even though radial drift stays small in Q30. +The CSV deliberately reports this effect instead of replacing the last computed +point. The 1024-step stress test can cover many full turns for coarse steps; +its displacement is not a typical error for a short 20-degree-step arc. + +For this i32 experiment, the public API is `Rotation::::with_precision(angle, precision)`, where the +parameter is the exponent: 3 means angle/8 and 4 means angle/16. It shares this +experiment's iteration selector. All `u32` exponents are accepted; >=32 uses +full precision. The conservative matrix-angle error bound in binary units is +`max(floor(|angle| / 2^precision), 5)`. + +The consumer obtains the **achieved angle** with `rotation.angle()` and uses +it for the number of rotations and the final remainder. CORDIC already tracks +the residual, so another atan2 is unnecessary. The returned angle differs from +the Q30 matrix angle by at most `Rotation::::ANGLE_MAX_ERROR` (2 binary units). +Tiny requests can report zero. A strict maximum-step policy also needs to +budget for rounding and bias the requested step inward; a consumer regression +test demonstrates this for i32 storage and exponents 3/4. Public +`Rotation::new` retains all 29 iterations; iOverlay is unchanged. + +## Explicit iteration and coefficient controls + +From the iFloat checkout: + +```sh +cargo run --release --example cordic_precision -- --iterations 24 --q 30 --step 1 --radius 65536 +cargo run --release --example cordic_precision -- --sweep > precision.csv +cargo run --release --example cordic_precision -- --sweep --reverse > precision-reverse.csv +cargo run --release --example cordic_precision -- --help +``` + +`--iterations` is a runtime integer in 1..29; it changes both vectoring and +rotation. `--q` is a runtime coefficient precision in 14..30. Neither requires +recompilation; after the initial build the executable in `target/release/examples/` +can be invoked directly with different values. The default is 29 iterations, +Q30, step 1 degree, radius 65536. Step accepts 0.3515625..45 degrees, radius +accepts integers 1..2^30. `--sweep` compares twelve profiles, varying iterations and +coefficient bits independently and together, plus two relative-error profiles; it cannot be combined with explicit +`--iterations`, `--q`, or `--relative`. Invalid options exit with status 2. `--reverse` reverses +the measurement order, including the fixed-library reference. + +The example includes `src/int/angle/cordic.rs` by source path. The measured +CORDIC loops and lookup tables are the same ones used by the library, with a +separate precomputed gain for each iteration count. Public `Angle` operations +and `Rotation::::new` retain constant 29/Q30 settings. `Rotation::::with_precision` +selects iterations while retaining Q30 coefficients. There is no user-selectable +normalization algorithm. `Angle::sin`, `cos`, `sin_cos`, and `atan2` expose the +fixed-precision operations themselves. + +The `library_fixed` reference calls the public `Angle::atan2` and `sin_cos`. +The `runtime` profiles pass settings through `black_box`, allowing the CSV to +show runtime-loop overhead separately from precision changes. The 29/Q30 runtime +profile is tested to match the public library bit for bit. The experiment's +matrix application is also checked against `Rotation::apply`. + +All vectors retain i32 Q30 storage. Coefficients quantized to Q14..Q30 are lifted +back to Q30 **without recovering lost bits** before multiplication. This gives +every profile the same i64 multiplication and constant division in the point +loop, isolating coefficient accuracy from runtime division and storage changes. +The library also has `Rotation` with Q14 coefficients and `Rotation` +with coefficients lifted exactly from Q30 to Q62. This experiment does not +benchmark i16/Q14 vector storage, SIMD, or +narrow multiply implementation. The computational kernel remains integer-only; +float references and reporting run outside the measured regions. + +## What the CSV measures + +- `setup_median_ns`: exact cross/dot, vectoring, step division, rotation CORDIC, + coefficient conversion. One setup for every input arc, including short arcs. +- `arc_median_ns`, `arc_min_ns`, `arc_max_ns`: setup plus matrix applications + and writing intermediate points into a warmed reusable buffer. Arcs with no + intermediate points skip setup using the precomputed reference segment count. +- `points_arc`: identical work across all configurations; `amortized_ns_point` + includes setup and is not a pure matrix-operation latency. +- `rotation_median_ns`: matrix coefficient construction alone, including adaptive + iteration selection when enabled, with step angles prepared outside timing. +- `api_rotation_median_ns`: public `Rotation::new` or `with_precision` construction, + including the achieved-angle result. Empty for explicit runtime profiles that + do not correspond to a public constructor. Older CSV snapshots lack this column. +- `rotation_iterations_min/max/mean`: iteration counts actually chosen for + rotation across the inputs, including short arcs in the setup measurement. +- `rotation_max_relative_error`: largest absolute matrix angle error divided + by its requested step angle, not accumulated arc error. +- `atan_max_error_bits`: maximum vectoring error, in 32-bit binary angle units. +- `arc_max_error_units`: maximum displacement over all computed arc points, + **including a computed endpoint**, relative to the ideal rotation of the + stored initial vector. It includes vectoring and matrix errors. +- `drift_1024_max_error_units` and `radial_1024_max_loss_units`: accumulated + displacement and radial loss over 1024 rotations, independently of vectoring. +- `max_step_excess_deg`, `nonmonotone_gaps`: observed step-limit excess and + wrong-order/zero angular gaps, including the gap to the exact input endpoint. + +Timing uses 2048 deterministic mixed arcs (0.1..179.9 degrees, both directions). +Input normalization and reference segment counts are outside timing. Counts +are computed by ceiling the accurate reference sweep divided by the requested +maximum step, and are held fixed for every profile. This prevents a low-precision +profile from appearing faster just because it underestimated the sweep and +emitted fewer points. It is a **controlled comparison**, not an end-to-end +measurement of adaptive subdivision. The existing `cordic_arc_bench` example +continues to measure subdivision as part of the build. + +Experimental profiles have **no maximum-step guarantee**: the usual final-gap +reserve is intentionally omitted to keep point counts identical. Violations +are reported, not hidden. The production example and its strict-gap tests still +include their reserve. Accuracy is relative to stored normalized inputs, before +integer grid rounding; initial normalization error is additional. The 1024-step +stress test covers 257 start/step combinations over the supported step range, +both orientations, every intermediate result, with no endpoint replacement. + +Timing: 40 ms warmup, seven batches of approximately 50 ms, median and range. +Errors are deterministic sample maxima, not exhaustive bounds. `--radius` scales +reported errors; it does not change the timed direction-generation workload. + +## Sample results + +The following explicit-profile measurements precede the relative-mode addition. +Mac ARM64, rustc 1.98.1, release defaults, step 1 degree, radius 65536. Each arc +has 88.342 intermediate points on average, for every profile. Ranges below span +the normal-order and reversed-order medians; see the raw CSVs for each run's +minimum and maximum. Values are local measurements, not portable speed promises. + +| Mode | Iterations | Q | Setup ns | Arc ns | Arc error | Error after 1024 | +|---|---:|---:|---:|---:|---:|---:| +| library_fixed | 29 | 30 | 63.5–65.8 | 250.1–252.9 | 0.065 | 0.308 | +| runtime | 29 | 30 | 68.1–71.3 | 255.2–255.3 | 0.065 | 0.308 | +| runtime | 24 | 30 | 50.9–52.8 | 239.7–241.4 | 1.359 | 7.967 | +| runtime | 20 | 30 | 38.7–41.9 | 231.8–232.8 | 20.994 | 127.902 | +| runtime | 15 | 30 | 26.6–26.8 | 224.0–224.3 | 699.105 | 4068.959 | +| runtime | 29 | 24 | 69.1–71.1 | 253.7–256.4 | 0.918 | 5.573 | +| runtime | 29 | 20 | 67.1–71.0 | 255.4–256.7 | 14.737 | 87.067 | +| runtime | 29 | 14 | 68.2–70.7 | 257.3–257.5 | 797.476 | 5208.792 | +| runtime | 24 | 24 | 49.5–52.8 | 238.6–240.7 | 2.054 | 11.169 | +| runtime | 20 | 20 | 38.2–38.3 | 232.7–233.8 | 25.246 | 177.744 | +| runtime | 15 | 14 | 27.0–27.3 | 222.1–224.2 | 1394.052 | 7545.789 | + +Reducing iteration count speeds up setup, but point generation dominates long +arcs. Lowering only coefficient precision barely changes speed in this storage +model and severely increases radial/angular drift. A dynamic 29/Q30 configuration +also costs more setup time than the constant public implementation; that overhead +should not be mistaken for an effect of numerical accuracy. + +Raw runs: [normal order](cordic/precision.csv), +[reverse order](cordic/precision-reverse.csv). + +Tests check public atan2 across full signed integer ranges (including i128::MIN), +public sin/cos consistency, runtime/default equivalence, invalid CLI arguments, +and coefficient non-expansion for every supported iteration/Q combination over +4097 full-turn samples. Existing accuracy and strict consumer-gap tests remain. + +## Relative-mode measurements + +Same Mac ARM64, release, 2048 mixed arcs, radius 65536, maximum step 20 degrees. +Counts are held at 3.934 interior points per arc. Ranges span the two run-order +medians, not different input sets. + +| Mode | Matrix setup ns | Full setup ns | Arc ns | Maximum single-step relative error | +|---|---:|---:|---:|---:| +| Fixed 29/Q30 | 21.67–21.68 | 62.25–65.79 | 55.65–58.65 | 0.000031% | +| Rotation step/16 | 6.37–6.39 | 50.05–50.58 | 44.48–45.11 | 5.880032% | +| Rotation step/8 | 5.61–5.63 | 46.43–52.48 | 39.22–43.04 | 11.811171% | + +At a nominal 20-degree step, the matrix setup alone is about 3–4x faster. +Full-arc speedup is smaller because atan2 stays accurate and output work remains. +Only the rotation iteration count changes; coefficient norm remains conservative +Q30 and radial drift after 1024 applications is still about 0.11 units at R65536. + +The unchanged-reference-count experiment also exposes why subdivision must use +the achieved matrix angle: at max step 20 and divisor 16, the worst final gap +exceeds the requested step by 6.991 degrees. At max step 1 and divisor 16, +1640 of 2048 arcs have a negative final remainder because the repeated rotations +have passed the input endpoint. This is systematic angular drift, not radial +contraction or vectoring error. Those configurations are a speed/accuracy +experiment, not a valid replacement arc builder with the current count policy. + +Raw runs: [step 20, divisor 16](cordic/relative16-20.csv), +[reversed](cordic/relative16-20-reverse.csv), +[step 20, divisor 8](cordic/relative8-20.csv), +[reversed](cordic/relative8-20-reverse.csv), +[step 1, divisor 16](cordic/relative16-1.csv). + +Additional tests verify the minimal sufficient iteration count, clockwise/CCW +symmetry, and the single-step angular budget on 65,536 angles per divisor and +five arc-step workloads. Atan2 is checked to retain its original accuracy. + +## Public constructor measurement + +After adding `with_precision` and the achieved-angle field, the same step-20 +workload measured public construction at 21.75–21.81 ns for `Rotation::new`, +7.45 ns for exponent 4, and 6.62 ns for exponent 3. These are local medians, +including the achieved-angle result; the old coefficient-only column remains +separate. The full-arc timing still uses fixed reference counts and is not a +validated approximate subdivision implementation. + +Raw runs: [exponent 4](cordic/rotation-api16-20.csv), +[exponent 3](cordic/rotation-api8-20.csv). diff --git a/docs/cordic.md b/docs/cordic.md new file mode 100644 index 0000000..8a9b54f --- /dev/null +++ b/docs/cordic.md @@ -0,0 +1,278 @@ +# Integer CORDIC for arc consumers + +The public implementation is `i_float::int::angle::{Angle, Rotation}`. It has +no allocations, float arithmetic, runtime trigonometry, physics dependencies, +serde support, coordinate/radius types, or arc topology. The crate remains +`no_std`. + +```rust +impl Angle { + pub const MAX_ERROR: u32 = 32; + pub const SIN_COS_SCALE: i32 = 1 << 30; + pub const fn from_bits(bits: u32) -> Self; + pub const fn bits(self) -> u32; + pub fn between(from: UnitIntVector, to: UnitIntVector) -> Self; + pub fn atan2(y: W, x: W) -> Option; + pub fn sin(self) -> i32; + pub fn cos(self) -> i32; + pub fn sin_cos(self) -> (i32, i32); +} +impl Rotation { + pub fn new(angle: Angle) -> Self; + pub fn with_precision(angle: Angle, precision: u32) -> Self; + pub const ANGLE_MAX_ERROR: u32; // 131072 for i16; 2 for i32/i64 + pub const MAX_ERROR: u32; // ANGLE_MAX_ERROR + 3 + pub const fn angle(self) -> Angle; + pub fn apply(self, vector: UnitIntVector) -> UnitIntVector; +} +``` + +`between(a,b)` measures counterclockwise sweep. For clockwise magnitude, use +`between(b,a)` and negate the eventual step with `wrapping_neg()`. Coincident +rays give zero even if their lengths differ. Opposite rays give exactly half a +turn. Axes are exact. Almost coincident rays retain the sign of the exact cross +product, so a tiny major arc cannot accidentally become an empty arc. There is +no dedicated full-circle API. + +`src/int/angle/angle.rs` contains an executable documentation example of the future +iOverlay loop. It clears and reuses a consumer-owned `Vec`, computes the sweep +once, divides it into segments, constructs one matrix, and applies that matrix +only to intermediate directions. It excludes both input endpoints. The +consumer retains the exact original contact points, center, radius, and mesh +topology. Replacing the endpoint does not correct the measured interior errors. + +The example uses a binary-angle step clamped to `[2^22, 2^29]`, exactly +`360/1024` through `45` degrees. These limits do **not** quantize the input rays. +The existing iOverlay `ArcStep` currently stores a squared chord. Integrating +this example will require adapting that representation or its conversion on +the iOverlay side; the current squared-chord constructor is not directly +accepted by this new API. No iOverlay files were modified. + +The module is split into `angle.rs` (public angle and trigonometry), +`rotation.rs` (matrix storage/application), and private `cordic.rs` (shared +integer kernels and tables), re-exported by `mod.rs`. `atan2(y,x)` accepts the +full built-in wide integer range, including MIN; only `(0,0)` returns `None`. +`sin_cos` returns `(sin,cos)` in Q30, with `SIN_COS_SCALE` representing one. +Separate `sin`/`cos` calls each run the kernel; use `sin_cos` when both are needed. +`Rotation::new` reuses the full 29-iteration kernel. For approximate steps, +`Rotation::::with_precision(angle, 3)` allows angle/8 error; exponent 4 allows +angle/16. The magnitude is the shortest signed angle, so clockwise steps have +the same budget. `Rotation` stores coefficients in `I`, with the same scale +as `UnitIntVector`: Q14/Q30/Q62. The precision parameter changes iteration +count, including a reserve for coefficient quantization in the selected type. +Every `u32` exponent is accepted, with >=32 selecting full precision. The +matrix-angle error bound in binary units is `max(floor(|angle| / 2^precision), Rotation::::MAX_ERROR)`. +Zero and cardinal rotations are exact. + +`rotation.angle()` returns the achieved angle from the CORDIC residual, without +another atan2. Its error relative to the matrix angle is at most +`Rotation::::ANGLE_MAX_ERROR` (131072 binary units for i16, 2 for i32/i64). Use this achieved step for counts +and remainders; tiny requests may report zero. A strict maximum-step policy +must also reserve room for rounding and bias the requested step inward. +The test `achieved_angle_drives_approximate_arc_counts` demonstrates those +margins for fresh i32 directions, both orientations, and exponents 3 and 4. + +For runtime precision experiments see [the precision benchmark](cordic-precision.md). + +## Numerical choices and contracts + +- Public angle: unsigned 32-bit binary turn, about `1.463e-9` radians per unit. + A 10-bit angle would lose direction information unrelated to the requested + subdivision resolution. Even Q24 angle rounding alone could accumulate about + 12.6 coordinate units at radius 65536 over 1024 steps. Q32 keeps this term + below 0.05 units for nearest rounding (below 0.1 for step truncation). +- `cordic::ITERATIONS` sets the full-precision count for both kernels and table + lengths. Its value 29 is an accuracy choice: after n iterations the residual + is below `2^-(n-1)` radians. At radius 65536 over 1024 applications, 29 + iterations contribute less than 0.25 coordinate units from the residual + alone; 28 would give a bound of 0.5. Coefficient truncation and coordinate + rounding add separate errors. This is not maximum attainable precision or + a count derived from the coordinate type; increasing it requires extending + the tables and revisiting the error contract. +- The atan table has 16 additional guard bits (48-bit binary turns internally). + This avoids accumulating the rounding errors of a table rounded directly to + Q32. It fits in `i64`, as do the internal rotation coordinates. +- Matrix coefficients are stored in `I`: Q14 for i16, Q30 for i32, Q62 for i64. + The kernel produces at most Q30 precision. i16 truncates toward zero to Q14; + i64 lifts Q30 to Q62 exactly, without gaining bits of accuracy. Coefficient + truncation contributes at most `sqrt(2)/2^min(I::BITS-2,30)` per application. +- Gain compensation: Q60 internally, with a 128-unit inward guard. The integer + shift error is less than `29 * sqrt(2) * 1.647 < 68` Q60 units, smaller than the + propagated gain guard. Truncating the resulting coordinates to Q14 or Q30 cannot + increase norm. In contrast, copying iPhysics's 128-unit **Q30** guard would + introduce roughly 13 coordinate units of contraction after 1024 steps at + radius 65536. +- Each matrix application truncates both output components **toward zero**. + Thus neither the matrix nor storage rounding can increase the input length. + There is no per-point CORDIC, square root, or normalization. All applications + use `I::Wide`: i32 for i16, i64 for i32, i128 for i64. With S = 2^(I::BITS-2), + each product is at most S² and the sum at most 2S², within the signed wide + range. No coordinate-type branch or conversion through usize is needed. +- Cross/dot products are computed at full input precision. With + `S=2^(I::BITS-2)`, each component product is at most `S²` and sums/differences + at most `2S²`; even the conservative i64 bound `2^125` fits signed i128. + Cross/dot are never squared. Only after preserving exact signs and axis cases + are their magnitudes reduced to 30 significant bits, then lifted for i64 + vectoring. The conversion also fits a 32-bit host's `usize`. + +`Angle::MAX_ERROR` bounds vectoring relative to the **stored** input directions, +not to an original vector before `fast_normalize`. All error bounds are for the +built-in i16/i32/i64 implementations of `IntNumber`. + +`Angle` remains a 32-bit binary turn independently of `I`. Coordinate storage, +angle representation, and the requested approximation budget are separate +choices. Making only `Angle` wider would not recover precision beyond the Q30 kernel; +making it narrower would introduce angle quantization even when the consumer +requests high accuracy. `with_precision` controls construction cost without +changing these representations. The 30 retained vectoring input bits and their +guard-bit lift likewise have separate constants, unrelated to iteration count. + +The supplied consumer example uses integer ceiling, with an upper sweep bound +`upper = sweep + Angle::MAX_ERROR`. A simple ceiling does not by itself account +for accumulated matrix rounding in the final gap. For fresh i32/i64 normalized +inputs and this step range, it increases `n` until +`upper + 8*n*n <= n*max_step`. Eight binary units cover step division and the +per-application angular error, which is below five units when the vector norm +is at least 0.99. Over at most 1030 applications, freshly normalized i32/i64 +vectors remain above that norm. This reserves room for the final exact endpoint +as well as intermediate gaps. Tests reached 1027 segments on almost full turns. +The margin is deliberately part of the consumer example, not an arc policy in +`Angle` or `Rotation`. It does not cover i16, heavily contracted input vectors, +or angular perturbations introduced when a consumer rounds scaled mesh points. + +## Accuracy measurements + +Run `cargo test --lib int::angle -- --nocapture`. The tests inspect the whole +sequence, including a computed endpoint, without replacing it with the input. + +Vectoring: 20,000 deterministic pairs per coordinate type, independent input +lengths, every quadrant, axes, same/opposite rays, and near-degenerate cases. +Maximum observed errors were 3.023 / 3.961 / 3.970 binary units for i16/i32/i64, +below the public conservative bound of 32 units. Extreme i64 axis and near-axis +cases exercise full-width i128 products. + +Rotation: 262,193 sampled angles across the full turn. Maximum observed matrix +coefficient error was `4.652e-9`, maximum norm loss `1.3134e-9`. Cardinal +rotations are exact, including after 1024 applications. + +Accumulation: 257 start directions and step values spanning the allowed range, +both orientations, every point up to 1024 applications. The large-step cases +intentionally run through multiple turns as a stress test. These are observed +maxima over a deterministic sample, not a claim of exhaustive coverage. + +| Storage | Added displacement, R=1024 | Added displacement, R=65536 | Radial loss, R=65536 | Angular error, radians | +|---|---:|---:|---:|---:| +| i32 | 0.004811 | 0.307934 | 0.125796 | 4.531e-6 | +| i64 | 0.004697 | 0.300580 | 0.085662 | 4.542e-6 | + +Relative to an ideal unit direction with the same initial angle, including +initial normalization contraction, observed maxima at R=65536 were 0.621145 +for i32 and 0.300580 for i64. After `scale(65536)` rounds to the coordinate grid, +added displacement relative to rotation of the stored initial vector reached +0.950276 / 0.925431, respectively. Errors before grid rounding scale linearly +with radius; these endpoints cover radii 2^10 through 2^16. Initial normalization +errors depend on input magnitude and direction and can exceed this sample's +values; the existing approximate normalization contract has not been tightened. + +With Q14 coefficients in `Rotation`, added displacement at R=1024 reached +**120.724** units, with angular error 0.060781 radians after 1024 applications. +Both coefficient quantization and per-step Q14 storage rounding contribute. +The earlier Q30-coefficient implementation measured 42.436 units; that result +no longer describes i16 rotation. The regression bound now sums the documented +coefficient and component-rounding errors over all applications. +This API supports i16 arithmetic but does not make fine, long i16 arcs accurate. +It preserves the length upper bound and documents accumulation explicitly on +`UnitIntVector`. No selectable normalization policy was introduced. + +A separate consumer test checks 64 start angles, both traversal directions, +three maximum steps, minor/major arcs, half turns, and near coincident/opposite +rays. With the final-gap reserve above, maximum added displacement at R=65536 +was 0.184918 and no angular gap exceeded the allowed maximum (test tolerance +`1e-12` radians). Before adding the reserve, a simple ceiling allowed a final-gap +excess of `3.700e-6` radians. That issue is covered by the regression test. + +## Reproducible performance comparison + +`examples/cordic_arc_bench.rs` adapts the provided temporary benchmark. Its +`examples/support/bisection.rs` is a benchmark-only snapshot of the current, +uncommitted iOverlay arc implementation, with imports adjusted and unused +API/docs removed. It is not another public arc implementation in iFloat. +Float uses the real `FloatNumber` methods (`libm::acos` and `libm::sincos`) and +ceiling for segment count; the old floor baseline is not used in the results. + +```sh +cargo run --release --example cordic_arc_bench > results.csv +cargo run --release --example cordic_arc_bench -- reverse > results-reverse.csv +``` + +Mac ARM64, rustc 1.98.1, release defaults, libm 0.2.16, 2026-09-13. Each group +contains the same 2048 deterministic arcs as the supplied benchmark: mixed +sweeps 0.1–179.9 degrees or short sweeps 0.1–10 degrees, both orientations. +Input normalization is outside timing. Each algorithm uses warmed, retained +buffers and emits intermediate directions. A 60 ms warmup precedes seven +batches of about 70 ms each; results are medians in ns/arc. Algorithms were +also measured in reversed order. Timing and counts include vectoring, matrix +setup, subdivision guards, and all emitted points. Center/radius/grid placement +is outside the timing for every implementation. These are local microbenchmarks, +not integrated mesh timings; major arcs are accuracy-tested, not benchmarked. + +Ranges below span the normal-order and reverse-order medians. Raw CSV files +include minima, maxima, f32 results, and short-arc groups. + +| Max step | Bisection i32 | CORDIC i32 | Bisection i64 | CORDIC i64 | f64 ceil | +|---|---:|---:|---:|---:|---:| +| 45° | 32.5–32.6 | 60.5–68.5 | 84.8–85.0 | 70.0–71.5 | 8.5–8.6 | +| 15° | 148.8–149.8 | 73.8–75.5 | 378.7–380.5 | 88.3–89.6 | 13.9–14.0 | +| 5° | 495.3–499.9 | 97.3–97.5 | 1292.6–1301.3 | 133.2–135.6 | 28.5–29.1 | +| 1° | 2428.5–2432.2 | 258.3–258.7 | 6456.8–6594.7 | 397.9–405.8 | 130.0–131.6 | +| 0.351562° | 5258.2–5332.4 | 624.9–635.3 | 15025.9–15208.0 | 1020.8–1029.9 | 381.3–389.2 | + +Counts must be considered alongside timing: + +| Max step | Bisection points/arc | CORDIC points/arc | f64 ceil points/arc | +|---|---:|---:|---:| +| 5° | 26.191 | 17.256 | 17.256 | +| 1° | 129.973 | 88.350 | 88.342 | +| 0.351562° | 316.454 | 252.383 | 252.225 | + +At 5 degrees CORDIC i32 is about 5x faster than bisection; at 1 degree about +9x faster. This improvement combines removal of per-point normalization with +fewer output points. CORDIC remains about 3.4x / 2x slower than f64 at those +steps. The conservative integer final-gap reserve accounts for the small extra +point count compared with f64 at fine steps. i64 CORDIC is about 1.4–1.6x slower +than i32 on these workloads, rather than the roughly 2.6–2.9x bisection gap. + +CORDIC loses on short/coarsely subdivided i32 arcs: the 45-degree mixed case +is roughly twice as slow as bisection; short arcs needing no intermediate points +still pay for vectoring. The float baseline can reject those with one dot +comparison. No early-rejection optimization or additional public API was added +just to improve that microbenchmark. + +Raw results: [normal order](cordic/results.csv), +[reverse order](cordic/results-reverse.csv). + +## Source and integration notes + +The worktree started at `1ff6354`, before `UnitIntVector` existed. The necessary +existing changes from the main iFloat checkout were brought forward without +commits: unsigned squared lengths (`1471fc2`, `5eada81`) and unit normalization +(`a3a24c5`). `IntPoint` squared-length/distance return types therefore also +reflect that already existing unsigned API. Unrelated later range-check and +adapter changes were not copied. The only new change to unit-vector behavior +is the explicitly documented rotation operation and its private invariant- +preserving component constructor. + +The source iPhysics implementation was verified in +`src/quantity/angle.rs`: Angle(u32), 31 iterations, Q30 output, guarded gain, +exact axes, and no vectoring/atan2. Only the integer rotation approach was +adapted. Its angular velocity, physical units, float conversions, and operators +were not transferred. Vectoring is new here. + +No commits, merges, or modifications of the source iPhysics/iOverlay checkouts +were performed. Integration into `outline`, `stroke`, and `variable_stroke`, +and changing iOverlay's `ArcStep` representation, remain separate work. + +Initial adaptation validation completed locally, offline: 96 tests in both debug and release with +all features, four doctests, formatting, Clippy with warnings denied, rustdoc +with warnings denied, all five CI feature configurations, Rust 1.85 MSRV, +wasm32-unknown-unknown, and `cargo package --allow-dirty` including verification. diff --git a/docs/cordic/precision-reverse.csv b/docs/cordic/precision-reverse.csv new file mode 100644 index 0000000..87e0a76 --- /dev/null +++ b/docs/cordic/precision-reverse.csv @@ -0,0 +1,12 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps +runtime,15,14,1.00000000,65536,26.993,222.128,220.249,225.876,2.514,88.342,41721.387676,1394.052284,7545.789141,5281.007402,1.162284220,0 +runtime,20,20,1.00000000,65536,38.184,232.685,228.834,233.377,2.634,88.342,1303.387676,25.245656,177.744099,84.689733,0.019341539,0 +runtime,24,24,1.00000000,65536,52.820,238.587,235.735,245.730,2.701,88.342,81.676073,2.053902,11.169138,5.339835,0.001157143,0 +runtime,29,14,1.00000000,65536,70.700,257.257,252.881,259.090,2.912,88.342,3.364758,797.476136,5208.791695,5045.142196,0.606978256,0 +runtime,29,20,1.00000000,65536,70.987,255.363,248.990,267.226,2.891,88.342,3.364758,14.737282,87.066685,85.405892,0.008799770,0 +runtime,29,24,1.00000000,65536,69.134,253.657,250.419,262.212,2.871,88.342,3.364758,0.917627,5.572877,5.214302,0.000564380,0 +runtime,15,30,1.00000000,65536,26.577,224.345,220.960,227.498,2.539,88.342,41721.387676,699.104760,4068.959059,0.124800,0.604112572,0 +runtime,20,30,1.00000000,65536,41.922,232.820,230.518,236.782,2.635,88.342,1303.387676,20.993921,127.901843,0.121945,0.016658852,0 +runtime,24,30,1.00000000,65536,52.846,239.659,238.845,245.253,2.713,88.342,81.676073,1.359282,7.966947,0.120423,0.000817472,0 +runtime,29,30,1.00000000,65536,68.068,255.196,254.568,256.329,2.889,88.342,3.364758,0.064596,0.307934,0.125796,0.000042079,0 +library_fixed,29,30,1.00000000,65536,65.844,252.871,247.594,258.935,2.862,88.342,3.364758,0.064596,0.307934,0.125796,0.000042079,0 diff --git a/docs/cordic/precision.csv b/docs/cordic/precision.csv new file mode 100644 index 0000000..24c48e5 --- /dev/null +++ b/docs/cordic/precision.csv @@ -0,0 +1,12 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps +library_fixed,29,30,1.00000000,65536,63.511,250.109,245.746,253.999,2.831,88.342,3.364758,0.064596,0.307934,0.125796,0.000042079,0 +runtime,29,30,1.00000000,65536,71.343,255.326,248.961,261.860,2.890,88.342,3.364758,0.064596,0.307934,0.125796,0.000042079,0 +runtime,24,30,1.00000000,65536,50.867,241.352,236.506,243.431,2.732,88.342,81.676073,1.359282,7.966947,0.120423,0.000817472,0 +runtime,20,30,1.00000000,65536,38.749,231.755,228.533,235.994,2.623,88.342,1303.387676,20.993921,127.901843,0.121945,0.016658852,0 +runtime,15,30,1.00000000,65536,26.805,224.038,219.899,227.915,2.536,88.342,41721.387676,699.104760,4068.959059,0.124800,0.604112572,0 +runtime,29,24,1.00000000,65536,71.099,256.440,253.265,264.243,2.903,88.342,3.364758,0.917627,5.572877,5.214302,0.000564380,0 +runtime,29,20,1.00000000,65536,67.051,256.658,253.016,377.983,2.905,88.342,3.364758,14.737282,87.066685,85.405892,0.008799770,0 +runtime,29,14,1.00000000,65536,68.195,257.492,249.132,265.256,2.915,88.342,3.364758,797.476136,5208.791695,5045.142196,0.606978256,0 +runtime,24,24,1.00000000,65536,49.549,240.658,238.428,244.889,2.724,88.342,81.676073,2.053902,11.169138,5.339835,0.001157143,0 +runtime,20,20,1.00000000,65536,38.348,233.811,232.822,238.085,2.647,88.342,1303.387676,25.245656,177.744099,84.689733,0.019341539,0 +runtime,15,14,1.00000000,65536,27.272,224.213,222.822,228.883,2.538,88.342,41721.387676,1394.052284,7545.789141,5281.007402,1.162284220,0 diff --git a/docs/cordic/relative16-1.csv b/docs/cordic/relative16-1.csv new file mode 100644 index 0000000..9f27d3a --- /dev/null +++ b/docs/cordic/relative16-1.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns +library_fixed,29,30,1.00000000,65536,73.916,261.055,256.384,262.270,2.955,88.342,3.364758,0.064596,0.307934,0.125796,0.000042079,0,0,29,29,29.000,0.000000315,21.844 +relative,29,30,1.00000000,65536,50.720,249.805,241.228,263.548,2.828,88.342,3.364758,4663.636074,131071.109761,0.105923,1.214668419,1640,16,11,14,11.030,0.058185043,9.838 diff --git a/docs/cordic/relative16-20-reverse.csv b/docs/cordic/relative16-20-reverse.csv new file mode 100644 index 0000000..7544fa2 --- /dev/null +++ b/docs/cordic/relative16-20-reverse.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns +relative,29,30,20.00000000,65536,50.050,45.109,40.803,53.648,11.468,3.934,3.364758,9175.182281,131071.109761,0.105923,6.990962194,0,16,7,14,7.235,0.058800323,6.387 +library_fixed,29,30,20.00000000,65536,62.251,55.647,53.039,64.945,14.146,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.672 diff --git a/docs/cordic/relative16-20.csv b/docs/cordic/relative16-20.csv new file mode 100644 index 0000000..29fc4ca --- /dev/null +++ b/docs/cordic/relative16-20.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns +library_fixed,29,30,20.00000000,65536,65.786,58.649,57.423,61.481,14.910,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.685 +relative,29,30,20.00000000,65536,50.576,44.475,38.270,51.442,11.307,3.934,3.364758,9175.182281,131071.109761,0.105923,6.990962194,0,16,7,14,7.235,0.058800323,6.373 diff --git a/docs/cordic/relative8-20-reverse.csv b/docs/cordic/relative8-20-reverse.csv new file mode 100644 index 0000000..8008d22 --- /dev/null +++ b/docs/cordic/relative8-20-reverse.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns +relative,29,30,20.00000000,65536,52.475,39.221,37.200,46.482,9.971,3.934,3.364758,15839.574841,131071.140520,0.106601,10.771968096,0,8,6,13,6.235,0.118111713,5.611 +library_fixed,29,30,20.00000000,65536,73.560,56.695,52.311,63.582,14.413,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.575 diff --git a/docs/cordic/relative8-20.csv b/docs/cordic/relative8-20.csv new file mode 100644 index 0000000..b0acea2 --- /dev/null +++ b/docs/cordic/relative8-20.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns +library_fixed,29,30,20.00000000,65536,64.692,60.358,55.498,64.499,15.344,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.934 +relative,29,30,20.00000000,65536,46.434,43.038,36.221,49.635,10.941,3.934,3.364758,15839.574841,131071.140520,0.106601,10.771968096,0,8,6,13,6.235,0.118111713,5.635 diff --git a/docs/cordic/results-reverse.csv b/docs/cordic/results-reverse.csv new file mode 100644 index 0000000..e851211 --- /dev/null +++ b/docs/cordic/results-reverse.csv @@ -0,0 +1,43 @@ +implementation,group,step_deg,median_ns_arc,min_ns_arc,max_ns_arc,points_arc +cordic_i64,mixed,45.00000000,71.453,69.837,85.124,1.476 +cordic_i32,mixed,45.00000000,60.463,54.470,70.953,1.476 +bisection_i64,mixed,45.00000000,84.836,83.851,85.842,1.725 +bisection_i32,mixed,45.00000000,32.571,31.753,33.108,1.725 +f64_ceil,mixed,45.00000000,8.463,8.381,8.627,1.476 +f32_ceil,mixed,45.00000000,6.313,6.264,6.381,1.476 +cordic_i64,mixed,15.00000000,89.600,85.702,102.941,5.424 +cordic_i32,mixed,15.00000000,75.483,68.361,81.817,5.424 +bisection_i64,mixed,15.00000000,378.680,373.670,381.611,7.642 +bisection_i32,mixed,15.00000000,148.803,147.420,150.558,7.642 +f64_ceil,mixed,15.00000000,13.996,13.767,14.102,5.424 +f32_ceil,mixed,15.00000000,12.020,11.857,12.076,5.424 +cordic_i64,mixed,5.00000000,135.581,130.453,153.080,17.256 +cordic_i32,mixed,5.00000000,97.309,93.952,103.188,17.256 +bisection_i64,mixed,5.00000000,1292.559,1278.702,1306.091,26.191 +bisection_i32,mixed,5.00000000,499.878,493.558,508.462,26.191 +f64_ceil,mixed,5.00000000,28.539,28.306,29.382,17.256 +f32_ceil,mixed,5.00000000,28.376,28.198,28.685,17.256 +cordic_i64,mixed,1.00000000,405.799,399.567,414.912,88.350 +cordic_i32,mixed,1.00000000,258.704,255.918,266.822,88.350 +bisection_i64,mixed,1.00000000,6456.816,6428.670,6511.353,129.973 +bisection_i32,mixed,1.00000000,2432.173,2409.604,2453.540,129.973 +f64_ceil,mixed,1.00000000,131.581,130.840,132.266,88.342 +f32_ceil,mixed,1.00000000,139.206,138.886,141.913,88.342 +cordic_i64,mixed,0.35156250,1029.875,1011.050,1037.369,252.383 +cordic_i32,mixed,0.35156250,635.333,621.489,642.739,252.383 +bisection_i64,mixed,0.35156250,15025.869,14801.676,15246.745,316.454 +bisection_i32,mixed,0.35156250,5258.236,5227.587,5342.167,316.454 +f64_ceil,mixed,0.35156250,381.285,376.908,394.934,252.225 +f32_ceil,mixed,0.35156250,392.241,391.842,395.054,252.225 +cordic_i64,short,45.00000000,36.111,33.760,51.006,0.000 +cordic_i32,short,45.00000000,35.011,30.560,42.867,0.000 +bisection_i64,short,45.00000000,5.398,5.397,5.416,0.000 +bisection_i32,short,45.00000000,4.994,4.986,5.003,0.000 +f64_ceil,short,45.00000000,0.519,0.518,0.520,0.000 +f32_ceil,short,45.00000000,0.540,0.538,0.551,0.000 +cordic_i64,short,5.00000000,54.889,48.893,64.514,0.498 +cordic_i32,short,5.00000000,55.645,47.241,68.083,0.498 +bisection_i64,short,5.00000000,29.957,29.642,30.231,0.498 +bisection_i32,short,5.00000000,13.729,13.664,13.836,0.498 +f64_ceil,short,5.00000000,5.685,5.657,5.732,0.498 +f32_ceil,short,5.00000000,4.274,4.219,4.297,0.498 diff --git a/docs/cordic/results.csv b/docs/cordic/results.csv new file mode 100644 index 0000000..7824aff --- /dev/null +++ b/docs/cordic/results.csv @@ -0,0 +1,43 @@ +implementation,group,step_deg,median_ns_arc,min_ns_arc,max_ns_arc,points_arc +f64_ceil,mixed,45.00000000,8.584,8.502,9.218,1.476 +f32_ceil,mixed,45.00000000,6.405,6.326,6.427,1.476 +bisection_i32,mixed,45.00000000,32.498,31.944,33.022,1.725 +bisection_i64,mixed,45.00000000,85.006,83.992,86.792,1.725 +cordic_i32,mixed,45.00000000,68.527,63.277,71.555,1.476 +cordic_i64,mixed,45.00000000,70.013,63.532,72.607,1.476 +f64_ceil,mixed,15.00000000,13.912,13.884,14.103,5.424 +f32_ceil,mixed,15.00000000,12.002,11.879,12.046,5.424 +bisection_i32,mixed,15.00000000,149.825,146.882,151.636,7.642 +bisection_i64,mixed,15.00000000,380.477,376.551,384.663,7.642 +cordic_i32,mixed,15.00000000,73.798,68.374,79.789,5.424 +cordic_i64,mixed,15.00000000,88.265,85.105,99.325,5.424 +f64_ceil,mixed,5.00000000,29.082,28.412,29.662,17.256 +f32_ceil,mixed,5.00000000,28.155,27.787,28.272,17.256 +bisection_i32,mixed,5.00000000,495.320,493.102,500.662,26.191 +bisection_i64,mixed,5.00000000,1301.338,1284.147,1307.833,26.191 +cordic_i32,mixed,5.00000000,97.497,94.887,108.585,17.256 +cordic_i64,mixed,5.00000000,133.247,129.309,140.947,17.256 +f64_ceil,mixed,1.00000000,130.041,129.783,131.390,88.342 +f32_ceil,mixed,1.00000000,138.672,138.258,139.180,88.342 +bisection_i32,mixed,1.00000000,2428.500,2405.158,2463.188,129.973 +bisection_i64,mixed,1.00000000,6594.698,6428.528,6643.624,129.973 +cordic_i32,mixed,1.00000000,258.281,253.074,265.723,88.350 +cordic_i64,mixed,1.00000000,397.890,395.291,408.068,88.350 +f64_ceil,mixed,0.35156250,389.199,378.943,395.330,252.225 +f32_ceil,mixed,0.35156250,393.231,392.274,405.471,252.225 +bisection_i32,mixed,0.35156250,5332.357,5266.110,7247.236,316.454 +bisection_i64,mixed,0.35156250,15208.028,15054.850,15374.044,316.454 +cordic_i32,mixed,0.35156250,624.898,614.672,629.674,252.383 +cordic_i64,mixed,0.35156250,1020.811,1010.478,1083.661,252.383 +f64_ceil,short,45.00000000,0.523,0.523,0.526,0.000 +f32_ceil,short,45.00000000,0.542,0.539,0.546,0.000 +bisection_i32,short,45.00000000,5.005,4.990,5.014,0.000 +bisection_i64,short,45.00000000,5.426,5.412,5.461,0.000 +cordic_i32,short,45.00000000,38.282,33.364,49.603,0.000 +cordic_i64,short,45.00000000,41.283,38.797,45.158,0.000 +f64_ceil,short,5.00000000,5.655,5.597,5.675,0.498 +f32_ceil,short,5.00000000,4.235,4.213,4.258,0.498 +bisection_i32,short,5.00000000,13.525,13.481,13.581,0.498 +bisection_i64,short,5.00000000,29.961,29.882,30.023,0.498 +cordic_i32,short,5.00000000,48.197,44.410,59.849,0.498 +cordic_i64,short,5.00000000,52.790,51.330,56.220,0.498 diff --git a/docs/cordic/rotation-api16-20.csv b/docs/cordic/rotation-api16-20.csv new file mode 100644 index 0000000..1df7ac2 --- /dev/null +++ b/docs/cordic/rotation-api16-20.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns,api_rotation_median_ns +library_fixed,29,30,20.00000000,65536,62.958,57.448,52.524,61.865,14.605,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.643,21.751 +relative,29,30,20.00000000,65536,45.162,40.989,36.154,50.591,10.420,3.934,3.364758,9175.182281,131071.109761,0.105923,6.990962194,0,16,7,14,7.235,0.058800323,6.444,7.453 diff --git a/docs/cordic/rotation-api8-20.csv b/docs/cordic/rotation-api8-20.csv new file mode 100644 index 0000000..f86a781 --- /dev/null +++ b/docs/cordic/rotation-api8-20.csv @@ -0,0 +1,3 @@ +mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns,api_rotation_median_ns +library_fixed,29,30,20.00000000,65536,62.539,58.471,51.463,66.889,14.864,3.934,3.364758,0.003103,0.307934,0.125796,0.000000000,0,0,29,29,29.000,0.000000315,21.715,21.811 +relative,29,30,20.00000000,65536,39.370,38.532,32.931,46.730,9.796,3.934,3.364758,15839.574841,131071.140520,0.106601,10.771968096,0,8,6,13,6.235,0.118111713,5.637,6.615 diff --git a/examples/cordic_arc_bench.rs b/examples/cordic_arc_bench.rs new file mode 100644 index 0000000..1ce6c8c --- /dev/null +++ b/examples/cordic_arc_bench.rs @@ -0,0 +1,252 @@ +// Adapted from /private/tmp/ioverlay-arc-bench.JE1wuF/src/main.rs. +// Run: cargo run --release --example cordic_arc_bench [-- reverse] + +use i_float::float::number::FloatNumber; +use i_float::int::number::{int::IntNumber, wide_int::WideIntNumber}; +use i_float::int::unit_vector::UnitIntVector; +use i_float::int::vector::IntVector; +#[path = "support/bisection.rs"] +mod bisection; +use bisection::{ArcBuilder, ArcDirection, ArcStep}; +use i_float::int::angle::{Angle, Rotation}; +use std::{ + hint::black_box, + time::{Duration, Instant}, +}; +const COUNT: usize = 2048; +#[derive(Clone, Copy)] +struct Case { + a: [f64; 2], + b: [f64; 2], + cw: bool, +} +fn cases(short: bool) -> Vec { + let mut rng = 0xd718af03_u64; + let mut next = || { + rng ^= rng << 13; + rng ^= rng >> 7; + rng ^= rng << 17; + (rng >> 11) as f64 / ((1u64 << 53) as f64) + }; + (0..COUNT) + .map(|_| { + let start = next() * std::f64::consts::TAU; + let sweep = (0.1 + next() * if short { 9.9 } else { 179.8 }).to_radians(); + let cw = next() < 0.5; + let end = start + if cw { -sweep } else { sweep }; + Case { + a: [start.cos(), start.sin()], + b: [end.cos(), end.sin()], + cw, + } + }) + .collect() +} +fn measure(mut f: impl FnMut() -> usize) -> (f64, f64, f64, f64) { + let points = f(); + let t = Instant::now(); + let mut loops = 0usize; + while t.elapsed() < Duration::from_millis(60) { + black_box(f()); + loops += 1; + } + let repeats = ((70_000_000.0 / (t.elapsed().as_nanos() as f64 / loops as f64)) as usize).max(1); + let mut times = Vec::new(); + for _ in 0..7 { + let t = Instant::now(); + for _ in 0..repeats { + black_box(f()); + } + times.push(t.elapsed().as_nanos() as f64 / (repeats * COUNT) as f64); + } + times.sort_by(f64::total_cmp); + (times[3], times[0], times[6], points as f64 / COUNT as f64) +} +fn report(name: &str, group: &str, deg: f64, result: (f64, f64, f64, f64)) { + println!( + "{name},{group},{deg:.8},{:.3},{:.3},{:.3},{:.3}", + result.0, result.1, result.2, result.3 + ); +} +fn bench_int(data: &[Case], group: &str, deg: f64, name: &str) { + let normalize = |p: [f64; 2]| { + IntVector::::new( + I::Wide::from_rounded_float(p[0] * 1_000_000.0), + I::Wide::from_rounded_float(p[1] * 1_000_000.0), + ) + .fast_normalize() + .unwrap() + }; + let input: Vec<_> = data + .iter() + .map(|c| { + ( + normalize(c.a), + normalize(c.b), + if c.cw { + ArcDirection::Clockwise + } else { + ArcDirection::Counterclockwise + }, + ) + }) + .collect(); + let s = UnitIntVector::::DENOMINATOR.to_f64(); + let chord = 4.0 * (deg.to_radians() / 2.0).sin().powi(2); + let step = ArcStep::from_squared_chord(I::Wide::from_rounded_float(chord * s * s).to_uint()); + let mut builder = ArcBuilder::::new(); + report( + name, + group, + deg, + measure(|| { + let mut count = 0; + for &(a, b, d) in black_box(&input) { + let points = builder.build(a, b, d, black_box(step)); + count += points.len(); + black_box(points); + } + count + }), + ); +} +// RoundJoinBuilder's angle/count/rotation kernel, with the same FloatNumber +// libm calls. Both variants emit intermediate unit directions into a reused +// buffer, matching ArcBuilder's output contract. ceil is the comparison with +// a true maximum step; floor reproduces the existing float implementation. +fn bench_float(data: &[Case], group: &str, deg: f64, ceil: bool, name: &str) { + let input: Vec<_> = data + .iter() + .map(|c| { + ( + [F::from_float(c.a[0]), F::from_float(c.a[1])], + [F::from_float(c.b[0]), F::from_float(c.b[1])], + c.cw, + ) + }) + .collect(); + let step = F::from_float(deg.to_radians()); + let inv = F::ONE / step; + let limit = step.cos(); + let mut output = Vec::<[F; 2]>::with_capacity(1024); + report( + name, + group, + deg, + measure(|| { + let mut count = 0; + for &(a, b, cw) in black_box(&input) { + output.clear(); + let dot = a[0] * b[0] + a[1] * b[1]; + if black_box(limit) < dot { + black_box(&output); + continue; + } + let angle = dot.max(-F::ONE).min(F::ONE).acos(); + let ratio = angle * black_box(inv); + let mut n = ratio.to_usize(); + if ceil && F::from_usize(n) < ratio { + n += 1; + } + n = n.max(1); + let delta = angle / F::from_usize(n); + let (sn, cs) = (if cw { -delta } else { delta }).sin_cos(); + let mut v = a; + for _ in 1..n { + v = [cs * v[0] - sn * v[1], sn * v[0] + cs * v[1]]; + output.push(v); + } + count += output.len(); + black_box(&output); + } + count + }), + ); +} +fn bench_cordic(data: &[Case], group: &str, deg: f64, name: &str) { + let normalize = |p: [f64; 2]| { + IntVector::::new( + I::Wide::from_rounded_float(p[0] * 1_000_000.0), + I::Wide::from_rounded_float(p[1] * 1_000_000.0), + ) + .fast_normalize() + .unwrap() + }; + let input: Vec<_> = data + .iter() + .map(|c| (normalize(c.a), normalize(c.b), c.cw)) + .collect(); + let max_step = ((deg / 360.0 * 4294967296.0) as u32).clamp(1 << 22, 1 << 29); + let mut output = Vec::with_capacity(1024); + report( + name, + group, + deg, + measure(|| { + let mut count = 0; + for &(a, b, cw) in black_box(&input) { + output.clear(); + let sweep = if cw { + Angle::between(b, a) + } else { + Angle::between(a, b) + }; + let raw = sweep.bits() as u64; + if raw == 0 { + black_box(&output); + continue; + } + // Conservative vectoring error margin; integer ceiling, never floor. + let upper = raw + Angle::MAX_ERROR as u64; + let max_step = black_box(max_step) as u64; + let mut n = upper.div_ceil(max_step); + // i32/i64 only: reserve accumulated error in the final gap. + while upper + 8 * n * n > n * max_step { + n += 1; + } + if n > 1 { + let step = (raw / n) as u32; + let matrix = Rotation::new(Angle::from_bits(if cw { step.wrapping_neg() } else { step })); + let mut v = a; + for _ in 1..n { + v = matrix.apply(v); + output.push(v); + } + } + count += output.len(); + black_box(&output); + } + count + }), + ); +} + +fn main() { + println!("implementation,group,step_deg,median_ns_arc,min_ns_arc,max_ns_arc,points_arc"); + let mixed = cases(false); + let short = cases(true); + let reverse = std::env::args().any(|a| a == "reverse"); + for (group, data, steps) in [ + ("mixed", &mixed[..], &[45.0, 15.0, 5.0, 1.0, 360.0 / 1024.0][..]), + ("short", &short[..], &[45.0, 5.0][..]), + ] { + for ° in steps { + if reverse { + bench_cordic::(data, group, deg, "cordic_i64"); + bench_cordic::(data, group, deg, "cordic_i32"); + bench_int::(data, group, deg, "bisection_i64"); + bench_int::(data, group, deg, "bisection_i32"); + } + + bench_float::(data, group, deg, true, "f64_ceil"); + + bench_float::(data, group, deg, true, "f32_ceil"); + if !reverse { + bench_int::(data, group, deg, "bisection_i32"); + bench_int::(data, group, deg, "bisection_i64"); + bench_cordic::(data, group, deg, "cordic_i32"); + bench_cordic::(data, group, deg, "cordic_i64"); + } + } + } +} diff --git a/examples/cordic_precision.rs b/examples/cordic_precision.rs new file mode 100644 index 0000000..3939365 --- /dev/null +++ b/examples/cordic_precision.rs @@ -0,0 +1,717 @@ +//! Controlled runtime precision experiment; not an arc builder or public policy. +//! The CORDIC kernel is shared by source with the library, never copied. +#[path = "../src/int/angle/cordic.rs"] +mod cordic; + +use i_float::int::{ + angle::{Angle, Rotation}, + vector::IntVector, +}; +use std::{ + f64::consts::TAU, + hint::black_box, + time::{Duration, Instant}, +}; + +const COUNT: usize = 2048; +const SCALE: f64 = (1u64 << 30) as f64; +const BINARY: f64 = (1u64 << 32) as f64; + +#[derive(Clone, Copy, Debug)] +struct Precision { + iterations: usize, + q: u32, + relative: u32, // 0: explicit iterations; 8/16: rotation error divisor +} +impl Precision { + fn rotation_iterations(self, bits: u32) -> usize { + if self.relative == 0 { + return self.iterations; + } + cordic::rotation_iterations( + bits, + self.relative.trailing_zeros(), + Rotation::::ANGLE_MAX_ERROR, + ) + } +} + +#[derive(Debug)] +struct Options { + precision: Precision, + step: f64, + radius: u32, + sweep: bool, + reverse: bool, +} + +fn options(args: impl IntoIterator) -> Result, String> { + let mut result = Options { + precision: Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 0, + }, + step: 1.0, + radius: 65536, + sweep: false, + reverse: false, + }; + let mut args = args.into_iter(); + let mut explicit = false; + while let Some(arg) = args.next() { + match arg.as_str() { + "--help" | "-h" => return Ok(None), + "--sweep" => result.sweep = true, + "--reverse" => result.reverse = true, + "--iterations" | "--q" | "--relative" | "--step" | "--radius" => { + let value = args.next().ok_or_else(|| format!("missing value for {arg}"))?; + let invalid = || format!("invalid value for {arg}: {value}"); + match arg.as_str() { + "--iterations" => { + result.precision.iterations = value.parse().map_err(|_| invalid())?; + explicit = true; + } + "--q" => { + result.precision.q = value.parse().map_err(|_| invalid())?; + explicit = true; + } + "--relative" => { + result.precision.relative = value.parse().map_err(|_| invalid())?; + if ![8, 16].contains(&result.precision.relative) { + return Err("--relative must be 8 or 16".into()); + } + } + "--step" => result.step = value.parse().map_err(|_| invalid())?, + _ => result.radius = value.parse().map_err(|_| invalid())?, + } + } + _ => return Err(format!("unknown argument: {arg}")), + } + } + if !(1..=cordic::ITERATIONS).contains(&result.precision.iterations) { + return Err(format!("--iterations must be 1..{}", cordic::ITERATIONS)); + } + if !(14..=30).contains(&result.precision.q) { + return Err("--q must be 14..30".into()); + } + if !result.step.is_finite() || !(360.0 / 1024.0..=45.0).contains(&result.step) { + return Err("--step must be 0.3515625..45 degrees".into()); + } + if result.radius == 0 || result.radius > 1 << 30 { + return Err("--radius must be 1..1073741824".into()); + } + if result.precision.relative != 0 && ![8, 16].contains(&result.precision.relative) { + return Err("--relative must be 8 or 16".into()); + } + if result.precision.relative != 0 && explicit { + return Err(format!( + "--relative keeps atan2 at {} iterations and coefficients Q30; do not combine with --iterations/--q", + cordic::ITERATIONS + )); + } + if result.sweep && (explicit || result.precision.relative != 0) { + return Err("use --sweep or individual --iterations/--q/--relative values".into()); + } + Ok(Some(result)) +} + +#[derive(Clone, Copy)] +struct Case { + a: [i32; 2], + b: [i32; 2], + cw: bool, + sweep: f64, + segments: usize, +} + +fn cases(step: f64) -> Vec { + let mut state = 0xd718af03u64; + let mut next = || { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + (state >> 11) as f64 / (1u64 << 53) as f64 + }; + let normalize = |angle: f64| { + let v = IntVector::::new( + (angle.cos() * 1e6).round() as i64, + (angle.sin() * 1e6).round() as i64, + ) + .fast_normalize() + .unwrap(); + [v.x(), v.y()] + }; + (0..COUNT) + .map(|_| { + let start = next() * TAU; + let sweep = (0.1 + next() * 179.8).to_radians(); + let cw = next() < 0.5; + let a = normalize(start); + let b = normalize(start + if cw { -sweep } else { sweep }); + let mut sweep = (b[1] as f64).atan2(b[0] as f64) - (a[1] as f64).atan2(a[0] as f64); + if cw { + sweep = -sweep; + } + if sweep < 0.0 { + sweep += TAU; + } + // Fixed reference work across configurations. No low-precision profile + // can appear faster merely by underestimating its segment count. + let segments = (sweep / step.to_radians()).ceil().max(1.0) as usize; + Case { + a, + b, + cw, + sweep, + segments, + } + }) + .collect() +} + +fn atan2(y: i64, x: i64, iterations: usize) -> u32 { + if y == 0 { + return if x < 0 { 1 << 31 } else { 0 }; + } + if x == 0 { + return if y > 0 { 1 << 30 } else { 3 << 30 }; + } + let (ax, ay) = (x.unsigned_abs(), y.unsigned_abs()); + let top = ax.max(ay).ilog2(); + let input_top = cordic::VECTOR_INPUT_BITS - 1; + let reduce = |v: u64| -> i64 { + (if top > input_top { + v >> (top - input_top) + } else { + v << (input_top - top) + }) as i64 + * (1i64 << cordic::VECTOR_GUARD_BITS) + }; + let mut z = cordic::vectoring(reduce(ax), reduce(ay), iterations); + if x < 0 { + z = (1 << 31) - z; + } + let magnitude = z.clamp(1, (1 << 31) - 1) as u32; + if y < 0 { + magnitude.wrapping_neg() + } else { + magnitude + } +} + +#[derive(Clone, Copy)] +struct Matrix { + sin: i64, + cos: i64, +} +impl Matrix { + fn from_coefficients((sin, cos): (i32, i32), q: u32) -> Self { + // Q14..Q30 quantization, all stored in Q30 for the same i64 multiply + // and constant division in every profile. This isolates precision + // from the cost of runtime division or changing vector storage. + Self { + sin: (sin as i64) << (30 - q), + cos: (cos as i64) << (30 - q), + } + } + #[inline] + fn apply(self, v: [i32; 2]) -> [i32; 2] { + let (x, y) = (v[0] as i64, v[1] as i64); + [ + ((self.cos * x - self.sin * y) / (1 << 30)) as i32, + ((self.sin * x + self.cos * y) / (1 << 30)) as i32, + ] + } +} + +#[inline] +fn setup(case: Case, p: Precision) -> (u32, Matrix) { + let (ax, ay, bx, by) = ( + case.a[0] as i64, + case.a[1] as i64, + case.b[0] as i64, + case.b[1] as i64, + ); + let cross = ax * by - ay * bx; + let cross = if case.cw { -cross } else { cross }; + let dot = ax * bx + ay * by; + let sweep = if FIXED || p.relative != 0 { + Angle::atan2(cross, dot).unwrap().bits() + } else { + atan2(cross, dot, p.iterations) + }; + let step = sweep / case.segments as u32; + let step = if case.cw { step.wrapping_neg() } else { step }; + let coefficients = if FIXED { + Angle::from_bits(step).sin_cos() + } else { + cordic::sin_cos(step, p.rotation_iterations(step), p.q) + }; + ( + sweep, + Matrix::from_coefficients(coefficients, if FIXED { 30 } else { p.q }), + ) +} + +fn measure(mut f: impl FnMut()) -> (f64, f64, f64) { + f(); + let start = Instant::now(); + let mut loops = 0usize; + while start.elapsed() < Duration::from_millis(40) { + f(); + loops += 1; + } + let repeats = (50_000_000.0 / (start.elapsed().as_nanos() as f64 / loops as f64)).max(1.0) as usize; + let mut times = [0.0; 7]; + for time in &mut times { + let start = Instant::now(); + for _ in 0..repeats { + f(); + } + *time = start.elapsed().as_nanos() as f64 / (repeats * COUNT) as f64; + } + times.sort_by(f64::total_cmp); + (times[3], times[0], times[6]) +} + +#[derive(Default, Debug)] +struct Error { + arc: f64, + atan_bits: f64, + drift: f64, + radial: f64, + step_excess: f64, + bad_order: usize, + rotation_iterations_min: usize, + rotation_iterations_max: usize, + rotation_iterations_sum: usize, + rotation_relative_error: f64, +} + +fn accuracy(data: &[Case], p: Precision, step: f64) -> Error { + let mut error = Error { + rotation_iterations_min: cordic::ITERATIONS, + ..Error::default() + }; + for &case in data { + let (sweep, matrix) = setup::(case, p); + let step_bits = sweep / case.segments as u32; + let n = if FIXED { + cordic::ITERATIONS + } else { + p.rotation_iterations(step_bits) + }; + error.rotation_iterations_min = error.rotation_iterations_min.min(n); + error.rotation_iterations_max = error.rotation_iterations_max.max(n); + error.rotation_iterations_sum += n; + let requested = step_bits as f64 * TAU / BINARY; + if requested > 0.0 { + let actual = (matrix.sin as f64).atan2(matrix.cos as f64) * if case.cw { -1.0 } else { 1.0 }; + error.rotation_relative_error = error + .rotation_relative_error + .max((actual - requested).abs() / requested); + } + error.atan_bits = error + .atan_bits + .max((sweep as f64 - case.sweep * BINARY / TAU).abs()); + let mut v = case.a; + let direction = if case.cw { -1.0 } else { 1.0 }; + let mut travelled = 0.0; + for k in 1..=case.segments { + let previous = v; + v = matrix.apply(v); + let angle = direction * case.sweep * k as f64 / case.segments as f64; + let (sin, cos) = angle.sin_cos(); + let expected = [ + (cos * case.a[0] as f64 - sin * case.a[1] as f64) / SCALE, + (sin * case.a[0] as f64 + cos * case.a[1] as f64) / SCALE, + ]; + error.arc = error + .arc + .max((v[0] as f64 / SCALE - expected[0]).hypot(v[1] as f64 / SCALE - expected[1])); + if k < case.segments { + let (x, y, px, py) = (v[0] as f64, v[1] as f64, previous[0] as f64, previous[1] as f64); + let delta = direction * (px * y - py * x).atan2(px * x + py * y); + if delta <= 0.0 { + error.bad_order += 1; + } + travelled += delta; + error.step_excess = error.step_excess.max(delta - step.to_radians()); + } + } + // Gap from the last emitted interior point to the exact input end. + let final_gap = case.sweep - travelled; + if final_gap < 0.0 { + error.bad_order += 1; + } + error.step_excess = error.step_excess.max(final_gap - step.to_radians()); + } + for sample in 0..257u32 { + let start = sample as f64 * TAU / 257.0; + let initial = IntVector::::new( + (start.cos() * 1e6).round() as i64, + (start.sin() * 1e6).round() as i64, + ) + .fast_normalize() + .unwrap(); + let initial = [initial.x(), initial.y()]; + let initial_norm = (initial[0] as f64).hypot(initial[1] as f64) / SCALE; + let step = (1 << 22) + ((sample as u64 * ((1 << 29) - (1 << 22))) / 256) as u32; + for cw in [false, true] { + let bits = if cw { step.wrapping_neg() } else { step }; + let coeff = if FIXED { + Angle::from_bits(bits).sin_cos() + } else { + cordic::sin_cos(bits, p.rotation_iterations(bits), p.q) + }; + let matrix = Matrix::from_coefficients(coeff, if FIXED { 30 } else { p.q }); + let mut v = initial; + for k in 1..=1024 { + v = matrix.apply(v); + let angle = (if cw { -(step as f64) } else { step as f64 }) * k as f64 * TAU / BINARY; + let (sin, cos) = angle.sin_cos(); + let expected = [ + (cos * initial[0] as f64 - sin * initial[1] as f64) / SCALE, + (sin * initial[0] as f64 + cos * initial[1] as f64) / SCALE, + ]; + error.drift = error + .drift + .max((v[0] as f64 / SCALE - expected[0]).hypot(v[1] as f64 / SCALE - expected[1])); + error.radial = error + .radial + .max(initial_norm - (v[0] as f64).hypot(v[1] as f64) / SCALE); + } + } + } + error +} + +fn bench(data: &[Case], p: Precision, options: &Options) { + let step_bits: Vec<_> = data + .iter() + .map(|&case| { + let (sweep, _) = setup::(case, p); + let bits = sweep / case.segments as u32; + if case.cw { bits.wrapping_neg() } else { bits } + }) + .collect(); + let api_rotation_time = if FIXED || p.relative != 0 { + Some( + measure(|| { + let p = black_box(p); + for &bits in black_box(&step_bits) { + let angle = Angle::from_bits(bits); + black_box(if FIXED { + Rotation::::new(angle) + } else { + Rotation::::with_precision(angle, p.relative.trailing_zeros()) + }); + } + }) + .0, + ) + } else { + None + }; + let rotation_time = measure(|| { + let p = black_box(p); + for &bits in black_box(&step_bits) { + black_box(if FIXED { + Angle::from_bits(bits).sin_cos() + } else { + cordic::sin_cos(bits, p.rotation_iterations(bits), p.q) + }); + } + }); + let setup_time = measure(|| { + let p = black_box(p); + for &case in black_box(data) { + black_box(setup::(case, p)); + } + }); + let mut output = Vec::with_capacity(1024); + let arc_time = measure(|| { + let p = black_box(p); + for &case in black_box(data) { + output.clear(); + if case.segments > 1 { + let (_, matrix) = setup::(case, p); + let mut v = case.a; + for _ in 1..case.segments { + v = matrix.apply(v); + output.push(v); + } + } + black_box(&output); + } + }); + let error = accuracy::(data, p, options.step); + let points = data.iter().map(|c| c.segments - 1).sum::() as f64 / data.len() as f64; + let mode = if FIXED { + "library_fixed" + } else if p.relative != 0 { + "relative" + } else { + "runtime" + }; + println!( + "{mode},{},{},{:.8},{},{:.3},{:.3},{:.3},{:.3},{:.3},{:.3},{:.6},{:.6},{:.6},{:.6},{:.9},{},{},{},{},{:.3},{:.9},{:.3},{}", + p.iterations, + p.q, + options.step, + options.radius, + setup_time.0, + arc_time.0, + arc_time.1, + arc_time.2, + arc_time.0 / points.max(1.0), + points, + error.atan_bits, + error.arc * options.radius as f64, + error.drift * options.radius as f64, + error.radial * options.radius as f64, + error.step_excess.to_degrees(), + error.bad_order, + p.relative, + error.rotation_iterations_min, + error.rotation_iterations_max, + error.rotation_iterations_sum as f64 / data.len() as f64, + error.rotation_relative_error, + rotation_time.0, + api_rotation_time.map(|v| format!("{v:.3}")).unwrap_or_default() + ); +} + +fn main() { + let options = match options(std::env::args().skip(1)) { + Ok(Some(v)) => v, + Ok(None) => { + println!( + "Usage: cordic_precision [--iterations 1..{iterations}] [--q 14..30] [--relative 8|16] [--step degrees] [--radius integer] [--sweep] [--reverse]\nDefaults: {iterations} iterations, Q30, step 1 degree, radius 65536.\n--relative 8/16 automatically reduces only rotation iterations for an angular error budget step/8 or step/16; atan2 stays accurate and coefficients stay Q30.\n--sweep compares explicit and relative precision profiles; cannot combine with --iterations/--q/--relative.\nCSV timings use fixed reference segment counts and exclude input normalization. Errors include computed endpoints, before grid rounding. Experimental profiles have no maximum-step guarantee.", + iterations = cordic::ITERATIONS + ); + return; + } + Err(e) => { + eprintln!("{e}; use --help"); + std::process::exit(2); + } + }; + let data = cases(options.step); + let mut profiles = if options.sweep { + vec![ + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 0, + }, + Precision { + iterations: 24, + q: 30, + relative: 0, + }, + Precision { + iterations: 20, + q: 30, + relative: 0, + }, + Precision { + iterations: 15, + q: 30, + relative: 0, + }, + Precision { + iterations: cordic::ITERATIONS, + q: 24, + relative: 0, + }, + Precision { + iterations: cordic::ITERATIONS, + q: 20, + relative: 0, + }, + Precision { + iterations: cordic::ITERATIONS, + q: 14, + relative: 0, + }, + Precision { + iterations: 24, + q: 24, + relative: 0, + }, + Precision { + iterations: 20, + q: 20, + relative: 0, + }, + Precision { + iterations: 15, + q: 14, + relative: 0, + }, + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 8, + }, + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 16, + }, + ] + } else { + vec![options.precision] + }; + println!( + "mode,iterations,q,step_deg,radius,setup_median_ns,arc_median_ns,arc_min_ns,arc_max_ns,amortized_ns_point,points_arc,atan_max_error_bits,arc_max_error_units,drift_1024_max_error_units,radial_1024_max_loss_units,max_step_excess_deg,nonmonotone_gaps,relative_divisor,rotation_iterations_min,rotation_iterations_max,rotation_iterations_mean,rotation_max_relative_error,rotation_median_ns,api_rotation_median_ns" + ); + if !options.reverse { + bench::( + &data, + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 0, + }, + &options, + ); + } + if options.reverse { + profiles.reverse(); + } + for p in profiles { + bench::(&data, p, &options); + } + if options.reverse { + bench::( + &data, + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 0, + }, + &options, + ); + } +} + +#[cfg(test)] +mod tests { + use super::*; + #[test] + fn default_profile_matches_public_kernel_and_rotation() { + use i_float::int::angle::Rotation; + let p = Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 0, + }; + for case in cases(1.0) { + let (a, ma) = setup::(case, p); + let (b, mb) = setup::(case, p); + assert_eq!(a, b); + assert_eq!((ma.sin, ma.cos), (mb.sin, mb.cos)); + let from = IntVector::::new(case.a[0] as i64, case.a[1] as i64) + .fast_normalize() + .unwrap(); + let bits = if case.cw { + (a / case.segments as u32).wrapping_neg() + } else { + a / case.segments as u32 + }; + let rotated = Rotation::new(Angle::from_bits(bits)).apply(from); + assert_eq!(ma.apply([from.x(), from.y()]), [rotated.x(), rotated.y()]); + } + } + #[test] + fn every_profile_keeps_coefficients_nonexpanding() { + for iterations in 1..=cordic::ITERATIONS { + for q in 14..=30 { + for bits in (0..=u32::MAX).step_by(1_048_573) { + let (sin, cos) = cordic::sin_cos(bits, iterations, q); + assert!(sin as i64 * sin as i64 + cos as i64 * cos as i64 <= 1i64 << (2 * q)); + } + } + } + } + #[test] + fn relative_budget_selects_minimum_count_and_preserves_direction() { + for divisor in [8, 16] { + let p = Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: divisor, + }; + for sample in 0..65536u32 { + let bits = (BINARY * (0.1 + 44.9 * sample as f64 / 65535.0) / 360.0) as u32; + let n = p.rotation_iterations(bits); + assert_eq!(n, p.rotation_iterations(bits.wrapping_neg())); + let budget = bits >> divisor.trailing_zeros(); + let bound = |n: usize| 683_565_276u32.div_ceil(1 << (n - 1)) + 2; + assert!(bound(n) <= budget); + if n > 1 { + assert!(bound(n - 1) > budget); + } + for signed in [bits, bits.wrapping_neg()] { + let (s, c) = cordic::sin_cos(signed, n, 30); + let expected = (signed as i32 as f64) * TAU / BINARY; + let actual = (s as f64).atan2(c as f64); + assert!((actual - expected).abs() <= expected.abs() / divisor as f64); + assert_eq!(s.signum(), (signed as i32).signum()); + let public = Rotation::with_precision(Angle::from_bits(signed), divisor.trailing_zeros()); + let axis = IntVector::::new(1, 0).fast_normalize().unwrap(); + let v = public.apply(axis); + assert_eq!([v.y(), v.x()], [s, c]); + } + } + for step in [0.3515625, 1.0, 5.0, 20.0, 45.0] { + let errors = accuracy::(&cases(step), p, step); + assert!(errors.rotation_relative_error <= 1.0 / divisor as f64); + assert!(errors.atan_bits <= Angle::MAX_ERROR as f64); + } + } + let angle20 = (BINARY / 18.0) as u32; + assert_eq!( + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 8 + } + .rotation_iterations(angle20), + 6 + ); + assert_eq!( + Precision { + iterations: cordic::ITERATIONS, + q: 30, + relative: 16 + } + .rotation_iterations(angle20), + 7 + ); + } + + #[test] + fn cli_rejects_invalid_ranges_and_conflicting_modes() { + for args in [ + vec!["--q", "31"], + vec!["--iterations", "0"], + vec!["--step", "NaN"], + vec!["--radius", "0"], + vec!["--sweep", "--q", "14"], + vec!["--q"], + vec!["--unknown"], + vec!["--relative", "0"], + vec!["--relative", "10"], + vec!["--relative", "16", "--iterations", "15"], + vec!["--relative", "8", "--q", "30"], + vec!["--relative", "16", "--sweep"], + ] { + assert!(options(args.into_iter().map(str::to_owned)).is_err()); + } + assert!(options(["--iterations", "15", "--q", "14"].into_iter().map(str::to_owned)).is_ok()); + } +} diff --git a/examples/support/bisection.rs b/examples/support/bisection.rs new file mode 100644 index 0000000..b322325 --- /dev/null +++ b/examples/support/bisection.rs @@ -0,0 +1,235 @@ +// Benchmark-only snapshot of iOverlay mesh/int/arc.rs on 2026-09-13. +//! Directed unit-circle arcs built by integer bisection. + +use i_float::int::number::int::IntNumber; +use i_float::int::number::uint::UIntNumber; +use i_float::int::number::wide_int::WideIntNumber; +use i_float::int::unit_vector::UnitIntVector; +use i_float::int::vector::IntVector; +use std::vec::Vec; + +/// Direction of traversal in Cartesian coordinates (y increases upward). +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ArcDirection { + Clockwise, + Counterclockwise, +} + +/// Maximum subdivision step expressed as a squared unit-circle chord. +/// +/// The integer scale is `UnitIntVector::::DENOMINATOR.pow(2)`: a value +/// of one in real units represents a chord as long as the radius. This is +/// independent of the radius used later to place the arc in a mesh. +/// +/// Limits correspond approximately to 1/1024 of a turn and 45 degrees. +/// Actual angular accuracy depends on [`IntVector::fast_normalize`] +/// (especially for `i16`); this is not an exact angular-error guarantee. +#[derive(Clone, Copy)] +pub struct ArcStep { + squared_chord: I::WideUInt, +} + +impl ArcStep { + /// Clamps a squared chord in the unit-vector scale to the supported limits. + /// Zero selects the finest subdivision; large values select the coarsest. + pub fn from_squared_chord(squared_chord: I::WideUInt) -> Self { + Self { + squared_chord: squared_chord.clamp(Self::min_squared_chord(), Self::max_squared_chord()), + } + } + + /// Squared chord for approximately 360/1024 degrees, rounded upward. + pub fn min_squared_chord() -> I::WideUInt { + // ceil(4 * sin(pi / 1024)^2 * 2^60). No runtime trigonometry. + Self::rescale(43_406_843_014_579, true) + } + + /// Squared chord for approximately 45 degrees, rounded downward. + pub fn max_squared_chord() -> I::WideUInt { + // floor((2 - sqrt(2)) * 2^60). + Self::rescale(675_365_781_047_096_175, false) + } + + fn rescale(value: u64, round_up: bool) -> I::WideUInt { + let bits = 2 * (I::BITS - 2); + if bits >= 60 { + I::WideUInt::from_u64(value) << (bits - 60) + } else { + let shift = 60 - bits; + let value = if round_up { + (value + (1u64 << shift) - 1) >> shift + } else { + value >> shift + }; + I::WideUInt::from_u64(value) + } + } +} + +impl Default for ArcStep { + /// Selects the coarsest supported step (approximately 45 degrees). + fn default() -> Self { + Self { + squared_chord: Self::max_squared_chord(), + } + } +} + +/// Reusable storage for constructing directed arcs without trigonometry. +/// +/// Output contains intermediate directions in traversal order, excluding both +/// input endpoints. The caller applies the center and radius and retains the +/// original contact points to avoid rounding seams. Equal directions describe +/// an empty arc; full turns are not supported. Opposite directions describe a +/// semicircle whose side is selected by `ArcDirection`. +/// +/// Arcs are first split into spans of at most 90 degrees, then bisected at +/// `normalize(from + to)`. Each span has at most eight bisection levels, so a +/// build produces at most 1279 intermediate directions. Subdivision also stops +/// when rounding prevents a midpoint strictly inside its span. Both buffers +/// retain their capacity between builds. +/// +pub struct ArcBuilder { + stack: Vec>, + directions: Vec>, +} + +struct ArcSpan { + from: UnitIntVector, + to: UnitIntVector, + depth: u8, +} + +impl ArcBuilder { + pub fn new() -> Self { + Self { + stack: Vec::new(), + directions: Vec::new(), + } + } + + /// Clears the previous result and returns the intermediate directions. + pub fn build( + &mut self, + from: UnitIntVector, + to: UnitIntVector, + direction: ArcDirection, + step: ArcStep, + ) -> &[UnitIntVector] { + self.stack.clear(); + self.directions.clear(); + let end = vector(to); + let start = vector(from); + if start.cross_product(end) == I::Wide::ZERO && start.dot_product(end) > I::Wide::ZERO { + return &self.directions; + } + + let mut from = from; + // Split at exact cardinal directions so approximate normalization + // cannot rotate a boundary past the endpoint and select a major arc. + // At most four cardinal boundaries precede the final span. + for _ in 0..4 { + let start = vector(from); + if directed_cross(start, end, direction) >= I::Wide::ZERO + && start.dot_product(end) >= I::Wide::ZERO + { + break; + } + let boundary = next_axis(start, direction); + self.build_span(from, boundary, direction, step); + from = boundary; + } + self.build_span(from, to, direction, step); + // Leaves emit their end, including quarter boundaries, exactly once. + // Only the final input endpoint is omitted from the returned buffer. + self.directions.pop(); + &self.directions + } + + fn build_span( + &mut self, + from: UnitIntVector, + to: UnitIntVector, + direction: ArcDirection, + step: ArcStep, + ) { + self.stack.push(ArcSpan { from, to, depth: 0 }); + while let Some(span) = self.stack.pop() { + let a = vector(span.from); + let b = vector(span.to); + let chord = IntVector::::new(b.x - a.x, b.y - a.y); + if span.depth < 8 && chord.sqr_length() > step.squared_chord { + if let Some(mid) = IntVector::::new(a.x + b.x, a.y + b.y).fast_normalize() { + let m = vector(mid); + // Compare directions, not just stored components: unequal + // approximate lengths can still represent the same ray. + if directed_cross(a, m, direction) > I::Wide::ZERO + && directed_cross(m, b, direction) > I::Wide::ZERO + { + let depth = span.depth + 1; + self.stack.push(ArcSpan { + from: mid, + to: span.to, + depth, + }); + self.stack.push(ArcSpan { + from: span.from, + to: mid, + depth, + }); + continue; + } + } + } + self.directions.push(span.to); + } + } +} + +impl Default for ArcBuilder { + fn default() -> Self { + Self::new() + } +} + +fn vector(unit: UnitIntVector) -> IntVector { + IntVector::new(unit.x().to_wide(), unit.y().to_wide()) +} + +fn directed_cross(a: IntVector, b: IntVector, direction: ArcDirection) -> I::Wide { + let cross = a.cross_product(b); + match direction { + ArcDirection::Counterclockwise => cross, + ArcDirection::Clockwise => -cross, + } +} + +fn next_axis(v: IntVector, direction: ArcDirection) -> UnitIntVector { + let zero = I::Wide::ZERO; + let one = I::Wide::ONE; + let axis = match direction { + ArcDirection::Counterclockwise => { + if v.x > zero && v.y >= zero { + IntVector::new(zero, one) + } else if v.x <= zero && v.y > zero { + IntVector::new(-one, zero) + } else if v.x < zero && v.y <= zero { + IntVector::new(zero, -one) + } else { + IntVector::new(one, zero) + } + } + ArcDirection::Clockwise => { + if v.x >= zero && v.y > zero { + IntVector::new(one, zero) + } else if v.x > zero && v.y <= zero { + IntVector::new(zero, -one) + } else if v.x <= zero && v.y < zero { + IntVector::new(-one, zero) + } else { + IntVector::new(zero, one) + } + } + }; + axis.fast_normalize().expect("an axis direction is nonzero") +} diff --git a/src/adapter.rs b/src/adapter.rs index 12e8436..234148f 100644 --- a/src/adapter.rs +++ b/src/adapter.rs @@ -1,91 +1,146 @@ use crate::float::compatible::FloatPointCompatible; use crate::float::number::FloatNumber; use crate::float::point::FloatPoint; -use crate::float::rect::FloatRect; +use crate::float::rect::{FloatRect, FloatRectError}; use crate::int::number::int::IntNumber; use crate::int::number::wide_int::WideIntNumber; use crate::int::point::IntPoint; -use core::marker::PhantomData; +use crate::int::rect::IntRect; #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum FloatPointAdapterScaleError { + /// Input bounds violate the floating-point coordinate contract. + InvalidRect(FloatRectError), /// Requested scale is larger than the safe adapter scale for the input bounds. ScaleTooLarge, + /// Scale is too small to have a finite reciprocal in the scalar type. + ScaleTooSmall, /// Requested scale is zero or negative. ScaleNonPositive, /// Requested scale is NaN or infinite. ScaleNotFinite, } +impl From for FloatPointAdapterScaleError { + fn from(error: FloatRectError) -> Self { + Self::InvalidRect(error) + } +} + #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum FloatPointAdapterRangeError { - /// Point is outside the adapter rectangle. + /// Point is outside the source rectangle or its enclosing integer grid. PointOutOfRange, } #[derive(Clone)] /// Maps a bounded floating-point coordinate space onto an integer grid. /// +/// Input coordinates and rectangle bounds must be finite, with absolute values +/// at most `2^60` for `f32` or `2^500` for `f64`. Bounds must satisfy `min <= max` +/// on each axis. See the [`crate::float`] coordinate-range contract. +/// Constructors validate these bounds, including rectangles assembled through +/// public fields. Fallible constructors return [`FloatPointAdapterScaleError::InvalidRect`]; +/// infallible constructors panic. Coordinates are never clamped. +/// The coordinate limits do not apply to the stored scales. +/// /// The adapter controls converted coordinate magnitude, but it cannot prove /// that every expression in a downstream integer algorithm is safe. For a /// conservative general-purpose budget covering point differences, dot and -/// cross products, and squared lengths, use [`Self::with_coordinate_bits`] -/// with `coordinate_bits <= I::BITS - 3`. A larger value is appropriate only -/// when the downstream algorithm has a stronger range analysis. +/// cross products, and squared lengths, use the conservative constructors such +/// as [`Self::new_conservative`] or [`Self::with_iter_conservative`]. They use +/// [`Self::CONSERVATIVE_COORDINATE_BITS`], reserving an extra bit for coordinate +/// rounding. Larger explicit budgets require a stronger downstream range analysis. pub struct FloatPointAdapter { dir_scale: P::Scalar, inv_scale: P::Scalar, offset: P, rect: FloatRect, - int: PhantomData, + int_rect: IntRect, } impl FloatPointAdapter { const SCALE_SAFETY_BITS: i32 = 3; + /// Conservative coordinate-bit budget for point differences and their products. + /// + /// Converted coordinates have magnitude at most `2^(I::BITS - 3)`, including + /// both endpoints. This reserves one extra bit for rounding inside the strict + /// arithmetic range `(-2^(I::BITS - 2), 2^(I::BITS - 2))` of [`IntPoint`]. + /// For `i16`, `i32`, and `i64`, the budgets are 13, 29, and 61 bits. + /// This does not guarantee that arbitrary downstream arithmetic fits. + pub const CONSERVATIVE_COORDINATE_BITS: u32 = I::BITS - 3; + + /// Creates an automatically scaled adapter with the conservative coordinate budget. + /// See [`Self::CONSERVATIVE_COORDINATE_BITS`] for the range and rounding margin. + /// + /// # Panics + /// Panics for invalid bounds, as in [`Self::with_coordinate_bits`]. + #[inline] + pub fn new_conservative(rect: FloatRect) -> Self { + Self::with_coordinate_bits(rect, Self::CONSERVATIVE_COORDINATE_BITS) + } + + /// Creates an automatically scaled adapter from points with the conservative budget. + /// Empty input uses a zero rectangle and scale one. + /// + /// # Panics + /// Panics for invalid coordinates, as in [`Self::with_iter_and_coordinate_bits`]. + #[inline] + pub fn with_iter_conservative<'a, Q>(iter: Q) -> Self + where + Q: Iterator, + P: 'a, + { + Self::with_iter_and_coordinate_bits(iter, Self::CONSERVATIVE_COORDINATE_BITS) + } + + /// Validates bounds and an explicit scale against the conservative coordinate budget. + /// Preserves the requested scale, including for single-point bounds. Returns the + /// same errors as [`Self::try_with_scale_and_coordinate_bits`]. + #[inline] + pub fn try_with_scale_conservative( + rect: FloatRect, + scale: P::Scalar, + ) -> Result { + Self::try_with_scale_and_coordinate_bits(rect, scale, Self::CONSERVATIVE_COORDINATE_BITS) + } + + /// Creates an adapter from points with an explicit scale and the conservative budget. + /// Empty input uses a zero rectangle. Preserves the requested scale, including for + /// empty or single-point input. Returns the same errors as + /// [`Self::try_with_iter_and_scale_and_coordinate_bits`]. + #[inline] + pub fn try_with_iter_and_scale_conservative<'a, Q>( + iter: Q, + scale: P::Scalar, + ) -> Result + where + Q: Iterator, + P: 'a, + { + Self::try_with_iter_and_scale_and_coordinate_bits(iter, scale, Self::CONSERVATIVE_COORDINATE_BITS) + } + /// Creates an adapter and automatically selects a power-of-two scale while /// reserving internal coordinate safety bits. /// /// Algorithms with a precise bit budget should prefer /// [`Self::with_coordinate_bits`]. + /// The scale is capped at the largest finite power of two in the scalar + /// type, which may reduce grid precision for very small bounds. + /// Within the supported coordinate range, automatic scales and their + /// reciprocals are finite. + /// + /// # Panics + /// Panics if rectangle bounds are non-finite, reversed, or exceed the + /// supported coordinate limits (`2^60` for `f32`, `2^500` for `f64`). #[inline] pub fn new(rect: FloatRect) -> Self { - let a = rect.width() * P::Scalar::HALF; - let b = rect.height() * P::Scalar::HALF; - - let x = rect.min_x + a; - let y = rect.min_y + b; - - let offset = P::from_xy(x, y); - - let max = a.max(b); - - // degenerate case - if max == P::Scalar::ZERO { - return Self { - dir_scale: P::Scalar::ONE, - inv_scale: P::Scalar::ONE, - offset, - rect, - int: PhantomData, - }; - } - - let log2 = max.log2().to_i32(); - let safe_bits = I::BITS as i32 - Self::SCALE_SAFETY_BITS; - let ie = safe_bits - log2; - let e = ie as f64; - - let dir_scale = FloatNumber::from_float(libm::exp2(e)); - let inv_scale = FloatNumber::from_float(libm::exp2(-e)); - - Self { - dir_scale, - inv_scale, - offset, - rect, - int: PhantomData, - } + rect.validate().expect("Invalid adapter bounds"); + let (offset, radius) = Self::center_and_radius(&rect); + let (dir_scale, inv_scale) = Self::automatic_scales(radius, None); + Self::from_transform(rect, offset, dir_scale, inv_scale) } /// Creates an adapter whose converted coordinates have magnitude at most @@ -96,137 +151,130 @@ impl FloatPointAdapter { /// magnitude but does not independently guarantee that every downstream /// product fits. `I::BITS - 3` is a conservative budget for the general /// point and vector operations provided by this crate. + /// The scale is capped at the largest finite power of two in the scalar type. + /// + /// # Panics + /// Panics for invalid rectangle bounds or if `coordinate_bits > I::BITS - 2`. #[inline] pub fn with_coordinate_bits(rect: FloatRect, coordinate_bits: u32) -> Self { + rect.validate().expect("Invalid adapter bounds"); assert!(coordinate_bits <= I::BITS - 2); - let a = rect.width() * P::Scalar::HALF; - let b = rect.height() * P::Scalar::HALF; - - let x = rect.min_x + a; - let y = rect.min_y + b; - - let offset = P::from_xy(x, y); - let max = a.max(b); - - // degenerate case - if max == P::Scalar::ZERO { - return Self { - dir_scale: P::Scalar::ONE, - inv_scale: P::Scalar::ONE, - offset, - rect, - int: PhantomData, - }; - } - - let log2 = max.log2().to_i32(); - let log2_scale = P::Scalar::from_float(libm::exp2(log2 as f64)); - let ceil_log2 = if log2_scale < max { log2 + 1 } else { log2 }; - let exponent = coordinate_bits as i32 - ceil_log2; - let e = exponent as f64; - - let dir_scale = P::Scalar::from_float(libm::exp2(e)); - let inv_scale = P::Scalar::from_float(libm::exp2(-e)); - - Self { - dir_scale, - inv_scale, - offset, - rect, - int: PhantomData, - } + let (offset, radius) = Self::center_and_radius(&rect); + let (dir_scale, inv_scale) = Self::automatic_scales(radius, Some(coordinate_bits)); + Self::from_transform(rect, offset, dir_scale, inv_scale) } + /// Creates an adapter with an explicit scale, without checking the coordinate + /// budget. A rectangle containing a single point uses scale one, as in + /// [`Self::new`]. + /// + /// # Panics + /// Panics for invalid rectangle bounds or if the supplied scale is + /// non-positive, non-finite, or has a non-finite reciprocal. + /// Use [`Self::try_with_scale`] to receive an error. #[inline] pub fn with_scale(rect: FloatRect, scale: P::Scalar) -> Self { - let a = rect.width() * P::Scalar::HALF; - let b = rect.height() * P::Scalar::HALF; - - let x = rect.min_x + a; - let y = rect.min_y + b; - - let offset = P::from_xy(x, y); - - let max = a.max(b); - - // degenerate case - if max == P::Scalar::ZERO { - return Self { - dir_scale: P::Scalar::ONE, - inv_scale: P::Scalar::ONE, - offset, - rect, - int: PhantomData, - }; - } - - let dir_scale = scale; - let inv_scale = P::Scalar::ONE / scale; - - Self { - dir_scale, - inv_scale, - offset, - rect, - int: PhantomData, + rect.validate().expect("Invalid adapter bounds"); + let inv_scale = Self::inverse_scale(scale).expect("Invalid adapter scale"); + let (offset, radius) = Self::center_and_radius(&rect); + if radius == P::Scalar::ZERO { + Self::from_transform(rect, offset, P::Scalar::ONE, P::Scalar::ONE) + } else { + Self::from_transform(rect, offset, scale, inv_scale) } } + /// Validates rectangle bounds and an explicit scale against the default budget. + /// Returns [`FloatPointAdapterScaleError::ScaleTooSmall`] if its reciprocal + /// is non-finite. The requested scale is retained, including for point bounds. #[inline] pub fn try_with_scale( rect: FloatRect, scale: P::Scalar, ) -> Result { - let scale = Self::validate_scale(scale)?; - let mut adapter = Self::new(rect); - - let is_degenerate = - adapter.rect.width() == P::Scalar::ZERO && adapter.rect.height() == P::Scalar::ZERO; - if !is_degenerate && adapter.dir_scale < scale { - return Err(FloatPointAdapterScaleError::ScaleTooLarge); - } - - adapter.dir_scale = scale; - adapter.inv_scale = P::Scalar::ONE / scale; - Ok(adapter) + Self::checked_with_scale(rect, scale, None) } /// Creates an adapter with an explicit scale, rejecting scales that can /// produce coordinates with magnitude greater than `2^coordinate_bits`. + /// The requested scale must also have a finite reciprocal. Rectangle bounds + /// are validated against the floating-point coordinate contract. #[inline] pub fn try_with_scale_and_coordinate_bits( rect: FloatRect, scale: P::Scalar, coordinate_bits: u32, ) -> Result { - let scale = Self::validate_scale(scale)?; - let mut adapter = Self::with_coordinate_bits(rect, coordinate_bits); - - let zero = P::Scalar::ZERO; - let is_degenerate = adapter.rect.width() == zero && adapter.rect.height() == zero; - if !is_degenerate && adapter.dir_scale < scale { - return Err(FloatPointAdapterScaleError::ScaleTooLarge); - } - - adapter.dir_scale = scale; - adapter.inv_scale = P::Scalar::ONE / scale; - Ok(adapter) + assert!(coordinate_bits <= I::BITS - 2); + Self::checked_with_scale(rect, scale, Some(coordinate_bits)) } #[inline] pub fn with_radius_and_scale(radius: P::Scalar, scale: P::Scalar) -> FloatPointAdapter { - let rect = FloatRect::new(-radius, radius, -radius, radius); + let rect = FloatRect::new(-radius, radius, -radius, radius).unwrap(); FloatPointAdapter::with_scale(rect, scale) } + /// Creates an adapter from input points. + /// An empty iterator uses a zero rectangle and scale one. + /// + /// # Panics + /// Panics if any coordinate is non-finite or exceeds the supported absolute + /// limit (`2^60` for `f32`, `2^500` for `f64`). #[inline] pub fn with_iter<'a, Q>(iter: Q) -> Self where Q: Iterator, P: 'a, { - Self::new(FloatRect::with_iter(iter).unwrap_or(FloatRect::zero())) + Self::new( + FloatRect::with_iter(iter) + .expect("Invalid adapter bounds") + .unwrap_or(FloatRect::zero()), + ) + } + + /// Creates an adapter from input points with an explicit coordinate-bit budget. + /// See [`Self::with_coordinate_bits`] for the scale and range guarantees. + /// An empty iterator uses a zero rectangle and scale one. + /// + /// # Panics + /// Panics if any coordinate violates the floating-point coordinate contract + /// or if `coordinate_bits > I::BITS - 2`. + #[inline] + pub fn with_iter_and_coordinate_bits<'a, Q>(iter: Q, coordinate_bits: u32) -> Self + where + Q: Iterator, + P: 'a, + { + let rect = FloatRect::with_iter(iter) + .expect("Invalid adapter bounds") + .unwrap_or(FloatRect::zero()); + Self::with_coordinate_bits(rect, coordinate_bits) + } + + /// Creates an adapter from input points with a checked explicit scale and + /// coordinate-bit budget. See [`Self::try_with_scale_and_coordinate_bits`] + /// for scale validation. Invalid input coordinates return + /// [`FloatPointAdapterScaleError::InvalidRect`]. + /// An empty iterator uses a zero rectangle. The requested scale is retained, + /// including for empty input or bounds containing a single point. + /// + /// # Panics + /// Panics if `coordinate_bits > I::BITS - 2`. + #[inline] + pub fn try_with_iter_and_scale_and_coordinate_bits<'a, Q>( + iter: Q, + scale: P::Scalar, + coordinate_bits: u32, + ) -> Result + where + Q: Iterator, + P: 'a, + { + let rect = FloatRect::with_iter(iter)?.unwrap_or(FloatRect::zero()); + Self::try_with_scale_and_coordinate_bits(rect, scale, coordinate_bits) } #[inline] @@ -238,18 +286,115 @@ impl FloatPointAdapter { Q: Iterator, P: 'a, { - Self::try_with_scale(FloatRect::with_iter(iter).unwrap_or(FloatRect::zero()), scale) + Self::try_with_scale(FloatRect::with_iter(iter)?.unwrap_or(FloatRect::zero()), scale) + } + + #[inline] + fn center_and_radius(rect: &FloatRect) -> (P, P::Scalar) { + let x = (rect.min_x + rect.max_x) * P::Scalar::HALF; + let y = (rect.min_y + rect.max_y) * P::Scalar::HALF; + // The rounded center can lie on a bound for adjacent floats. Measure + // from it rather than halving the span, which can underestimate radius + // or underflow to zero for a one-subnormal span. + let radius = (x - rect.min_x) + .max(rect.max_x - x) + .max(y - rect.min_y) + .max(rect.max_y - y); + (P::from_xy(x, y), radius) } #[inline] - fn validate_scale(scale: P::Scalar) -> Result { + fn scale_exponent(radius: P::Scalar, coordinate_bits: Option) -> i32 { + let log2 = radius.log2().to_i32(); + if let Some(bits) = coordinate_bits { + let power = P::Scalar::from_float(libm::exp2(log2 as f64)); + let ceil_log2 = if power < radius { log2 + 1 } else { log2 }; + bits as i32 - ceil_log2 + } else { + I::BITS as i32 - Self::SCALE_SAFETY_BITS - log2 + } + } + + #[inline] + fn automatic_scales(radius: P::Scalar, coordinate_bits: Option) -> (P::Scalar, P::Scalar) { + if radius == P::Scalar::ZERO { + return (P::Scalar::ONE, P::Scalar::ONE); + } + + // Validated coordinate bounds keep the reciprocal finite, even with + // a zero-bit budget. Only tiny bounds need the upper scale cap. + let max_exponent = P::Scalar::MAX_EXP - 1; + let exponent = Self::scale_exponent(radius, coordinate_bits).min(max_exponent); + let scale = P::Scalar::from_float(libm::exp2(exponent as f64)); + (scale, P::Scalar::ONE / scale) + } + + #[inline] + fn checked_with_scale( + rect: FloatRect, + scale: P::Scalar, + coordinate_bits: Option, + ) -> Result { + rect.validate()?; + let inv_scale = Self::inverse_scale(scale)?; + let (offset, radius) = Self::center_and_radius(&rect); + if radius != P::Scalar::ZERO { + let exponent = Self::scale_exponent(radius, coordinate_bits); + // This comparison limit may exceed the scalar range. Unlike the + // stored automatic scale, it must not cap a valid explicit scale. + let limit = P::Scalar::from_float(libm::exp2(exponent as f64)); + if limit < scale { + return Err(FloatPointAdapterScaleError::ScaleTooLarge); + } + } + Ok(Self::from_transform(rect, offset, scale, inv_scale)) + } + + #[inline] + fn from_transform( + rect: FloatRect, + offset: P, + dir_scale: P::Scalar, + inv_scale: P::Scalar, + ) -> Self { + Self { + dir_scale, + inv_scale, + offset, + rect, + int_rect: Self::grid_rect(&rect, &offset, dir_scale), + } + } + + #[inline] + fn grid_rect(rect: &FloatRect, offset: &P, scale: P::Scalar) -> IntRect { + // Enclose the source bounds in the existing grid without changing its + // scale or origin. Use the same scalar arithmetic as float_to_int. + let min_x = ((rect.min_x - offset.x()) * scale).to_f64(); + let max_x = ((rect.max_x - offset.x()) * scale).to_f64(); + let min_y = ((rect.min_y - offset.y()) * scale).to_f64(); + let max_y = ((rect.max_y - offset.y()) * scale).to_f64(); + IntRect::new( + I::from_float(libm::floor(min_x)), + I::from_float(libm::ceil(max_x)), + I::from_float(libm::floor(min_y)), + I::from_float(libm::ceil(max_y)), + ) + } + + #[inline] + fn inverse_scale(scale: P::Scalar) -> Result { if !scale.is_finite() { return Err(FloatPointAdapterScaleError::ScaleNotFinite); } if scale <= P::Scalar::ZERO { return Err(FloatPointAdapterScaleError::ScaleNonPositive); } - Ok(scale) + let inverse = P::Scalar::ONE / scale; + if !inverse.is_finite() { + return Err(FloatPointAdapterScaleError::ScaleTooSmall); + } + Ok(inverse) } #[inline(always)] @@ -267,45 +412,39 @@ impl FloatPointAdapter { self.offset } + /// Returns the original floating-point bounds used to validate input points. + /// Snapped points may lie outside these bounds, on the enclosing integer grid. #[inline(always)] pub fn rect(&self) -> &FloatRect { &self.rect } + /// Converts a point from the enclosing integer grid back to floating point. + /// Debug builds assert that the point is inside that grid's bounds. #[inline(always)] pub fn int_to_float(&self, point: &IntPoint) -> P { + debug_assert!( + self.int_rect.contains(*point), + "Integer point [{}, {}] is outside the adapter grid", + point.x, + point.y + ); let fx: P::Scalar = FloatNumber::from_int(point.x); let fy: P::Scalar = FloatNumber::from_int(point.y); let x = fx * self.inv_scale + self.offset.x(); let y = fy * self.inv_scale + self.offset.y(); - let float = P::from_xy(x, y); - - if cfg!(debug_assertions) { - let radius = self.rect.height().max(self.rect.width()) * P::Scalar::from_float(0.01); - if !self.rect.contains_with_radius(&float, radius) { - panic!( - "You are trying to convert a point[{}, {}] which is out of rect: {}", - x, y, self.rect - ); - } - } - - float + P::from_xy(x, y) } + /// Validates the integer point against the enclosing grid before conversion. + /// The result may lie outside [`Self::rect`]. #[inline(always)] pub fn try_int_to_float(&self, point: &IntPoint) -> Result { - let fx: P::Scalar = FloatNumber::from_int(point.x); - let fy: P::Scalar = FloatNumber::from_int(point.y); - let x = fx * self.inv_scale + self.offset.x(); - let y = fy * self.inv_scale + self.offset.y(); - let float = P::from_xy(x, y); - - if self.rect.contains(&float) { - Ok(float) - } else { - Err(FloatPointAdapterRangeError::PointOutOfRange) + if !self.int_rect.contains(*point) { + return Err(FloatPointAdapterRangeError::PointOutOfRange); } + + Ok(self.int_to_float(point)) } #[inline(always)] @@ -332,7 +471,7 @@ impl FloatPointAdapter { #[inline(always)] pub fn try_float_to_int(&self, point: &P) -> Result, FloatPointAdapterRangeError> { - if !point.is_finite() || !self.rect.contains(point) { + if !self.rect.contains(point) { return Err(FloatPointAdapterRangeError::PointOutOfRange); } @@ -344,22 +483,26 @@ impl FloatPointAdapter { Ok(IntPoint { x, y }) } + /// Rounds to the nearest grid point, which may lie outside [`Self::rect`]. #[inline(always)] pub fn snap_to_grid(&self, point: &P) -> P { self.int_to_float(&self.float_to_int(point)) } + /// Checks the input against [`Self::rect`] and rounds it to the nearest grid + /// point. Rounding may move the result outside the original float bounds. #[inline(always)] pub fn try_snap_to_grid(&self, point: &P) -> Result { self.try_float_to_int(point) - .map(|point| self.int_to_float(&point)) + .and_then(|point| self.try_int_to_float(&point)) } #[inline(always)] pub fn round_sqr_len_to_int(&self, value: P::Scalar) -> I::Wide { let scale = self.dir_scale; - let sqr_scale = scale * scale; - I::Wide::from_rounded_float(sqr_scale * value) + // Multiply the area first: scale * scale can overflow even when + // the scaled area fits, especially with f32 and the i64 engine. + I::Wide::from_rounded_float((value * scale) * scale) } #[inline(always)] @@ -380,7 +523,7 @@ impl FloatPointAdapter { inv_scale: self.inv_scale, offset: FloatPoint::from_point(self.offset), rect: self.rect, - int: self.int, + int_rect: self.int_rect, } } } @@ -389,10 +532,34 @@ impl FloatPointAdapter { mod tests { use crate::adapter::{FloatPointAdapter, FloatPointAdapterRangeError, FloatPointAdapterScaleError}; use crate::float::compatible::FloatPointCompatible; + use crate::float::number::FloatNumber; use crate::float::point::FloatPoint; use crate::float::rect::FloatRect; + use crate::int::number::int::IntNumber; use crate::int::point::IntPoint; + #[test] + fn round_sqr_len_to_int_avoids_intermediate_overflow() { + let points = [[0.0_f32, 0.0], [0.01, 0.01]]; + let adapter = FloatPointAdapter::<[f32; 2], i64>::with_iter(points.iter()); + let scale = adapter.dir_scale(); + assert!(scale.is_finite()); + assert!((scale * scale).is_infinite()); + + let value = 1e-6_f32; + let wide_scale = f64::from(scale); + let expected = (wide_scale * f64::from(value) * wide_scale) as i128; + assert!(expected > 0 && expected < i128::MAX); + assert_eq!(adapter.round_sqr_len_to_int(value), expected); + assert_eq!(adapter.round_sqr_len_to_int(0.0), 0); + + // Round the final area, rather than a length that is then squared. + for (area, expected) in [(1.5_f64, 2_i128), (2.5, 3)] { + let value = (area / wide_scale / wide_scale) as f32; + assert_eq!(adapter.round_sqr_len_to_int(value), expected); + } + } + #[test] fn test_0() { let rect = FloatRect { @@ -462,7 +629,7 @@ mod tests { #[test] fn coordinate_bits_bound_integer_magnitude() { let adapter = FloatPointAdapter::<[f64; 2], i32>::with_coordinate_bits( - FloatRect::new(-3.0, 3.0, -0.25, 0.25), + FloatRect::new(-3.0, 3.0, -0.25, 0.25).unwrap(), 10, ); @@ -474,7 +641,7 @@ mod tests { #[test] fn coordinate_bits_include_exact_power_of_two_boundary() { let adapter = FloatPointAdapter::<[f64; 2], i32>::with_coordinate_bits( - FloatRect::new(-4.0, 4.0, -0.25, 0.25), + FloatRect::new(-4.0, 4.0, -0.25, 0.25).unwrap(), 10, ); @@ -485,8 +652,10 @@ mod tests { #[test] fn test_round_point_round_length() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ); assert_eq!(adapter.float_to_int(&[0.16, -0.16]), IntPoint::new(2, -2)); assert_eq!(adapter.round_len_to_int(0.16), 2); @@ -494,8 +663,10 @@ mod tests { #[test] fn test_snap_to_grid() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ); assert_eq!(adapter.snap_to_grid(&[0.16, -0.16]), [0.2, -0.2]); assert_eq!(adapter.snap_to_grid(&[0.14, -0.14]), [0.1, -0.1]); @@ -503,8 +674,10 @@ mod tests { #[test] fn test_try_snap_to_grid_range() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ); assert_eq!(adapter.try_snap_to_grid(&[1.0, -1.0]).unwrap(), [1.0, -1.0]); assert_eq!( @@ -513,11 +686,129 @@ mod tests { ); } + #[test] + fn test_try_snap_to_grid_accepts_rounded_point_outside_float_rect() { + let rect = FloatRect::new(-0.6, 0.6, -0.6, 0.6).unwrap(); + let adapter = FloatPointAdapter::<[f64; 2], i32>::try_with_scale(rect, 1.0).unwrap(); + + assert_eq!(adapter.try_snap_to_grid(&[0.4, -0.4]), Ok([0.0, 0.0])); + + for (point, snapped) in [ + ([0.6, 0.0], [1.0, 0.0]), + ([-0.6, 0.0], [-1.0, 0.0]), + ([0.0, 0.6], [0.0, 1.0]), + ([0.0, -0.6], [0.0, -1.0]), + ] { + assert!(rect.contains(&point)); + assert!(!adapter.rect().contains(&snapped)); + assert_eq!(adapter.try_snap_to_grid(&point), Ok(snapped)); + assert_eq!(adapter.snap_to_grid(&point), snapped); + } + // Expanding the grid must not expand the accepted float inputs. + assert_eq!( + adapter.try_snap_to_grid(&[0.61, 0.0]), + Err(FloatPointAdapterRangeError::PointOutOfRange) + ); + assert_eq!( + adapter.try_int_to_float(&IntPoint::new(2, 0)), + Err(FloatPointAdapterRangeError::PointOutOfRange) + ); + } + + fn assert_grid_bounds() { + let f = F::from_float::; + for rect in [ + FloatRect::new(f(9.375), f(10.625), f(-2.625), f(-1.375)).unwrap(), + FloatRect::new(f(10.0), f(10.0), f(-2.625), f(-1.375)).unwrap(), + FloatRect::new(f(9.375), f(10.625), f(-2.0), f(-2.0)).unwrap(), + FloatRect::new(f(10.0), f(10.0), f(-2.0), f(-2.0)).unwrap(), + ] { + for adapter in [ + FloatPointAdapter::<[F; 2], I>::new(rect), + FloatPointAdapter::with_coordinate_bits(rect, 2), + FloatPointAdapter::with_scale(rect, F::ONE), + FloatPointAdapter::try_with_scale(rect, F::ONE).unwrap(), + FloatPointAdapter::try_with_scale_and_coordinate_bits(rect, F::ONE, 2).unwrap(), + ] { + let copied = adapter.clone().to_float_point_adapter(); + for x in [rect.min_x, rect.max_x] { + for y in [rect.min_y, rect.max_y] { + let integer = adapter.try_float_to_int(&[x, y]).unwrap(); + let restored = adapter.try_int_to_float(&integer).unwrap(); + let converted = copied.try_int_to_float(&integer).unwrap(); + assert!(restored == adapter.int_to_float(&integer)); + assert!(restored == adapter.try_snap_to_grid(&[x, y]).unwrap()); + assert!(restored == [converted.x, converted.y]); + } + } + // Reject the next grid cell on every side, including after a + // checked constructor changes its automatically chosen scale. + let min = adapter.try_float_to_int(&[rect.min_x, rect.min_y]).unwrap(); + let max = adapter.try_float_to_int(&[rect.max_x, rect.max_y]).unwrap(); + for point in [ + IntPoint::new(min.x - I::ONE, I::ZERO), + IntPoint::new(max.x + I::ONE, I::ZERO), + IntPoint::new(I::ZERO, min.y - I::ONE), + IntPoint::new(I::ZERO, max.y + I::ONE), + ] { + assert!(adapter.try_int_to_float(&point).is_err()); + assert!(copied.try_int_to_float(&point).is_err()); + } + } + } + } + + #[test] + fn grid_bounds_cover_all_coordinate_types_and_constructors() { + assert_grid_bounds::(); + assert_grid_bounds::(); + assert_grid_bounds::(); + assert_grid_bounds::(); + assert_grid_bounds::(); + assert_grid_bounds::(); + } + + #[test] + fn integer_range_check_preserves_precision_beyond_f64() { + let limit = 1_i64 << 60; + let adapter = FloatPointAdapter::<[f64; 2], i64>::try_with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + limit as f64, + ) + .unwrap(); + assert_eq!(limit as f64, (limit + 1) as f64); + for sign in [-1, 1] { + assert!(adapter.try_int_to_float(&IntPoint::new(sign * limit, 0)).is_ok()); + assert_eq!( + adapter.try_int_to_float(&IntPoint::new(sign * (limit + 1), 0)), + Err(FloatPointAdapterRangeError::PointOutOfRange) + ); + } + } + + #[test] + fn grid_bounds_enclose_even_inward_rounded_float_bounds() { + let adapter = FloatPointAdapter::<[f64; 2], i32>::try_with_scale( + FloatRect::new(-0.25, 0.25, -0.25, 0.25).unwrap(), + 1.0, + ) + .unwrap(); + assert_eq!(adapter.try_snap_to_grid(&[0.25, -0.25]), Ok([0.0, 0.0])); + assert_eq!(adapter.try_int_to_float(&IntPoint::new(-1, 1)), Ok([-1.0, 1.0])); + assert_eq!(adapter.try_int_to_float(&IntPoint::new(1, -1)), Ok([1.0, -1.0])); + assert_eq!( + adapter.try_int_to_float(&IntPoint::new(2, 0)), + Err(FloatPointAdapterRangeError::PointOutOfRange) + ); + } + #[test] fn test_try_with_scale() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::try_with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0) - .unwrap(); + let adapter = FloatPointAdapter::<[f64; 2], i32>::try_with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ) + .unwrap(); assert_eq!(adapter.dir_scale(), 10.0); assert_eq!(adapter.inv_scale(), 0.1); @@ -526,7 +817,7 @@ mod tests { #[test] fn test_try_with_scale_and_coordinate_bits() { - let rect = FloatRect::new(-4.0, 4.0, -1.0, 1.0); + let rect = FloatRect::new(-4.0, 4.0, -1.0, 1.0).unwrap(); let adapter = FloatPointAdapter::<[f64; 2], i32>::try_with_scale_and_coordinate_bits(rect.clone(), 128.0, 10) .unwrap(); @@ -542,7 +833,7 @@ mod tests { #[test] fn test_try_with_scale_errors() { - let rect = FloatRect::new(-1.0, 1.0, -1.0, 1.0); + let rect = FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(); assert_eq!( FloatPointAdapter::<[f64; 2], i32>::try_with_scale(rect.clone(), 0.0) @@ -577,9 +868,11 @@ mod tests { #[test] fn test_degenerate_try_with_scale_keeps_requested_scale() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::try_with_scale(FloatRect::new(1.0, 1.0, -2.0, -2.0), 1000.0) - .unwrap(); + let adapter = FloatPointAdapter::<[f64; 2], i32>::try_with_scale( + FloatRect::new(1.0, 1.0, -2.0, -2.0).unwrap(), + 1000.0, + ) + .unwrap(); assert_eq!(adapter.dir_scale(), 1000.0); assert_eq!(adapter.float_to_int(&[1.0, -2.0]), IntPoint::ZERO); @@ -587,8 +880,10 @@ mod tests { #[test] fn test_try_float_to_int_range() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ); assert_eq!( adapter.try_float_to_int(&[1.0, -1.0]).unwrap(), @@ -602,23 +897,29 @@ mod tests { #[test] fn test_try_float_to_int_not_finite() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); - - assert_eq!( - adapter.try_float_to_int(&[f64::NAN, 0.0]).err().unwrap(), - FloatPointAdapterRangeError::PointOutOfRange - ); - assert_eq!( - adapter.try_float_to_int(&[0.0, f64::INFINITY]).err().unwrap(), - FloatPointAdapterRangeError::PointOutOfRange - ); + fn check(nan: F, infinity: F) { + let limit = F::MAX_COORDINATE; + let rect = FloatRect::new(-limit, limit, -limit, limit).unwrap(); + let adapter = FloatPointAdapter::<[F; 2], i32>::with_coordinate_bits(rect, 0); + for value in [nan, infinity, -infinity] { + for point in [[value, F::ZERO], [F::ZERO, value]] { + assert_eq!( + adapter.try_float_to_int(&point), + Err(FloatPointAdapterRangeError::PointOutOfRange), + ); + } + } + } + check(f32::NAN, f32::INFINITY); + check(f64::NAN, f64::INFINITY); } #[test] fn test_try_int_to_float_range() { - let adapter = - FloatPointAdapter::<[f64; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0), 10.0); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_scale( + FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), + 10.0, + ); assert_eq!( adapter.try_int_to_float(&IntPoint::new(10, -10)).unwrap(), diff --git a/src/float/compatible.rs b/src/float/compatible.rs index 1889060..cadb8d9 100644 --- a/src/float/compatible.rs +++ b/src/float/compatible.rs @@ -10,6 +10,13 @@ pub trait FloatPointCompatible: Copy { fn is_finite(&self) -> bool { self.x().is_finite() && self.y().is_finite() } + + /// Whether both coordinates satisfy the inclusive [`crate::float`] input + /// limits. NaN and infinity are rejected. + #[inline(always)] + fn is_in_safe_range(&self) -> bool { + self.x().is_in_safe_range() && self.y().is_in_safe_range() + } } impl FloatPointCompatible for [T; 2] { diff --git a/src/float/mod.rs b/src/float/mod.rs index a83892d..0bd5d35 100644 --- a/src/float/mod.rs +++ b/src/float/mod.rs @@ -1,3 +1,36 @@ +//! Floating-point geometry with a conservative coordinate range. +//! +//! Supported input coordinates and rectangle bounds are finite and satisfy: +//! +//! | Scalar | Maximum absolute coordinate | +//! | --- | --- | +//! | `f32` | `2^60` (approximately `1.15e18`) | +//! | `f64` | `2^500` (approximately `3.27e150`) | +//! +//! These inclusive bounds leave room for point differences, their dot and cross +//! products, squared lengths, and midpoints without floating-point overflow. +//! They do not guarantee exact arithmetic: rounding, cancellation, and underflow +//! still follow the scalar's floating-point behavior. Arbitrary scaling and +//! repeated operations require their own range analysis. +//! +//! Normalization additionally requires a positive, finite, normal squared +//! length (at least `f32::MIN_POSITIVE` or `f64::MIN_POSITIVE`). A nonzero vector +//! alone is insufficient: squaring tiny components can underflow. +//! +//! [`rect::FloatRect`] constructors and checked mutators validate coordinates +//! and bound ordering in every build, returning [`rect::FloatRectError`] for +//! invalid input. Public field writes and [`rect::FloatRect::unsafe_add_point`] +//! remain the caller's responsibility; the latter checks only in debug builds. +//! +//! [`point::FloatPoint::new`] and [`point::FloatPoint::from_point`] assert input +//! coordinate limits in debug builds. Both normalization implementations also +//! assert the squared-length precondition in debug builds. Arithmetic results +//! are not checked against the input limits, since point differences can be +//! larger. Coordinates are never clamped; other preconditions remain the +//! caller's responsibility. Behavior outside this contract is not guaranteed. +//! The coordinate limits do not restrict scalar constants, adapter scales, or +//! every use of [`number::FloatNumber`]. + pub mod compatible; pub mod number; diff --git a/src/float/number.rs b/src/float/number.rs index 6320397..d242f77 100644 --- a/src/float/number.rs +++ b/src/float/number.rs @@ -3,6 +3,10 @@ use crate::int::number::wide_int::WideIntNumber; use core::fmt::Display; use core::ops::{Add, Div, Mul, Neg, Sub}; +/// Scalar operations used by floating-point geometry. +/// +/// The [`crate::float`] coordinate limits apply to geometry inputs, not to this +/// trait's full scalar range or to adapter scales. pub trait FloatNumber where Self: Copy @@ -16,7 +20,15 @@ where { const MAX: Self; const MIN: Self; + /// Inclusive maximum absolute input coordinate for floating-point geometry. + /// This does not limit adapter scales or the scalar's representable range. + const MAX_COORDINATE: Self; + /// Smallest positive normal scalar, used as the minimum squared length + /// supported by floating-point normalization. + const MIN_POSITIVE: Self; const BITS: u32; + /// One greater than the largest exponent of a finite power of two. + const MAX_EXP: i32; const ZERO: Self; const ONE: Self; const TWO: Self; @@ -42,8 +54,19 @@ where fn sin_cos(self) -> (Self, Self); fn acos(self) -> Self; fn asin(self) -> Self; + /// Returns the signed angle in radians for the vector (x, self). + fn atan2(self, x: Self) -> Self; fn signum(self) -> Self; fn is_finite(self) -> bool; + + /// Whether this value is within the inclusive geometry input coordinate + /// limits defined by [`Self::MAX_COORDINATE`]. NaN and infinity are rejected. + /// This check does not apply to adapter scales or intermediate results. + #[inline(always)] + fn is_in_safe_range(self) -> bool { + self.abs() <= Self::MAX_COORDINATE + } + // Truncating casts. fn to_i16(self) -> i16; fn to_i32(self) -> i32; @@ -63,7 +86,10 @@ where impl FloatNumber for f32 { const MAX: Self = f32::MAX; const MIN: Self = f32::MIN; + const MAX_COORDINATE: Self = f32::from_bits((127 + 60) << 23); + const MIN_POSITIVE: Self = f32::MIN_POSITIVE; const BITS: u32 = 32; + const MAX_EXP: i32 = f32::MAX_EXP; const ZERO: Self = 0.0; const ONE: Self = 1.0; const TWO: Self = 2.0; @@ -148,6 +174,11 @@ impl FloatNumber for f32 { libm::asinf(self) } + #[inline(always)] + fn atan2(self, x: Self) -> Self { + libm::atan2f(self, x) + } + #[inline(always)] fn signum(self) -> Self { self.signum() @@ -192,34 +223,37 @@ impl FloatNumber for f32 { // Rounding casts. #[inline(always)] fn to_round_i16(self) -> i16 { - (self + Self::HALF.copysign(self)) as i16 + libm::roundf(self) as i16 } #[inline(always)] fn to_round_i32(self) -> i32 { - (self + Self::HALF.copysign(self)) as i32 + libm::roundf(self) as i32 } #[inline(always)] fn to_round_i64(self) -> i64 { - (self + Self::HALF.copysign(self)) as i64 + libm::roundf(self) as i64 } #[inline(always)] fn to_round_i128(self) -> i128 { - (self + Self::HALF.copysign(self)) as i128 + libm::roundf(self) as i128 } #[inline(always)] fn to_round_usize(self) -> usize { - (self + Self::HALF) as usize + libm::roundf(self) as usize } } impl FloatNumber for f64 { const MAX: Self = f64::MAX; const MIN: Self = f64::MIN; + const MAX_COORDINATE: Self = f64::from_bits((1023 + 500) << 52); + const MIN_POSITIVE: Self = f64::MIN_POSITIVE; const BITS: u32 = 64; + const MAX_EXP: i32 = f64::MAX_EXP; const ZERO: Self = 0.0; const ONE: Self = 1.0; const TWO: Self = 2.0; @@ -303,6 +337,11 @@ impl FloatNumber for f64 { libm::asin(self) } + #[inline(always)] + fn atan2(self, x: Self) -> Self { + libm::atan2(self, x) + } + #[inline(always)] fn signum(self) -> Self { self.signum() @@ -347,26 +386,26 @@ impl FloatNumber for f64 { // Rounding casts. #[inline(always)] fn to_round_i16(self) -> i16 { - (self + Self::HALF.copysign(self)) as i16 + libm::round(self) as i16 } #[inline(always)] fn to_round_i32(self) -> i32 { - (self + Self::HALF.copysign(self)) as i32 + libm::round(self) as i32 } #[inline(always)] fn to_round_i64(self) -> i64 { - (self + Self::HALF.copysign(self)) as i64 + libm::round(self) as i64 } #[inline(always)] fn to_round_i128(self) -> i128 { - (self + Self::HALF.copysign(self)) as i128 + libm::round(self) as i128 } #[inline(always)] fn to_round_usize(self) -> usize { - (self + Self::HALF) as usize + libm::round(self) as usize } } diff --git a/src/float/point.rs b/src/float/point.rs index bf9abec..3e5a78f 100644 --- a/src/float/point.rs +++ b/src/float/point.rs @@ -3,6 +3,13 @@ use crate::float::number::FloatNumber; use core::fmt; use core::ops::{Add, AddAssign, Mul, Neg, Sub}; +/// A floating-point point or vector subject to the [`crate::float`] +/// coordinate-range contract. Input coordinates must be finite, with absolute +/// values at most `2^60` for `f32` or `2^500` for `f64`. +/// [`Self::new`] and [`Self::from_point`] check these limits in debug builds. +/// Arithmetic results and construction through public fields or +/// [`FloatPointCompatible::from_xy`] are not range-checked: intermediate +/// vectors such as point differences may exceed the input coordinate limits. #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] #[derive(Debug, Clone, Copy)] pub struct FloatPoint { @@ -19,14 +26,17 @@ impl FloatPoint { } } + /// Creates a point, asserting the coordinate limits in debug builds. #[inline(always)] pub fn new(x: T, y: T) -> Self { - Self { x, y } + let point = Self { x, y }; + debug_assert!(point.is_in_safe_range(), "FloatPoint coordinates out of range"); + point } #[inline(always)] pub fn from_point>(p: P) -> Self { - Self { x: p.x(), y: p.y() } + Self::new(p.x(), p.y()) } #[inline(always)] @@ -54,15 +64,29 @@ impl FloatPoint { self.sqr_length().sqrt() } + /// Returns an approximately unit-length vector. + /// + /// Requires a positive, finite, normal [`Self::sqr_length`]: at least + /// `f32::MIN_POSITIVE` or `f64::MIN_POSITIVE`. Rescale vectors with smaller + /// squared lengths before normalizing. Zero vectors and vectors outside + /// this contract are unsupported and may produce zero, infinity, or NaN. + /// Debug builds assert these preconditions. #[inline(always)] pub fn normalize(&self) -> Self { - let l = self.length(); + let sqr_length = self.sqr_length(); + debug_assert!( + sqr_length.is_finite() && sqr_length >= T::MIN_POSITIVE, + "Normalization requires a positive finite normal squared length" + ); + let l = sqr_length.sqrt(); Self { x: self.x / l, y: self.y / l, } } + /// Returns the midpoint using ordinary floating-point rounding. + /// Both points must satisfy the [`crate::float`] coordinate-range contract. #[inline(always)] pub fn midpoint(self, other: Self) -> Self { (self + other) * T::HALF diff --git a/src/float/rect.rs b/src/float/rect.rs index 34a9fdd..b07da8c 100644 --- a/src/float/rect.rs +++ b/src/float/rect.rs @@ -2,6 +2,32 @@ use crate::float::compatible::FloatPointCompatible; use crate::float::number::FloatNumber; use core::fmt; +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum FloatRectError { + /// A coordinate is non-finite or exceeds the supported absolute limit. + CoordinatesOutOfRange, + /// A minimum bound is greater than its maximum. + InvalidBounds, +} + +impl fmt::Display for FloatRectError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::CoordinatesOutOfRange => { + f.write_str("A coordinate is non-finite or exceeds the supported absolute limit") + } + Self::InvalidBounds => f.write_str("A minimum bound is greater than its maximum"), + } + } +} + +impl core::error::Error for FloatRectError {} + +/// Floating-point bounds under the [`crate::float`] coordinate-range contract. +/// Bounds must be finite, ordered (`min <= max`), and have absolute values at +/// most `2^60` for `f32` or `2^500` for `f64`. Constructors and checked mutators +/// validate them in every build. Public field writes and [`Self::unsafe_add_point`] +/// remain the caller's responsibility. #[derive(Debug, Copy, Clone)] pub struct FloatRect { pub min_x: T, @@ -11,6 +37,29 @@ pub struct FloatRect { } impl FloatRect { + /// Whether the bounds are ordered, finite, and within the inclusive + /// [`FloatNumber::MAX_COORDINATE`] limits on both axes. + #[inline(always)] + pub fn is_in_safe_range(&self) -> bool { + self.validate().is_ok() + } + + /// Validates public bounds, including rectangles assembled through fields. + #[inline] + pub fn validate(&self) -> Result<(), FloatRectError> { + if !(self.min_x.is_in_safe_range() + && self.max_x.is_in_safe_range() + && self.min_y.is_in_safe_range() + && self.max_y.is_in_safe_range()) + { + return Err(FloatRectError::CoordinatesOutOfRange); + } + if self.min_x > self.max_x || self.min_y > self.max_y { + return Err(FloatRectError::InvalidBounds); + } + Ok(()) + } + #[inline(always)] pub fn width(&self) -> T { self.max_x - self.min_x @@ -21,14 +70,17 @@ impl FloatRect { self.max_y - self.min_y } + /// Creates bounds, rejecting out-of-range coordinates and reversed axes. #[inline(always)] - pub fn new(min_x: T, max_x: T, min_y: T, max_y: T) -> Self { - Self { + pub fn new(min_x: T, max_x: T, min_y: T, max_y: T) -> Result { + let rect = Self { min_x, max_x, min_y, max_y, - } + }; + rect.validate()?; + Ok(rect) } #[inline(always)] @@ -42,48 +94,47 @@ impl FloatRect { } #[inline] - pub fn with_point>(point: P) -> Self { - Self { - min_x: point.x(), - max_x: point.x(), - min_y: point.y(), - max_y: point.y(), - } + pub fn with_point>(point: P) -> Result { + Self::new(point.x(), point.x(), point.y(), point.y()) } + /// Returns `Ok(None)` for empty input, or an error for any invalid point. #[inline] - pub fn with_points

(points: &[P]) -> Option + pub fn with_points

(points: &[P]) -> Result, FloatRectError> where P: FloatPointCompatible, { Self::with_iter(points.iter()) } + /// Returns `Ok(None)` for empty input, or an error for any invalid point. #[inline] - pub fn with_iter<'a, I, P>(iter: I) -> Option + pub fn with_iter<'a, I, P>(iter: I) -> Result, FloatRectError> where I: Iterator, P: FloatPointCompatible + 'a, T: FloatNumber, { let mut iter = iter; - let first_point = iter.next()?; - let mut rect = Self { - min_x: first_point.x(), - max_x: first_point.x(), - min_y: first_point.y(), - max_y: first_point.y(), + let Some(first_point) = iter.next() else { + return Ok(None); }; + let mut rect = Self::with_point(*first_point)?; for p in iter { + if !p.is_in_safe_range() { + return Err(FloatRectError::CoordinatesOutOfRange); + } rect.unsafe_add_point(p); } - Some(rect) + Ok(Some(rect)) } #[inline] - pub fn with_rects(rect_0: Self, rect_1: Self) -> Self { + pub fn with_rects(rect_0: Self, rect_1: Self) -> Result { + rect_0.validate()?; + rect_1.validate()?; let min_x = rect_0.min_x.min(rect_1.min_x); let max_x = rect_0.max_x.max(rect_1.max_x); let min_y = rect_0.min_y.min(rect_1.min_y); @@ -92,42 +143,56 @@ impl FloatRect { } #[inline] - pub fn with_optional_rects(rect_0: Option, rect_1: Option) -> Option { + pub fn with_optional_rects( + rect_0: Option, + rect_1: Option, + ) -> Result, FloatRectError> { match (rect_0, rect_1) { - (Some(r0), Some(r1)) => Some(Self::with_rects(r0, r1)), - (Some(r0), None) => Some(r0), - (None, Some(r1)) => Some(r1), - (None, None) => None, + (Some(r0), Some(r1)) => Self::with_rects(r0, r1).map(Some), + (Some(rect), None) | (None, Some(rect)) => { + rect.validate()?; + Ok(Some(rect)) + } + (None, None) => Ok(None), } } + /// Expands the bounds to include a valid point. On error, leaves `self` unchanged. #[inline] - pub fn add_point>(&mut self, point: &P) { - if self.min_x > point.x() { - self.min_x = point.x() - } - if self.max_x < point.x() { - self.max_x = point.x() - } - - if self.min_y > point.y() { - self.min_y = point.y() - } - if self.max_y < point.y() { - self.max_y = point.y() + pub fn add_point>( + &mut self, + point: &P, + ) -> Result<(), FloatRectError> { + self.validate()?; + if !point.is_in_safe_range() { + return Err(FloatRectError::CoordinatesOutOfRange); } + self.unsafe_add_point(point); + Ok(()) } + /// Expands each side by `offset`; negative values shrink the bounds. + /// Returns an error without changing `self` if the resulting bounds are invalid. #[inline] - pub fn add_offset(&mut self, offset: T) { - self.max_x = self.max_x + offset; - self.max_y = self.max_y + offset; - self.min_x = self.min_x - offset; - self.min_y = self.min_y - offset; + pub fn add_offset(&mut self, offset: T) -> Result<(), FloatRectError> { + self.validate()?; + let rect = Self::new( + self.min_x - offset, + self.max_x + offset, + self.min_y - offset, + self.max_y + offset, + )?; + *self = rect; + Ok(()) } + /// Expands valid bounds with a point known to satisfy the coordinate limits. + /// Only debug builds check these preconditions; use [`Self::add_point`] for + /// validation in every build. #[inline] pub fn unsafe_add_point>(&mut self, point: &P) { + debug_assert!(self.is_in_safe_range(), "FloatRect bounds out of range"); + debug_assert!(point.is_in_safe_range(), "FloatPoint coordinates out of range"); if self.min_x > point.x() { self.min_x = point.x() } else if self.max_x < point.x() { @@ -142,10 +207,16 @@ impl FloatRect { } #[inline] - pub fn optional_add_point>(rect: &mut Option, point: &P) { + pub fn optional_add_point>( + rect: &mut Option, + point: &P, + ) -> Result<(), FloatRectError> { match rect { - Some(rect) => rect.unsafe_add_point(point), - None => *rect = Some(FloatRect::with_point(*point)), + Some(rect) => rect.add_point(point), + None => { + *rect = Some(FloatRect::with_point(*point)?); + Ok(()) + } } } @@ -192,7 +263,7 @@ mod tests { fn test_0() { let points = [[-2.0, -4.0], [-2.0, 3.0], [5.0, 3.0], [5.0, -4.0]]; - let rect: FloatRect = FloatRect::with_iter(points.iter()).unwrap(); + let rect: FloatRect = FloatRect::with_iter(points.iter()).unwrap().unwrap(); assert_eq!((rect.max_x - 5.0).abs() < 0.000_0001, true); assert_eq!((rect.min_x + 2.0).abs() < 0.000_0001, true); @@ -202,30 +273,34 @@ mod tests { #[test] fn test_1() { - let r0 = Some(FloatRect::new(-2.0, 2.0, -2.0, 2.0)); - let r1 = Some(FloatRect::new(-4.0, 4.0, -4.0, 4.0)); - let rr = FloatRect::with_optional_rects(r0, r1).unwrap(); + let r0 = Some(FloatRect::new(-2.0, 2.0, -2.0, 2.0).unwrap()); + let r1 = Some(FloatRect::new(-4.0, 4.0, -4.0, 4.0).unwrap()); + let rr = FloatRect::with_optional_rects(r0, r1).unwrap().unwrap(); assert_eq!(-4.0, rr.min_x); assert_eq!(-4.0, rr.min_y); assert_eq!(4.0, rr.max_x); assert_eq!(4.0, rr.max_y); - assert!(FloatRect::with_optional_rects(r0, None).is_some()); - assert!(FloatRect::with_optional_rects(None, r1).is_some()); - assert!(FloatRect::::with_optional_rects(None, None).is_none()); + assert!(FloatRect::with_optional_rects(r0, None).unwrap().is_some()); + assert!(FloatRect::with_optional_rects(None, r1).unwrap().is_some()); + assert!( + FloatRect::::with_optional_rects(None, None) + .unwrap() + .is_none() + ); } #[test] fn is_intersect_with_padding() { - let rect = FloatRect::new(0.0, 10.0, 0.0, 10.0); - - assert!(rect.is_intersect_with_padding(&FloatRect::new(11.0, 12.0, 0.0, 10.0), 2.0)); - assert!(rect.is_intersect_with_padding(&FloatRect::new(-3.0, -1.0, 0.0, 10.0), 2.0)); - assert!(rect.is_intersect_with_padding(&FloatRect::new(9.0, 11.0, 9.0, 11.0), 0.0)); - assert!(!rect.is_intersect_with_padding(&FloatRect::new(10.0, 12.0, 0.0, 10.0), 0.0)); - assert!(rect.is_intersect_with_padding(&FloatRect::new(10.5, 12.0, 0.0, 10.0), 1.0)); - assert!(!rect.is_intersect_with_padding(&FloatRect::new(11.0, 12.0, 0.0, 10.0), 1.0)); - assert!(!rect.is_intersect_with_padding(&FloatRect::new(12.0, 13.0, 0.0, 10.0), 1.0)); - assert!(!rect.is_intersect_with_padding(&FloatRect::new(9.0, 11.0, 11.0, 12.0), 1.0)); + let rect = FloatRect::new(0.0, 10.0, 0.0, 10.0).unwrap(); + + assert!(rect.is_intersect_with_padding(&FloatRect::new(11.0, 12.0, 0.0, 10.0).unwrap(), 2.0)); + assert!(rect.is_intersect_with_padding(&FloatRect::new(-3.0, -1.0, 0.0, 10.0).unwrap(), 2.0)); + assert!(rect.is_intersect_with_padding(&FloatRect::new(9.0, 11.0, 9.0, 11.0).unwrap(), 0.0)); + assert!(!rect.is_intersect_with_padding(&FloatRect::new(10.0, 12.0, 0.0, 10.0).unwrap(), 0.0)); + assert!(rect.is_intersect_with_padding(&FloatRect::new(10.5, 12.0, 0.0, 10.0).unwrap(), 1.0)); + assert!(!rect.is_intersect_with_padding(&FloatRect::new(11.0, 12.0, 0.0, 10.0).unwrap(), 1.0)); + assert!(!rect.is_intersect_with_padding(&FloatRect::new(12.0, 13.0, 0.0, 10.0).unwrap(), 1.0)); + assert!(!rect.is_intersect_with_padding(&FloatRect::new(9.0, 11.0, 11.0, 12.0).unwrap(), 1.0)); } } diff --git a/src/float/vector.rs b/src/float/vector.rs index 6a2ba37..331bafc 100644 --- a/src/float/vector.rs +++ b/src/float/vector.rs @@ -1,6 +1,9 @@ use crate::float::compatible::FloatPointCompatible; use crate::float::number::FloatNumber; +/// Arithmetic for compatible points under the [`crate::float`] +/// coordinate-range contract. Normalization checks its squared length in debug +/// builds; other arithmetic does not perform range checks. pub struct FloatPointMath

{ _phantom: core::marker::PhantomData

, } @@ -31,9 +34,21 @@ impl FloatPointMath

{ Self::sqr_length(p).sqrt() } + /// Returns an approximately unit-length vector. + /// + /// Requires a positive, finite, normal [`Self::sqr_length`]: at least + /// `f32::MIN_POSITIVE` or `f64::MIN_POSITIVE`. Rescale vectors with smaller + /// squared lengths before normalizing. Zero vectors and vectors outside + /// this contract are unsupported and may produce zero, infinity, or NaN. + /// Debug builds assert these preconditions. #[inline(always)] pub fn normalize(p: &P) -> P { - let inv_len = P::Scalar::ONE / Self::length(p); + let sqr_length = Self::sqr_length(p); + debug_assert!( + sqr_length.is_finite() && sqr_length >= P::Scalar::MIN_POSITIVE, + "Normalization requires a positive finite normal squared length" + ); + let inv_len = P::Scalar::ONE / sqr_length.sqrt(); P::from_xy(p.x() * inv_len, p.y() * inv_len) } diff --git a/src/int/angle/angle.rs b/src/int/angle/angle.rs new file mode 100644 index 0000000..4f1310c --- /dev/null +++ b/src/int/angle/angle.rs @@ -0,0 +1,319 @@ +use super::cordic; +use crate::float::number::FloatNumber; +use crate::int::number::{int::IntNumber, uint::UIntNumber, wide_int::WideIntNumber}; +use crate::int::unit_vector::UnitIntVector; +use core::f64::consts::TAU; + +const ANGLE_SCALE: f64 = (1_u64 << 32) as f64; + +/// Counterclockwise binary angle: one turn is `2^32` units. +/// +/// This resolution is independent of an arc's minimum subdivision step. +/// Radians can be converted with [`Self::from_radians`]. Arc traversal policies +/// live outside this type. +/// +/// Example of the future consumer's loop (the reusable buffer and subdivision +/// policy live outside this library): +/// ``` +/// use i_float::int::{angle::{Angle, Rotation}, vector::IntVector}; +/// let from = IntVector::::new(1, 0).fast_normalize().unwrap(); +/// let to = IntVector::::new(0, -1).fast_normalize().unwrap(); +/// let clockwise = false; // a 270-degree counterclockwise arc +/// let max_step = (1u32 << 27).clamp(1 << 22, 1 << 29); +/// let sweep = if clockwise { Angle::between(to, from) } +/// else { Angle::between(from, to) }; +/// let mut directions = Vec::new(); // retain this buffer between builds +/// directions.clear(); +/// if sweep.bits() != 0 { +/// // Include vectoring uncertainty before rounding the count upward. +/// let upper = sweep.bits() as u64 + Angle::MAX_ERROR as u64; +/// let mut segments = upper.div_ceil(max_step as u64); +/// // Budget the final gap too: see Rotation's error contract. For fresh +/// // i32/i64 normalized inputs, 8 binary units per step cover matrix +/// // error, component truncation, and the integer step division below. +/// while upper + 8 * segments * segments > segments * max_step as u64 { +/// segments += 1; +/// } +/// if segments > 1 { +/// let step = (sweep.bits() as u64 / segments) as u32; +/// let rotation = Rotation::new(Angle::from_bits( +/// if clockwise { step.wrapping_neg() } else { step })); +/// let mut v = from; +/// for _ in 1..segments { +/// v = rotation.apply(v); +/// directions.push(v); +/// } +/// } +/// } +/// // Endpoints are excluded; the consumer uses its original exact contacts. +/// assert!(!directions.is_empty()); +/// ``` +/// This example's error budget applies to fresh `i32`/`i64` normalized inputs, +/// and at most 1030 steps. It does not apply to `i16` or already heavily +/// contracted directions. Subdivision and its policy remain consumer code. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)] +pub struct Angle(u32); + +impl Angle { + pub const ZERO: Self = Self(0); + pub const QUARTER_TURN: Self = Self(1 << 30); + pub const HALF_TURN: Self = Self(1 << 31); + pub const THREE_QUARTER_TURN: Self = Self(3 << 30); + + /// Conservative absolute error of [`Self::between`], in binary angle units + /// (about 4.7e-8 radians), relative to the *stored* input directions. Also applies to [`Self::atan2`]. + /// Input normalization error is additional. Exact axes have zero error. + pub const MAX_ERROR: u32 = 32; + + #[inline] + pub const fn from_bits(bits: u32) -> Self { + Self(bits) + } + + /// Converts radians modulo one turn to binary angle units. + /// Rounds to the nearest integer, with ties away from zero. + /// Returns `None` for NaN or infinity. + #[inline] + pub fn from_radians(radians: T) -> Option { + let radians = radians.to_f64(); + if !radians.is_finite() { + return None; + } + + let scaled = (radians % TAU) * (ANGLE_SCALE / TAU); + let truncated = scaled as i64; + let fraction = scaled - truncated as f64; + let rounded = if fraction >= 0.5 { + truncated + 1 + } else if fraction <= -0.5 { + truncated - 1 + } else { + truncated + }; + + Some(Self(rounded as u32)) + } + + #[inline] + pub const fn bits(self) -> u32 { + self.0 + } + + /// Returns the angle in radians in `[0, 2π)` before rounding to `T`. + /// For `f32`, values near the upper boundary may round to `2π`. + #[inline(always)] + pub fn to_radians(self) -> T { + T::from_float(self.0 as f64 * (TAU / ANGLE_SCALE)) + } + + /// Returns the equivalent angle in radians in `[-π, π)` before rounding to `T`. + /// For `f32`, values near the upper boundary may round to `π`. + #[inline(always)] + pub fn to_signed_radians(self) -> T { + T::from_float(self.0 as i32 as f64 * (TAU / ANGLE_SCALE)) + } + + /// Adds a signed angle delta. Overflow performs the required full-turn wrap. + #[inline(always)] + pub const fn wrapping_add(self, delta: AngleDelta) -> Self { + Self(self.0.wrapping_add(delta.0 as u32)) + } + + /// Returns the shortest signed delta from `self` to `target`. + /// Exactly opposite angles yield a negative half-turn. + #[inline(always)] + pub const fn delta_to(self, target: Self) -> AngleDelta { + AngleDelta(target.0.wrapping_sub(self.0) as i32) + } + + /// Counterclockwise sweep from `from` to `to`, in `[0, one turn)`. + /// Swap the arguments to obtain the clockwise sweep magnitude. + /// + /// Equal rays return zero regardless of their approximate lengths; + /// opposite rays return exactly `2^31`. A nonzero cross product is never + /// rounded to an empty arc, even arbitrarily close to a full turn. + /// Products fit in the associated wide type, including `i128` for `i64`: + /// components are approximately bounded by S = 2^(I::BITS-2), so products + /// and sums/differences fit with room for float rounding near unit length. + /// No products of cross/dot are formed. + pub fn between(from: UnitIntVector, to: UnitIntVector) -> Self { + let (ax, ay) = (from.x().to_wide(), from.y().to_wide()); + let (bx, by) = (to.x().to_wide(), to.y().to_wide()); + let cross = ax * by - ay * bx; + let dot = ax * bx + ay * by; + Self::atan2(cross, dot).unwrap_or(Self::from_bits(0)) + } + + /// Floating-point alternative to [`Self::between`], using [`Self::atan2_with_float`]. + /// Returns the counterclockwise sweep in `[0, one turn)`; swap the arguments + /// for the clockwise sweep. Equal rays return zero and opposite rays return + /// exactly half a turn. Cross/dot products are computed in wide integer + /// arithmetic before conversion, preserving the side of near-parallel turns. + #[inline] + pub fn between_with_float(from: UnitIntVector, to: UnitIntVector) -> Self { + let (ax, ay) = (from.x().to_wide(), from.y().to_wide()); + let (bx, by) = (to.x().to_wide(), to.y().to_wide()); + let cross = ax * by - ay * bx; + let dot = ax * bx + ay * by; + Self::atan2_with_float(cross, dot).unwrap_or(Self::from_bits(0)) + } + + /// Returns atan2(y, x) as a counterclockwise angle in [0, one turn). + /// Every value of the built-in wide integer types is valid, including MIN. + /// Returns None only for (0, 0); no normalization is required. + pub fn atan2(y: W, x: W) -> Option { + if y == W::ZERO { + return if x == W::ZERO { + None + } else { + Some(Self(if x < W::ZERO { 1 << 31 } else { 0 })) + }; + } + if x == W::ZERO { + return Some(Self(if y > W::ZERO { 1 << 30 } else { 3 << 30 })); + } + let (negative_x, negative_y) = (x < W::ZERO, y < W::ZERO); + let (x, y) = (x.unsigned_abs(), y.unsigned_abs()); + let top = x.max(y).ilog2(); + // Preserve signs before reducing magnitudes; this conversion fits on + // 32-bit hosts too. Unsigned abs also handles i128::MIN. + let input_top = cordic::VECTOR_INPUT_BITS - 1; + let reduce = |v: W::UInt| -> i64 { + let v = if top > input_top { + v >> (top - input_top) + } else { + v << (input_top - top) + }; + (v.to_usize() as i64) << cordic::VECTOR_GUARD_BITS + }; + let mut z = cordic::vectoring(reduce(x), reduce(y), cordic::ITERATIONS); + if negative_x { + z = (1 << 31) - z; + } + let magnitude = z.clamp(1, (1 << 31) - 1) as u32; + Some(Self(if negative_y { + magnitude.wrapping_neg() + } else { + magnitude + })) + } + + /// Floating-point alternative to [`Self::atan2`], with the same integer inputs + /// and binary-angle output. Uses f64 internally, including for i128 inputs. + /// Returns None only for (0, 0); axes are exact. For nonzero y, rounding + /// preserves the open upper/lower half-turn, including near zero and pi. + #[inline] + pub fn atan2_with_float(y: W, x: W) -> Option { + if y == W::ZERO { + return if x == W::ZERO { + None + } else { + Some(Self(if x < W::ZERO { 1 << 31 } else { 0 })) + }; + } + if x == W::ZERO { + return Some(Self(if y > W::ZERO { 1 << 30 } else { 3 << 30 })); + } + let radians = FloatNumber::atan2(y.to_f64().abs(), x.to_f64()); + let magnitude = Self::from_radians(radians)?.bits().clamp(1, (1 << 31) - 1); + Some(Self(if y < W::ZERO { + magnitude.wrapping_neg() + } else { + magnitude + })) + } + + /// Integer value representing one in sin/cos results (Q30). + pub const SIN_COS_SCALE: i32 = 1 << 30; + + /// Returns (sine, cosine) in Q30, using a single CORDIC rotation. + /// Axes are exact; the coefficient pair has length at most [`Self::SIN_COS_SCALE`]. + /// See [`super::Rotation`] for error bounds. The result uses the same scale for + /// every coordinate type. + #[inline] + pub fn sin_cos(self) -> (i32, i32) { + cordic::sin_cos(self.0, cordic::ITERATIONS, cordic::COEFFICIENT_BITS) + } + + /// Returns (sine, cosine) in Q30, using f64 trigonometry internally. + /// Like [`Self::sin_cos`], axes are exact and the coefficient pair has length + /// at most [`Self::SIN_COS_SCALE`]. Components are rounded to the nearest + /// integer, then contracted toward zero if needed to preserve that bound. + #[inline] + pub fn sin_cos_with_float(self) -> (i32, i32) { + let scale = Self::SIN_COS_SCALE; + match self.0 { + 0 => return (0, scale), + 0x4000_0000 => return (scale, 0), + 0x8000_0000 => return (0, -scale), + 0xc000_0000 => return (-scale, 0), + _ => {} + } + let radians = self.0 as f64 * (TAU / 4294967296.0); + let (sin, cos) = FloatNumber::sin_cos(radians); + let (mut sin, mut cos) = ( + (sin * scale as f64).to_round_i32(), + (cos * scale as f64).to_round_i32(), + ); + // Check the rounded norm exactly. Reducing the larger component removes + // the most excess squared length per integer unit of adjustment. + while sin as i64 * sin as i64 + cos as i64 * cos as i64 > scale as i64 * scale as i64 { + if sin.abs() >= cos.abs() { + sin -= sin.signum(); + } else { + cos -= cos.signum(); + } + } + (sin, cos) + } + + /// Returns the Q30 sine. Use sin_cos when both components are needed. + #[inline] + pub fn sin(self) -> i32 { + self.sin_cos().0 + } + + /// Returns the Q30 cosine. Use sin_cos when both components are needed. + #[inline] + pub fn cos(self) -> i32 { + self.sin_cos().1 + } +} + +impl core::ops::Add for Angle { + type Output = Angle; + + /// Adds angles modulo one full turn. + #[inline(always)] + fn add(self, rhs: Self) -> Self::Output { + Self::from_bits(self.bits().wrapping_add(rhs.bits())) + } +} + +/// Signed binary angle difference in the range of one half-turn. +/// +/// The raw range maps to `[-π, π)`, with the same approximately `1.46e-9 rad` +/// resolution as [`Angle`]. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)] +pub struct AngleDelta(i32); + +impl AngleDelta { + pub const ZERO: Self = Self(0); + + #[inline(always)] + pub const fn from_raw(raw: i32) -> Self { + Self(raw) + } + + #[inline(always)] + pub const fn raw(self) -> i32 { + self.0 + } + + /// Returns the signed delta in radians in `[-π, π)` before rounding to `T`. + /// For `f32`, values near the upper boundary may round to `π`. + #[inline(always)] + pub fn to_radians(self) -> T { + T::from_float(self.0 as f64 * (TAU / ANGLE_SCALE)) + } +} diff --git a/src/int/angle/cordic.rs b/src/int/angle/cordic.rs new file mode 100644 index 0000000..b71490e --- /dev/null +++ b/src/int/angle/cordic.rs @@ -0,0 +1,178 @@ +//! Shared integer CORDIC kernel. Public coefficients stay Q30; approximate +//! Rotation construction selects iterations from a relative angular budget. +// Accuracy budget, independent of coordinate storage: n rotations leave a +// residual below 2^-(n-1) radians. At n=29, 1024 steps at radius 2^16 +// contribute <0.25 coordinate units from that residual alone (n=28: <0.5). +// Q30 coefficient truncation and application rounding add separate errors. +pub(crate) const ITERATIONS: usize = 29; +// Vectoring retains 30 significant input bits, then lifts the top bit to 58. +// These widths are independent of the iteration count and fit 32-bit usize +// conversion and i64 CORDIC growth, respectively. +pub(crate) const VECTOR_INPUT_BITS: u32 = 30; +pub(crate) const VECTOR_GUARD_BITS: u32 = 59 - VECTOR_INPUT_BITS; +pub(crate) const COEFFICIENT_BITS: u32 = 30; +// atan(2^-i), 48-bit binary turns (16 guard bits beyond the public angle). +const ATAN: [i64; ITERATIONS] = [ + 35184372088832, + 20770547670515, + 10974586953444, + 5570871696862, + 2796246208089, + 1399486241028, + 699913886760, + 349978300884, + 174991820497, + 87496244017, + 43748163730, + 21874087080, + 10937044192, + 5468522177, + 2734261099, + 1367130551, + 683565276, + 341782638, + 170891319, + 85445659, + 42722830, + 21361415, + 10680707, + 5340354, + 2670177, + 1335088, + 667544, + 333772, + 166886, +]; + +// floor(2^60 / product(sqrt(1 + 2^-2i), i=0..n-1)) - 128. +// A gain for each iteration count avoids conflating early termination with +// using the wrong gain. The inward guard dominates shift rounding (<68 units). +const GAINS: [i64; ITERATIONS] = [ + 815238614083298760, + 729171583589189357, + 707400343138147019, + 701937710475640438, + 700570741874588230, + 700228916656934686, + 700143455142409185, + 700122089437857531, + 700116747991345093, + 700115412628443506, + 700115078787638515, + 700114995327432293, + 700114974462380426, + 700114969246117440, + 700114967942051693, + 700114967616035256, + 700114967534531146, + 700114967514155119, + 700114967509061112, + 700114967507787611, + 700114967507469235, + 700114967507389641, + 700114967507369743, + 700114967507364768, + 700114967507363525, + 700114967507363214, + 700114967507363136, + 700114967507363116, + 700114967507363112, +]; + +/// Smallest iteration count covering |signed angle| / 2^precision, down to +/// the full-precision floor. The rounding budget covers coefficient and +/// table rounding in the selected output format. Any exponent is valid, including >= 32 (full precision). +#[inline] +pub(crate) fn rotation_iterations(bits: u32, precision: u32, rounding_budget: u32) -> usize { + let magnitude = (bits as i32).unsigned_abs(); + let budget = magnitude.checked_shr(precision).unwrap_or(0); + if budget <= rounding_budget { + return ITERATIONS; + } + let residual_budget = budget - rounding_budget; + if residual_budget >= 683_565_276 { + return 1; + } + // The residual bound is ceil((2^32/tau) / 2^(n-1)). Start from the next + // power of two and check one fewer iteration: no division or float log. + let mut n = 31 - residual_budget.ilog2(); + let shift = n - 2; + let previous_bound = (683_565_276u32 + (1 << shift) - 1) >> shift; + if previous_bound <= residual_budget { + n -= 1; + } + (n as usize).min(ITERATIONS) +} + +/// First-quadrant vectoring, positive coordinates scaled to top bit 58. +/// The caller handles exact axes and folds the result into the original quadrant. +#[inline] +pub(crate) fn vectoring(mut x: i64, mut y: i64, iterations: usize) -> i64 { + debug_assert!((1..=ITERATIONS).contains(&iterations)); + let mut z = 0i64; + for (shift, angle) in ATAN[..iterations].iter().enumerate() { + let old_x = x; + if y > 0 { + x += y >> shift; + y -= old_x >> shift; + z += angle; + } else if y < 0 { + x -= y >> shift; + y += old_x >> shift; + z -= angle; + } + } + (z + (1 << 15)) >> 16 +} + +#[inline] +pub(crate) fn sin_cos(bits: u32, iterations: usize, coefficient_bits: u32) -> (i32, i32) { + let (sin, cos, _) = rotation(bits, iterations, coefficient_bits); + (sin, cos) +} + +/// Returns sine, cosine and achieved binary angle from the final residual. +/// The angle precedes coefficient quantization; the Q30 angular discrepancy, +/// including rounding the angle to 32 bits, is bounded by two binary units. +#[inline] +pub(crate) fn rotation(bits: u32, iterations: usize, coefficient_bits: u32) -> (i32, i32, u32) { + debug_assert!((1..=ITERATIONS).contains(&iterations)); + debug_assert!((14..=COEFFICIENT_BITS).contains(&coefficient_bits)); + let one = 1 << coefficient_bits; + match bits { + 0 => return (0, one, bits), + 0x4000_0000 => return (one, 0, bits), + 0x8000_0000 => return (0, -one, bits), + 0xc000_0000 => return (-one, 0, bits), + _ => {} + } + let mut z = (bits as i32 as i64) << 16; + let mut sign = 1; + if z > 1 << 46 { + z -= 1 << 47; + sign = -1; + } else if z < -(1 << 46) { + z += 1 << 47; + sign = -1; + } + let (mut x, mut y) = (GAINS[iterations - 1], 0i64); + for (shift, angle) in ATAN[..iterations].iter().enumerate() { + let old_x = x; + if z >= 0 { + x -= y >> shift; + y += old_x >> shift; + z -= angle; + } else { + x += y >> shift; + y -= old_x >> shift; + z += angle; + } + } + let scale = 1i64 << (60 - coefficient_bits); + let (sin, cos) = ((sign * (y / scale)) as i32, (sign * (x / scale)) as i32); + debug_assert!(sin as i64 * sin as i64 + cos as i64 * cos as i64 <= 1i64 << (2 * coefficient_bits)); + // Keep the original unsigned angle here: folding above only changes the + // CORDIC convergence interval. Casting back to u32 wraps across zero. + let achieved = ((((bits as i64) << 16) - z + (1 << 15)) >> 16) as u32; + (sin, cos, achieved) +} diff --git a/src/int/angle/mod.rs b/src/int/angle/mod.rs new file mode 100644 index 0000000..155e08d --- /dev/null +++ b/src/int/angle/mod.rs @@ -0,0 +1,16 @@ +//! Integer binary angles and reusable rotation matrices, without allocation. +//! Angle and rotation operations use integer arithmetic by default. +//! [`Angle::from_radians`], [`Angle::sin_cos_with_float`] and [`Angle::atan2_with_float`] +//! use floating-point arithmetic; [`Angle::between_with_float`] combines exact integer +//! products with floating-point angle evaluation. +// The private implementation file mirrors the public type name. +#[allow(clippy::module_inception)] +mod angle; +mod cordic; +mod rotation; + +pub use angle::{Angle, AngleDelta}; +pub use rotation::Rotation; + +#[cfg(test)] +mod tests; diff --git a/src/int/angle/rotation.rs b/src/int/angle/rotation.rs new file mode 100644 index 0000000..e6db5f3 --- /dev/null +++ b/src/int/angle/rotation.rs @@ -0,0 +1,125 @@ +use super::{Angle, cordic}; +use crate::int::number::{fixed_scale::FixedScale, int::IntNumber, wide_int::WideIntNumber}; +use crate::int::unit_vector::UnitIntVector; + +/// A reusable integer rotation matrix with coefficients stored in `I`. +/// Coefficients share the vector scale: Q14 for i16, Q30 for i32, Q62 for i64. +/// +/// Construct once per arc step, then call [`Self::apply`] for each intermediate +/// direction. Coefficients have norm at most one; component truncation toward +/// zero also cannot increase length. Cardinal rotations are exact. +/// +/// Each non-cardinal application can lose up to `sqrt(2)` storage units due +/// to truncation, in addition to coefficient error. Thus repeated rotations +/// accumulate error; in particular, `i16` storage is unsuitable for long fine +/// arcs. `i64` coefficients are lifted from Q30 into Q62 without gaining precision. +/// +/// For [`Self::new`], coefficient error is at most `2^-28 + sqrt(2)*2^-q + 4e-13` in Euclidean +/// norm, relative to the requested angle, with q = min(I::BITS - 2, 30). Application adds at most +/// `sqrt(2)/DENOMINATOR` in component norm. For input length at least 0.99, +/// `i32`/`i64` angular error is less than 5 binary angle units per application; +/// subdivision must also account for step division and accumulated endpoint +/// error (see the example on [`Angle`]). These bounds concern directions, +/// before a consumer rounds scaled coordinates to the integer grid. +#[derive(Debug, Clone, Copy)] +pub struct Rotation { + sin: I, + cos: I, + angle: Angle, +} + +impl Rotation { + /// Builds a counterclockwise matrix. Negate the binary angle with + /// `bits.wrapping_neg()` for clockwise rotation. + pub fn new(angle: Angle) -> Self { + Self::with_iterations(angle, cordic::ITERATIONS) + } + + /// Builds a matrix with an angular error budget of `|angle| / 2^precision`. + /// `precision` is an exponent: 3 allows 1/8 (12.5%), 4 allows 1/16 (6.25%). + /// The magnitude is that of the shortest signed angle in [-half turn, half + /// turn], so clockwise steps use the same precision as counterclockwise. + /// The coefficient format is determined by `I`; this parameter changes only + /// the number of CORDIC iterations, accounting for coefficient quantization. + /// + /// All exponents are valid. Larger exponents request more accuracy, up to + /// the precision of [`Self::new`]. When the requested budget is below that + /// floor, the constructor uses the full iteration count. A conservative error + /// bound in binary units is `max(floor(|angle| / 2^precision), Self::MAX_ERROR)`; + /// exponents >= 32 select full precision. Zero and cardinal angles are exact. + /// + /// This bounds the matrix angle, not accumulated errors from repeated + /// application or integer coordinate rounding. Use [`Self::angle`] to + /// calculate the number of repeated rotations and the final remainder; + /// an approximate matrix may rotate farther than the requested angle. + /// + /// ``` + /// use i_float::int::angle::{Angle, Rotation}; + /// let requested = Angle::from_bits(((1u64 << 32) / 18) as u32); // ~20 degrees + /// let rotation = Rotation::::with_precision(requested, 4); // angle / 16 + /// let achieved = rotation.angle(); + /// assert!(achieved.bits() > 0); + /// let clockwise = Rotation::::with_precision( + /// Angle::from_bits(requested.bits().wrapping_neg()), 4); + /// assert!(clockwise.angle().bits() > 1 << 31); + /// ``` + #[inline] + pub fn with_precision(angle: Angle, precision: u32) -> Self { + Self::with_iterations( + angle, + cordic::rotation_iterations(angle.bits(), precision, Self::ANGLE_MAX_ERROR), + ) + } + + /// Maximum difference between [`Self::angle`] and the stored matrix's actual + /// angle, in 32-bit binary angle units, for either constructor. + /// Cardinal rotations have zero error. Component rounding in [`Self::apply`] + /// adds its own error, dependent on coordinate type and input length. + /// The bound is 131072 units for i16 (about 0.011 degrees), and 2 for i32/i64. + pub const ANGLE_MAX_ERROR: u32 = 2 << cordic::COEFFICIENT_BITS.saturating_sub(FixedScale::::SHIFT); + + /// Conservative matrix-angle error of [`Self::new`], in binary angle units. + /// Also the accuracy floor for [`Self::with_precision`]. The extra three + /// units cover the full-iteration CORDIC residual, below 2^-28 radians. + pub const MAX_ERROR: u32 = Self::ANGLE_MAX_ERROR + 3; + + /// Returns the achieved rotation angle, estimated from the CORDIC residual. + /// This can differ from the requested angle, particularly with + /// [`Self::with_precision`]. No additional atan2/CORDIC is performed here. + /// The matrix angle is within [`Self::ANGLE_MAX_ERROR`] binary units of it. + /// Tiny nonzero requests may report zero; check before dividing by it. + #[inline] + pub const fn angle(self) -> Angle { + self.angle + } + + #[inline] + fn with_iterations(angle: Angle, iterations: usize) -> Self { + let q = FixedScale::::SHIFT.min(cordic::COEFFICIENT_BITS); + let (sin, cos, bits) = cordic::rotation(angle.bits(), iterations, q); + let coefficient = |value: i32| { + let magnitude = I::Wide::from_u32(value.unsigned_abs()) << (FixedScale::::SHIFT - q); + I::from_wide(if value < 0 { -magnitude } else { magnitude }) + }; + Self { + sin: coefficient(sin), + cos: coefficient(cos), + angle: Angle::from_bits(bits), + } + } + + /// Applies the matrix without CORDIC, square roots, or normalization. + /// Does not increase the input length; an approximate float input may + /// remain slightly longer than one. + #[inline] + pub fn apply(self, vector: UnitIntVector) -> UnitIntVector { + let (sin, cos) = (self.sin.to_wide(), self.cos.to_wide()); + let (x, y) = (vector.x().to_wide(), vector.y().to_wide()); + let scale = FixedScale::::DENOMINATOR; + // Products and sums fit I::Wide for approximately unit inputs. + UnitIntVector::from_components( + I::from_wide((cos * x - sin * y) / scale), + I::from_wide((sin * x + cos * y) / scale), + ) + } +} diff --git a/src/int/angle/tests.rs b/src/int/angle/tests.rs new file mode 100644 index 0000000..b4f72fb --- /dev/null +++ b/src/int/angle/tests.rs @@ -0,0 +1,556 @@ +extern crate std; + +use super::*; +use crate::int::number::{int::IntNumber, uint::UIntNumber, wide_int::WideIntNumber}; +use crate::int::unit_vector::UnitIntVector; +use crate::int::vector::IntVector; +use core::f64::consts::TAU; + +fn unit(x: i64, y: i64) -> UnitIntVector { + let component = |v: i64| { + let n = I::Wide::from_uint(I::WideUInt::from_u64(v.unsigned_abs())); + if v < 0 { -n } else { n } + }; + IntVector::::new(component(x), component(y)) + .fast_normalize() + .unwrap() +} + +fn check_angles() { + let axes = [ + unit::(1, 0), + unit::(0, 1), + unit::(-1, 0), + unit::(0, -1), + ]; + for i in 0..4 { + for j in 0..4 { + assert_eq!( + Angle::between(axes[i], axes[j]).bits(), + ((j as u32).wrapping_sub(i as u32)).wrapping_mul(1 << 30) + ); + } + } + let a = unit::(3, 4); + assert_eq!(Angle::between(a, unit::(6, 8)).bits(), 0); + assert_eq!(Angle::between(a, unit::(-3, -4)).bits(), 1 << 31); + let denominator = UnitIntVector::::DENOMINATOR; + let near = UnitIntVector::::from_components(I::from_wide(denominator - I::Wide::ONE), I::ONE); + assert!(Angle::between(axes[0], near).bits() > 0); + assert!(Angle::between(near, axes[0]).bits() > 1 << 31); + let near_opposite = UnitIntVector::::from_components(-near.x(), near.y()); + assert!(Angle::between(axes[0], near_opposite).bits() < 1 << 31); + assert!(Angle::between(near_opposite, axes[0]).bits() > 1 << 31); + + let mut random = 0x876abcd987u64; + let mut next = || { + random ^= random << 13; + random ^= random >> 7; + random ^= random << 17; + random + }; + let mut max_error = 0.0f64; + for _ in 0..20_000 { + // Includes wide i128 products around 2^124 for i64, with independent + // input lengths. Components chosen inside the unit disk without sqrt. + let make = |raw: u64| { + let x = I::Wide::from_u32(raw as u32 & 0xffff); + let y = I::Wide::from_u32((raw >> 32) as u32 & 0xffff); + let shift = I::BITS.saturating_sub(19); + let (x, y) = (x << shift, y << shift); + // i16 needs small values too; construction always nonzero. + let x = if I::BITS == 16 { + (x >> 3) + I::Wide::ONE + } else { + x + I::Wide::ONE + }; + let y = if I::BITS == 16 { y >> 3 } else { y }; + UnitIntVector::::from_components( + I::from_wide(if raw & 1 == 0 { x } else { -x }), + I::from_wide(if raw & 2 == 0 { y } else { -y }), + ) + }; + let (a, b) = (make(next()), make(next())); + let expected = + libm::atan2(b.y().to_f64(), b.x().to_f64()) - libm::atan2(a.y().to_f64(), a.x().to_f64()); + let measured = Angle::between(a, b).bits() as f64 * (TAU / 4294967296.0); + let error = libm::atan2(libm::sin(measured - expected), libm::cos(measured - expected)).abs() + * 4294967296.0 + / TAU; + max_error = max_error.max(error); + assert!(error <= Angle::MAX_ERROR as f64, "{}: {error}", I::BITS); + let forward = Angle::between(a, b).bits(); + let backward = Angle::between(b, a).bits(); + assert_eq!(forward.wrapping_add(backward), 0); + } + std::println!("i{} vectoring max error: {max_error:.3} binary units", I::BITS); +} + +#[test] +fn signs_quadrants_and_wide_products() { + check_angles::(); + check_angles::(); + check_angles::(); +} + +#[test] +fn coefficient_norm_and_accuracy() { + let mut max_error = 0.0f64; + let mut min_norm = 1.0f64; + for bits in (0..=u32::MAX).step_by(16_381) { + let (sin, cos) = Angle::from_bits(bits).sin_cos(); + let norm2 = sin as i64 * sin as i64 + cos as i64 * cos as i64; + assert!(norm2 <= 1 << 60); + let angle = bits as f64 * TAU / 4294967296.0; + let (s, c) = (sin as f64 / 1073741824.0, cos as f64 / 1073741824.0); + let error = libm::hypot(s - libm::sin(angle), c - libm::cos(angle)); + max_error = max_error.max(error); + min_norm = min_norm.min(libm::hypot(s, c)); + let bound = 1.0 / (1u64 << 28) as f64 + core::f64::consts::SQRT_2 / (1u64 << 30) as f64 + 4e-13; + assert!(error <= bound, "{bits}: {error}"); + } + std::println!( + "rotation coefficient error {max_error:e}, max norm loss {:e}", + 1.0 - min_norm + ); +} + +fn drift() { + let s = UnitIntVector::::DENOMINATOR.to_f64(); + let mut angular = 0.0f64; + let mut radial = 0.0f64; + let mut displacement = 0.0f64; + let mut total = 0.0f64; + let mut scaled = 0.0f64; + // All points, including the computed endpoint: no endpoint replacement. + for sample in 0..257u32 { + let start_angle = sample as f64 * TAU / 257.0; + let magnitude = if I::BITS == 16 { 10_000.0 } else { 1_000_000.0 }; + let initial = unit::( + libm::round(libm::cos(start_angle) * magnitude) as i64, + libm::round(libm::sin(start_angle) * magnitude) as i64, + ); + let ix = initial.x().to_f64() / s; + let iy = initial.y().to_f64() / s; + let initial_length = libm::hypot(ix, iy); + // 1/1024 turn through 1/8 turn, both orientations. Not just table angles. + let step = (1 << 22) + ((sample as u64 * ((1 << 29) - (1 << 22))) / 256) as u32; + for cw in [false, true] { + let bits = if cw { step.wrapping_neg() } else { step }; + let r = Rotation::new(Angle::from_bits(bits)); + let mut v = initial; + for index in 1..=1024 { + let previous = v; + v = r.apply(v); + let length2 = |v: UnitIntVector| { + let x = v.x().to_wide(); + let y = v.y().to_wide(); + x * x + y * y + }; + assert!(length2(v) <= length2(previous)); + let angle = + (if cw { -(step as f64) } else { step as f64 }) * index as f64 * TAU / 4294967296.0; + let (sn, cs) = (libm::sin(angle), libm::cos(angle)); + let (ex, ey) = (cs * ix - sn * iy, sn * ix + cs * iy); + let (x, y) = (v.x().to_f64() / s, v.y().to_f64() / s); + let length = libm::hypot(x, y); + displacement = displacement.max(libm::hypot(x - ex, y - ey)); + radial = radial.max(initial_length - length); + angular = angular.max(libm::atan2(x * ey - y * ex, x * ex + y * ey).abs()); + total = total.max(libm::hypot(x - ex / initial_length, y - ey / initial_length)); + if I::BITS != 16 { + let p = v.scale(I::from_u32(65536)); + scaled = scaled.max(libm::hypot( + p.x.to_f64() - 65536.0 * ex, + p.y.to_f64() - 65536.0 * ey, + )); + } + } + } + } + std::println!( + "i{} /1024 steps: angular={angular:e} rad; at R=1024 / 65536: added={:.6} / {:.6}, radial={:.6} / {:.6}, incl input={:.6} / {:.6}; rounded R65536={scaled:.6}", + I::BITS, + displacement * 1024.0, + displacement * 65536.0, + radial * 1024.0, + radial * 65536.0, + total * 1024.0, + total * 65536.0 + ); + if I::BITS >= 32 { + assert!(displacement * 65536.0 < 1.0); + assert!(radial * 65536.0 < 0.3); + assert!(scaled < 1.7); + } else { + // Q14 loses bits in both the coefficients and every stored result. + // Matrix norm <=1 makes the sum of the per-step errors a valid bound. + let per_step = 1.0 / (1u64 << 28) as f64 + 2.0 * libm::sqrt(2.0) / s + 4e-13; + assert!(displacement <= 1024.0 * per_step); + } +} + +#[test] +fn repeated_rotations_report_accumulation() { + drift::(); + drift::(); + drift::(); +} + +#[test] +fn exact_cardinal_rotations_all_types() { + fn check() { + let v = unit::(3, 4); + let quarter = Rotation::new(Angle::from_bits(1 << 30)); + let rotated = quarter.apply(v); + assert!(rotated.x() == -v.y() && rotated.y() == v.x()); + let mut p = v; + for _ in 0..1024 { + p = quarter.apply(p); + } + assert!(p.x() == v.x() && p.y() == v.y()); + let p = Rotation::new(Angle::from_bits(0)).apply(v); + assert!(p.x() == v.x() && p.y() == v.y()); + let p = Rotation::new(Angle::from_bits(1 << 31)).apply(v); + assert!(p.x() == -v.x() && p.y() == -v.y()); + } + check::(); + check::(); + check::(); +} + +#[test] +fn directed_arc_subdivision_and_endpoint_gap() { + let mut worst_gap_excess = 0.0f64; + let mut worst_error = 0.0f64; + let mut max_points = 0u64; + for start in 0..64 { + let a = start as f64 * TAU / 64.0; + for sweep in [ + 0.0, + 1e-7, + 0.003, + TAU / 8.0, + TAU / 2.0 - 1e-7, + TAU / 2.0, + TAU / 2.0 + 1e-7, + TAU * 0.75, + TAU - 1e-7, + ] { + for clockwise in [false, true] { + let sign = if clockwise { -1.0 } else { 1.0 }; + let b = a + sign * sweep; + let make = |v: f64| { + unit::( + libm::round(libm::cos(v) * 1e9) as i64, + libm::round(libm::sin(v) * 1e9) as i64, + ) + }; + let from = make(a); + let to = make(b); + let angle = if clockwise { + Angle::between(to, from) + } else { + Angle::between(from, to) + }; + if sweep == 0.0 { + assert_eq!(angle.bits(), 0); + continue; + } + let get_angle = |v: UnitIntVector| libm::atan2(v.y() as f64, v.x() as f64); + let mut exact = sign * (get_angle(to) - get_angle(from)); + if exact < 0.0 { + exact += TAU; + } + for max_step in [1u32 << 22, (4294967296u64 / 360) as u32, 1 << 29] { + let upper = angle.bits() as u64 + Angle::MAX_ERROR as u64; + let mut n = upper.div_ceil(max_step as u64); + while upper + 8 * n * n > n * max_step as u64 { + n += 1; + } + assert!(n >= libm::ceil(exact / (max_step as f64 * TAU / 4294967296.0)) as u64); + assert!(n <= 1030); + max_points = max_points.max(n - 1); + let step = (angle.bits() as u64 / n) as u32; + let rotation = Rotation::new(Angle::from_bits(if clockwise { + step.wrapping_neg() + } else { + step + })); + let mut v = from; + let mut previous = 0.0; + let norm = libm::hypot(from.x() as f64, from.y() as f64) / 1073741824.0; + for k in 1..n { + v = rotation.apply(v); + let mut travelled = sign * (get_angle(v) - get_angle(from)); + if travelled < 0.0 { + travelled += TAU; + } + assert!(travelled > previous && travelled < exact); + worst_gap_excess = + worst_gap_excess.max(travelled - previous - max_step as f64 * TAU / 4294967296.0); + previous = travelled; + let expected = get_angle(from) + sign * exact * k as f64 / n as f64; + let error = libm::hypot( + v.x() as f64 / 1073741824.0 - norm * libm::cos(expected), + v.y() as f64 / 1073741824.0 - norm * libm::sin(expected), + ); + worst_error = worst_error.max(error * 65536.0); + } + // Includes the gap to the exact stored endpoint, not an + // artificially replaced computed point. + worst_gap_excess = + worst_gap_excess.max(exact - previous - max_step as f64 * TAU / 4294967296.0); + } + } + } + } + std::println!( + "directed arcs: max {max_points} intermediate points; added error R65536={worst_error:.6}; max step excess from accumulated arithmetic={worst_gap_excess:e} rad" + ); + assert!(worst_error < 0.5); + assert!(worst_gap_excess < 1e-12); +} + +#[test] +fn public_trigonometry_and_atan2_cover_full_integer_range() { + fn check() { + assert_eq!(Angle::atan2(W::ZERO, W::ZERO), None); + for y in [W::MIN, W::MIN + W::ONE, -W::ONE, W::ZERO, W::ONE, W::MAX] { + for x in [W::MIN, W::MIN + W::ONE, -W::ONE, W::ZERO, W::ONE, W::MAX] { + let Some(angle) = Angle::atan2(y, x) else { + continue; + }; + let expected = libm::atan2(y.to_f64(), x.to_f64()); + let actual = angle.bits() as f64 * TAU / 4294967296.0; + let difference = + libm::atan2(libm::sin(actual - expected), libm::cos(actual - expected)).abs(); + assert!(difference * 4294967296.0 / TAU <= Angle::MAX_ERROR as f64); + } + } + } + check::(); + check::(); + check::(); + for (bits, expected) in [ + (0, (0, 1 << 30)), + (1 << 30, (1 << 30, 0)), + (1 << 31, (0, -(1 << 30))), + (3 << 30, (-(1 << 30), 0)), + ] { + let angle = Angle::from_bits(bits); + assert_eq!(angle.sin_cos(), expected); + assert_eq!(angle.sin(), expected.0); + assert_eq!(angle.cos(), expected.1); + } + for bits in (0..=u32::MAX).step_by(65537) { + let angle = Angle::from_bits(bits); + assert_eq!(angle.sin_cos(), (angle.sin(), angle.cos())); + } +} + +#[test] +fn rotation_uses_coordinate_scale() { + macro_rules! check { + ($t:ty, $q:expr) => { + for bits in (0..=u32::MAX).step_by(1_048_573) { + let angle = Angle::from_bits(bits); + let rotation = Rotation::<$t>::new(angle); + let (sin, cos) = angle.sin_cos(); + // Independent rescaling of the public Q30 coefficients. + let convert = |v: i32| (v as i128 * (1i128 << $q)) / (1i128 << 30); + let (sin, cos) = (convert(sin), convert(cos)); + for (x, y) in [(1, 0), (0, -1), (3, 4), (-5, 2), (-1, -7)] { + let from = unit::<$t>(x, y); + let to = rotation.apply(from); + let (x, y) = (from.x() as i128, from.y() as i128); + assert_eq!(to.x() as i128, (cos * x - sin * y) / (1i128 << $q)); + assert_eq!(to.y() as i128, (sin * x + cos * y) / (1i128 << $q)); + } + } + }; + } + check!(i16, 14); + check!(i32, 30); + check!(i64, 62); +} + +#[test] +fn approximate_rotation_reports_achieved_angle_and_respects_budget() { + fn check_type() { + let axis = unit::(1, 0); + let check = |bits: u32, precision: u32| { + let requested = Angle::from_bits(bits); + let rotation = Rotation::with_precision(requested, precision); + // Applying an exact axis exposes the stored coefficients without adding + // component rounding, independently of the residual-based estimate. + let v = rotation.apply(axis); + let actual = libm::atan2(v.y().to_f64(), v.x().to_f64()); + let to_radians = |a: Angle| a.bits() as f64 * TAU / 4294967296.0; + let difference = + |a: f64, b: f64| libm::atan2(libm::sin(a - b), libm::cos(a - b)).abs() * 4294967296.0 / TAU; + assert!( + difference(actual, to_radians(rotation.angle())) <= Rotation::::ANGLE_MAX_ERROR as f64, + "achieved angle: bits={bits}, precision={precision}" + ); + let budget = (bits as i32) + .unsigned_abs() + .checked_shr(precision) + .unwrap_or(0) + .max(Rotation::::MAX_ERROR); + assert!( + difference(actual, to_radians(requested)) <= budget as f64, + "requested angle: bits={bits}, precision={precision}" + ); + let (x, y) = (v.x().to_wide(), v.y().to_wide()); + let scale = UnitIntVector::::DENOMINATOR; + assert!(x * x + y * y <= scale * scale); + if precision >= 32 { + let precise = Rotation::new(requested); + assert_eq!(rotation.angle(), precise.angle()); + assert!(v == precise.apply(axis)); + } + if bits & ((1 << 30) - 1) == 0 { + assert_eq!(rotation.angle(), requested); + assert!(v == Rotation::new(requested).apply(axis)); + } + }; + for precision in [0, 1, 2, 3, 4, 8, 16, 29, 31, 32, u32::MAX] { + for bits in (0..=u32::MAX).step_by(65537) { + check(bits, precision); + } + for axis in [0u32, 1 << 30, 1 << 31, 3 << 30] { + for delta in [0u32, 1, 2, 3, 4, u32::MAX, u32::MAX - 1, u32::MAX - 2] { + check(axis.wrapping_add(delta), precision); + } + } + } + } + check_type::(); + check_type::(); + check_type::(); +} + +#[test] +fn achieved_angle_drives_approximate_arc_counts() { + let mut max_points = 0; + // Test-only consumer policy: leave margin for one-step error and the final + // gap, then count by the achieved matrix angle, not the requested one. + // i32, fresh normalization, relative exponents 3/4, at most 1400 rotations. + const LIMIT: u64 = 1400; + const ERROR: u64 = 4; // two metadata units plus component rounding + for precision in [3u32, 4] { + for max_step in [ + 1u32 << 22, + (4294967296u64 / 360) as u32, + (4294967296u64 / 18) as u32, + 1 << 29, + ] { + let divisor = 1u64 << precision; + let reserve = 2 * ERROR * (LIMIT + 1) + 2 * Angle::MAX_ERROR as u64; + let requested = ((max_step as u64 - reserve) * divisor / (divisor + 1)) as u32; + for clockwise in [false, true] { + let rotation = Rotation::with_precision( + Angle::from_bits(if clockwise { + requested.wrapping_neg() + } else { + requested + }), + precision, + ); + let step = (rotation.angle().bits() as i32).unsigned_abs() as u64; + assert!(step > ERROR); + for start in 0..16 { + for sweep in [0.001, 0.3, TAU / 2.0, 3.0 * TAU / 4.0, TAU - 0.0001] { + let initial = start as f64 * TAU / 16.0; + let sign = if clockwise { -1.0 } else { 1.0 }; + let make = + |a: f64| unit::((libm::cos(a) * 1e6) as i64, (libm::sin(a) * 1e6) as i64); + let from = make(initial); + let to = make(initial + sign * sweep); + let measured = if clockwise { + Angle::between(to, from) + } else { + Angle::between(from, to) + }; + let lower = measured.bits().saturating_sub(Angle::MAX_ERROR) as u64; + let count = lower.saturating_sub(1) / (step + ERROR); + assert!(count < LIMIT); + max_points = max_points.max(count); + let angle = |v: UnitIntVector| libm::atan2(v.y() as f64, v.x() as f64); + let mut exact = sign * (angle(to) - angle(from)); + if exact < 0.0 { + exact += TAU; + } + let mut previous = 0.0; + let mut v = from; + for _ in 0..count { + v = rotation.apply(v); + let mut travelled = sign * (angle(v) - angle(from)); + if travelled < 0.0 { + travelled += TAU; + } + assert!(travelled > previous && travelled < exact); + assert!(travelled - previous <= max_step as f64 * TAU / 4294967296.0 + 1e-12); + previous = travelled; + } + assert!(exact - previous <= max_step as f64 * TAU / 4294967296.0 + 1e-12); + } + } + } + } + } + std::println!( + "approximate arcs using achieved angle: max {max_points} intermediate points; all points ordered, all gaps within maximum" + ); +} + +#[test] +fn float_between_preserves_exact_products_and_turn_side() { + fn check() { + let axes = [ + unit::(1, 0), + unit::(0, 1), + unit::(-1, 0), + unit::(0, -1), + ]; + for i in 0..4 { + for j in 0..4 { + assert_eq!( + Angle::between_with_float(axes[i], axes[j]).bits(), + ((j as u32).wrapping_sub(i as u32)).wrapping_mul(1 << 30) + ); + } + } + let n = UnitIntVector::::DENOMINATOR >> 1; + let make = |x, y| UnitIntVector::::from_components(I::from_wide(x), I::from_wide(y)); + let a = make(n, n - I::Wide::ONE); + let b = make(n - I::Wide::ONE, n - I::Wide::TWO); + // The exact cross product is -1. Converting i64 components before + // multiplication would erase it and incorrectly produce an empty arc. + if I::BITS == 64 { + assert_eq!( + a.x().to_f64() * b.y().to_f64() - a.y().to_f64() * b.x().to_f64(), + 0.0 + ); + } + assert!(Angle::between_with_float(a, b).bits() > 1 << 31); + assert!(Angle::between_with_float(b, a).bits() > 0); + assert!(Angle::between_with_float(b, a).bits() < 1 << 31); + assert_eq!( + Angle::between_with_float(a, b) + .bits() + .wrapping_add(Angle::between_with_float(b, a).bits()), + 0 + ); + let opposite = make(-b.x().to_wide(), -b.y().to_wide()); + assert!(Angle::between_with_float(a, opposite).bits() < 1 << 31); + assert!(Angle::between_with_float(opposite, a).bits() > 1 << 31); + let ray = make(n, n); + assert_eq!(Angle::between_with_float(ray, make(n >> 1, n >> 1)).bits(), 0); + assert_eq!(Angle::between_with_float(ray, make(-n, -n)).bits(), 1 << 31); + } + check::(); + check::(); + check::(); +} diff --git a/src/int/mod.rs b/src/int/mod.rs index ddacebd..b46de56 100644 --- a/src/int/mod.rs +++ b/src/int/mod.rs @@ -1,4 +1,6 @@ +pub mod angle; pub mod number; pub mod point; pub mod rect; +pub mod unit_vector; pub mod vector; diff --git a/src/int/number/fixed_scale.rs b/src/int/number/fixed_scale.rs index 6855677..9bcf9f5 100644 --- a/src/int/number/fixed_scale.rs +++ b/src/int/number/fixed_scale.rs @@ -20,6 +20,14 @@ impl FixedScale { I::from_wide(Self::div_round(scaled, Self::DENOMINATOR)) } + /// Returns `numerator * DENOMINATOR / denominator`, rounded to the nearest + /// integer with midpoint values away from zero. + /// + /// # Preconditions + /// The caller must ensure that `denominator` is neither zero nor + /// `I::Wide::MIN`, and that the rounded scaled result fits in `I::Wide`. + /// The intermediate product uses extended-width arithmetic. + /// Inputs outside this contract may panic or produce an incorrect result. #[inline(always)] pub fn div_to_scaled_round(numerator: I::Wide, denominator: I::Wide) -> I::Wide { debug_assert!(denominator != I::Wide::ZERO); @@ -34,6 +42,13 @@ impl FixedScale { Self::from_unsigned_abs(quotient, negative) } + /// Returns `numerator / denominator`, rounded to the nearest integer with + /// midpoint values away from zero. + /// + /// # Preconditions + /// The caller must ensure that `denominator` is neither zero nor + /// `I::Wide::MIN`, and that the rounded result fits in `I::Wide`. + /// Inputs outside this contract may panic or produce an incorrect result. #[inline(always)] pub fn div_round(numerator: I::Wide, denominator: I::Wide) -> I::Wide { debug_assert!(denominator != I::Wide::ZERO); diff --git a/src/int/number/uint.rs b/src/int/number/uint.rs index 20d8cc1..19d9cb7 100644 --- a/src/int/number/uint.rs +++ b/src/int/number/uint.rs @@ -28,6 +28,8 @@ pub trait UIntNumber: const LAST_BIT_INDEX: u32; const ZERO: Self; const ONE: Self; + const TWO: Self; + const FOUR: Self; const MAX: Self; const LAST_BIT: Self; @@ -49,6 +51,8 @@ impl UIntNumber for u32 { const LAST_BIT_INDEX: u32 = 31; const ZERO: Self = 0; const ONE: Self = 1; + const TWO: Self = 2; + const FOUR: Self = 4; const MAX: Self = Self::MAX; const LAST_BIT: Self = Self::ONE << Self::LAST_BIT_INDEX; @@ -99,6 +103,8 @@ impl UIntNumber for u64 { const LAST_BIT_INDEX: u32 = 63; const ZERO: Self = 0; const ONE: Self = 1; + const TWO: Self = 2; + const FOUR: Self = 4; const MAX: Self = Self::MAX; const LAST_BIT: Self = Self::ONE << Self::LAST_BIT_INDEX; @@ -150,6 +156,8 @@ impl UIntNumber for u128 { const LAST_BIT_INDEX: u32 = 127; const ZERO: Self = 0; const ONE: Self = 1; + const TWO: Self = 2; + const FOUR: Self = 4; const MAX: Self = Self::MAX; const LAST_BIT: Self = Self::ONE << Self::LAST_BIT_INDEX; #[inline(always)] diff --git a/src/int/point.rs b/src/int/point.rs index 39d7f90..f27a075 100644 --- a/src/int/point.rs +++ b/src/int/point.rs @@ -1,4 +1,5 @@ use crate::int::number::int::IntNumber; +use crate::int::number::wide_int::WideIntNumber; use crate::int::vector::IntVector; use core::cmp::Ordering; use core::{fmt, ops}; @@ -19,6 +20,8 @@ use core::{fmt, ops}; /// -2^(T::BITS - 2) < coordinate < 2^(T::BITS - 2) /// ``` /// +/// Use [`Self::is_in_safe_range`] to check this precondition. +/// /// This guarantees enough headroom for a point difference and for the sum or /// difference of two products. Arithmetic is unchecked beyond Rust's normal /// debug overflow checks. Callers may use a wider range only when they prove @@ -40,6 +43,17 @@ impl IntPoint { Self { x, y } } + /// Returns whether both coordinates are in the conservative arithmetic range. + /// + /// Each coordinate must be strictly between `-2^(T::BITS - 2)` and + /// `2^(T::BITS - 2)`. See the type's arithmetic range documentation. + #[inline(always)] + pub fn is_in_safe_range(&self) -> bool { + let limit = T::ONE << (T::BITS - 2); + let min = -limit; + self.x > min && self.x < limit && self.y > min && self.y < limit + } + #[inline(always)] pub fn cross_product(self, v: Self) -> T::Wide { let a = self.x.to_wide() * v.y.to_wide(); @@ -56,14 +70,16 @@ impl IntPoint { } #[inline(always)] - pub fn sqr_length(self) -> T::Wide { + pub fn sqr_length(self) -> T::WideUInt { let x = self.x.to_wide(); let y = self.y.to_wide(); - x * x + y * y + let xx = x * x; + let yy = y * y; + xx.to_uint() + yy.to_uint() } #[inline(always)] - pub fn sqr_distance(self, other: Self) -> T::Wide { + pub fn sqr_distance(self, other: Self) -> T::WideUInt { (self - other).sqr_length() } } @@ -146,14 +162,46 @@ impl From> for IntPoint { #[macro_export] macro_rules! int_pnt { ($x:expr, $y:expr) => { - IntPoint::new($x, $y) + $crate::int::point::IntPoint::new($x, $y) }; } #[cfg(test)] mod tests { + use crate::int::number::int::IntNumber; use crate::int::point::IntPoint; + fn assert_safe_range(limit: T) { + assert!(IntPoint::::ZERO.is_in_safe_range()); + assert!(!IntPoint::::EMPTY.is_in_safe_range()); + + for x in [-limit + T::ONE, T::ZERO, limit - T::ONE] { + for y in [-limit + T::ONE, T::ZERO, limit - T::ONE] { + assert!(IntPoint::new(x, y).is_in_safe_range()); + } + } + + for value in [T::MIN, -limit - T::ONE, -limit, limit, limit + T::ONE, T::MAX] { + assert!(!IntPoint::new(value, T::ZERO).is_in_safe_range()); + assert!(!IntPoint::new(T::ZERO, value).is_in_safe_range()); + } + } + + #[test] + fn test_safe_range_i16() { + assert_safe_range::(16_384); + } + + #[test] + fn test_safe_range_i32() { + assert_safe_range::(1_073_741_824); + } + + #[test] + fn test_safe_range_i64() { + assert_safe_range::(4_611_686_018_427_387_904); + } + #[test] fn test_0() { let p: IntPoint = (1, 2).into(); diff --git a/src/int/rect.rs b/src/int/rect.rs index 10b25b8..75f43fe 100644 --- a/src/int/rect.rs +++ b/src/int/rect.rs @@ -11,6 +11,18 @@ pub struct IntRect { } impl IntRect { + /// Returns whether all four bounds are in the conservative arithmetic range. + /// + /// Each bound must be strictly between `-2^(T::BITS - 2)` and + /// `2^(T::BITS - 2)`, as in [`IntPoint::is_in_safe_range`]. + /// This checks coordinate range only, not whether each minimum is at most + /// its corresponding maximum. + #[inline(always)] + pub fn is_in_safe_range(&self) -> bool { + IntPoint::new(self.min_x, self.min_y).is_in_safe_range() + && IntPoint::new(self.max_x, self.max_y).is_in_safe_range() + } + #[inline(always)] pub fn width(&self) -> T { self.max_x - self.min_x @@ -183,9 +195,43 @@ impl From<[IntPoint; 2]> for IntRect { #[cfg(test)] mod tests { + use crate::int::number::int::IntNumber; use crate::int::point::IntPoint; use crate::int::rect::IntRect; + fn assert_safe_range(limit: T) { + let min = -limit + T::ONE; + let max = limit - T::ONE; + assert!(IntRect::new(min, max, min, max).is_in_safe_range()); + assert!(IntRect::new(max, min, max, min).is_in_safe_range()); + + for value in [min, T::ZERO, max] { + assert!(IntRect::with_point(IntPoint::new(value, value)).is_in_safe_range()); + } + + for value in [T::MIN, -limit - T::ONE, -limit, limit, limit + T::ONE, T::MAX] { + assert!(!IntRect::new(value, max, min, max).is_in_safe_range()); + assert!(!IntRect::new(min, value, min, max).is_in_safe_range()); + assert!(!IntRect::new(min, max, value, max).is_in_safe_range()); + assert!(!IntRect::new(min, max, min, value).is_in_safe_range()); + } + } + + #[test] + fn test_safe_range_i16() { + assert_safe_range::(16_384); + } + + #[test] + fn test_safe_range_i32() { + assert_safe_range::(1_073_741_824); + } + + #[test] + fn test_safe_range_i64() { + assert_safe_range::(4_611_686_018_427_387_904); + } + #[test] fn test_0() { let rect = if let Some(rect) = diff --git a/src/int/unit_vector.rs b/src/int/unit_vector.rs new file mode 100644 index 0000000..151e388 --- /dev/null +++ b/src/int/unit_vector.rs @@ -0,0 +1,260 @@ +use crate::float::number::FloatNumber; +use crate::int::number::fixed_scale::FixedScale; +use crate::int::number::int::IntNumber; +use crate::int::number::uint::UIntNumber; +use crate::int::number::wide_int::WideIntNumber; +use crate::int::vector::IntVector; +use core::ops::Mul; + +/// An approximate unit direction stored as fixed-scale integer components. +/// +/// Obtain one with [`IntVector::fast_normalize`] or [`Self::normalize_with_float`]. +/// The represented components are `x() / DENOMINATOR` and `y() / DENOMINATOR`. +/// Floating-point construction is approximate: rounding may make the length +/// slightly greater than one. Integer normalization via `fast_normalize` and +/// checked conversion via [`Self::try_from_float`] keep the length at most one. +/// +/// Integer normalization via `fast_normalize` favors speed over precision: +/// it keeps about 6, 14, or 30 bits +/// of direction precision for `i16`, `i32`, or `i64`, respectively. +/// The storage scale does not imply that all stored bits are accurate. +/// [`Rotation::apply`](crate::int::angle::Rotation::apply) does not increase +/// length, but repeated rotations accumulate contraction and angular +/// error; the normalization precision above is not a bound on that accumulation. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub struct UnitIntVector { + x: T, + y: T, +} + +impl UnitIntVector { + // The caller supplies an approximately unit or shorter nonzero vector. + #[inline(always)] + pub(crate) fn from_components(x: T, y: T) -> Self { + Self { x, y } + } + /// Converts floating-point components to the fixed-scale representation. + /// + /// Components are truncated toward zero. This does not normalize the input: + /// callers supplying a direction should normalize it first. Returns `None` + /// for non-finite components, components outside [-1, 1], a quantized zero + /// vector, or an integer representation whose squared length exceeds one. + /// Floating-point normalization/rotation may overshoot one slightly; callers + /// must account for that error before conversion. The final length check is + /// exact integer arithmetic, including for i64. + #[inline] + pub fn try_from_float(x: F, y: F) -> Option { + let (x, y) = (x.to_f64(), y.to_f64()); + if !x.is_in_safe_range() || !y.is_in_safe_range() || x.abs() > 1.0 || y.abs() > 1.0 { + return None; + } + let unit = Self::from_float_unchecked(x, y); + let (wx, wy) = (unit.x.to_wide(), unit.y.to_wide()); + let squared = wx * wx + wy * wy; + if squared == T::Wide::ZERO || squared > Self::DENOMINATOR * Self::DENOMINATOR { + return None; + } + Some(unit) + } + + /// Converts already prepared floating-point components without validation + /// or normalization. Components are scaled and truncated toward zero. + /// + /// The caller supplies finite components of an approximately unit or shorter + /// vector that remains nonzero after quantization. Small floating-point + /// overshoots above unit length are accepted. Use [`Self::try_from_float`] + /// when the fixed-scale length must be checked against one. + #[inline] + pub fn from_float_unchecked(x: F, y: F) -> Self { + let scale = Self::DENOMINATOR.to_f64(); + Self { + x: T::from_float(x.to_f64() * scale), + y: T::from_float(y.to_f64() * scale), + } + } + + /// Normalizes an integer vector using f64 arithmetic, returning `None` only + /// for the zero vector. Supports the full built-in wide component range, + /// including i128::MIN; no integer products of the input are formed. + /// + /// This is an approximate, unchecked normalization: rounding may make the + /// output length slightly greater than one, including for axis directions. + /// Precision is limited by f64 arithmetic and the output component scale. + /// When normalizing a point difference, subtract the integer points first + /// to retain small differences between large coordinates. + #[inline] + pub fn normalize_with_float(vector: IntVector) -> Option { + if vector.x == T::Wide::ZERO && vector.y == T::Wide::ZERO { + return None; + } + let x = vector.x.to_f64(); + let y = vector.y.to_f64(); + // Squaring and adding even two i128 components fits in f64. + let inv_length = 1.0 / FloatNumber::sqrt(x * x + y * y); + Some(Self::from_float_unchecked(x * inv_length, y * inv_length)) + } + + /// The stored integer value representing one: 2^14, 2^30, or 2^62. + pub const DENOMINATOR: T::Wide = FixedScale::::DENOMINATOR; + + /// Returns the stored, scaled x component. + #[inline(always)] + pub fn x(self) -> T { + self.x + } + + /// Returns the stored, scaled y component. + #[inline(always)] + pub fn y(self) -> T { + self.y + } + + /// Scales this direction into an integer vector, rounding to the nearest + /// integer with midpoint values away from zero. Every value of `T` is valid. + /// Direction approximation error grows with the scalar's magnitude. + #[inline(always)] + pub fn scale(self, scalar: T) -> IntVector { + let scalar = scalar.to_wide(); + IntVector { + x: (self.x.to_wide() * scalar).shr_round(FixedScale::::SHIFT), + y: (self.y.to_wide() * scalar).shr_round(FixedScale::::SHIFT), + } + } + + #[inline(always)] + pub(crate) fn with_vector(vector: IntVector) -> Option { + let sqr_length = vector.sqr_length(); + if sqr_length == T::WideUInt::ZERO { + return None; + } + + let precision = T::HALF_POWER_OF_TWO - 1; + let shift = (sqr_length.ilog2().saturating_sub(2 * precision) + 1) >> 1; + // Round the reduced squared length upward so the reciprocal never + // makes the represented direction longer than one. This also works + // with shift == 0 and avoids an overflowing addition before shifting. + let one = T::WideUInt::ONE; + let reduced_length = ((sqr_length - one) >> (2 * shift)) + one; + let sqr_unit = one << (2 * FixedScale::::SHIFT); + let scale = (sqr_unit / reduced_length).isqrt(); + + // scale represents S * 2^shift / length. Apply it to the original + // components, then discard the extra fractional bits toward zero. + let x = T::from_uint((scale * vector.x.unsigned_abs()) >> shift); + let y = T::from_uint((scale * vector.y.unsigned_abs()) >> shift); + Some(Self { + x: if vector.x < T::Wide::ZERO { -x } else { x }, + y: if vector.y < T::Wide::ZERO { -y } else { y }, + }) + } +} + +impl Mul for UnitIntVector { + type Output = IntVector; + + #[inline(always)] + fn mul(self, scalar: T) -> Self::Output { + self.scale(scalar) + } +} + +#[cfg(test)] +mod tests { + use super::UnitIntVector; + use crate::int::point::IntPoint; + use crate::int::vector::IntVector; + + #[test] + fn scales_direction_into_integer_vector() { + let direction = IntVector::::new(3, 4).fast_normalize().unwrap(); + assert_eq!(direction * 10, IntVector::new(6, 8)); + assert_eq!(direction * -10, IntVector::new(-6, -8)); + assert_eq!(direction * 0, IntVector::new(0, 0)); + } + + macro_rules! check_type { + ($name:ident, $t:ty, $wide:ty, $tolerance:expr) => { + #[test] + fn $name() { + assert!(IntVector::<$t>::new(0, 0).fast_normalize().is_none()); + let denominator = UnitIntVector::<$t>::DENOMINATOR; + let axis = IntVector::<$t>::new(-1, 0).fast_normalize().unwrap(); + // Negating the minimum scalar must produce a positive wide result. + assert_eq!((axis * <$t>::MIN).x, -(<$t>::MIN as $wide)); + assert_eq!((axis * <$t>::MIN).y, 0); + + let check = |x: $wide, y: $wide| { + let vector = IntVector::<$t>::new(x, y); + let Some(unit) = vector.fast_normalize() else { + assert_eq!((x, y), (0, 0)); + return; + }; + let (a, b) = ( + unit.x() as f64 / denominator as f64, + unit.y() as f64 / denominator as f64, + ); + assert!(a.abs() <= 1.0 && b.abs() <= 1.0); + let (raw_x, raw_y) = (unit.x() as $wide, unit.y() as $wide); + assert!(raw_x * raw_x + raw_y * raw_y <= denominator * denominator); + let length = libm::hypot(x as f64, y as f64); + assert!( + (a - x as f64 / length).abs() < $tolerance, + "x: ({x}, {y}) -> ({a}, {b})" + ); + assert!( + (b - y as f64 / length).abs() < $tolerance, + "y: ({x}, {y}) -> ({a}, {b})" + ); + if x != <$wide>::MIN && y != <$wide>::MIN { + let opposite = IntVector::<$t>::new(-x, -y).fast_normalize().unwrap(); + assert_eq!((opposite.x(), opposite.y()), (-unit.x(), -unit.y())); + } + for scalar in [<$t>::MIN, <$t>::MAX, -3, -1, 0, 1, 3] { + let actual = unit * scalar; + for (component, result) in [(unit.x(), actual.x), (unit.y(), actual.y)] { + let product = component as i128 * scalar as i128; + let divisor = denominator as i128; + let quotient = product / divisor; + let remainder = product % divisor; + let rounded = quotient + + if remainder.abs() * 2 >= divisor { + product.signum() + } else { + 0 + }; + assert_eq!(result as i128, rounded); + } + } + }; + + let limit = 2 * denominator - 2; + for x in [-limit, -limit + 1, -3, -1, 0, 1, 3, limit] { + for y in [-limit, -4, -1, 0, 1, 4, limit] { + check(x, y); + } + } + // A valid point difference longer than S must not collapse to zero. + let coordinate = (3 * (denominator / 4)) as $t; + let long = IntPoint::new(coordinate, 0) - IntPoint::new(-coordinate, 0); + check(long.x, long.y); + let direction = long.fast_normalize().unwrap(); + assert!(direction.x() as f64 / denominator as f64 > 1.0 - $tolerance); + assert_eq!(direction.y(), 0); + // Exercise zero and nonzero reciprocal shifts near bit boundaries. + for power in 0..(<$t>::BITS - 1) { + let magnitude = (1 as $wide) << power; + for x in [magnitude - 1, magnitude, magnitude.saturating_add(1)] { + for y in [1, x / 3, x / 2, x] { + check(x, y); + check(-x, y); + } + } + } + } + }; + } + + check_type!(normalize_and_scale_i16, i16, i32, 0.04); + check_type!(normalize_and_scale_i32, i32, i64, 0.00015); + check_type!(normalize_and_scale_i64, i64, i128, 1e-9); +} diff --git a/src/int/vector.rs b/src/int/vector.rs index 8f9b743..1699e77 100644 --- a/src/int/vector.rs +++ b/src/int/vector.rs @@ -1,4 +1,6 @@ use crate::int::number::int::IntNumber; +use crate::int::number::wide_int::WideIntNumber; +use crate::int::unit_vector::UnitIntVector; use core::fmt; #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] @@ -23,6 +25,21 @@ impl IntVector { Self { x, y } } + /// Returns an approximate unit direction, or `None` for the zero vector. + /// + /// Uses integer arithmetic with about 6, 14, or 30 bits of direction precision + /// for `i16`, `i32`, or `i64`, respectively. The squared length is shifted + /// to retain fractional precision in the reciprocal; the original vector + /// components are preserved. The resulting direction has length at most one. + /// + /// Requires the arithmetic range of [`Self::sqr_length`], as guaranteed for + /// point differences by the coordinate range documented on + /// [`IntPoint`](crate::int::point::IntPoint). + #[inline(always)] + pub fn fast_normalize(self) -> Option> { + UnitIntVector::with_vector(self) + } + #[inline(always)] pub fn cross_product(self, v: Self) -> T::Wide { let a = self.x * v.y; @@ -39,10 +56,12 @@ impl IntVector { } #[inline(always)] - pub fn sqr_length(self) -> T::Wide { + pub fn sqr_length(self) -> T::WideUInt { let x = self.x; let y = self.y; - x * x + y * y + let xx = x * x; + let yy = y * y; + xx.to_uint() + yy.to_uint() } } impl fmt::Display for IntVector { diff --git a/src/lib.rs b/src/lib.rs index bf5aae3..6851603 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,9 @@ #![no_std] -#![doc = include_str!("../README.md")] +#![cfg_attr(feature = "core", doc = include_str!("../README.md"))] +#![cfg_attr( + not(feature = "core"), + doc = "Enable the `core` feature (enabled by default) to use the numeric and geometry API." +)] extern crate alloc; #[cfg(feature = "core")] diff --git a/src/triangle.rs b/src/triangle.rs index 688ddbd..d8921c4 100644 --- a/src/triangle.rs +++ b/src/triangle.rs @@ -42,6 +42,12 @@ impl Triangle { Self::area_two(p0, p2, p1).signum() } + /// Tests whether the point is inside the triangle, including its border. + /// + /// # Preconditions + /// The triangle must be non-degenerate: `area_two(p0, p1, p2) != 0`. + /// Collinear or coincident vertices are unsupported and may produce + /// incorrect containment results. Use [`Self::is_not_line`] to check first. #[inline] pub fn is_contain( p: IntPoint, @@ -59,6 +65,12 @@ impl Triangle { !(has_neg && has_pos) } + /// Tests whether the point is outside the triangle's strict interior, + /// including points on its border. + /// + /// # Preconditions + /// The triangle must be non-degenerate. Collinear or coincident vertices + /// are unsupported; see [`Self::is_contain`] for the containment contract. #[inline] pub fn is_not_contain( p: IntPoint, @@ -75,6 +87,11 @@ impl Triangle { has_neg && has_pos } + /// Tests whether the point is strictly inside the triangle, excluding its border. + /// + /// # Preconditions + /// The triangle must be non-degenerate. Collinear or coincident vertices + /// are unsupported; see [`Self::is_contain`] for the containment contract. #[inline] pub fn is_contain_exclude_borders( p: IntPoint, diff --git a/tests/adapter_conservative_tests.rs b/tests/adapter_conservative_tests.rs new file mode 100644 index 0000000..daa7eaa --- /dev/null +++ b/tests/adapter_conservative_tests.rs @@ -0,0 +1,50 @@ +#![cfg(feature = "core")] + +use i_float::adapter::{FloatPointAdapter, FloatPointAdapterScaleError}; +use i_float::float::{number::FloatNumber, rect::FloatRect}; +use i_float::int::number::int::IntNumber; + +fn check_rounding_margin() { + let limit = I::ONE << (I::BITS - 3); + for radius in [1.0, 1.99999, 3.0] { + let r = F::from_float(radius); + let points = [[-r, -r], [r, r]]; + let rect = FloatRect::new(-r, r, -r, r).unwrap(); + let automatic = FloatPointAdapter::<[F; 2], I>::new_conservative(rect); + let scale = automatic.dir_scale(); + let adapters = [ + automatic, + FloatPointAdapter::with_iter_conservative(points.iter()), + FloatPointAdapter::try_with_scale_conservative(rect, scale).unwrap(), + FloatPointAdapter::try_with_iter_and_scale_conservative(points.iter(), scale).unwrap(), + ]; + for adapter in adapters { + for point in &points { + let integer = adapter.try_float_to_int(point).unwrap(); + assert!(integer.is_in_safe_range()); + assert!(-limit <= integer.x && integer.x <= limit); + assert!(-limit <= integer.y && integer.y <= limit); + } + } + let too_large = scale * F::from_float(2.0); + assert_eq!( + FloatPointAdapter::<[F; 2], I>::try_with_scale_conservative(rect, too_large).err(), + Some(FloatPointAdapterScaleError::ScaleTooLarge), + ); + assert_eq!( + FloatPointAdapter::<[F; 2], I>::try_with_iter_and_scale_conservative(points.iter(), too_large) + .err(), + Some(FloatPointAdapterScaleError::ScaleTooLarge), + ); + } +} + +#[test] +fn conservative_constructors_reserve_rounding_margin() { + check_rounding_margin::(); + check_rounding_margin::(); + check_rounding_margin::(); + check_rounding_margin::(); + check_rounding_margin::(); + check_rounding_margin::(); +} diff --git a/tests/adapter_extreme_scale_tests.rs b/tests/adapter_extreme_scale_tests.rs new file mode 100644 index 0000000..c1a470f --- /dev/null +++ b/tests/adapter_extreme_scale_tests.rs @@ -0,0 +1,270 @@ +#![cfg(feature = "core")] + +use i_float::adapter::{FloatPointAdapter, FloatPointAdapterScaleError}; +use i_float::float::number::FloatNumber; +use i_float::float::rect::FloatRect; +use i_float::int::number::int::IntNumber; +use i_float::int::point::IntPoint; + +#[test] +fn automatic_constructors_reject_invalid_bounds() { + // Public fields can bypass FloatRect construction. The adapter must still + // reject these bounds and invalid points supplied through an iterator. + let rect = FloatRect { + min_x: 0.0_f64, + max_x: f64::INFINITY, + min_y: 0.0, + max_y: 0.0, + }; + assert!(std::panic::catch_unwind(|| FloatPointAdapter::<[f64; 2], i32>::new(rect)).is_err()); + let points = [[0.0_f64, 0.0], [f64::INFINITY, 0.0]]; + assert!( + std::panic::catch_unwind(|| { FloatPointAdapter::<[f64; 2], i32>::with_iter(points.iter()) }) + .is_err() + ); +} + +#[test] +fn empty_iterator_uses_zero_bounds_and_unit_scale() { + let points: [[f64; 2]; 0] = []; + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_iter(points.iter()); + assert_eq!(adapter.offset(), [0.0, 0.0]); + assert_eq!(adapter.dir_scale(), 1.0); + assert_eq!(adapter.inv_scale(), 1.0); + assert_eq!(adapter.try_float_to_int(&[0.0, 0.0]), Ok(IntPoint::ZERO)); + assert!(adapter.try_float_to_int(&[1.0, 0.0]).is_err()); +} + +#[test] +fn tiny_finite_bounds_keep_both_scales_finite() { + let rect = FloatRect::new(-1e-35_f32, 1e-35, -1e-35, 1e-35).unwrap(); + for (name, adapter) in [ + ("new", FloatPointAdapter::<[f32; 2], i32>::new(rect)), + ( + "with_coordinate_bits", + FloatPointAdapter::with_coordinate_bits(rect, 29), + ), + ] { + let scale = adapter.dir_scale(); + let inverse = adapter.inv_scale(); + assert!(scale.is_finite() && scale > 0.0, "{name}: dir_scale={scale}"); + assert!( + inverse.is_finite() && inverse > 0.0, + "{name}: inv_scale={inverse}" + ); + + let source = [rect.max_x, 0.0]; + let point = adapter.try_float_to_int(&source).unwrap(); + assert!(point.is_in_safe_range(), "{name}: point={point}"); + let restored = adapter.try_int_to_float(&point).unwrap(); + assert!( + (restored[0] - source[0]).abs() <= inverse, + "{name}: restored={restored:?}" + ); + assert_eq!(restored[1], 0.0); + } +} + +#[test] +fn checked_constructors_reject_scale_with_infinite_reciprocal() { + let rect = FloatRect::new(-1.0_f32, 1.0, -1.0, 1.0).unwrap(); + let scale = 1e-40_f32; + assert!(scale.is_finite() && scale > 0.0); + assert!((1.0 / scale).is_infinite()); + + for (name, result) in [ + ( + "try_with_scale", + FloatPointAdapter::<[f32; 2], i32>::try_with_scale(rect, scale), + ), + ( + "try_with_scale_and_coordinate_bits", + FloatPointAdapter::try_with_scale_and_coordinate_bits(rect, scale, 29), + ), + ] { + assert_eq!( + result.err(), + Some(FloatPointAdapterScaleError::ScaleTooSmall), + "{name}" + ); + } +} + +#[test] +fn supported_boundary_keeps_center_and_scales_finite() { + let radius = f64::MAX_COORDINATE; + let rect = FloatRect::new(-radius, radius, -1.0, 1.0).unwrap(); + assert!(rect.width().is_finite()); + for (name, adapter) in [ + ("new", FloatPointAdapter::<[f64; 2], i32>::new(rect)), + ( + "with_coordinate_bits", + FloatPointAdapter::with_coordinate_bits(rect, 29), + ), + ("with_scale", FloatPointAdapter::with_scale(rect, 1e-145)), + ( + "try_with_scale", + FloatPointAdapter::try_with_scale(rect, 1e-145).unwrap(), + ), + ( + "try_with_scale_and_coordinate_bits", + FloatPointAdapter::try_with_scale_and_coordinate_bits(rect, 1e-145, 29).unwrap(), + ), + ] { + assert_eq!(adapter.offset(), [0.0, 0.0], "{name}"); + let scale = adapter.dir_scale(); + let inverse = adapter.inv_scale(); + assert!(scale.is_finite() && scale > 0.0, "{name}: dir_scale={scale}"); + assert!( + inverse.is_finite() && inverse > 0.0, + "{name}: inv_scale={inverse}" + ); + assert_eq!( + adapter.try_float_to_int(&[0.0, 0.0]), + Ok(IntPoint::ZERO), + "{name}" + ); + assert_eq!( + adapter.try_int_to_float(&IntPoint::ZERO), + Ok([0.0, 0.0]), + "{name}" + ); + } +} + +fn check_extremes(small: F, large: F) { + for radius in [small, large] { + let rect = FloatRect::new(-radius, radius, -radius, radius).unwrap(); + for adapter in [ + FloatPointAdapter::<[F; 2], I>::new(rect), + FloatPointAdapter::with_coordinate_bits(rect, 0), + FloatPointAdapter::with_coordinate_bits(rect, I::BITS - 3), + ] { + assert!(adapter.dir_scale().is_finite() && adapter.dir_scale() > F::ZERO); + assert!(adapter.inv_scale().is_finite() && adapter.inv_scale() > F::ZERO); + assert!(adapter.offset() == [F::ZERO, F::ZERO]); + for source in [[-radius, radius], [F::ZERO, F::ZERO], [radius, -radius]] { + let integer = adapter.try_float_to_int(&source).unwrap(); + assert!(integer.is_in_safe_range()); + let restored = adapter.try_int_to_float(&integer).unwrap(); + assert!(restored[0].is_finite() && restored[1].is_finite()); + assert!((restored[0] - source[0]).abs() <= adapter.inv_scale()); + assert!((restored[1] - source[1]).abs() <= adapter.inv_scale()); + } + } + } + // Halving this span underflows to zero, but the bounds are not a point. + let rect = FloatRect::new(F::ZERO, small, F::ZERO, F::ZERO).unwrap(); + let adapter = FloatPointAdapter::<[F; 2], I>::new(rect); + assert!(adapter.dir_scale() > F::ONE); + let explicit = FloatPointAdapter::<[F; 2], I>::with_scale(rect, F::TWO); + assert!(explicit.dir_scale() == F::TWO); +} + +#[test] +fn rounded_center_preserves_coordinate_budget_for_adjacent_floats() { + fn check(next: F) { + let rect = FloatRect::new(F::ONE, next, -next, -F::ONE).unwrap(); + let adapter = FloatPointAdapter::<[F; 2], i32>::with_coordinate_bits(rect, 0); + // The midpoint rounds to an endpoint even around ordinary values like + // 1.0. Using half the span as radius would exceed the zero-bit budget. + assert!(adapter.offset() == [F::ONE, -F::ONE]); + let expected_scale = F::ONE / (next - F::ONE); + assert!(adapter.dir_scale() == expected_scale); + for (source, expected) in [ + ([F::ONE, -next], IntPoint::new(0, -1)), + ([next, -F::ONE], IntPoint::new(1, 0)), + ] { + let integer = adapter.try_float_to_int(&source).unwrap(); + assert_eq!(integer, expected); + assert!(adapter.try_int_to_float(&integer).unwrap() == source); + } + assert_eq!( + FloatPointAdapter::<[F; 2], i32>::try_with_scale_and_coordinate_bits( + rect, + expected_scale * F::TWO, + 0, + ) + .err(), + Some(FloatPointAdapterScaleError::ScaleTooLarge), + ); + } + check(f32::from_bits(1.0_f32.to_bits() + 1)); + check(f64::from_bits(1.0_f64.to_bits() + 1)); +} + +#[test] +fn finite_scales_cover_subnormals_and_all_coordinate_types() { + check_extremes::(f32::from_bits(1), f32::MAX_COORDINATE); + check_extremes::(f32::from_bits(1), f32::MAX_COORDINATE); + check_extremes::(f32::from_bits(1), f32::MAX_COORDINATE); + check_extremes::(f64::from_bits(1), f64::MAX_COORDINATE); + check_extremes::(f64::from_bits(1), f64::MAX_COORDINATE); + check_extremes::(f64::from_bits(1), f64::MAX_COORDINATE); +} + +fn check_explicit_scale_limits(small: F) { + for rect in [ + FloatRect::zero(), + FloatRect::new(-F::ONE, F::ONE, -F::ONE, F::ONE).unwrap(), + ] { + assert_eq!( + FloatPointAdapter::<[F; 2], i32>::try_with_scale(rect, small).err(), + Some(FloatPointAdapterScaleError::ScaleTooSmall) + ); + assert_eq!( + FloatPointAdapter::<[F; 2], i32>::try_with_scale_and_coordinate_bits(rect, small, 29).err(), + Some(FloatPointAdapterScaleError::ScaleTooSmall) + ); + } + // An explicit finite scale above the automatic power-of-two cap is valid + // when the bounds permit it. Do not silently replace it with that cap. + let rect = FloatRect::new(-small, small, -small, small).unwrap(); + for adapter in [ + FloatPointAdapter::<[F; 2], i32>::with_scale(rect, F::MAX), + FloatPointAdapter::try_with_scale(rect, F::MAX).unwrap(), + FloatPointAdapter::try_with_scale_and_coordinate_bits(rect, F::MAX, 0).unwrap(), + ] { + assert!(adapter.dir_scale() == F::MAX); + assert!(adapter.inv_scale().is_finite() && adapter.inv_scale() > F::ZERO); + } +} + +#[test] +fn explicit_scale_limits_apply_to_both_float_types_and_point_bounds() { + check_explicit_scale_limits(f32::from_bits(1)); + check_explicit_scale_limits(f64::from_bits(1)); +} + +#[test] +#[should_panic(expected = "Invalid adapter scale: ScaleTooSmall")] +fn infallible_explicit_constructor_rejects_infinite_reciprocal() { + FloatPointAdapter::<[f32; 2], i32>::with_scale(FloatRect::new(-1.0, 1.0, -1.0, 1.0).unwrap(), 1e-40); +} + +#[test] +#[should_panic(expected = "Invalid adapter bounds")] +fn infallible_constructor_rejects_bounds_outside_the_contract() { + let rect = FloatRect { + min_x: -1e308, + max_x: 1e308, + min_y: -1.0, + max_y: 1.0, + }; + FloatPointAdapter::<[f64; 2], i32>::with_coordinate_bits(rect, 0); +} + +#[test] +fn checked_constructor_respects_budget_at_coordinate_limit() { + let radius = f64::MAX_COORDINATE; + let rect = FloatRect::new(-radius, radius, -1.0, 1.0).unwrap(); + let scale = libm::exp2(-499.0); + assert!((1.0 / scale).is_finite()); + assert_eq!( + FloatPointAdapter::<[f64; 2], i32>::try_with_scale_and_coordinate_bits(rect, scale, 0).err(), + Some(FloatPointAdapterScaleError::ScaleTooLarge) + ); + let adapter = FloatPointAdapter::<[f64; 2], i32>::with_coordinate_bits(rect, 1); + assert_eq!(adapter.dir_scale(), scale); + assert!(adapter.inv_scale().is_finite()); +} diff --git a/tests/adapter_iter_tests.rs b/tests/adapter_iter_tests.rs new file mode 100644 index 0000000..121de40 --- /dev/null +++ b/tests/adapter_iter_tests.rs @@ -0,0 +1,113 @@ +#![cfg(feature = "core")] + +use i_float::adapter::{FloatPointAdapter, FloatPointAdapterScaleError}; +use i_float::float::number::FloatNumber; +use i_float::float::rect::FloatRectError; +use i_float::int::number::int::IntNumber; +use i_float::int::point::IntPoint; + +fn check_budget() { + let three = F::from_float(3.0); + let points = [[-three, -F::ONE], [three, F::ONE]]; + for bits in [0, 2, I::BITS - 3] { + let adapter = FloatPointAdapter::<[F; 2], I>::with_iter_and_coordinate_bits(points.iter(), bits); + let limit = I::ONE << bits; + for point in &points { + let integer = adapter.try_float_to_int(point).unwrap(); + assert!(-limit <= integer.x && integer.x <= limit); + assert!(-limit <= integer.y && integer.y <= limit); + } + let scale = adapter.dir_scale(); + assert!( + FloatPointAdapter::<[F; 2], I>::try_with_iter_and_scale_and_coordinate_bits( + points.iter(), + scale, + bits, + ) + .is_ok() + ); + assert_eq!( + FloatPointAdapter::<[F; 2], I>::try_with_iter_and_scale_and_coordinate_bits( + points.iter(), + scale * F::from_float(2.0), + bits, + ) + .err(), + Some(FloatPointAdapterScaleError::ScaleTooLarge) + ); + } +} + +#[test] +fn iterator_constructors_enforce_selected_budget() { + check_budget::(); + check_budget::(); + check_budget::(); + check_budget::(); + check_budget::(); + check_budget::(); +} + +type Adapter = FloatPointAdapter<[f64; 2], i32>; + +#[test] +fn empty_and_single_point_inputs_preserve_scale_policy() { + for points in [&[][..], &[[7.0, -3.0]][..]] { + let origin = points.first().copied().unwrap_or([0.0, 0.0]); + let automatic = Adapter::with_iter_and_coordinate_bits(points.iter(), 29); + let fixed = Adapter::try_with_iter_and_scale_and_coordinate_bits(points.iter(), 100.0, 29).unwrap(); + assert_eq!(automatic.dir_scale(), 1.0); + assert_eq!(fixed.dir_scale(), 100.0); + for adapter in [automatic, fixed] { + assert_eq!(adapter.offset(), origin); + assert_eq!(adapter.try_float_to_int(&origin), Ok(IntPoint::ZERO)); + assert_eq!(adapter.try_int_to_float(&IntPoint::ZERO), Ok(origin)); + } + for (scale, error) in [ + (0.0, FloatPointAdapterScaleError::ScaleNonPositive), + (-1.0, FloatPointAdapterScaleError::ScaleNonPositive), + (f64::NAN, FloatPointAdapterScaleError::ScaleNotFinite), + (f64::INFINITY, FloatPointAdapterScaleError::ScaleNotFinite), + (f64::from_bits(1), FloatPointAdapterScaleError::ScaleTooSmall), + ] { + assert_eq!( + Adapter::try_with_iter_and_scale_and_coordinate_bits(points.iter(), scale, 29).err(), + Some(error) + ); + } + } +} + +#[test] +fn invalid_coordinates_are_rejected_anywhere_in_iterator() { + for invalid in [f64::NAN, f64::INFINITY, 2.0 * f64::MAX_COORDINATE] { + for points in [[[invalid, 0.0], [0.0, 0.0]], [[0.0, 0.0], [0.0, invalid]]] { + assert_eq!( + Adapter::try_with_iter_and_scale_and_coordinate_bits(points.iter(), 1.0, 29).err(), + Some(FloatPointAdapterScaleError::InvalidRect( + FloatRectError::CoordinatesOutOfRange + )) + ); + assert!( + std::panic::catch_unwind(|| { Adapter::with_iter_and_coordinate_bits(points.iter(), 29) }) + .is_err() + ); + } + } +} + +#[test] +fn invalid_budget_panics_even_for_empty_input() { + assert!( + std::panic::catch_unwind(|| { + Adapter::with_iter_and_coordinate_bits(core::iter::empty(), i32::BITS - 1) + }) + .is_err() + ); + assert!( + std::panic::catch_unwind(|| { + Adapter::try_with_iter_and_scale_and_coordinate_bits(core::iter::empty(), 1.0, i32::BITS - 1) + }) + .is_err() + ); +} diff --git a/tests/angle_float_tests.rs b/tests/angle_float_tests.rs new file mode 100644 index 0000000..ae6c551 --- /dev/null +++ b/tests/angle_float_tests.rs @@ -0,0 +1,200 @@ +#![cfg(feature = "core")] + +use i_float::float::number::FloatNumber; +use i_float::int::angle::{Angle, AngleDelta}; +use i_float::int::number::wide_int::WideIntNumber; +use std::f64::consts::TAU; + +#[test] +fn angle_radians_support_both_float_types() { + for (angle, unsigned, signed) in [ + (Angle::ZERO, 0.0, 0.0), + (Angle::QUARTER_TURN, TAU / 4.0, TAU / 4.0), + (Angle::HALF_TURN, TAU / 2.0, -TAU / 2.0), + (Angle::THREE_QUARTER_TURN, 3.0 * TAU / 4.0, -TAU / 4.0), + ] { + assert_eq!(angle.to_radians::(), unsigned); + assert_eq!(angle.to_signed_radians::(), signed); + assert_eq!(angle.to_radians::(), unsigned as f32); + assert_eq!(angle.to_signed_radians::(), signed as f32); + let delta = Angle::ZERO.delta_to(angle); + assert_eq!(delta.to_radians::(), signed); + assert_eq!(delta.to_radians::(), signed as f32); + } + let last = Angle::from_bits(u32::MAX); + assert!(last.to_radians::() < TAU); + assert!(last.to_signed_radians::() < 0.0); + assert_eq!(last.to_radians::(), TAU as f32); + let before_half = AngleDelta::from_raw(i32::MAX); + assert!(before_half.to_radians::() < TAU / 2.0); + assert_eq!(before_half.to_radians::(), (TAU / 2.0) as f32); +} + +#[test] +fn angle_arithmetic_wraps_and_chooses_shortest_delta() { + const BACKWARD: AngleDelta = Angle::ZERO.delta_to(Angle::THREE_QUARTER_TURN); + const WRAPPED: Angle = Angle::ZERO.wrapping_add(BACKWARD); + assert_eq!(BACKWARD.raw(), -(1 << 30)); + assert_eq!(WRAPPED, Angle::THREE_QUARTER_TURN); + assert_eq!(Angle::default(), Angle::ZERO); + assert_eq!(AngleDelta::default(), AngleDelta::ZERO); + assert_eq!(Angle::HALF_TURN + Angle::HALF_TURN, Angle::ZERO); + assert_eq!(Angle::THREE_QUARTER_TURN + Angle::HALF_TURN, Angle::QUARTER_TURN); + assert_eq!(Angle::ZERO.delta_to(Angle::HALF_TURN).raw(), i32::MIN); + assert_eq!(Angle::HALF_TURN.delta_to(Angle::ZERO).raw(), i32::MIN); + assert_eq!(Angle::from_bits(u32::MAX).delta_to(Angle::ZERO).raw(), 1); + assert_eq!(Angle::ZERO.delta_to(Angle::from_bits(u32::MAX)).raw(), -1); + let boundaries = [0, 1, (1 << 31) - 1, 1 << 31, (1 << 31) + 1, u32::MAX]; + for from in boundaries.map(Angle::from_bits) { + for target in boundaries.map(Angle::from_bits) { + assert_eq!(from.wrapping_add(from.delta_to(target)), target); + } + } +} + +#[test] +fn radians_reject_non_finite_inputs() { + for radians in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert_eq!(Angle::from_radians(radians), None); + assert_eq!(Angle::from_radians(radians as f32), None); + } + assert_eq!(Angle::from_radians(0.0f32).unwrap().bits(), 0); + assert_eq!(Angle::from_radians(1.0f32), Angle::from_radians(1.0f64)); +} + +#[test] +fn radians_wrap_turns_and_preserve_axes() { + for (radians, bits) in [ + (0.0, 0), + (-0.0, 0), + (TAU, 0), + (-TAU, 0), + (TAU / 4.0, 1 << 30), + (-TAU / 4.0, 3 << 30), + (TAU / 2.0, 1 << 31), + (-TAU / 2.0, 1 << 31), + (3.0 * TAU / 4.0, 3 << 30), + (5.0 * TAU / 4.0, 1 << 30), + (-5.0 * TAU / 4.0, 3 << 30), + ] { + assert_eq!(Angle::from_radians(radians).unwrap().bits(), bits); + } +} + +#[test] +fn radians_round_half_units_away_from_zero() { + let half_unit = TAU / (1_u64 << 33) as f64; + let below_half = f64::from_bits(half_unit.to_bits() - 1); + let above_half = f64::from_bits(half_unit.to_bits() + 1); + for (radians, bits) in [ + (below_half, 0), + (-below_half, 0), + (half_unit, 1), + (-half_unit, u32::MAX), + (above_half, 1), + (-above_half, u32::MAX), + (TAU - half_unit / 2.0, 0), + (-TAU + half_unit / 2.0, 0), + ] { + assert_eq!(Angle::from_radians(radians).unwrap().bits(), bits); + } +} + +#[test] +fn float_coefficients_preserve_axes_norm_and_accuracy() { + let scale = Angle::SIN_COS_SCALE; + for (bits, expected) in [ + (0, (0, scale)), + (1 << 30, (scale, 0)), + (1 << 31, (0, -scale)), + (3 << 30, (-scale, 0)), + ] { + assert_eq!(Angle::from_bits(bits).sin_cos_with_float(), expected); + } + let near_axes = [0u32, 1 << 30, 1 << 31, 3 << 30] + .into_iter() + .flat_map(|axis| [0u32, 1, 2, 3, u32::MAX, u32::MAX - 1].map(|delta| axis.wrapping_add(delta))); + for bits in (0..=u32::MAX).step_by(16_381).chain(near_axes) { + let (sin, cos) = Angle::from_bits(bits).sin_cos_with_float(); + let norm2 = sin as i64 * sin as i64 + cos as i64 * cos as i64; + assert!(norm2 <= scale as i64 * scale as i64, "bits={bits}"); + let radians = bits as f64 * TAU / 4294967296.0; + let (expected_sin, expected_cos) = radians.sin_cos(); + // Independent native-float reference; allow rounding and contraction. + let error = + (sin as f64 - expected_sin * scale as f64).hypot(cos as f64 - expected_cos * scale as f64); + assert!(error < 2.0, "bits={bits}, error={error}"); + } +} + +fn check_atan2() { + assert_eq!(Angle::atan2_with_float(W::ZERO, W::ZERO), None); + for (y, x, bits) in [ + (W::ZERO, W::ONE, 0), + (W::ONE, W::ZERO, 1 << 30), + (W::ZERO, W::MIN, 1 << 31), + (W::MIN, W::ZERO, 3 << 30), + ] { + assert_eq!(Angle::atan2_with_float(y, x).unwrap().bits(), bits); + } + let check = |y: W, x: W| { + let Some(angle) = Angle::atan2_with_float(y, x) else { + assert!(y == W::ZERO && x == W::ZERO); + return; + }; + let bits = angle.bits(); + if y > W::ZERO { + assert!(bits > 0 && bits < 1 << 31); + } else if y < W::ZERO { + assert!(bits > 1 << 31); + } + let expected = y.to_f64().atan2(x.to_f64()); + let actual = bits as f64 * TAU / 4294967296.0; + let difference = (actual - expected).sin().atan2((actual - expected).cos()).abs(); + // Rounding is within half a binary unit, except for the one-unit + // adjustment that preserves the side of an almost parallel ray. + assert!(difference * 4294967296.0 / TAU < 1.00001, "y={y}, x={x}"); + }; + let values = [W::MIN, W::MIN + W::ONE, -W::ONE, W::ZERO, W::ONE, W::MAX]; + for y in values { + for x in values { + check(y, x); + } + } + // Unequal magnitudes at every power of two, in all quadrants. + let mut magnitude = W::MAX; + while magnitude > W::ZERO { + for x in [magnitude, -magnitude] { + for y in [W::ONE, -W::ONE, magnitude, -magnitude] { + check(y, x); + check(x, y); + assert_eq!( + Angle::atan2_with_float(y, x) + .unwrap() + .bits() + .wrapping_add(Angle::atan2_with_float(-y, x).unwrap().bits()), + 0 + ); + } + } + magnitude = magnitude >> 1; + } +} + +#[test] +fn float_atan2_handles_full_wide_range_and_preserves_turn_side() { + check_atan2::(); + check_atan2::(); + check_atan2::(); +} + +#[test] +fn float_number_atan2_returns_signed_radians() { + for (y, x) in [(1.0f64, 1.0), (-1.0, 1.0), (1.0, -1.0), (-1.0, -1.0)] { + assert!((FloatNumber::atan2(y, x) - y.atan2(x)).abs() < 1e-15); + let (y, x) = (y as f32, x as f32); + assert!((FloatNumber::atan2(y, x) - y.atan2(x)).abs() < 1e-6); + } + assert!(FloatNumber::atan2(-0.0f64, 1.0).is_sign_negative()); + assert!(FloatNumber::atan2(-0.0f32, 1.0).is_sign_negative()); +} diff --git a/tests/float_coordinate_range_tests.rs b/tests/float_coordinate_range_tests.rs new file mode 100644 index 0000000..70ea3c3 --- /dev/null +++ b/tests/float_coordinate_range_tests.rs @@ -0,0 +1,239 @@ +#![cfg(feature = "core")] + +use i_float::adapter::{FloatPointAdapter, FloatPointAdapterScaleError}; +use i_float::float::compatible::FloatPointCompatible; +use i_float::float::number::FloatNumber; +use i_float::float::point::FloatPoint; +use i_float::float::rect::{FloatRect, FloatRectError}; +use i_float::float::vector::FloatPointMath; + +fn check_coordinate_limits(limit: F, outside: F, nan: F, infinity: F) { + assert!(F::MAX_COORDINATE == limit); + for value in [-F::MIN_POSITIVE, -F::ZERO, F::MIN_POSITIVE] { + assert!(value.is_in_safe_range()); + } + for x in [-limit, F::ZERO, limit] { + assert!(x.is_in_safe_range()); + for y in [-limit, F::ZERO, limit] { + assert!([x, y].is_in_safe_range()); + assert!(FloatPoint::new(x, y).is_in_safe_range()); + assert!(FloatRect::with_point([x, y]).unwrap().is_in_safe_range()); + } + } + assert!( + FloatRect::new(-limit, limit, -limit, limit) + .unwrap() + .is_in_safe_range() + ); + for invalid in [outside, -outside, nan, infinity, -infinity] { + assert!(!invalid.is_in_safe_range()); + assert!(![invalid, F::ZERO].is_in_safe_range()); + assert!(![F::ZERO, invalid].is_in_safe_range()); + for bounds in [ + [invalid, limit, -limit, limit], + [-limit, invalid, -limit, limit], + [-limit, limit, invalid, limit], + [-limit, limit, -limit, invalid], + ] { + let [min_x, max_x, min_y, max_y] = bounds; + assert_eq!( + FloatRect::new(min_x, max_x, min_y, max_y).err(), + Some(FloatRectError::CoordinatesOutOfRange) + ); + let raw = FloatRect { + min_x, + max_x, + min_y, + max_y, + }; + assert!(!raw.is_in_safe_range()); + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + FloatPointAdapter::<[F; 2], i32>::new(raw) + })) + .is_err() + ); + for result in [ + FloatPointAdapter::<[F; 2], i32>::try_with_scale(raw, F::ONE), + FloatPointAdapter::try_with_scale_and_coordinate_bits(raw, F::ONE, 29), + ] { + assert_eq!( + result.err(), + Some(FloatPointAdapterScaleError::InvalidRect( + FloatRectError::CoordinatesOutOfRange + )) + ); + } + } + for points in [ + [[invalid, F::ZERO], [F::ZERO, F::ZERO]], + [[F::ZERO, F::ZERO], [F::ZERO, invalid]], + ] { + assert_eq!( + FloatRect::with_points(&points).err(), + Some(FloatRectError::CoordinatesOutOfRange) + ); + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + FloatPointAdapter::<[F; 2], i32>::with_iter(points.iter()) + })) + .is_err() + ); + assert_eq!( + FloatPointAdapter::<[F; 2], i32>::with_iter_and_scale_checked(points.iter(), F::ONE).err(), + Some(FloatPointAdapterScaleError::InvalidRect( + FloatRectError::CoordinatesOutOfRange + )) + ); + } + #[cfg(debug_assertions)] + { + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + FloatPoint::new(invalid, F::ZERO) + })) + .is_err() + ); + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + FloatPoint::from_point([F::ZERO, invalid]) + })) + .is_err() + ); + } + } + // Differences can be larger than input coordinates while their products fit. + let a = FloatPoint::new(-limit, -limit); + let b = FloatPoint::new(limit, limit); + let difference = b - a; + assert!(difference.sqr_length().is_finite()); + assert!(difference.dot_product(difference).is_finite()); + assert!(difference.cross_product(-difference).is_finite()); + assert!(a.midpoint(b).x == F::ZERO); + assert!(b.midpoint(b).x == limit); +} + +#[test] +fn rectangle_builders_and_mutators_preserve_valid_bounds() { + fn check() { + let limit = F::MAX_COORDINATE; + assert_eq!( + FloatRect::new(F::ONE, F::ZERO, F::ZERO, F::ONE).err(), + Some(FloatRectError::InvalidBounds) + ); + assert_eq!( + FloatRect::new(F::ZERO, F::ONE, F::ONE, F::ZERO).err(), + Some(FloatRectError::InvalidBounds) + ); + let invalid = FloatRect { + min_x: F::ONE, + max_x: F::ZERO, + min_y: F::ZERO, + max_y: F::ONE, + }; + assert!(!invalid.is_in_safe_range()); + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { + FloatPointAdapter::<[F; 2], i32>::new(invalid) + })) + .is_err() + ); + assert_eq!( + FloatRect::with_rects(invalid, FloatRect::zero()).err(), + Some(FloatRectError::InvalidBounds) + ); + assert_eq!( + FloatRect::with_optional_rects(None, Some(invalid)).err(), + Some(FloatRectError::InvalidBounds) + ); + + let empty: [[F; 2]; 0] = []; + assert!(FloatRect::with_points(&empty).unwrap().is_none()); + let mut optional = None; + FloatRect::optional_add_point(&mut optional, &[F::ZERO, F::ONE]).unwrap(); + FloatRect::optional_add_point(&mut optional, &[-limit, limit]).unwrap(); + let mut rect = optional.unwrap(); + assert!(rect.min_x == -limit && rect.max_y == limit); + let before = rect; + assert_eq!(rect.add_offset(limit), Err(FloatRectError::CoordinatesOutOfRange)); + assert_eq!(rect.add_offset(-limit), Err(FloatRectError::InvalidBounds)); + assert_eq!( + rect.add_point(&[limit * F::TWO, F::ZERO]), + Err(FloatRectError::CoordinatesOutOfRange) + ); + assert!( + rect.min_x == before.min_x + && rect.max_x == before.max_x + && rect.min_y == before.min_y + && rect.max_y == before.max_y + ); + rect.add_point(&[limit, -limit]).unwrap(); + assert!(rect.min_x == -limit && rect.max_x == limit && rect.min_y == -limit && rect.max_y == limit); + rect.add_offset(-limit).unwrap(); + assert!( + rect.min_x == F::ZERO && rect.max_x == F::ZERO && rect.min_y == F::ZERO && rect.max_y == F::ZERO + ); + } + check::(); + check::(); +} + +#[test] +fn f32_coordinate_limits_are_inclusive() { + let limit = libm::exp2f(60.0); + check_coordinate_limits( + limit, + f32::from_bits(limit.to_bits() + 1), + f32::NAN, + f32::INFINITY, + ); +} + +#[test] +fn f64_coordinate_limits_are_inclusive() { + let limit = libm::exp2(500.0); + check_coordinate_limits( + limit, + f64::from_bits(limit.to_bits() + 1), + f64::NAN, + f64::INFINITY, + ); +} + +fn check_normalization() { + for x in [F::MIN_POSITIVE.sqrt(), F::MAX_COORDINATE] { + let point = FloatPoint::new(x, F::ZERO); + let direction = point.normalize(); + assert!(direction.x == F::ONE && direction.y == F::ZERO); + let direction = FloatPointMath::normalize(&[x, F::ZERO]); + assert!(direction[0] == F::ONE && direction[1] == F::ZERO); + } +} + +#[test] +fn normalization_supports_the_documented_boundaries() { + check_normalization::(); + check_normalization::(); +} + +#[cfg(debug_assertions)] +#[test] +fn normalization_debug_checks_reject_unsupported_squared_lengths() { + fn check() { + for x in [F::ZERO, F::MIN_POSITIVE, F::MIN_POSITIVE.sqrt() * F::HALF, F::MAX] { + // Public fields allow constructing a vector to test normalization's + // own check independently of the point constructor. + let point = FloatPoint { x, y: F::ZERO }; + assert!(std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| point.normalize())).is_err()); + assert!( + std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| FloatPointMath::normalize(&[ + x, + F::ZERO + ]))) + .is_err() + ); + } + } + check::(); + check::(); +} diff --git a/tests/float_normalize_tests.rs b/tests/float_normalize_tests.rs new file mode 100644 index 0000000..4aa993f --- /dev/null +++ b/tests/float_normalize_tests.rs @@ -0,0 +1,72 @@ +#![cfg(feature = "core")] +use i_float::int::{unit_vector::UnitIntVector, vector::IntVector}; + +macro_rules! check_normalize { + ($name:ident, $int:ty, $wide:ty) => { + #[test] + fn $name() { + let scale = UnitIntVector::<$int>::DENOMINATOR; + let check = |x: $wide, y: $wide| { + let result = UnitIntVector::<$int>::normalize_with_float(IntVector::new(x, y)); + if x == 0 && y == 0 { + assert!(result.is_none()); + return; + } + let unit = result.expect("every nonzero wide vector can be normalized"); + let (a, b) = (unit.x() as $wide, unit.y() as $wide); + assert!(a != 0 || b != 0); + let length = (x as f64).hypot(y as f64); + let tolerance = 1.0 / scale as f64 + 8.0 * f64::EPSILON; + let norm = (a as f64 / scale as f64).hypot(b as f64 / scale as f64); + assert!((norm - 1.0).abs() <= 2.0 * tolerance); + assert!((a as f64 / scale as f64 - x as f64 / length).abs() <= tolerance); + assert!((b as f64 / scale as f64 - y as f64 / length).abs() <= tolerance); + if x != <$wide>::MIN && y != <$wide>::MIN { + let opposite = + UnitIntVector::<$int>::normalize_with_float(IntVector::new(-x, -y)).unwrap(); + assert_eq!((opposite.x(), opposite.y()), (-unit.x(), -unit.y())); + } + }; + let values = [ + <$wide>::MIN, + <$wide>::MIN + 1, + -scale, + -4, + -1, + 0, + 1, + 3, + scale, + <$wide>::MAX, + ]; + for x in values { + for y in values { + check(x, y); + } + } + for power in 0..(<$wide>::BITS - 1) { + let n = (1 as $wide) << power; + for x in [n - 1, n, n + 1] { + for y in [1, n / 3, n - 1, n] { + check(x, y); + check(-x, y); + } + } + } + } + }; +} +check_normalize!(normalize_i16, i16, i32); +check_normalize!(normalize_i32, i32, i64); +check_normalize!(normalize_i64, i64, i128); + +#[test] +fn near_axis_roundoff_does_not_reject_nonzero_vector() { + let unit = UnitIntVector::::normalize_with_float(IntVector::new(1 << 30, 1)).unwrap(); + let scale = UnitIntVector::::DENOMINATOR; + let (x, y) = (unit.x() as i128, unit.y() as i128); + // f64 loses the small squared component. A slight overshoot is intentional. + assert!(x * x + y * y > scale * scale); + assert!((x as f64 / scale as f64 - 1.0).abs() <= f64::EPSILON); + assert!(y > 0); +} diff --git a/tests/float_rect_error_tests.rs b/tests/float_rect_error_tests.rs new file mode 100644 index 0000000..82e6df1 --- /dev/null +++ b/tests/float_rect_error_tests.rs @@ -0,0 +1,30 @@ +#![cfg(feature = "core")] + +use core::error::Error; +use i_float::float::rect::{FloatRect, FloatRectError}; + +#[test] +fn rectangle_errors_support_display_and_error_from_public_api() { + let cases = [ + ( + FloatRect::new(f64::NAN, 1.0, 0.0, 1.0).unwrap_err(), + FloatRectError::CoordinatesOutOfRange, + "A coordinate is non-finite or exceeds the supported absolute limit", + ), + ( + FloatRect::new(1.0, 0.0, 0.0, 1.0).unwrap_err(), + FloatRectError::InvalidBounds, + "A minimum bound is greater than its maximum", + ), + ]; + + for (error, expected_variant, expected_message) in cases { + assert_eq!(error, expected_variant); + assert_eq!(error.to_string(), expected_message); + + let source: &(dyn Error + 'static) = &error; + assert_eq!(source.to_string(), expected_message); + assert!(source.source().is_none()); + assert_eq!(source.downcast_ref::(), Some(&error)); + } +} diff --git a/tests/float_rounding_tests.rs b/tests/float_rounding_tests.rs new file mode 100644 index 0000000..f50fd18 --- /dev/null +++ b/tests/float_rounding_tests.rs @@ -0,0 +1,74 @@ +#![cfg(feature = "core")] + +use i_float::float::number::FloatNumber; + +fn assert_small_round(value: F, expected: i16) { + assert_eq!(value.to_round_i16(), expected); + assert_eq!(value.to_round_i32(), i32::from(expected)); + assert_eq!(value.to_round_i64(), i64::from(expected)); + assert_eq!(value.to_round_i128(), i128::from(expected)); + assert_eq!(value.to_round_usize(), expected.max(0) as usize); +} + +macro_rules! rounding_tests { + ($module:ident, $float:ty, $boundary:expr) => { + mod $module { + use super::*; + + #[test] + fn rounds_half_and_neighbors_away_from_zero() { + for (half, below, rounded) in [(0.5, 0, 1), (1.5, 1, 2), (2.5, 2, 3)] { + let half = half as $float; + for (value, expected) in [ + (<$float>::from_bits(half.to_bits() - 1), below), + (half, rounded), + (<$float>::from_bits(half.to_bits() + 1), rounded), + ] { + assert_small_round(value, expected); + assert_small_round(-value, -expected); + } + } + assert_small_round(0.0 as $float, 0); + assert_small_round(-0.0 as $float, 0); + } + + #[test] + fn preserves_integers_at_precision_boundary() { + for expected in ($boundary - 2)..=($boundary + 3) { + for signed in [expected, -expected] { + let value = signed as $float; + assert_eq!(value.to_round_i64(), signed); + assert_eq!(value.to_round_i128(), i128::from(signed)); + if let Ok(expected_i32) = i32::try_from(signed) { + assert_eq!(value.to_round_i32(), expected_i32); + } + if let Ok(expected_usize) = usize::try_from(signed) { + assert_eq!(value.to_round_usize(), expected_usize); + } + } + } + } + + #[test] + fn preserves_saturating_cast_behavior() { + for value in [<$float>::INFINITY, <$float>::MAX] { + assert_eq!(value.to_round_i16(), i16::MAX); + assert_eq!(value.to_round_i32(), i32::MAX); + assert_eq!(value.to_round_i64(), i64::MAX); + assert_eq!(value.to_round_i128(), i128::MAX); + assert_eq!(value.to_round_usize(), usize::MAX); + + assert_eq!((-value).to_round_i16(), i16::MIN); + assert_eq!((-value).to_round_i32(), i32::MIN); + assert_eq!((-value).to_round_i64(), i64::MIN); + assert_eq!((-value).to_round_i128(), i128::MIN); + assert_eq!((-value).to_round_usize(), 0); + } + assert_small_round(<$float>::NAN, 0); + } + } + }; +} + +rounding_tests!(f32_rounding, f32, 1_i64 << 23); +rounding_tests!(f64_rounding, f64, 1_i64 << 52); diff --git a/tests/float_unit_vector_tests.rs b/tests/float_unit_vector_tests.rs new file mode 100644 index 0000000..d60f970 --- /dev/null +++ b/tests/float_unit_vector_tests.rs @@ -0,0 +1,68 @@ +#![cfg(feature = "core")] +use i_float::int::unit_vector::UnitIntVector; + +macro_rules! float_conversion { + ($name:ident, $int:ty, $wide:ty) => { + #[test] + fn $name() { + let s = UnitIntVector::<$int>::DENOMINATOR; + for (x, y) in [ + (1.0, 0.0), + (0.0, -1.0), + (0.3, 0.4), + (-0.3, -0.4), + (0.5, 0.5), + ] { + let unit = UnitIntVector::<$int>::try_from_float(x, y).unwrap(); + assert_eq!(unit, UnitIntVector::<$int>::from_float_unchecked(x, y)); + let (a, b) = (unit.x() as $wide, unit.y() as $wide); + assert!(a * a + b * b <= s * s); + assert!((a as f64 / s as f64 - x).abs() <= 1.0 / s as f64); + assert!((b as f64 / s as f64 - y).abs() <= 1.0 / s as f64); + let opposite = UnitIntVector::<$int>::try_from_float(-x, -y).unwrap(); + assert_eq!((opposite.x(), opposite.y()), (-unit.x(), -unit.y())); + } + for step in 0..10000 { + let angle = step as f64 * core::f64::consts::TAU / 10000.0; + let (sin, cos) = libm::sincos(angle); + let inward = 1.0 - 8.0 * f64::EPSILON; + let unit = UnitIntVector::<$int>::try_from_float(cos * inward, sin * inward).unwrap(); + assert_eq!( + unit, + UnitIntVector::<$int>::from_float_unchecked(cos * inward, sin * inward) + ); + let (a, b) = (unit.x() as $wide, unit.y() as $wide); + assert!(a * a + b * b <= s * s); + } + let axis = UnitIntVector::<$int>::try_from_float(1.0_f32, 0.0).unwrap(); + assert_eq!(axis.x() as $wide, s); + assert_eq!(axis.y(), 0); + for (x, y) in [ + (0.0, 0.0), + (1.0, 1.0), + (1.01, 0.0), + (f64::NAN, 0.0), + (0.0, f64::INFINITY), + ] { + assert!(UnitIntVector::<$int>::try_from_float(x, y).is_none()); + } + for non_finite in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + for finite in [0.0, 0.5, -1.0] { + for (x, y) in [(non_finite, finite), (finite, non_finite)] { + assert!(UnitIntVector::<$int>::try_from_float(x, y).is_none()); + assert!(UnitIntVector::<$int>::try_from_float(x as f32, y as f32).is_none()); + } + } + } + } + }; +} +float_conversion!(float_i16, i16, i32); +float_conversion!(float_i32, i32, i64); +float_conversion!(float_i64, i64, i128); + +#[test] +fn i64_checks_length_below_float_resolution() { + // In f64, 1 + 2^-60 rounds to 1. The exact fixed-scale norm is larger. + assert!(UnitIntVector::::try_from_float(1.0, 1.0 / (1_u64 << 30) as f64).is_none()); +} diff --git a/tests/int_point_macro_tests.rs b/tests/int_point_macro_tests.rs new file mode 100644 index 0000000..de4a336 --- /dev/null +++ b/tests/int_point_macro_tests.rs @@ -0,0 +1,7 @@ +#![cfg(feature = "core")] + +#[test] +fn exported_macro_does_not_require_an_int_point_import() { + let point = i_float::int_pnt!(3_i64, -4_i64); + assert_eq!((point.x, point.y), (3_i64, -4_i64)); +}