cranpose-ui 0.1.90

//! Wear's `setLuminance`: one colour taken to another lightness.
//!
//! Wear Material 3 builds several of its roles by moving one scheme colour to a
//! fixed L\*: the scroll indicator's thumb is `onBackground` at L\* 80 and its
//! track is the same colour at L\* 20 (`ScrollIndicatorDefaults`, which calls
//! `ColorScheme.getDefaultScrollIndicatorColorsCached` and, on a miss, builds
//! `ScrollIndicatorColors(setLuminance(fromToken(OnBackground), 80f),
//! setLuminance(fromToken(OnBackground), 20f))`).
//!
//! **The operation is not a CIE L\*a\*b\* round trip.** That is the guessable
//! answer and the wrong one: substituting L\* in Lab keeps Lab's own a\* and
//! b\*, and Wear keeps CAM16's hue and chroma instead. `setLuminance` is
//!
//! ```text
//! val cam = Cam.fromInt(color.toArgb())
//! Cam.getInt(cam.hue, cam.chroma, lstar)
//! ```
//!
//! — the CAM16 appearance model forward, then the HCT solver backwards. The two
//! agree at some lightnesses and not at others, and the difference is a whole
//! level of one channel where it bites: `#DFF6FF` at L\* 20 is `(30, 51, 58)`
//! through CAM16 and `(31, 51, 58)` through Lab, while at L\* 80 both give
//! `(180, 202, 211)`. A port that checks only the lighter of the two constants
//! concludes its arithmetic is right.
//!
//! # Precision
//!
//! The two halves work at different widths, and this matters at the byte:
//! `Cam.fromInt` is `float` throughout (its matrices are `float[][]`, its
//! transcendentals are `(float) Math.pow((double) x, (double) y)`), and its
//! `hue` and `chroma` are widened to `double` only when they are handed to
//! `HctSolver`, which is `double` throughout. This module keeps that split —
//! `f32` in [`cam16_hue_chroma`], `f64` in [`hct_solve`] — because collapsing it
//! to one width moves the answer by a level on some inputs.
//!
//! Validated against the shipped `compose-material3-1.6.2.aar`: 3024
//! `(source, L*)` pairs spanning both of the solver's paths and both
//! short-circuits, generated by calling the AAR's own `Cam` and `HctSolver` on a
//! JVM. A representative subset is asserted in this module's tests.

use crate::modifier::Color;

/// A colour moved to a given CIE L\*, keeping its CAM16 hue and chroma.
///
/// This is `androidx.wear.compose.material3.setLuminance`. `l_star` is in
/// 0..=100; outside a hair of that range Wear returns the neutral grey of that
/// lightness rather than solving, which is what keeps pure black and pure white
/// from being asked for a hue they do not have.
///
/// The result is **opaque**, as the platform's is: `setLuminance` goes through
/// an ARGB int whose alpha byte the solver writes as `0xFF`.
pub fn set_luminance(color: Color, l_star: f32) -> Color {
    if !SOLVED_RANGE.contains(&(l_star as f64)) {
        return color_from_argb(argb_from_lstar(l_star as f64));
    }
    let (hue, chroma) = cam16_hue_chroma(argb_from_color(color));
    color_from_argb(hct_solve(hue as f64, chroma as f64, l_star as f64))
}

/// A colour's CAM16 hue in degrees and its chroma, in the default viewing
/// conditions — `Cam.fromInt`, which is `Cam.fromIntInFrame(argb,
/// Frame.Default)`.
///
/// Public because the pair is the whole of what `setLuminance` carries across,
/// and a caller building its own Wear role wants to see it rather than infer it.
pub fn cam16_hue_chroma(argb: u32) -> (f32, f32) {
    let xyz = xyz_from_argb(argb);
    let r_t = xyz[0] * XYZ_TO_CAM16RGB[0][0]
        + xyz[1] * XYZ_TO_CAM16RGB[0][1]
        + xyz[2] * XYZ_TO_CAM16RGB[0][2];
    let g_t = xyz[0] * XYZ_TO_CAM16RGB[1][0]
        + xyz[1] * XYZ_TO_CAM16RGB[1][1]
        + xyz[2] * XYZ_TO_CAM16RGB[1][2];
    let b_t = xyz[0] * XYZ_TO_CAM16RGB[2][0]
        + xyz[1] * XYZ_TO_CAM16RGB[2][1]
        + xyz[2] * XYZ_TO_CAM16RGB[2][2];

    let r_d = FRAME.rgb_d[0] * r_t;
    let g_d = FRAME.rgb_d[1] * g_t;
    let b_d = FRAME.rgb_d[2] * b_t;

    let r_af = powf32(FRAME.fl * r_d.abs() / 100.0, 0.42);
    let g_af = powf32(FRAME.fl * g_d.abs() / 100.0, 0.42);
    let b_af = powf32(FRAME.fl * b_d.abs() / 100.0, 0.42);

    let r_a = signum_f32(r_d) * 400.0 * r_af / (r_af + 27.13);
    let g_a = signum_f32(g_d) * 400.0 * g_af / (g_af + 27.13);
    let b_a = signum_f32(b_d) * 400.0 * b_af / (b_af + 27.13);

    let a = (11.0 * r_a + -12.0 * g_a + b_a) / 11.0;
    let b = (r_a + g_a - 2.0 * b_a) / 9.0;
    let u = (20.0 * r_a + 20.0 * g_a + 21.0 * b_a) / 20.0;
    let p2 = (40.0 * r_a + 20.0 * g_a + b_a) / 20.0;

    let atan_degrees = (b as f64).atan2(a as f64) as f32 * 180.0 / PI_F32;
    let hue = if atan_degrees < 0.0 {
        atan_degrees + 360.0
    } else if atan_degrees >= 360.0 {
        atan_degrees - 360.0
    } else {
        atan_degrees
    };

    let ac = p2 * FRAME.nbb;
    let j = 100.0 * powf32(ac / FRAME.aw, FRAME.c * FRAME.z);

    // The eccentricity term is stated on a hue wrapped past its own origin, and
    // the comparison is made after widening to double, as `Cam.kt` does.
    let hue_prime = if (hue as f64) < 20.14 {
        hue + 360.0
    } else {
        hue
    };
    let e_hue = 0.25 * (((hue_prime * PI_F32 / 180.0 + 2.0) as f64).cos() as f32 + 3.8);
    let p1 = 3846.1538 * e_hue * FRAME.nc * FRAME.ncb;
    let t = p1 * ((a * a + b * b) as f64).sqrt() as f32 / (u + 0.305);
    let alpha = powf32(t, 0.9) * powf32(1.64 - powf32(0.29, FRAME.n), 0.73);
    let chroma = alpha * ((j / 100.0) as f64).sqrt() as f32;
    (hue, chroma)
}

/// `HctSolver.solveToInt`: the sRGB colour with this CAM16 hue and chroma at
/// this L\*, or the closest one the cube holds.
///
/// Two paths, and a Wear palette reaches both — 41% of a 3024-row sample fell
/// through to the second. Newton's method on lightness answers whenever the
/// requested chroma is inside the gamut at that lightness; when it is not, the
/// answer is on the cube's surface and is found by bisecting first to a segment
/// between two vertices and then along that segment, one 8-bit plane at a time.
pub fn hct_solve(hue_degrees: f64, chroma: f64, l_star: f64) -> u32 {
    if chroma < 1.0e-4 || !SOLVED_RANGE.contains(&l_star) {
        return argb_from_lstar(l_star);
    }
    let hue_degrees = sanitize_degrees(hue_degrees);
    let hue_radians = hue_degrees / 180.0 * std::f64::consts::PI;
    let y = y_from_lstar(l_star);
    if let Some(exact) = find_result_by_j(hue_radians, chroma, y) {
        return exact;
    }
    let linrgb = bisect_to_limit(y, hue_radians);
    argb_from_linrgb(linrgb)
}

// -- the appearance model's viewing conditions -------------------------------

/// `Frame.Default`, which is the only frame `Cam.fromInt` and `HctSolver` ever
/// use.
///
/// AOSP derives these from the D65 white point at an adapting luminance of
/// `(200 / PI) * yFromLstar(50) / 100`, a background L\* of 50, an average
/// surround and no illuminant discounting. They are held here as the exact
/// `f32` values the shipped `compose-material3-1.6.2` object carries — read out
/// of the AAR by reflection — rather than re-derived, so that a rounding
/// difference in the derivation cannot move a channel.
struct ViewingConditions {
    n: f32,
    aw: f32,
    nbb: f32,
    ncb: f32,
    c: f32,
    nc: f32,
    rgb_d: [f32; 3],
    fl: f32,
    z: f32,
}

const FRAME: ViewingConditions = ViewingConditions {
    n: 0.184_186_52,
    aw: 29.981_003,
    nbb: 1.016_919_3,
    ncb: 1.016_919_3,
    c: 0.690_000_06,
    nc: 1.0,
    rgb_d: [1.021_177_8, 0.986_307_7, 0.933_960_5],
    fl: 0.388_481_47,
    z: 1.909_169_6,
};

/// `CamUtils.XYZ_TO_CAM16RGB`.
const XYZ_TO_CAM16RGB: [[f32; 3]; 3] = [
    [0.401_288, 0.650_173, -0.051_461],
    [-0.250_268, 1.204_414, 0.045_854],
    [-0.002_079, 0.048_952, 0.953_127],
];

/// `CamUtils.SRGB_TO_XYZ`. Double, even though what it multiplies and what it
/// produces are floats.
const SRGB_TO_XYZ: [[f64; 3]; 3] = [
    [0.412_338_95, 0.357_620_64, 0.180_510_42],
    [0.2126, 0.7152, 0.0722],
    [0.019_321_41, 0.119_163_82, 0.950_344_78],
];

/// `CamUtils.XYZ_TO_SRGB`.
const XYZ_TO_SRGB: [[f64; 3]; 3] = [
    [
        3.241_377_479_238_868_5,
        -1.537_665_240_285_185_1,
        -0.498_853_668_462_680_53,
    ],
    [
        -0.969_145_251_300_532_1,
        1.875_885_345_106_787_2,
        0.041_565_856_169_120_61,
    ],
    [
        0.055_620_936_896_913_05,
        -0.203_955_245_647_421_23,
        1.057_179_911_122_033_5,
    ],
];

/// `CamUtils.WHITE_POINT_D65`.
const WHITE_POINT_D65: [f32; 3] = [95.047, 100.0, 108.883];

/// `HctSolver.Y_FROM_LINRGB`.
const Y_FROM_LINRGB: [f64; 3] = [0.2126, 0.7152, 0.0722];

/// `HctSolver.SCALED_DISCOUNT_FROM_LINRGB`.
const SCALED_DISCOUNT_FROM_LINRGB: [[f64; 3]; 3] = [
    [
        0.001_200_833_568_784_504,
        0.002_389_694_492_170_889,
        2.795_742_885_861_124e-4,
    ],
    [
        5.891_086_651_375_999e-4,
        0.002_978_550_257_343_875_8,
        3.270_666_104_008_398e-4,
    ],
    [
        1.014_669_249_164_057_2e-4,
        5.364_214_359_186_694e-4,
        0.003_297_940_177_071_207_6,
    ],
];

/// `HctSolver.LINRGB_FROM_SCALED_DISCOUNT`.
const LINRGB_FROM_SCALED_DISCOUNT: [[f64; 3]; 3] = [
    [
        1373.2198709594231,
        -1100.4251190754821,
        -7.278_681_089_101_213,
    ],
    [
        -271.815_969_077_903,
        559.658_046_594_073_3,
        -32.460_474_827_911_94,
    ],
    [
        1.962_289_959_966_566_6,
        -57.173_814_538_844_006,
        308.723_319_781_238_5,
    ],
];

/// The linear-RGB coordinate of every boundary between two 8-bit channel
/// values: `CRITICAL_PLANES[i]` delinearizes to exactly `i + 0.5` of 255.
///
/// `HctSolver` bisects along a segment by *plane index* rather than by
/// coordinate, so that each step of the search lands on a value the framebuffer
/// can actually hold. The table is `HctSolver.CRITICAL_PLANES` verbatim;
/// `the_critical_planes_are_the_eight_bit_boundaries` re-derives every entry
/// from [`true_delinearized`] so the numbers are not opaque.
const CRITICAL_PLANES: [f64; 255] = [
    0.015176349177441876,
    0.045529047532325624,
    0.07588174588720938,
    0.10623444424209313,
    0.13658714259697685,
    0.16693984095186062,
    0.19729253930674434,
    0.2276452376616281,
    0.2579979360165119,
    0.28835063437139563,
    0.3188300904430532,
    0.350925934958123,
    0.3848314933096426,
    0.42057480301049466,
    0.458183274052838,
    0.4976837250274023,
    0.5391024159806381,
    0.5824650784040898,
    0.6277969426914107,
    0.6751227633498623,
    0.7244668422128921,
    0.775853049866786,
    0.829304845476233,
    0.8848452951698498,
    0.942497089126609,
    1.0022825574869039,
    1.0642236851973577,
    1.1283421258858297,
    1.1946592148522128,
    1.2631959812511864,
    1.3339731595349034,
    1.407011200216447,
    1.4823302800086415,
    1.5599503113873272,
    1.6398909516233677,
    1.7221716113234105,
    1.8068114625156377,
    1.8938294463134073,
    1.9832442801866852,
    2.075074464868551,
    2.1693382909216234,
    2.2660538449872063,
    2.36523901573795,
    2.4669114995532007,
    2.5710888059345764,
    2.6777882626779785,
    2.7870270208169257,
    2.898822059350997,
    3.0131901897720907,
    3.1301480604002863,
    3.2497121605402226,
    3.3718988244681087,
    3.4967242352587946,
    3.624204428461639,
    3.754355295633311,
    3.887192587735158,
    4.022731918402185,
    4.160988767090289,
    4.301978482107941,
    4.445716283538092,
    4.592217266055746,
    4.741496401646282,
    4.893568542229298,
    5.048448422192488,
    5.20615066083972,
    5.3666897647573375,
    5.5300801301023865,
    5.696336044816294,
    5.865471690767354,
    6.037501145825082,
    6.212438385869475,
    6.390297286737924,
    6.571091626112461,
    6.7548350853498045,
    6.941541251256611,
    7.131223617812143,
    7.323895587840543,
    7.5195704746346665,
    7.7182615035334345,
    7.919981813454504,
    8.124744458384042,
    8.332562408825165,
    8.543448553206703,
    8.757415699253682,
    8.974476575321063,
    9.194643831691977,
    9.417930041841839,
    9.644347703669503,
    9.873909240696694,
    10.106627003236781,
    10.342513269534024,
    10.58158024687427,
    10.8238400726681,
    11.069304815507364,
    11.317986476196008,
    11.569896988756009,
    11.825048221409341,
    12.083451977536606,
    12.345119996613247,
    12.610063955123938,
    12.878295467455942,
    13.149826086772048,
    13.42466730586372,
    13.702830557985108,
    13.984327217668513,
    14.269168601521828,
    14.55736596900856,
    14.848930523210871,
    15.143873411576273,
    15.44220572664832,
    15.743938506781891,
    16.04908273684337,
    16.35764934889634,
    16.66964922287304,
    16.985093187232053,
    17.30399201960269,
    17.62635644741625,
    17.95219714852476,
    18.281524751807332,
    18.614349837764564,
    18.95068293910138,
    19.290534541298456,
    19.633915083172692,
    19.98083495742689,
    20.331304511189067,
    20.685334046541502,
    21.042933821039977,
    21.404114048223256,
    21.76888489811322,
    22.137256497705877,
    22.50923893145328,
    22.884842241736916,
    23.264076429332462,
    23.6469514538663,
    24.033477234264016,
    24.42366364919083,
    24.817520537484558,
    25.21505769858089,
    25.61628489293138,
    26.021211842414342,
    26.429848230738664,
    26.842203703840827,
    27.258287870275353,
    27.678110301598522,
    28.10168053274597,
    28.529008062403893,
    28.96010235337422,
    29.39497283293396,
    29.83362889318845,
    30.276079891419332,
    30.722335150426627,
    31.172403958865512,
    31.62629557157785,
    32.08401920991837,
    32.54558406207592,
    33.010999283389665,
    33.4802739966603,
    33.953417292456834,
    34.430438229418264,
    34.911345834551085,
    35.39614910352207,
    35.88485700094671,
    36.37747846067349,
    36.87402238606382,
    37.37449765026789,
    37.87891309649659,
    38.38727753828926,
    38.89959975977785,
    39.41588851594697,
    39.93615253289054,
    40.460400508064545,
    40.98864111053629,
    41.520882981230194,
    42.05713473317016,
    42.597404951718396,
    43.141702194811224,
    43.6900349931913,
    44.24241185063697,
    44.798841244188324,
    45.35933162437017,
    45.92389141541209,
    46.49252901546552,
    47.065252796817916,
    47.64207110610409,
    48.22299226451468,
    48.808024568002054,
    49.3971762874833,
    49.9904556690408,
    50.587870934119984,
    51.189430279724725,
    51.79514187861014,
    52.40501387947288,
    53.0190544071392,
    53.637271562750364,
    54.259673423945976,
    54.88626804504493,
    55.517063457223934,
    56.15206766869424,
    56.79128866487574,
    57.43473440856916,
    58.08241284012621,
    58.734331877617365,
    59.39049941699807,
    60.05092333227251,
    60.715611475655585,
    61.38457167773311,
    62.057811747619894,
    62.7353394731159,
    63.417162620860914,
    64.10328893648692,
    64.79372614476921,
    65.48848194977529,
    66.18756403501224,
    66.89098006357258,
    67.59873767827808,
    68.31084450182222,
    69.02730813691093,
    69.74813616640164,
    70.47333615344107,
    71.20291564160104,
    71.93688215501312,
    72.67524319850172,
    73.41800625771542,
    74.16517879925733,
    74.9167682708136,
    75.67278210128072,
    76.43322770089146,
    77.1981124613393,
    77.96744375590167,
    78.74122893956174,
    79.51947534912904,
    80.30219030335869,
    81.08938110306934,
    81.88105503125999,
    82.67721935322541,
    83.4778813166706,
    84.28304815182372,
    85.09272707154808,
    85.90692527145302,
    86.72564993000343,
    87.54890820862819,
    88.3767072518277,
    89.2090541872801,
    90.04595612594655,
    90.88742016217518,
    91.73345337380438,
    92.58406282226491,
    93.43925555268066,
    94.29903859396902,
    95.16341895893969,
    96.03240364439274,
    96.9059996312159,
    97.78421388448044,
    98.6670533535366,
    99.55452497210776,
];

/// The lightness band Wear actually solves in. Outside it `setLuminance` and
/// `HctSolver.solveToInt` both return the neutral grey of that lightness
/// rather than asking for a hue that black and white do not have.
const SOLVED_RANGE: std::ops::RangeInclusive<f64> = 1.0e-4..=99.9999;

/// `3.1415927f`, the literal `Cam.kt` divides degrees by. Named rather than
/// spelled so it cannot drift from `std::f32::consts::PI`, which is the same
/// value.
const PI_F32: f32 = std::f32::consts::PI;

// -- CAM16 forward helpers ----------------------------------------------------

/// `(float) Math.pow((double) x, (double) y)`: the exponent is evaluated at
/// double width and the result narrowed, which is not the same as `f32::powf`.
fn powf32(base: f32, exponent: f32) -> f32 {
    (base as f64).powf(exponent as f64) as f32
}

/// `Math.signum(float)`, which is `0.0` at zero rather than `1.0`.
fn signum_f32(value: f32) -> f32 {
    if value > 0.0 {
        1.0
    } else if value < 0.0 {
        -1.0
    } else {
        value
    }
}

/// `CamUtils.linearized`: an 8-bit channel as linear light on a 0..=100 scale.
fn linearized(channel: u32) -> f32 {
    let normalized = channel as f32 / 255.0;
    if normalized <= 0.04045 {
        normalized / 12.92 * 100.0
    } else {
        (((normalized + 0.055) / 1.055) as f64).powf(2.4_f32 as f64) as f32 * 100.0
    }
}

/// `CamUtils.xyzFromInt`.
fn xyz_from_argb(argb: u32) -> [f32; 3] {
    let r = linearized((argb >> 16) & 0xFF);
    let g = linearized((argb >> 8) & 0xFF);
    let b = linearized(argb & 0xFF);
    let mut xyz = [0.0f32; 3];
    for (row, out) in SRGB_TO_XYZ.iter().zip(xyz.iter_mut()) {
        *out = (r as f64 * row[0] + g as f64 * row[1] + b as f64 * row[2]) as f32;
    }
    xyz
}

// -- HCT solver ---------------------------------------------------------------

fn sanitize_degrees(degrees: f64) -> f64 {
    let wrapped = degrees % 360.0;
    if wrapped < 0.0 {
        wrapped + 360.0
    } else {
        wrapped
    }
}

fn sanitize_radians(angle: f64) -> f64 {
    (angle + std::f64::consts::PI * 8.0) % (std::f64::consts::PI * 2.0)
}

/// `CamUtils.yFromLstar`: relative luminance on a 0..=100 scale.
///
/// The cube is `Math.pow(x, 3.0)` in the original and stays a `powf` here.
/// Writing it as `x * x * x` is the same number to within an ulp and the
/// bisection that reads it is sensitive enough for that ulp to move a channel.
fn y_from_lstar(l_star: f64) -> f64 {
    if l_star > 8.0 {
        ((l_star + 16.0) / 116.0).powf(3.0) * 100.0
    } else {
        l_star / KAPPA * 100.0
    }
}

/// The CIE standard's `kappa`, spelled as the decimal the platform holds.
const KAPPA: f64 = 903.296_296_296_296_3;
const EPSILON: f64 = 0.008_856_451_679_035_631;
/// `1 / 2.4` as the platform spells it. The quotient is a *different* double
/// from this literal, and both `delinearized` and `trueDelinearized` use the
/// literal.
const INVERSE_GAMMA: f64 = 0.416_666_666_666_666_7;

/// `CamUtils.argbFromLstar`: the neutral grey of a given lightness.
fn argb_from_lstar(l_star: f64) -> u32 {
    let component = (l_star + 16.0) / 116.0;
    // Here the cube really is written out, unlike in `y_from_lstar`.
    let cube = component * component * component;
    let y = if l_star > 8.0 { cube } else { l_star / KAPPA };
    let xz = if cube > EPSILON { cube } else { l_star / KAPPA };
    argb_from_xyz(
        xz * WHITE_POINT_D65[0] as f64,
        y * WHITE_POINT_D65[1] as f64,
        xz * WHITE_POINT_D65[2] as f64,
    )
}

fn argb_from_xyz(x: f64, y: f64, z: f64) -> u32 {
    let linear = |row: &[f64; 3]| row[0] * x + row[1] * y + row[2] * z;
    argb_from_rgb(
        delinearized(linear(&XYZ_TO_SRGB[0])),
        delinearized(linear(&XYZ_TO_SRGB[1])),
        delinearized(linear(&XYZ_TO_SRGB[2])),
    )
}

/// `CamUtils.delinearized`: linear light on 0..=100 to an 8-bit channel.
fn delinearized(component: f64) -> u32 {
    let normalized = component / 100.0;
    let delinearized = if normalized <= 0.0031308 {
        normalized * 12.92
    } else {
        1.055 * normalized.powf(INVERSE_GAMMA) - 0.055
    };
    // `Math.round(double)` is `floor(x + 0.5)`, which differs from Rust's
    // half-away-from-zero on the negatives the clamp then discards anyway.
    (((delinearized * 255.0) + 0.5).floor() as i64).clamp(0, 255) as u32
}

/// `HctSolver.trueDelinearized`: the same map without the rounding, which is
/// what a plane index is read off.
fn true_delinearized(rgb_component: f64) -> f64 {
    let normalized = rgb_component / 100.0;
    let delinearized = if normalized <= 0.0031308 {
        normalized * 12.92
    } else {
        1.055 * normalized.powf(INVERSE_GAMMA) - 0.055
    };
    delinearized * 255.0
}

fn argb_from_rgb(r: u32, g: u32, b: u32) -> u32 {
    0xFF00_0000 | ((r & 0xFF) << 16) | ((g & 0xFF) << 8) | (b & 0xFF)
}

fn argb_from_linrgb(linrgb: [f64; 3]) -> u32 {
    argb_from_rgb(
        delinearized(linrgb[0]),
        delinearized(linrgb[1]),
        delinearized(linrgb[2]),
    )
}

/// `CamUtils.signum`, which is **zero at zero** — `f64::signum` is not, and the
/// difference reaches a channel through the adaptation below.
fn signum(value: f64) -> f64 {
    if value < 0.0 {
        -1.0
    } else if value == 0.0 {
        0.0
    } else {
        1.0
    }
}

fn chromatic_adaptation(component: f64) -> f64 {
    let af = component.abs().powf(0.42);
    signum(component) * 400.0 * af / (af + 27.13)
}

/// `1 / 0.42`, again as the literal the platform holds rather than the
/// quotient.
const INVERSE_ADAPTATION_EXPONENT: f64 = 2.380_952_380_952_381;

fn inverse_chromatic_adaptation(adapted: f64) -> f64 {
    let adapted_abs = adapted.abs();
    let base = (0.0f64).max(27.13 * adapted_abs / (400.0 - adapted_abs));
    signum(adapted) * base.powf(INVERSE_ADAPTATION_EXPONENT)
}

fn matrix_multiply(row: [f64; 3], matrix: &[[f64; 3]; 3]) -> [f64; 3] {
    [
        row[0] * matrix[0][0] + row[1] * matrix[0][1] + row[2] * matrix[0][2],
        row[0] * matrix[1][0] + row[1] * matrix[1][1] + row[2] * matrix[1][2],
        row[0] * matrix[2][0] + row[1] * matrix[2][1] + row[2] * matrix[2][2],
    ]
}

/// The CAM16 hue of a linear-RGB point, in radians — `HctSolver.hueOf`.
fn hue_of(linrgb: [f64; 3]) -> f64 {
    let scaled = matrix_multiply(linrgb, &SCALED_DISCOUNT_FROM_LINRGB);
    let r_a = chromatic_adaptation(scaled[0]);
    let g_a = chromatic_adaptation(scaled[1]);
    let b_a = chromatic_adaptation(scaled[2]);
    let a = (11.0 * r_a + -12.0 * g_a + b_a) / 11.0;
    let b = (r_a + g_a - 2.0 * b_a) / 9.0;
    b.atan2(a)
}

fn are_in_cyclic_order(a: f64, b: f64, c: f64) -> bool {
    sanitize_radians(b - a) < sanitize_radians(c - a)
}

fn intercept(source: f64, mid: f64, target: f64) -> f64 {
    (mid - source) / (target - source)
}

fn lerp_point(source: [f64; 3], t: f64, target: [f64; 3]) -> [f64; 3] {
    [
        source[0] + (target[0] - source[0]) * t,
        source[1] + (target[1] - source[1]) * t,
        source[2] + (target[2] - source[2]) * t,
    ]
}

fn set_coordinate(source: [f64; 3], coordinate: f64, target: [f64; 3], axis: usize) -> [f64; 3] {
    lerp_point(
        source,
        intercept(source[axis], coordinate, target[axis]),
        target,
    )
}

fn is_bounded(value: f64) -> bool {
    (0.0..=100.0).contains(&value)
}

/// The `n`th place where a plane of constant Y meets an edge of the RGB cube —
/// `HctSolver.nthVertex`. A miss is reported as a negative first coordinate,
/// which is how the caller skips it.
fn nth_vertex(y: f64, n: i32) -> [f64; 3] {
    let [k_r, k_g, k_b] = Y_FROM_LINRGB;
    let coord_a = if n % 4 <= 1 { 0.0 } else { 100.0 };
    let coord_b = if n % 2 == 0 { 0.0 } else { 100.0 };
    const MISS: [f64; 3] = [-1.0, -1.0, -1.0];
    if n < 4 {
        let (g, b) = (coord_a, coord_b);
        let r = (y - g * k_g - b * k_b) / k_r;
        if is_bounded(r) {
            [r, g, b]
        } else {
            MISS
        }
    } else if n < 8 {
        let (b, r) = (coord_a, coord_b);
        let g = (y - r * k_r - b * k_b) / k_g;
        if is_bounded(g) {
            [r, g, b]
        } else {
            MISS
        }
    } else {
        let (r, g) = (coord_a, coord_b);
        let b = (y - r * k_r - g * k_g) / k_b;
        if is_bounded(b) {
            [r, g, b]
        } else {
            MISS
        }
    }
}

/// The two cube vertices the target hue lies between, on the plane of constant
/// Y — `HctSolver.bisectToSegment`.
fn bisect_to_segment(y: f64, target_hue: f64) -> ([f64; 3], [f64; 3]) {
    let mut left = [-1.0f64; 3];
    let mut right = left;
    let mut left_hue = 0.0;
    let mut right_hue = 0.0;
    let mut initialized = false;
    let mut uncut = true;
    for n in 0..12 {
        let mid = nth_vertex(y, n);
        if mid[0] < 0.0 {
            continue;
        }
        let mid_hue = hue_of(mid);
        if !initialized {
            left = mid;
            right = mid;
            left_hue = mid_hue;
            right_hue = mid_hue;
            initialized = true;
            continue;
        }
        if uncut || are_in_cyclic_order(left_hue, mid_hue, right_hue) {
            uncut = false;
            if are_in_cyclic_order(left_hue, target_hue, mid_hue) {
                right = mid;
                right_hue = mid_hue;
            } else {
                left = mid;
                left_hue = mid_hue;
            }
        }
    }
    (left, right)
}

fn critical_plane_below(x: f64) -> i32 {
    (x - 0.5).floor() as i32
}

fn critical_plane_above(x: f64) -> i32 {
    (x - 0.5).ceil() as i32
}

/// The surface point of the requested hue at the requested lightness —
/// `HctSolver.bisectToLimit`.
fn bisect_to_limit(y: f64, target_hue: f64) -> [f64; 3] {
    let (mut left, mut right) = bisect_to_segment(y, target_hue);
    let mut left_hue = hue_of(left);
    for axis in 0..3 {
        if left[axis] == right[axis] {
            continue;
        }
        let (mut l_plane, mut r_plane) = if left[axis] < right[axis] {
            (
                critical_plane_below(true_delinearized(left[axis])),
                critical_plane_above(true_delinearized(right[axis])),
            )
        } else {
            (
                critical_plane_above(true_delinearized(left[axis])),
                critical_plane_below(true_delinearized(right[axis])),
            )
        };
        for _ in 0..8 {
            if (r_plane - l_plane).abs() <= 1 {
                break;
            }
            let m_plane = ((l_plane + r_plane) as f64 / 2.0).floor() as i32;
            let Some(&coordinate) = CRITICAL_PLANES.get(m_plane as usize) else {
                break;
            };
            let mid = set_coordinate(left, coordinate, right, axis);
            let mid_hue = hue_of(mid);
            if are_in_cyclic_order(left_hue, target_hue, mid_hue) {
                right = mid;
                r_plane = m_plane;
            } else {
                left = mid;
                left_hue = mid_hue;
                l_plane = m_plane;
            }
        }
    }
    [
        (left[0] + right[0]) / 2.0,
        (left[1] + right[1]) / 2.0,
        (left[2] + right[2]) / 2.0,
    ]
}

/// Newton's method on lightness — `HctSolver.findResultByJ`. `None` where the
/// requested chroma is outside the gamut at this lightness, which is what sends
/// the solver to the cube's surface.
fn find_result_by_j(hue_radians: f64, chroma: f64, y: f64) -> Option<u32> {
    let mut j = y.sqrt() * 11.0;
    let n = FRAME.n as f64;
    let t_inner_coeff = 1.0 / (1.64 - 0.29f64.powf(n)).powf(0.73);
    let e_hue = 0.25 * ((hue_radians + 2.0).cos() + 3.8);
    let p1 = e_hue * 3846.153846153846 * FRAME.nc as f64 * FRAME.ncb as f64;
    let h_sin = hue_radians.sin();
    let h_cos = hue_radians.cos();
    for round in 0..5 {
        let j_normalized = j / 100.0;
        let alpha = if chroma == 0.0 || j == 0.0 {
            0.0
        } else {
            chroma / j_normalized.sqrt()
        };
        let t = (alpha * t_inner_coeff).powf(1.111_111_111_111_111_2);
        let ac = FRAME.aw as f64 * j_normalized.powf(1.0 / FRAME.c as f64 / FRAME.z as f64);
        let p2 = ac / FRAME.nbb as f64;
        let gamma = 23.0 * (p2 + 0.305) * t / (23.0 * p1 + 11.0 * t * h_cos + 108.0 * t * h_sin);
        let a = gamma * h_cos;
        let b = gamma * h_sin;
        let r_a = (460.0 * p2 + 451.0 * a + 288.0 * b) / 1403.0;
        let g_a = (460.0 * p2 - 891.0 * a - 261.0 * b) / 1403.0;
        let b_a = (460.0 * p2 - 220.0 * a - 6300.0 * b) / 1403.0;
        let linrgb = matrix_multiply(
            [
                inverse_chromatic_adaptation(r_a),
                inverse_chromatic_adaptation(g_a),
                inverse_chromatic_adaptation(b_a),
            ],
            &LINRGB_FROM_SCALED_DISCOUNT,
        );
        if linrgb[0] < 0.0 || linrgb[1] < 0.0 || linrgb[2] < 0.0 {
            return None;
        }
        let fnj = Y_FROM_LINRGB[0] * linrgb[0]
            + Y_FROM_LINRGB[1] * linrgb[1]
            + Y_FROM_LINRGB[2] * linrgb[2];
        if fnj <= 0.0 {
            return None;
        }
        if round == 4 || (fnj - y).abs() < 0.002 {
            if linrgb[0] > 100.01 || linrgb[1] > 100.01 || linrgb[2] > 100.01 {
                return None;
            }
            return Some(argb_from_linrgb(linrgb));
        }
        j -= (fnj - y) * j / (2.0 * fnj);
    }
    None
}

// -- colour plumbing ----------------------------------------------------------

/// `Color.toArgb()`, which is a read of the packed value the colour already
/// holds — see [`Color::srgb_8bit`] for why a Compose colour is eight bits
/// before anything paints with it.
fn argb_from_color(color: Color) -> u32 {
    let channel = |value: f32| (value.clamp(0.0, 1.0) * 255.0).round() as u32;
    argb_from_rgb(channel(color.0), channel(color.1), channel(color.2))
}

fn color_from_argb(argb: u32) -> Color {
    Color::from_rgb_u8(
        ((argb >> 16) & 0xFF) as u8,
        ((argb >> 8) & 0xFF) as u8,
        (argb & 0xFF) as u8,
    )
}

#[cfg(test)]
mod tests {
    use super::*;

    fn rgb(color: Color) -> (u8, u8, u8) {
        (
            (color.0 * 255.0).round() as u8,
            (color.1 * 255.0).round() as u8,
            (color.2 * 255.0).round() as u8,
        )
    }

    #[test]
    fn the_wear_scroll_indicator_colours_are_the_ones_compose_draws() {
        // Measured on the shipping Compose build: an `onBackground` of
        // #DFF6FF gives a thumb of (180, 202, 211) and a track of (30, 51, 58).
        let on_background = Color::from_rgb_u8(0xDF, 0xF6, 0xFF);
        assert_eq!(rgb(set_luminance(on_background, 80.0)), (180, 202, 211));
        assert_eq!(rgb(set_luminance(on_background, 20.0)), (30, 51, 58));
    }

    #[test]
    fn a_lab_round_trip_would_have_given_a_different_track() {
        // The pair the two models disagree on. Substituting L* in Lab keeps
        // Lab's a*/b* and lands on 31 of red; CAM16 keeps its own hue and
        // chroma and lands on 30. Both agree at L* 80, which is why checking
        // only the thumb clears an implementation that is wrong.
        let track = set_luminance(Color::from_rgb_u8(0xDF, 0xF6, 0xFF), 20.0);
        assert_eq!(rgb(track).0, 30, "not the Lab answer of 31");
    }

    #[test]
    fn set_luminance_matches_the_shipped_aar() {
        // (source ARGB, L*, expected ARGB), generated by calling
        // `androidx.wear.compose.material3`'s own `Cam` and `HctSolver` out of
        // compose-material3-1.6.2.aar on a JVM. Both solver paths and both
        // short-circuits are represented.
        const GOLDEN: &[(u32, f32, u32)] = &[
            (0xFFDF_F6FF, 0.0, 0xFF00_0000),
            (0xFFDF_F6FF, 1.0, 0xFF00_0407),
            (0xFFDF_F6FF, 5.0, 0xFF00_131A),
            (0xFFDF_F6FF, 20.0, 0xFF1E_333A),
            (0xFFDF_F6FF, 33.3, 0xFF3D_5259),
            (0xFFDF_F6FF, 50.0, 0xFF65_7A82),
            (0xFFDF_F6FF, 66.7, 0xFF90_A6AE),
            (0xFFDF_F6FF, 80.0, 0xFFB4_CAD3),
            (0xFFDF_F6FF, 95.0, 0xFFDE_F5FE),
            (0xFFDF_F6FF, 100.0, 0xFFFF_FFFF),
            (0xFFE3_E3E3, 20.0, 0xFF2F_3131),
            (0xFFE3_E3E3, 80.0, 0xFFC6_C6C6),
            (0xFFB9_F2FF, 20.0, 0xFF00_363E),
            (0xFFB9_F2FF, 80.0, 0xFF97_D0DC),
            (0xFF5E_7E93, 20.0, 0xFF10_3446),
            (0xFF5E_7E93, 80.0, 0xFFAA_CBE2),
            (0xFFFF_FFFF, 20.0, 0xFF2F_3131),
            (0xFF00_0000, 50.0, 0xFF77_7777),
            (0xFFFF_4B32, 20.0, 0xFF67_0500),
            (0xFFFF_4B32, 80.0, 0xFFFF_B4A7),
            (0xFF2F_A8F5, 20.0, 0xFF00_3351),
            (0xFF66_D9FF, 80.0, 0xFF61_D4FA),
        ];
        for &(source, l_star, expected) in GOLDEN {
            let got = argb_from_color(set_luminance(color_from_argb(source), l_star));
            assert_eq!(
                got, expected,
                "setLuminance(#{source:08X}, {l_star}) = #{got:08X}, want #{expected:08X}"
            );
        }
    }

    #[test]
    fn the_cam16_forward_reports_the_hue_and_chroma_the_platform_reports() {
        // `Cam.fromInt(0xFFDFF6FF)` on the AAR: hue 220.671203613,
        // chroma 15.368903160. The two are floats there and floats here.
        let (hue, chroma) = cam16_hue_chroma(0xFFDF_F6FF);
        assert!((hue - 220.671_2).abs() < 1e-3, "hue {hue}");
        assert!((chroma - 15.368_903).abs() < 1e-4, "chroma {chroma}");
    }

    #[test]
    fn the_ends_of_the_range_are_neutral_rather_than_solved() {
        // Wear short-circuits outside 0.0001..99.9999 rather than asking the
        // solver for a hue that black and white do not have.
        assert_eq!(
            rgb(set_luminance(Color::from_rgb_u8(255, 0, 0), 0.0)),
            (0, 0, 0)
        );
        assert_eq!(
            rgb(set_luminance(Color::from_rgb_u8(255, 0, 0), 100.0)),
            (255, 255, 255)
        );
    }

    #[test]
    fn a_grey_source_stays_grey_at_any_lightness() {
        for l_star in [10.0, 20.0, 50.0, 80.0, 90.0] {
            let (r, g, b) = rgb(set_luminance(Color::from_rgb_u8(0x80, 0x80, 0x80), l_star));
            assert!(
                r.abs_diff(g) <= 1 && g.abs_diff(b) <= 1,
                "L* {l_star} gave ({r}, {g}, {b})"
            );
        }
    }

    #[test]
    fn the_critical_planes_are_the_eight_bit_boundaries() {
        // Every entry is the linear-RGB coordinate that delinearizes to exactly
        // half a level, which is what makes bisecting by plane index land on
        // values the framebuffer can hold.
        for (index, &plane) in CRITICAL_PLANES.iter().enumerate() {
            let boundary = index as f64 + 0.5;
            assert!(
                (true_delinearized(plane) - boundary).abs() < 1e-9,
                "plane {index} is {plane}, which delinearizes to {}",
                true_delinearized(plane)
            );
        }
    }
}