Port osu!taiko updates since f08134f

This commit is contained in:
MaxOhn
2024-10-11 00:35:31 +02:00
committed by tsunyoku
parent 0c79c4073e
commit 43c5788ffb
11 changed files with 483 additions and 92 deletions
+4 -2
View File
@@ -65,8 +65,10 @@ impl ManiaGradualDifficulty {
let clock_rate = difficulty.get_clock_rate();
let mut params = ObjectParams::new(converted);
let HitWindows { od: hit_window, .. } =
converted.attributes().difficulty(&difficulty).hit_windows();
let HitWindows {
od_great: hit_window,
..
} = converted.attributes().difficulty(&difficulty).hit_windows();
let mania_objects = converted
.hit_objects
+1 -1
View File
@@ -28,7 +28,7 @@ pub fn difficulty(
.attributes()
.difficulty(difficulty)
.hit_windows()
.od;
.od_great;
ManiaDifficultyAttributes {
stars: values.strain.difficulty_value() * DIFFICULTY_MULTIPLIER,
+67 -30
View File
@@ -27,7 +27,11 @@ pub struct HitWindows {
/// Hit window for approach rate i.e. `TimePreempt` in milliseconds.
pub ar: f64,
/// Hit window for overall difficulty i.e. time to hit a 300 ("Great") in milliseconds.
pub od: f64,
pub od_great: f64,
/// Hit window for overall difficulty i.e. time to hit a 100 ("Ok") in milliseconds.
///
/// `None` for osu!mania.
pub od_ok: Option<f64>,
}
/// A builder for [`BeatmapAttributes`] and [`HitWindows`].
@@ -44,15 +48,43 @@ pub struct BeatmapAttributesBuilder {
clock_rate: Option<f64>,
}
struct GameModeHitWindows {
min: f64,
avg: f64,
max: f64,
}
const OSU_GREAT: GameModeHitWindows = GameModeHitWindows {
min: 80.0,
avg: 50.0,
max: 20.0,
};
const OSU_OK: GameModeHitWindows = GameModeHitWindows {
min: 140.0,
avg: 100.0,
max: 60.0,
};
const TAIKO_GREAT: GameModeHitWindows = GameModeHitWindows {
min: 50.0,
avg: 35.0,
max: 20.0,
};
const TAIKO_OK: GameModeHitWindows = GameModeHitWindows {
min: 120.0,
avg: 80.0,
max: 50.0,
};
const AR_WINDOWS: GameModeHitWindows = GameModeHitWindows {
min: 1800.0,
avg: 1200.0,
max: 450.0,
};
impl BeatmapAttributesBuilder {
const OSU_MIN: f64 = 80.0;
const OSU_AVG: f64 = 50.0;
const OSU_MAX: f64 = 20.0;
const TAIKO_MIN: f64 = 50.0;
const TAIKO_AVG: f64 = 35.0;
const TAIKO_MAX: f64 = 20.0;
/// Create a new [`BeatmapAttributesBuilder`].
///
/// The mode will be `GameMode::Osu` and attributes are set to `5.0`.
@@ -219,10 +251,10 @@ impl BeatmapAttributesBuilder {
mod_mult(self.ar.value(mods, GameMods::ar))
};
let preempt = difficulty_range(f64::from(raw_ar), 1800.0, 1200.0, 450.0) / ar_clock_rate;
let preempt = difficulty_range(f64::from(raw_ar), AR_WINDOWS) / ar_clock_rate;
// OD
let hit_window = match self.mode {
let (great, ok) = match self.mode {
GameMode::Osu | GameMode::Catch => {
let raw_od = if self.od.with_mods() {
self.od.value(mods, GameMods::od)
@@ -230,12 +262,10 @@ impl BeatmapAttributesBuilder {
mod_mult(self.od.value(mods, GameMods::od))
};
difficulty_range(
f64::from(raw_od),
Self::OSU_MIN,
Self::OSU_AVG,
Self::OSU_MAX,
) / od_clock_rate
let great = difficulty_range(f64::from(raw_od), OSU_GREAT) / od_clock_rate;
let ok = difficulty_range(f64::from(raw_od), OSU_OK) / od_clock_rate;
(great, Some(ok))
}
GameMode::Taiko => {
let raw_od = if self.od.with_mods() {
@@ -244,14 +274,10 @@ impl BeatmapAttributesBuilder {
mod_mult(self.od.value(mods, GameMods::od))
};
let diff_range = difficulty_range(
f64::from(raw_od),
Self::TAIKO_MIN,
Self::TAIKO_AVG,
Self::TAIKO_MAX,
);
let great = difficulty_range(f64::from(raw_od), TAIKO_GREAT) / od_clock_rate;
let ok = difficulty_range(f64::from(raw_od), TAIKO_OK) / od_clock_rate;
diff_range / od_clock_rate
(great, Some(ok))
}
GameMode::Mania => {
let mut value = if !self.is_convert {
@@ -270,13 +296,16 @@ impl BeatmapAttributesBuilder {
}
}
((f64::from(value) * od_clock_rate).floor() / od_clock_rate).ceil()
let great = ((f64::from(value) * od_clock_rate).floor() / od_clock_rate).ceil();
(great, None)
}
};
HitWindows {
ar: preempt,
od: hit_window,
od_great: great,
od_ok: ok,
}
}
@@ -308,7 +337,11 @@ impl BeatmapAttributesBuilder {
}
let hit_windows = self.hit_windows();
let HitWindows { ar, od } = hit_windows;
let HitWindows {
ar,
od_great,
od_ok: _,
} = hit_windows;
// AR
let ar = if ar > 1200.0 {
@@ -319,8 +352,10 @@ impl BeatmapAttributesBuilder {
// OD
let od = match self.mode {
GameMode::Osu => (Self::OSU_MIN - od) / 6.0,
GameMode::Taiko => (Self::TAIKO_MIN - od) / (Self::TAIKO_MIN - Self::TAIKO_AVG) * 5.0,
GameMode::Osu => (OSU_GREAT.min - od_great) / 6.0,
GameMode::Taiko => {
(TAIKO_GREAT.min - od_great) / (TAIKO_GREAT.min - TAIKO_GREAT.avg) * 5.0
}
GameMode::Catch | GameMode::Mania => f64::from(self.od.value(mods, GameMods::od)),
};
@@ -347,7 +382,9 @@ impl<M> From<&Converted<'_, M>> for BeatmapAttributesBuilder {
}
}
fn difficulty_range(difficulty: f64, min: f64, mid: f64, max: f64) -> f64 {
fn difficulty_range(difficulty: f64, windows: GameModeHitWindows) -> f64 {
let GameModeHitWindows { min, avg: mid, max } = windows;
if difficulty > 5.0 {
mid + (max - mid) * (difficulty - 5.0) / 5.0
} else if difficulty < 5.0 {
+1 -1
View File
@@ -26,7 +26,7 @@ impl OsuSkills {
map_attrs: &BeatmapAttributes,
time_preempt: f64,
) -> Self {
let hit_window = 2.0 * map_attrs.hit_windows.od;
let hit_window = 2.0 * map_attrs.hit_windows.od_great;
// * Preempt time can go below 450ms. Normally, this is achieved via the DT mod
// * which uniformly speeds up all animations game wide regardless of AR.
+5 -1
View File
@@ -12,7 +12,9 @@ pub struct TaikoDifficultyAttributes {
/// The difficulty of the hardest parts of the map.
pub peak: f64,
/// The perceived hit window for an n300 inclusive of rate-adjusting mods (DT/HT/etc)
pub hit_window: f64,
pub great_hit_window: f64,
/// The perceived hit window for an n100 inclusive of rate-adjusting mods (DT/HT/etc)
pub ok_hit_window: f64,
/// The final star rating.
pub stars: f64,
/// The maximum combo.
@@ -55,6 +57,8 @@ pub struct TaikoPerformanceAttributes {
pub pp_difficulty: f64,
/// Scaled miss count based on total hits.
pub effective_miss_count: f64,
/// Upper bound on the player's tap deviation.
pub estimated_unstable_rate: Option<f64>,
}
impl TaikoPerformanceAttributes {
+26 -19
View File
@@ -1,9 +1,6 @@
use std::cmp;
use rosu_map::{
section::{general::GameMode, hit_objects::CurveBuffers},
util::Pos,
};
use rosu_map::{section::general::GameMode, util::Pos};
use crate::{
model::{
@@ -20,7 +17,7 @@ use super::Taiko;
/// A [`Beatmap`] for [`Taiko`] calculations.
pub type TaikoBeatmap<'a> = Converted<'a, Taiko>;
const LEGACY_TAIKO_VELOCITY_MULTIPLIER: f32 = 1.4;
const VELOCITY_MULTIPLIER: f32 = 1.4;
const OSU_BASE_SCORING_DIST: f32 = 100.0;
pub const fn check_convert(map: &Beatmap) -> ConvertStatus {
@@ -44,8 +41,6 @@ pub fn try_convert(map: &mut Beatmap) -> ConvertStatus {
}
fn convert(map: &mut Beatmap) {
map.slider_multiplier *= f64::from(LEGACY_TAIKO_VELOCITY_MULTIPLIER);
let mut new_objects = Vec::new();
let mut new_sounds = Vec::new();
@@ -125,32 +120,46 @@ fn convert(map: &mut Beatmap) {
fn should_convert_slider_to_taiko_hits(map: &Beatmap, params: &mut SliderParams<'_>) -> bool {
let SliderParams {
slider,
bufs,
duration,
start_time,
tick_spacing,
} = params;
let curve = slider.curve(GameMode::Taiko, bufs);
// * The true distance, accounting for any repeats. This ends up being the drum roll distance later
let spans = slider.span_count() as f64;
let dist = curve.dist() * spans * f64::from(LEGACY_TAIKO_VELOCITY_MULTIPLIER);
let mut dist = slider.expected_dist.unwrap_or(0.0);
// * Do not combine the following two lines!
dist *= f64::from(VELOCITY_MULTIPLIER);
dist *= spans;
let timing_beat_len = map
.timing_point_at(*start_time)
.map_or(TimingPoint::DEFAULT_BEAT_LEN, |point| point.beat_len);
let bpm_multiplier = map
let slider_velocity = map
.difficulty_point_at(*start_time)
.map_or(DifficultyPoint::DEFAULT_BPM_MULTIPLIER, |point| {
point.bpm_multiplier
.map_or(DifficultyPoint::DEFAULT_SLIDER_VELOCITY, |point| {
point.slider_velocity
});
let mut beat_len = timing_beat_len * bpm_multiplier;
fn get_precision_adjusted_beat_len(slider_velocity_multiplier: f64, beat_len: f64) -> f64 {
let slider_velocity_as_beat_len = -100.0 / slider_velocity_multiplier;
let slider_scoring_point_dist =
f64::from(OSU_BASE_SCORING_DIST) * map.slider_multiplier / map.slider_tick_rate;
let bpm_multiplier = if slider_velocity_as_beat_len < 0.0 {
f64::from(((-slider_velocity_as_beat_len) as f32).clamp(10.0, 10_000.0)) / 100.0
} else {
1.0
};
beat_len * bpm_multiplier
}
let mut beat_len = get_precision_adjusted_beat_len(slider_velocity, timing_beat_len);
let slider_scoring_point_dist = f64::from(OSU_BASE_SCORING_DIST)
* (map.slider_multiplier * f64::from(VELOCITY_MULTIPLIER))
/ map.slider_tick_rate;
// * The velocity and duration of the taiko hit object - calculated as the velocity of a drum roll.
let taiko_vel = slider_scoring_point_dist * map.slider_tick_rate;
@@ -171,7 +180,6 @@ fn should_convert_slider_to_taiko_hits(map: &Beatmap, params: &mut SliderParams<
struct SliderParams<'c> {
slider: &'c Slider,
bufs: CurveBuffers,
duration: u32,
start_time: f64,
tick_spacing: f64,
@@ -181,7 +189,6 @@ impl<'c> SliderParams<'c> {
fn new(start_time: f64, slider: &'c Slider) -> Self {
Self {
slider,
bufs: CurveBuffers::default(),
start_time,
duration: 0,
tick_spacing: 0.0,
+7 -3
View File
@@ -83,8 +83,11 @@ impl TaikoGradualDifficulty {
(Some(true), Some(true)) => FirstTwoCombos::Both,
};
let HitWindows { od: hit_window, .. } =
converted.attributes().difficulty(&difficulty).hit_windows();
let HitWindows {
od_great,
od_ok,
ar: _,
} = converted.attributes().difficulty(&difficulty).hit_windows();
let mut n_diff_objects = 0;
let mut max_combo = 0;
@@ -100,7 +103,8 @@ impl TaikoGradualDifficulty {
let skills = TaikoSkills::new();
let attrs = TaikoDifficultyAttributes {
hit_window,
great_hit_window: od_great,
ok_hit_window: od_ok.unwrap_or(0.0),
is_convert: converted.is_convert,
..Default::default()
};
+8 -6
View File
@@ -2,6 +2,7 @@ use std::cmp;
use crate::{
any::difficulty::skills::Skill,
model::beatmap::HitWindows,
taiko::{
difficulty::{
color::preprocessor::ColorDifficultyPreprocessor,
@@ -31,16 +32,17 @@ pub fn difficulty(
difficulty: &Difficulty,
converted: &TaikoBeatmap<'_>,
) -> TaikoDifficultyAttributes {
let hit_window = converted
.attributes()
.difficulty(difficulty)
.hit_windows()
.od;
let HitWindows {
od_great,
od_ok,
ar: _,
} = converted.attributes().difficulty(difficulty).hit_windows();
let DifficultyValues { skills, max_combo } = DifficultyValues::calculate(difficulty, converted);
let mut attrs = TaikoDifficultyAttributes {
hit_window,
great_hit_window: od_great,
ok_hit_window: od_ok.unwrap_or(0.0),
max_combo,
is_convert: converted.is_convert,
..Default::default()
+86 -29
View File
@@ -4,7 +4,7 @@ use crate::{
any::{Difficulty, HitResultPriority, IntoModePerformance, IntoPerformance},
model::mods::GameMods,
osu::OsuPerformance,
util::map_or_attrs::MapOrAttrs,
util::{map_or_attrs::MapOrAttrs, special_functions},
Performance,
};
@@ -403,6 +403,10 @@ impl TaikoPerformanceInner<'_> {
// * and increasing the miss penalty for shorter object counts lower than 1000.
let total_successful_hits = self.total_successful_hits();
let estimated_unstable_rate = self
.compute_deviation_upper_bound(total_successful_hits)
.map(|v| v * 10.0);
let effective_miss_count = if total_successful_hits > 0 {
(1000.0 / f64::from(total_successful_hits)).max(1.0) * f64::from(self.state.misses)
} else {
@@ -419,8 +423,9 @@ impl TaikoPerformanceInner<'_> {
multiplier *= 0.975;
}
let diff_value = self.compute_difficulty_value(effective_miss_count);
let acc_value = self.compute_accuracy_value();
let diff_value =
self.compute_difficulty_value(effective_miss_count, estimated_unstable_rate);
let acc_value = self.compute_accuracy_value(estimated_unstable_rate);
let pp = (diff_value.powf(1.1) + acc_value.powf(1.1)).powf(1.0 / 1.1) * multiplier;
@@ -430,10 +435,19 @@ impl TaikoPerformanceInner<'_> {
pp_acc: acc_value,
pp_difficulty: diff_value,
effective_miss_count,
estimated_unstable_rate,
}
}
fn compute_difficulty_value(&self, effective_miss_count: f64) -> f64 {
fn compute_difficulty_value(
&self,
effective_miss_count: f64,
estimated_unstable_rate: Option<f64>,
) -> f64 {
let Some(estimated_unstable_rate) = estimated_unstable_rate else {
return 0.0;
};
let attrs = &self.attrs;
let exp_base = 5.0 * (attrs.stars / 0.115).max(1.0) - 4.0;
let mut diff_value = exp_base.powf(2.25) / 1150.0;
@@ -452,39 +466,95 @@ impl TaikoPerformanceInner<'_> {
}
if self.mods.hr() {
diff_value *= 1.05;
diff_value *= 1.10;
}
if self.mods.fl() {
diff_value *= 1.05 * len_bonus;
}
let acc = self.custom_accuracy();
diff_value * acc.powf(2.0)
diff_value
* (special_functions::erf(400.0 / (2.0_f64.sqrt() * estimated_unstable_rate))).powf(2.0)
}
fn compute_accuracy_value(&self) -> f64 {
if self.attrs.hit_window <= 0.0 {
fn compute_accuracy_value(&self, estimated_unstable_rate: Option<f64>) -> f64 {
if self.attrs.great_hit_window <= 0.0 {
return 0.0;
}
let mut acc_value = (60.0 / self.attrs.hit_window).powf(1.1)
* self.custom_accuracy().powf(8.0)
* self.attrs.stars.powf(0.4)
* 27.0;
let Some(estimated_unstable_rate) = estimated_unstable_rate else {
return 0.0;
};
let mut acc_value =
(70.0 / estimated_unstable_rate).powf(1.1) * self.attrs.stars.powf(0.4) * 100.0;
let len_bonus = (self.total_hits() / 1500.0).powf(0.3).min(1.15);
acc_value *= len_bonus;
// * Slight HDFL Bonus for accuracy. A clamp is used to prevent against negative values.
if self.mods.hd() && self.mods.fl() && !self.attrs.is_convert {
acc_value *= (1.075 * len_bonus).max(1.05);
acc_value *= (1.05 * len_bonus).max(1.0);
}
acc_value
}
// * Computes an upper bound on the player's tap deviation based on the OD, number of circles and sliders,
// * and the hit judgements, assuming the player's mean hit error is 0. The estimation is consistent in that
// * two SS scores on the same map with the same settings will always return the same deviation.
fn compute_deviation_upper_bound(&self, total_successful_hits: u32) -> Option<f64> {
if total_successful_hits == 0 || self.attrs.great_hit_window <= 0.0 {
return None;
}
let h300 = self.attrs.great_hit_window;
let h100 = self.attrs.ok_hit_window;
let n = self.total_hits();
// * 99% critical value for the normal distribution (one-tailed).
const Z: f64 = 2.32634787404;
// * The upper bound on deviation, calculated with the ratio of 300s to objects, and the great hit window.
let calc_deviation_great_window = || {
if self.state.n300 == 0 {
return None;
}
// * Proportion of greats hit.
let p = f64::from(self.state.n300) / n;
// * We can be 99% confident that p is at least this value.
let p_lower_bound = (n * p + Z * Z / 2.0) / (n + Z * Z)
- Z / (n + Z * Z) * (n * p * (1.0 - p) + Z * Z / 4.0).sqrt();
// * We can be 99% confident that the deviation is not higher than:
Some(h300 / (2.0_f64.sqrt() * special_functions::erf_inv(p_lower_bound)))
};
// * The upper bound on deviation, calculated with the ratio of 300s + 100s to objects, and the good hit window.
// * This will return a lower value than the first method when the number of 100s is high, but the miss count is low.
let calc_deviation_good_window = || {
// * Proportion of greats + goods hit.
let p = f64::from(total_successful_hits) / n;
// * We can be 99% confident that p is at least this value.
let p_lower_bound = (n * p + Z * Z / 2.0) / (n + Z * Z)
- Z / (n + Z * Z) * (n * p * (1.0 - p) + Z * Z / 4.0).sqrt();
// * We can be 99% confident that the deviation is not higher than:
h100 / (2.0_f64.sqrt() * special_functions::erf_inv(p_lower_bound))
};
let deviation_great_window = calc_deviation_great_window();
let deviation_good_window = calc_deviation_good_window();
let Some(deviation_great_window) = deviation_great_window else {
return Some(deviation_good_window);
};
Some(deviation_great_window.min(deviation_good_window))
}
const fn total_hits(&self) -> f64 {
self.state.total_hits() as f64
}
@@ -492,19 +562,6 @@ impl TaikoPerformanceInner<'_> {
const fn total_successful_hits(&self) -> u32 {
self.state.n300 + self.state.n100
}
fn custom_accuracy(&self) -> f64 {
let total_hits = self.state.total_hits();
if total_hits == 0 {
return 0.0;
}
let numerator = self.state.n300 * 300 + self.state.n100 * 150;
let denominator = total_hits * 300;
f64::from(numerator) / f64::from(denominator)
}
}
fn accuracy(n300: u32, n100: u32, misses: u32) -> f64 {
+1
View File
@@ -4,5 +4,6 @@ pub mod limited_queue;
pub mod map_or_attrs;
pub mod random;
pub mod sort;
pub mod special_functions;
pub mod strains_vec;
pub mod sync;
+277
View File
@@ -0,0 +1,277 @@
#[rustfmt::skip]
mod consts {
pub const ERF_IMP_AN: &[f64] = &[ 0.00337916709551257388990745, -0.00073695653048167948530905, -0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468, 0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225 ];
pub const ERF_IMP_AD: &[f64] = &[ 1.0, -0.218088218087924645390535, 0.412542972725442099083918, -0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266, 0.00408165558926174048329689, -0.000615900721557769691924509 ];
pub const ERF_IMP_BN: &[f64] = &[ -0.0361790390718262471360258, 0.292251883444882683221149, 0.281447041797604512774415, 0.125610208862766947294894, 0.0274135028268930549240776, 0.00250839672168065762786937 ];
pub const ERF_IMP_BD: &[f64] = &[ 1.0, 1.8545005897903486499845, 1.43575803037831418074962, 0.582827658753036572454135, 0.124810476932949746447682, 0.0113724176546353285778481 ];
pub const ERF_IMP_CN: &[f64] = &[ -0.0397876892611136856954425, 0.153165212467878293257683, 0.191260295600936245503129, 0.10276327061989304213645, 0.029637090615738836726027, 0.0046093486780275489468812, 0.000307607820348680180548455 ];
pub const ERF_IMP_CD: &[f64] = &[ 1.0, 1.95520072987627704987886, 1.64762317199384860109595, 0.768238607022126250082483, 0.209793185936509782784315, 0.0319569316899913392596356, 0.00213363160895785378615014 ];
pub const ERF_IMP_DN: &[f64] = &[ -0.0300838560557949717328341, 0.0538578829844454508530552, 0.0726211541651914182692959, 0.0367628469888049348429018, 0.00964629015572527529605267, 0.00133453480075291076745275, 0.778087599782504251917881e-4 ];
pub const ERF_IMP_DD: &[f64] = &[ 1.0, 1.75967098147167528287343, 1.32883571437961120556307, 0.552528596508757581287907, 0.133793056941332861912279, 0.0179509645176280768640766, 0.00104712440019937356634038, -0.106640381820357337177643e-7 ];
pub const ERF_IMP_EN: &[f64] = &[ -0.0117907570137227847827732, 0.014262132090538809896674, 0.0202234435902960820020765, 0.00930668299990432009042239, 0.00213357802422065994322516, 0.00025022987386460102395382, 0.120534912219588189822126e-4 ];
pub const ERF_IMP_ED: &[f64] = &[ 1.0, 1.50376225203620482047419, 0.965397786204462896346934, 0.339265230476796681555511, 0.0689740649541569716897427, 0.00771060262491768307365526, 0.000371421101531069302990367 ];
pub const ERF_IMP_FN: &[f64] = &[ -0.00546954795538729307482955, 0.00404190278731707110245394, 0.0054963369553161170521356, 0.00212616472603945399437862, 0.000394984014495083900689956, 0.365565477064442377259271e-4, 0.135485897109932323253786e-5 ];
pub const ERF_IMP_FD: &[f64] = &[ 1.0, 1.21019697773630784832251, 0.620914668221143886601045, 0.173038430661142762569515, 0.0276550813773432047594539, 0.00240625974424309709745382, 0.891811817251336577241006e-4, -0.465528836283382684461025e-11 ];
pub const ERF_IMP_GN: &[f64] = &[ -0.00270722535905778347999196, 0.0013187563425029400461378, 0.00119925933261002333923989, 0.00027849619811344664248235, 0.267822988218331849989363e-4, 0.923043672315028197865066e-6 ];
pub const ERF_IMP_GD: &[f64] = &[ 1.0, 0.814632808543141591118279, 0.268901665856299542168425, 0.0449877216103041118694989, 0.00381759663320248459168994, 0.000131571897888596914350697, 0.404815359675764138445257e-11 ];
pub const ERF_IMP_HN: &[f64] = &[ -0.00109946720691742196814323, 0.000406425442750422675169153, 0.000274499489416900707787024, 0.465293770646659383436343e-4, 0.320955425395767463401993e-5, 0.778286018145020892261936e-7 ];
pub const ERF_IMP_HD: &[f64] = &[ 1.0, 0.588173710611846046373373, 0.139363331289409746077541, 0.0166329340417083678763028, 0.00100023921310234908642639, 0.24254837521587225125068e-4 ];
pub const ERF_IMP_IN: &[f64] = &[ -0.00056907993601094962855594, 0.000169498540373762264416984, 0.518472354581100890120501e-4, 0.382819312231928859704678e-5, 0.824989931281894431781794e-7 ];
pub const ERF_IMP_ID: &[f64] = &[ 1.0, 0.339637250051139347430323, 0.043472647870310663055044, 0.00248549335224637114641629, 0.535633305337152900549536e-4, -0.117490944405459578783846e-12 ];
pub const ERF_IMP_JN: &[f64] = &[ -0.000241313599483991337479091, 0.574224975202501512365975e-4, 0.115998962927383778460557e-4, 0.581762134402593739370875e-6, 0.853971555085673614607418e-8 ];
pub const ERF_IMP_JD: &[f64] = &[ 1.0, 0.233044138299687841018015, 0.0204186940546440312625597, 0.000797185647564398289151125, 0.117019281670172327758019e-4 ];
pub const ERF_IMP_KN: &[f64] = &[ -0.000146674699277760365803642, 0.162666552112280519955647e-4, 0.269116248509165239294897e-5, 0.979584479468091935086972e-7, 0.101994647625723465722285e-8 ];
pub const ERF_IMP_KD: &[f64] = &[ 1.0, 0.165907812944847226546036, 0.0103361716191505884359634, 0.000286593026373868366935721, 0.298401570840900340874568e-5 ];
pub const ERF_IMP_LN: &[f64] = &[ -0.583905797629771786720406e-4, 0.412510325105496173512992e-5, 0.431790922420250949096906e-6, 0.993365155590013193345569e-8, 0.653480510020104699270084e-10 ];
pub const ERF_IMP_LD: &[f64] = &[ 1.0, 0.105077086072039915406159, 0.00414278428675475620830226, 0.726338754644523769144108e-4, 0.477818471047398785369849e-6 ];
pub const ERF_IMP_MN: &[f64] = &[ -0.196457797609229579459841e-4, 0.157243887666800692441195e-5, 0.543902511192700878690335e-7, 0.317472492369117710852685e-9 ];
pub const ERF_IMP_MD: &[f64] = &[ 1.0, 0.052803989240957632204885, 0.000926876069151753290378112, 0.541011723226630257077328e-5, 0.535093845803642394908747e-15 ];
pub const ERF_IMP_NN: &[f64] = &[ -0.789224703978722689089794e-5, 0.622088451660986955124162e-6, 0.145728445676882396797184e-7, 0.603715505542715364529243e-10 ];
pub const ERF_IMP_ND: &[f64] = &[ 1.0, 0.0375328846356293715248719, 0.000467919535974625308126054, 0.193847039275845656900547e-5 ];
pub const ERV_INV_IMP_AN: &[f64] = &[ -0.000508781949658280665617, -0.00836874819741736770379, 0.0334806625409744615033, -0.0126926147662974029034, -0.0365637971411762664006, 0.0219878681111168899165, 0.00822687874676915743155, -0.00538772965071242932965 ];
pub const ERV_INV_IMP_AD: &[f64] = &[ 1.0, -0.970005043303290640362, -1.56574558234175846809, 1.56221558398423026363, 0.662328840472002992063, -0.71228902341542847553, -0.0527396382340099713954, 0.0795283687341571680018, -0.00233393759374190016776, 0.000886216390456424707504 ];
pub const ERV_INV_IMP_BN: &[f64] = &[ -0.202433508355938759655, 0.105264680699391713268, 8.37050328343119927838, 17.6447298408374015486, -18.8510648058714251895, -44.6382324441786960818, 17.445385985570866523, 21.1294655448340526258, -3.67192254707729348546 ];
pub const ERV_INV_IMP_BD: &[f64] = &[ 1.0, 6.24264124854247537712, 3.9713437953343869095, -28.6608180499800029974, -20.1432634680485188801, 48.5609213108739935468, 10.8268667355460159008, -22.6436933413139721736, 1.72114765761200282724 ];
pub const ERV_INV_IMP_CN: &[f64] = &[ -0.131102781679951906451, -0.163794047193317060787, 0.117030156341995252019, 0.387079738972604337464, 0.337785538912035898924, 0.142869534408157156766, 0.0290157910005329060432, 0.00214558995388805277169, -0.679465575181126350155e-6, 0.285225331782217055858e-7, -0.681149956853776992068e-9 ];
pub const ERV_INV_IMP_CD: &[f64] = &[ 1.0, 3.46625407242567245975, 5.38168345707006855425, 4.77846592945843778382, 2.59301921623620271374, 0.848854343457902036425, 0.152264338295331783612, 0.01105924229346489121 ];
pub const ERV_INV_IMP_DN: &[f64] = &[ -0.0350353787183177984712, -0.00222426529213447927281, 0.0185573306514231072324, 0.00950804701325919603619, 0.00187123492819559223345, 0.000157544617424960554631, 0.460469890584317994083e-5, -0.230404776911882601748e-9, 0.266339227425782031962e-11 ];
pub const ERV_INV_IMP_DD: &[f64] = &[ 1.0, 1.3653349817554063097, 0.762059164553623404043, 0.220091105764131249824, 0.0341589143670947727934, 0.00263861676657015992959, 0.764675292302794483503e-4 ];
pub const ERV_INV_IMP_EN: &[f64] = &[ -0.0167431005076633737133, -0.00112951438745580278863, 0.00105628862152492910091, 0.000209386317487588078668, 0.149624783758342370182e-4, 0.449696789927706453732e-6, 0.462596163522878599135e-8, -0.281128735628831791805e-13, 0.99055709973310326855e-16 ];
pub const ERV_INV_IMP_ED: &[f64] = &[ 1.0, 0.591429344886417493481, 0.138151865749083321638, 0.0160746087093676504695, 0.000964011807005165528527, 0.275335474764726041141e-4, 0.282243172016108031869e-6 ];
pub const ERV_INV_IMP_FN: &[f64] = &[ -0.0024978212791898131227, -0.779190719229053954292e-5, 0.254723037413027451751e-4, 0.162397777342510920873e-5, 0.396341011304801168516e-7, 0.411632831190944208473e-9, 0.145596286718675035587e-11, -0.116765012397184275695e-17 ];
pub const ERV_INV_IMP_FD: &[f64] = &[ 1.0, 0.207123112214422517181, 0.0169410838120975906478, 0.000690538265622684595676, 0.145007359818232637924e-4, 0.144437756628144157666e-6, 0.509761276599778486139e-9 ];
pub const ERV_INV_IMP_GN: &[f64] = &[ -0.000539042911019078575891, -0.28398759004727721098e-6, 0.899465114892291446442e-6, 0.229345859265920864296e-7, 0.225561444863500149219e-9, 0.947846627503022684216e-12, 0.135880130108924861008e-14, -0.348890393399948882918e-21 ];
pub const ERV_INV_IMP_GD: &[f64] = &[ 1.0, 0.0845746234001899436914, 0.00282092984726264681981, 0.468292921940894236786e-4, 0.399968812193862100054e-6, 0.161809290887904476097e-8, 0.231558608310259605225e-11 ];
}
use consts::*;
pub fn erf(x: f64) -> f64 {
if x.abs() < f64::EPSILON {
return 0.0;
}
if x == f64::INFINITY {
return 1.0;
}
if x == f64::NEG_INFINITY {
return -1.0;
}
if x.is_nan() {
return f64::NAN;
}
erf_imp(x, false)
}
pub fn erf_inv(z: f64) -> f64 {
if z.abs() < f64::EPSILON {
return 0.0;
}
if z >= 1.0 {
return f64::INFINITY;
}
if z <= -1.0 {
return f64::NEG_INFINITY;
}
if z < 0.0 {
erf_inv_impl(-z, 1.0 - (-z), -1.0)
} else {
erf_inv_impl(z, 1.0 - z, 1.0)
}
}
fn erf_imp(z: f64, mut invert: bool) -> f64 {
if z < 0.0 {
if !invert {
return -erf_imp(-z, false);
}
if z < -0.5 {
return 2.0 - erf_imp(-z, true);
}
return 1.0 + erf_imp(-z, false);
}
let result = if z < 0.5 {
if z < 1e-10 {
(z * 1.125) + (z * 0.003379167095512573896158903121545171688)
} else {
(z * 1.125)
+ (z * evaluate_polynomial(z, ERF_IMP_AN) / evaluate_polynomial(z, ERF_IMP_AD))
}
} else if z < 110.0 {
invert = !invert;
let (r, b) = if z < 0.75 {
(
evaluate_polynomial(z - 0.5, ERF_IMP_BN) / evaluate_polynomial(z - 0.5, ERF_IMP_BD),
f64::from(0.3440242112_f32),
)
} else if z < 1.25 {
(
evaluate_polynomial(z - 0.75, ERF_IMP_CN)
/ evaluate_polynomial(z - 0.75, ERF_IMP_CD),
f64::from(0.419990927_f32),
)
} else if z < 2.25 {
(
evaluate_polynomial(z - 1.25, ERF_IMP_DN)
/ evaluate_polynomial(z - 1.25, ERF_IMP_DD),
f64::from(0.4898625016_f32),
)
} else if z < 3.5 {
(
evaluate_polynomial(z - 2.25, ERF_IMP_EN)
/ evaluate_polynomial(z - 2.25, ERF_IMP_ED),
f64::from(0.5317370892_f32),
)
} else if z < 5.25 {
(
evaluate_polynomial(z - 3.5, ERF_IMP_FN) / evaluate_polynomial(z - 3.5, ERF_IMP_FD),
f64::from(0.5489973426_f32),
)
} else if z < 8.0 {
(
evaluate_polynomial(z - 5.25, ERF_IMP_GN)
/ evaluate_polynomial(z - 5.25, ERF_IMP_GD),
f64::from(0.5571740866_f32),
)
} else if z < 11.5 {
(
evaluate_polynomial(z - 8.0, ERF_IMP_HN) / evaluate_polynomial(z - 8.0, ERF_IMP_HD),
f64::from(0.5609807968_f32),
)
} else if z < 17.0 {
(
evaluate_polynomial(z - 11.5, ERF_IMP_IN)
/ evaluate_polynomial(z - 11.5, ERF_IMP_ID),
f64::from(0.5626493692_f32),
)
} else if z < 24.0 {
(
evaluate_polynomial(z - 17.0, ERF_IMP_JN)
/ evaluate_polynomial(z - 17.0, ERF_IMP_JD),
f64::from(0.5634598136_f32),
)
} else if z < 38.0 {
(
evaluate_polynomial(z - 24.0, ERF_IMP_KN)
/ evaluate_polynomial(z - 24.0, ERF_IMP_KD),
f64::from(0.5638477802_f32),
)
} else if z < 60.0 {
(
evaluate_polynomial(z - 38.0, ERF_IMP_LN)
/ evaluate_polynomial(z - 38.0, ERF_IMP_LD),
f64::from(0.5640528202_f32),
)
} else if z < 85.0 {
(
evaluate_polynomial(z - 60.0, ERF_IMP_MN)
/ evaluate_polynomial(z - 60.0, ERF_IMP_MD),
f64::from(0.5641309023_f32),
)
} else {
(
evaluate_polynomial(z - 85.0, ERF_IMP_NN)
/ evaluate_polynomial(z - 85.0, ERF_IMP_ND),
f64::from(0.5641584396_f32),
)
};
let g = (-z * z).exp() / z;
(g * b) + (g * r)
} else {
invert = !invert;
0.0
};
if invert {
1.0 - result
} else {
result
}
}
fn erf_inv_impl(p: f64, q: f64, s: f64) -> f64 {
let result = if p <= 0.5 {
const Y: f32 = 0.0891314744949340820313;
let g = p * (p + 10.0);
let r = evaluate_polynomial(p, ERV_INV_IMP_AN) / evaluate_polynomial(p, ERV_INV_IMP_AD);
(g * f64::from(Y)) + (g * r)
} else if q >= 0.25 {
const Y: f32 = 2.249481201171875;
let g = (-2.0 * q.ln()).sqrt();
let xs = q - 0.25;
let r = evaluate_polynomial(xs, ERV_INV_IMP_BN) / evaluate_polynomial(xs, ERV_INV_IMP_BD);
g / (f64::from(Y) + r)
} else {
let x = (-q.ln()).sqrt();
if x < 3.0 {
const Y: f32 = 0.807220458984375;
let xs = x - 1.125;
let r =
evaluate_polynomial(xs, ERV_INV_IMP_CN) / evaluate_polynomial(xs, ERV_INV_IMP_CD);
(f64::from(Y) * x) + (r * x)
} else if x < 6.0 {
const Y: f32 = 0.93995571136474609375;
let xs = x - 3.0;
let r =
evaluate_polynomial(xs, ERV_INV_IMP_DN) / evaluate_polynomial(xs, ERV_INV_IMP_DD);
(f64::from(Y) * x) + (r * x)
} else if x < 18.0 {
const Y: f32 = 0.98362827301025390625;
let xs = x - 6.0;
let r =
evaluate_polynomial(xs, ERV_INV_IMP_EN) / evaluate_polynomial(xs, ERV_INV_IMP_ED);
(f64::from(Y) * x) + (r * x)
} else if x < 44.0 {
const Y: f32 = 0.99714565277099609375;
let xs = x - 18.0;
let r =
evaluate_polynomial(xs, ERV_INV_IMP_FN) / evaluate_polynomial(xs, ERV_INV_IMP_FD);
(f64::from(Y) * x) + (r * x)
} else {
const Y: f32 = 0.99941349029541015625;
let xs = x - 44.0;
let r =
evaluate_polynomial(xs, ERV_INV_IMP_GN) / evaluate_polynomial(xs, ERV_INV_IMP_GD);
(f64::from(Y) * x) + (r * x)
}
};
result * s
}
fn evaluate_polynomial(z: f64, coefficients: &[f64]) -> f64 {
let mut coefficients = coefficients.iter().copied().rev();
let Some(last) = coefficients.next() else {
return 0.0;
};
coefficients.fold(last, |sum, coefficient| (sum * z) + coefficient)
}