h1or4dme8/osu-lazer
1
0
Fork 0
mirror of https://github.com/ppy/osu.git synced 2024-12-05 09:42:54 +08:00
Code Issues Releases Wiki Activity

Implement aim probability

Browse Source
This commit is contained in:
Nathen 2024-03-08 12:35:11 -05:00
parent 9018579ecd
commit 8a582eef35
6 changed files with 958 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
osu.Game.Rulesets.Osu/Difficulty/OsuPerformanceCalculator.cs
View File

@ -90,7 +90,6 @@ namespace osu.Game.Rulesets.Osu.Difficulty
double lengthBonus = 0.95 + 0.4 * Math.Min(1.0, totalHits / 2000.0) +
(totalHits > 2000 ? Math.Log10(totalHits / 2000.0) * 0.5 : 0.0);
aimValue *= lengthBonus;
// Penalize misses by assessing # of misses relative to the total # of objects. Default a 3% reduction for any # of misses.
if (effectiveMissCount > 0)

6
osu.Game.Rulesets.Osu/Difficulty/Skills/Aim.cs
View File

@ -11,7 +11,7 @@ namespace osu.Game.Rulesets.Osu.Difficulty.Skills
/// <summary>
/// Represents the skill required to correctly aim at every object in the map with a uniform CircleSize and normalized distances.
/// </summary>
public class Aim : OsuStrainSkill
public class Aim : OsuProbSkill
{
public Aim(Mod[] mods, bool withSliders)
: base(mods)
@ -23,13 +23,11 @@ namespace osu.Game.Rulesets.Osu.Difficulty.Skills
private double currentStrain;
private double skillMultiplier => 23.55;
private double skillMultiplier => 125;
private double strainDecayBase => 0.15;
private double strainDecay(double ms) => Math.Pow(strainDecayBase, ms / 1000);
protected override double CalculateInitialStrain(double time, DifficultyHitObject current) => currentStrain * strainDecay(time - current.Previous(0).StartTime);
protected override double StrainValueAt(DifficultyHitObject current)
{
currentStrain *= strainDecay(current.DeltaTime);

143
osu.Game.Rulesets.Osu/Difficulty/Skills/OsuProbSkill.cs Normal file
View File

@ -0,0 +1,143 @@
// Copyright (c) ppy Pty Ltd <contact@ppy.sh>. Licensed under the MIT Licence.
// See the LICENCE file in the repository root for full licence text.
using System;
using System.Collections.Generic;
using System.Linq;
using osu.Game.Rulesets.Difficulty.Preprocessing;
using osu.Game.Rulesets.Difficulty.Skills;
using osu.Game.Rulesets.Mods;
using osu.Game.Rulesets.Osu.Difficulty.Utils;
namespace osu.Game.Rulesets.Osu.Difficulty.Skills
{
public abstract class OsuProbSkill : Skill
{
protected OsuProbSkill(Mod[] mods)
: base(mods)
{
}
private const double fc_probability = 0.02;
private const int bin_count = 32;
private readonly List<double> difficulties = new List<double>();
/// <summary>
/// Returns the strain value at <see cref="DifficultyHitObject"/>. This value is calculated with or without respect to previous objects.
/// </summary>
protected abstract double StrainValueAt(DifficultyHitObject current);
public override void Process(DifficultyHitObject current)
{
difficulties.Add(StrainValueAt(current));
}
private static double hitProbability(double skill, double difficulty)
{
if (skill <= 0) return 0;
if (difficulty <= 0) return 1;
return SpecialFunctions.Erf(skill / (Math.Sqrt(2) * difficulty));
}
private static double fcProbabilityAtSkill(double skill, IEnumerable<Bin> bins)
{
if (skill <= 0) return 0;
return bins.Aggregate(1.0, (current, bin) => current * bin.FcProbability(skill));
}
private double fcProbabilityAtSkill(double skill)
{
if (skill <= 0) return 0;
return difficulties.Aggregate<double, double>(1, (current, d) => current * hitProbability(skill, d));
}
/// <summary>
/// Create an array of equally spaced bins. Count is linearly interpolated into each bin.
/// For example, if we have bins with values [1,2,3,4,5] and want to insert the value 3.2,
/// we will add 0.8 to the count of 3's and 0.2 to the count of 4's
/// </summary>
private Bin[] createBins(double maxDifficulty)
{
var bins = new Bin[bin_count];
for (int i = 0; i < bin_count; i++)
{
bins[i].Difficulty = maxDifficulty * (i + 1) / bin_count;
}
foreach (double d in difficulties)
{
double binIndex = bin_count * (d / maxDifficulty) - 1;
int lowerBound = (int)Math.Floor(binIndex);
double t = binIndex - lowerBound;
//This can be -1, corresponding to the zero difficulty bucket.
//We don't store that since it doesn't contribute to difficulty
if (lowerBound >= 0)
{
bins[lowerBound].Count += (1 - t);
}
int upperBound = lowerBound + 1;
// this can be == bin_count for the maximum difficulty object, in which case t will be 0 anyway
if (upperBound < bin_count)
{
bins[upperBound].Count += t;
}
}
return bins;
}
private double difficultyValueBinned()
{
double maxDiff = difficulties.Max();
if (maxDiff <= 1e-10) return 0;
var bins = createBins(maxDiff);
double lowerBoundEstimate = 0.5 * maxDiff;
double upperBoundEstimate = 3.0 * maxDiff;
double skill = Chandrupatla.FindRootExpand(
skill => fcProbabilityAtSkill(skill, bins) - fc_probability,
lowerBoundEstimate,
upperBoundEstimate,
accuracy: 1e-4);
return skill;
}
private double difficultyValueExact()
{
double maxDiff = difficulties.Max();
if (maxDiff <= 1e-10) return 0;
double lowerBoundEstimate = 0.5 * maxDiff;
double upperBoundEstimate = 3.0 * maxDiff;
double skill = Chandrupatla.FindRootExpand(
skill => fcProbabilityAtSkill(skill) - fc_probability,
lowerBoundEstimate,
upperBoundEstimate,
accuracy: 1e-4);
return skill;
}
public override double DifficultyValue()
{
if (difficulties.Count == 0)
return 0;
return difficulties.Count < 2 * bin_count ? difficultyValueExact() : difficultyValueBinned();
}
}
}

17
osu.Game.Rulesets.Osu/Difficulty/Utils/Bin.cs Normal file
View File

@ -0,0 +1,17 @@
// Copyright (c) ppy Pty Ltd <contact@ppy.sh>. Licensed under the MIT Licence.
// See the LICENCE file in the repository root for full licence text.
using System;
namespace osu.Game.Rulesets.Osu.Difficulty.Utils
{
public struct Bin
{
public double Difficulty;
public double Count;
public double HitProbability(double skill) => SpecialFunctions.Erf(skill / (Math.Sqrt(2) * Difficulty));
public double FcProbability(double skill) => Math.Pow(HitProbability(skill), Count);
}
}

102
osu.Game.Rulesets.Osu/Difficulty/Utils/Chandrupatla.cs Normal file
View File

@ -0,0 +1,102 @@
// Copyright (c) ppy Pty Ltd <contact@ppy.sh>. Licensed under the MIT Licence.
// See the LICENCE file in the repository root for full licence text.
using System;
namespace osu.Game.Rulesets.Osu.Difficulty.Utils
{
public static class Chandrupatla
{
/// <summary>
/// Finds the root of a <paramref name="function"/>, expanding the bounds if the root is not located within.
/// Expansion only occurs for the upward bound, as this function is optimized for functions of range [0, x),
/// which is useful for finding skill level (skill can never be below 0).
/// </summary>
/// <param name="function">The function of which to find the root.</param>
/// <param name="guessLowerBound">The lower bound of the function inputs.</param>
/// <param name="guessUpperBound">The upper bound of the function inputs.</param>
/// <param name="maxIterations">The maximum number of iterations before the function throws an error.</param>
/// <param name="accuracy">The desired precision in which the root is returned.</param>
/// <param name="expansionFactor">The multiplier on the upper bound when no root is found within the provided bounds.</param>
/// <returns></returns>
public static double FindRootExpand(Func<double, double> function, double guessLowerBound, double guessUpperBound, int maxIterations = 25, double accuracy = 1e-6D, double expansionFactor = 2)
{
double a = guessLowerBound;
double b = guessUpperBound;
double fa = function(a);
double fb = function(b);
while (fa * fb > 0)
{
a = b;
b *= expansionFactor;
fa = function(a);
fb = function(b);
}
double t = 0.5;
for (int i = 0; i < maxIterations; i++)
{
double xt = a + t * (b - a);
double ft = function(xt);
double c;
double fc;
if (Math.Sign(ft) == Math.Sign(fa))
{
c = a;
fc = fa;
}
else
{
c = b;
b = a;
fc = fb;
fb = fa;
}
a = xt;
fa = ft;
double xm, fm;
if (Math.Abs(fa) < Math.Abs(fb))
{
xm = a;
fm = fa;
}
else
{
xm = b;
fm = fb;
}
if (fm == 0)
return xm;
double tol = 2 * accuracy * Math.Abs(xm) + 2 * accuracy;
double tlim = tol / Math.Abs(b - c);
if (tlim > 0.5)
{
return xm;
}
double chi = (a - b) / (c - b);
double phi = (fa - fb) / (fc - fb);
bool iqi = phi * phi < chi && (1 - phi) * (1 - phi) < chi;
if (iqi)
t = fa / (fb - fa) * fc / (fb - fc) + (c - a) / (b - a) * fa / (fc - fa) * fb / (fc - fb);
else
t = 0.5;
t = Math.Min(1 - tlim, Math.Max(tlim, t));
}
return 0;
}
}
}

694
osu.Game.Rulesets.Osu/Difficulty/Utils/SpecialFunctions.cs Normal file
View File

@ -0,0 +1,694 @@
// Copyright (c) ppy Pty Ltd <contact@ppy.sh>. Licensed under the MIT Licence.
// See the LICENCE file in the repository root for full licence text.
// All code is referenced from the following:
// https://github.com/mathnet/mathnet-numerics/blob/master/src/Numerics/SpecialFunctions/Erf.cs
// https://github.com/mathnet/mathnet-numerics/blob/master/src/Numerics/Optimization/NelderMeadSimplex.cs
/*
Copyright (c) 2002-2022 Math.NET
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
namespace osu.Game.Rulesets.Osu.Difficulty.Utils
{
public class SpecialFunctions
{
private const double sqrt2_pi = 2.5066282746310005024157652848110452530069867406099d;
/// <summary>
/// **************************************
/// COEFFICIENTS FOR METHOD ErfImp *
/// **************************************
/// </summary>
/// <summary> Polynomial coefficients for a numerator of ErfImp
/// calculation for Erf(x) in the interval [1e-10, 0.5].
/// </summary>
private static readonly double[] erf_imp_an = { 0.00337916709551257388990745, -0.00073695653048167948530905, -0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468, 0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225 };
/// <summary> Polynomial coefficients for a denominator of ErfImp
/// calculation for Erf(x) in the interval [1e-10, 0.5].
/// </summary>
private static readonly double[] erf_imp_ad = { 1, -0.218088218087924645390535, 0.412542972725442099083918, -0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266, 0.00408165558926174048329689, -0.000615900721557769691924509 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [0.5, 0.75].
/// </summary>
private static readonly double[] erf_imp_bn = { -0.0361790390718262471360258, 0.292251883444882683221149, 0.281447041797604512774415, 0.125610208862766947294894, 0.0274135028268930549240776, 0.00250839672168065762786937 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [0.5, 0.75].
/// </summary>
private static readonly double[] erf_imp_bd = { 1, 1.8545005897903486499845, 1.43575803037831418074962, 0.582827658753036572454135, 0.124810476932949746447682, 0.0113724176546353285778481 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [0.75, 1.25].
/// </summary>
private static readonly double[] erf_imp_cn = { -0.0397876892611136856954425, 0.153165212467878293257683, 0.191260295600936245503129, 0.10276327061989304213645, 0.029637090615738836726027, 0.0046093486780275489468812, 0.000307607820348680180548455 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [0.75, 1.25].
/// </summary>
private static readonly double[] erf_imp_cd = { 1, 1.95520072987627704987886, 1.64762317199384860109595, 0.768238607022126250082483, 0.209793185936509782784315, 0.0319569316899913392596356, 0.00213363160895785378615014 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [1.25, 2.25].
/// </summary>
private static readonly double[] erf_imp_dn = { -0.0300838560557949717328341, 0.0538578829844454508530552, 0.0726211541651914182692959, 0.0367628469888049348429018, 0.00964629015572527529605267, 0.00133453480075291076745275, 0.778087599782504251917881e-4 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [1.25, 2.25].
/// </summary>
private static readonly double[] erf_imp_dd = { 1, 1.75967098147167528287343, 1.32883571437961120556307, 0.552528596508757581287907, 0.133793056941332861912279, 0.0179509645176280768640766, 0.00104712440019937356634038, -0.106640381820357337177643e-7 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [2.25, 3.5].
/// </summary>
private static readonly double[] erf_imp_en = { -0.0117907570137227847827732, 0.014262132090538809896674, 0.0202234435902960820020765, 0.00930668299990432009042239, 0.00213357802422065994322516, 0.00025022987386460102395382, 0.120534912219588189822126e-4 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [2.25, 3.5].
/// </summary>
private static readonly double[] erf_imp_ed = { 1, 1.50376225203620482047419, 0.965397786204462896346934, 0.339265230476796681555511, 0.0689740649541569716897427, 0.00771060262491768307365526, 0.000371421101531069302990367 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [3.5, 5.25].
/// </summary>
private static readonly double[] erf_imp_fn = { -0.00546954795538729307482955, 0.00404190278731707110245394, 0.0054963369553161170521356, 0.00212616472603945399437862, 0.000394984014495083900689956, 0.365565477064442377259271e-4, 0.135485897109932323253786e-5 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [3.5, 5.25].
/// </summary>
private static readonly double[] erf_imp_fd = { 1, 1.21019697773630784832251, 0.620914668221143886601045, 0.173038430661142762569515, 0.0276550813773432047594539, 0.00240625974424309709745382, 0.891811817251336577241006e-4, -0.465528836283382684461025e-11 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [5.25, 8].
/// </summary>
private static readonly double[] erf_imp_gn = { -0.00270722535905778347999196, 0.0013187563425029400461378, 0.00119925933261002333923989, 0.00027849619811344664248235, 0.267822988218331849989363e-4, 0.923043672315028197865066e-6 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [5.25, 8].
/// </summary>
private static readonly double[] erf_imp_gd = { 1, 0.814632808543141591118279, 0.268901665856299542168425, 0.0449877216103041118694989, 0.00381759663320248459168994, 0.000131571897888596914350697, 0.404815359675764138445257e-11 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [8, 11.5].
/// </summary>
private static readonly double[] erf_imp_hn = { -0.00109946720691742196814323, 0.000406425442750422675169153, 0.000274499489416900707787024, 0.465293770646659383436343e-4, 0.320955425395767463401993e-5, 0.778286018145020892261936e-7 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [8, 11.5].
/// </summary>
private static readonly double[] erf_imp_hd = { 1, 0.588173710611846046373373, 0.139363331289409746077541, 0.0166329340417083678763028, 0.00100023921310234908642639, 0.24254837521587225125068e-4 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [11.5, 17].
/// </summary>
private static readonly double[] erf_imp_in = { -0.00056907993601094962855594, 0.000169498540373762264416984, 0.518472354581100890120501e-4, 0.382819312231928859704678e-5, 0.824989931281894431781794e-7 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [11.5, 17].
/// </summary>
private static readonly double[] erf_imp_id = { 1, 0.339637250051139347430323, 0.043472647870310663055044, 0.00248549335224637114641629, 0.535633305337152900549536e-4, -0.117490944405459578783846e-12 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [17, 24].
/// </summary>
private static readonly double[] erf_imp_jn = { -0.000241313599483991337479091, 0.574224975202501512365975e-4, 0.115998962927383778460557e-4, 0.581762134402593739370875e-6, 0.853971555085673614607418e-8 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [17, 24].
/// </summary>
private static readonly double[] erf_imp_jd = { 1, 0.233044138299687841018015, 0.0204186940546440312625597, 0.000797185647564398289151125, 0.117019281670172327758019e-4 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [24, 38].
/// </summary>
private static readonly double[] erf_imp_kn = { -0.000146674699277760365803642, 0.162666552112280519955647e-4, 0.269116248509165239294897e-5, 0.979584479468091935086972e-7, 0.101994647625723465722285e-8 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [24, 38].
/// </summary>
private static readonly double[] erf_imp_kd = { 1, 0.165907812944847226546036, 0.0103361716191505884359634, 0.000286593026373868366935721, 0.298401570840900340874568e-5 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [38, 60].
/// </summary>
private static readonly double[] erf_imp_ln = { -0.583905797629771786720406e-4, 0.412510325105496173512992e-5, 0.431790922420250949096906e-6, 0.993365155590013193345569e-8, 0.653480510020104699270084e-10 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [38, 60].
/// </summary>
private static readonly double[] erf_imp_ld = { 1, 0.105077086072039915406159, 0.00414278428675475620830226, 0.726338754644523769144108e-4, 0.477818471047398785369849e-6 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [60, 85].
/// </summary>
private static readonly double[] erf_imp_mn = { -0.196457797609229579459841e-4, 0.157243887666800692441195e-5, 0.543902511192700878690335e-7, 0.317472492369117710852685e-9 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [60, 85].
/// </summary>
private static readonly double[] erf_imp_md = { 1, 0.052803989240957632204885, 0.000926876069151753290378112, 0.541011723226630257077328e-5, 0.535093845803642394908747e-15 };
/// <summary> Polynomial coefficients for a numerator in ErfImp
/// calculation for Erfc(x) in the interval [85, 110].
/// </summary>
private static readonly double[] erf_imp_nn = { -0.789224703978722689089794e-5, 0.622088451660986955124162e-6, 0.145728445676882396797184e-7, 0.603715505542715364529243e-10 };
/// <summary> Polynomial coefficients for a denominator in ErfImp
/// calculation for Erfc(x) in the interval [85, 110].
/// </summary>
private static readonly double[] erf_imp_nd = { 1, 0.0375328846356293715248719, 0.000467919535974625308126054, 0.193847039275845656900547e-5 };
/// <summary>
/// **************************************
/// COEFFICIENTS FOR METHOD ErfInvImp *
/// **************************************
/// </summary>
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0, 0.5].
/// </summary>
private static readonly double[] erv_inv_imp_an = { -0.000508781949658280665617, -0.00836874819741736770379, 0.0334806625409744615033, -0.0126926147662974029034, -0.0365637971411762664006, 0.0219878681111168899165, 0.00822687874676915743155, -0.00538772965071242932965 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0, 0.5].
/// </summary>
private static readonly double[] erv_inv_imp_ad = { 1, -0.970005043303290640362, -1.56574558234175846809, 1.56221558398423026363, 0.662328840472002992063, -0.71228902341542847553, -0.0527396382340099713954, 0.0795283687341571680018, -0.00233393759374190016776, 0.000886216390456424707504 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.5, 0.75].
/// </summary>
private static readonly double[] erv_inv_imp_bn = { -0.202433508355938759655, 0.105264680699391713268, 8.37050328343119927838, 17.6447298408374015486, -18.8510648058714251895, -44.6382324441786960818, 17.445385985570866523, 21.1294655448340526258, -3.67192254707729348546 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.5, 0.75].
/// </summary>
private static readonly double[] erv_inv_imp_bd = { 1, 6.24264124854247537712, 3.9713437953343869095, -28.6608180499800029974, -20.1432634680485188801, 48.5609213108739935468, 10.8268667355460159008, -22.6436933413139721736, 1.72114765761200282724 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x less than 3.
/// </summary>
private static readonly double[] erv_inv_imp_cn = { -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 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x less than 3.
/// </summary>
private static readonly double[] erv_inv_imp_cd = { 1, 3.46625407242567245975, 5.38168345707006855425, 4.77846592945843778382, 2.59301921623620271374, 0.848854343457902036425, 0.152264338295331783612, 0.01105924229346489121 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 3 and 6.
/// </summary>
private static readonly double[] erv_inv_imp_dn = { -0.0350353787183177984712, -0.00222426529213447927281, 0.0185573306514231072324, 0.00950804701325919603619, 0.00187123492819559223345, 0.000157544617424960554631, 0.460469890584317994083e-5, -0.230404776911882601748e-9, 0.266339227425782031962e-11 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 3 and 6.
/// </summary>
private static readonly double[] erv_inv_imp_dd = { 1, 1.3653349817554063097, 0.762059164553623404043, 0.220091105764131249824, 0.0341589143670947727934, 0.00263861676657015992959, 0.764675292302794483503e-4 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 6 and 18.
/// </summary>
private static readonly double[] erv_inv_imp_en = { -0.0167431005076633737133, -0.00112951438745580278863, 0.00105628862152492910091, 0.000209386317487588078668, 0.149624783758342370182e-4, 0.449696789927706453732e-6, 0.462596163522878599135e-8, -0.281128735628831791805e-13, 0.99055709973310326855e-16 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 6 and 18.
/// </summary>
private static readonly double[] erv_inv_imp_ed = { 1, 0.591429344886417493481, 0.138151865749083321638, 0.0160746087093676504695, 0.000964011807005165528527, 0.275335474764726041141e-4, 0.282243172016108031869e-6 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 18 and 44.
/// </summary>
private static readonly double[] erv_inv_imp_fn = { -0.0024978212791898131227, -0.779190719229053954292e-5, 0.254723037413027451751e-4, 0.162397777342510920873e-5, 0.396341011304801168516e-7, 0.411632831190944208473e-9, 0.145596286718675035587e-11, -0.116765012397184275695e-17 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x between 18 and 44.
/// </summary>
private static readonly double[] erv_inv_imp_fd = { 1, 0.207123112214422517181, 0.0169410838120975906478, 0.000690538265622684595676, 0.145007359818232637924e-4, 0.144437756628144157666e-6, 0.509761276599778486139e-9 };
/// <summary> Polynomial coefficients for a numerator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x greater than 44.
/// </summary>
private static readonly double[] erv_inv_imp_gn = { -0.000539042911019078575891, -0.28398759004727721098e-6, 0.899465114892291446442e-6, 0.229345859265920864296e-7, 0.225561444863500149219e-9, 0.947846627503022684216e-12, 0.135880130108924861008e-14, -0.348890393399948882918e-21 };
/// <summary> Polynomial coefficients for a denominator of ErfInvImp
/// calculation for Erf^-1(z) in the interval [0.75, 1] with x greater than 44.
/// </summary>
private static readonly double[] erv_inv_imp_gd = { 1, 0.0845746234001899436914, 0.00282092984726264681981, 0.468292921940894236786e-4, 0.399968812193862100054e-6, 0.161809290887904476097e-8, 0.231558608310259605225e-11 };
/// <summary>Calculates the error function.</summary>
/// <param name="x">The value to evaluate.</param>
/// <returns>the error function evaluated at given value.</returns>
/// <remarks>
/// <list type="bullet">
/// <item>returns 1 if <c>x == double.PositiveInfinity</c>.</item>
/// <item>returns -1 if <c>x == double.NegativeInfinity</c>.</item>
/// </list>
/// </remarks>
public static double Erf(double x)
{
if (x == 0)
{
return 0;
}
if (double.IsPositiveInfinity(x))
{
return 1;
}
if (double.IsNegativeInfinity(x))
{
return -1;
}
if (double.IsNaN(x))
{
return double.NaN;
}
return erfImp(x, false);
}
/// <summary>Calculates the complementary error function.</summary>
/// <param name="x">The value to evaluate.</param>
/// <returns>the complementary error function evaluated at given value.</returns>
/// <remarks>
/// <list type="bullet">
/// <item>returns 0 if <c>x == double.PositiveInfinity</c>.</item>
/// <item>returns 2 if <c>x == double.NegativeInfinity</c>.</item>
/// </list>
/// </remarks>
public static double Erfc(double x)
{
if (x == 0)
{
return 1;
}
if (double.IsPositiveInfinity(x))
{
return 0;
}
if (double.IsNegativeInfinity(x))
{
return 2;
}
if (double.IsNaN(x))
{
return double.NaN;
}
return erfImp(x, true);
}
/// <summary>Calculates the inverse error function evaluated at z.</summary>
/// <returns>The inverse error function evaluated at given value.</returns>
/// <remarks>
/// <list type="bullet">
/// <item>returns double.PositiveInfinity if <c>z &gt;= 1.0</c>.</item>
/// <item>returns double.NegativeInfinity if <c>z &lt;= -1.0</c>.</item>
/// </list>
/// </remarks>
/// <summary>Calculates the inverse error function evaluated at z.</summary>
/// <param name="z">value to evaluate.</param>
/// <returns>the inverse error function evaluated at Z.</returns>
public static double ErfInv(double z)
{
if (z == 0.0)
{
return 0.0;
}
if (z >= 1.0)
{
return double.PositiveInfinity;
}
if (z <= -1.0)
{
return double.NegativeInfinity;
}
double p, q, s;
if (z < 0)
{
p = -z;
q = 1 - p;
s = -1;
}
else
{
p = z;
q = 1 - z;
s = 1;
}
return erfInvImpl(p, q, s);
}
/// <summary>
/// Implementation of the error function.
/// </summary>
/// <param name="z">Where to evaluate the error function.</param>
/// <param name="invert">Whether to compute 1 - the error function.</param>
/// <returns>the error function.</returns>
private static double erfImp(double z, bool invert)
{
if (z < 0)
{
if (!invert)
{
return -erfImp(-z, false);
}
if (z < -0.5)
{
return 2 - erfImp(-z, true);
}
return 1 + erfImp(-z, false);
}
double result;
// Big bunch of selection statements now to pick which
// implementation to use, try to put most likely options
// first:
if (z < 0.5)
{
// We're going to calculate erf:
if (z < 1e-10)
{
result = (z * 1.125) + (z * 0.003379167095512573896158903121545171688);
}
else
{
// Worst case absolute error found: 6.688618532e-21
result = (z * 1.125) + (z * evaluatePolynomial(z, erf_imp_an) / evaluatePolynomial(z, erf_imp_ad));
}
}
else if (z < 110)
{
// We'll be calculating erfc:
invert = !invert;
double r, b;
if (z < 0.75)
{
// Worst case absolute error found: 5.582813374e-21
r = evaluatePolynomial(z - 0.5, erf_imp_bn) / evaluatePolynomial(z - 0.5, erf_imp_bd);
b = 0.3440242112F;
}
else if (z < 1.25)
{
// Worst case absolute error found: 4.01854729e-21
r = evaluatePolynomial(z - 0.75, erf_imp_cn) / evaluatePolynomial(z - 0.75, erf_imp_cd);
b = 0.419990927F;
}
else if (z < 2.25)
{
// Worst case absolute error found: 2.866005373e-21
r = evaluatePolynomial(z - 1.25, erf_imp_dn) / evaluatePolynomial(z - 1.25, erf_imp_dd);
b = 0.4898625016F;
}
else if (z < 3.5)
{
// Worst case absolute error found: 1.045355789e-21
r = evaluatePolynomial(z - 2.25, erf_imp_en) / evaluatePolynomial(z - 2.25, erf_imp_ed);
b = 0.5317370892F;
}
else if (z < 5.25)
{
// Worst case absolute error found: 8.300028706e-22
r = evaluatePolynomial(z - 3.5, erf_imp_fn) / evaluatePolynomial(z - 3.5, erf_imp_fd);
b = 0.5489973426F;
}
else if (z < 8)
{
// Worst case absolute error found: 1.700157534e-21
r = evaluatePolynomial(z - 5.25, erf_imp_gn) / evaluatePolynomial(z - 5.25, erf_imp_gd);
b = 0.5571740866F;
}
else if (z < 11.5)
{
// Worst case absolute error found: 3.002278011e-22
r = evaluatePolynomial(z - 8, erf_imp_hn) / evaluatePolynomial(z - 8, erf_imp_hd);
b = 0.5609807968F;
}
else if (z < 17)
{
// Worst case absolute error found: 6.741114695e-21
r = evaluatePolynomial(z - 11.5, erf_imp_in) / evaluatePolynomial(z - 11.5, erf_imp_id);
b = 0.5626493692F;
}
else if (z < 24)
{
// Worst case absolute error found: 7.802346984e-22
r = evaluatePolynomial(z - 17, erf_imp_jn) / evaluatePolynomial(z - 17, erf_imp_jd);
b = 0.5634598136F;
}
else if (z < 38)
{
// Worst case absolute error found: 2.414228989e-22
r = evaluatePolynomial(z - 24, erf_imp_kn) / evaluatePolynomial(z - 24, erf_imp_kd);
b = 0.5638477802F;
}
else if (z < 60)
{
// Worst case absolute error found: 5.896543869e-24
r = evaluatePolynomial(z - 38, erf_imp_ln) / evaluatePolynomial(z - 38, erf_imp_ld);
b = 0.5640528202F;
}
else if (z < 85)
{
// Worst case absolute error found: 3.080612264e-21
r = evaluatePolynomial(z - 60, erf_imp_mn) / evaluatePolynomial(z - 60, erf_imp_md);
b = 0.5641309023F;
}
else
{
// Worst case absolute error found: 8.094633491e-22
r = evaluatePolynomial(z - 85, erf_imp_nn) / evaluatePolynomial(z - 85, erf_imp_nd);
b = 0.5641584396F;
}
double g = Math.Exp(-z * z) / z;
result = (g * b) + (g * r);
}
else
{
// Any value of z larger than 28 will underflow to zero:
result = 0;
invert = !invert;
}
if (invert)
{
result = 1 - result;
}
return result;
}
/// <summary>Calculates the complementary inverse error function evaluated at z.</summary>
/// <returns>The complementary inverse error function evaluated at given value.</returns>
/// <remarks> We have tested this implementation against the arbitrary precision mpmath library
/// and found cases where we can only guarantee 9 significant figures correct.
/// <list type="bullet">
/// <item>returns double.PositiveInfinity if <c>z &lt;= 0.0</c>.</item>
/// <item>returns double.NegativeInfinity if <c>z &gt;= 2.0</c>.</item>
/// </list>
/// </remarks>
/// <summary>calculates the complementary inverse error function evaluated at z.</summary>
/// <param name="z">value to evaluate.</param>
/// <returns>the complementary inverse error function evaluated at Z.</returns>
public static double ErfcInv(double z)
{
if (z <= 0.0)
{
return double.PositiveInfinity;
}
if (z >= 2.0)
{
return double.NegativeInfinity;
}
double p, q, s;
if (z > 1)
{
q = 2 - z;
p = 1 - q;
s = -1;
}
else
{
p = 1 - z;
q = z;
s = 1;
}
return erfInvImpl(p, q, s);
}
/// <summary>
/// The implementation of the inverse error function.
/// </summary>
/// <param name="p">First intermediate parameter.</param>
/// <param name="q">Second intermediate parameter.</param>
/// <param name="s">Third intermediate parameter.</param>
/// <returns>the inverse error function.</returns>
private static double erfInvImpl(double p, double q, double s)
{
double result;
if (p <= 0.5)
{
// Evaluate inverse erf using the rational approximation:
//
// x = p(p+10)(Y+R(p))
//
// Where Y is a constant, and R(p) is optimized for a low
// absolute error compared to |Y|.
//
// double: Max error found: 2.001849e-18
// long double: Max error found: 1.017064e-20
// Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
const float y = 0.0891314744949340820313f;
double g = p * (p + 10);
double r = evaluatePolynomial(p, erv_inv_imp_an) / evaluatePolynomial(p, erv_inv_imp_ad);
result = (g * y) + (g * r);
}
else if (q >= 0.25)
{
// Rational approximation for 0.5 > q >= 0.25
//
// x = sqrt(-2*log(q)) / (Y + R(q))
//
// Where Y is a constant, and R(q) is optimized for a low
// absolute error compared to Y.
//
// double : Max error found: 7.403372e-17
// long double : Max error found: 6.084616e-20
// Maximum Deviation Found (error term) 4.811e-20
const float y = 2.249481201171875f;
double g = Math.Sqrt(-2 * Math.Log(q));
double xs = q - 0.25;
double r = evaluatePolynomial(xs, erv_inv_imp_bn) / evaluatePolynomial(xs, erv_inv_imp_bd);
result = g / (y + r);
}
else
{
// For q < 0.25 we have a series of rational approximations all
// of the general form:
//
// let: x = sqrt(-log(q))
//
// Then the result is given by:
//
// x(Y+R(x-B))
//
// where Y is a constant, B is the lowest value of x for which
// the approximation is valid, and R(x-B) is optimized for a low
// absolute error compared to Y.
//
// Note that almost all code will really go through the first
// or maybe second approximation. After than we're dealing with very
// small input values indeed: 80 and 128 bit long double's go all the
// way down to ~ 1e-5000 so the "tail" is rather long...
double x = Math.Sqrt(-Math.Log(q));
if (x < 3)
{
// Max error found: 1.089051e-20
const float y = 0.807220458984375f;
double xs = x - 1.125;
double r = evaluatePolynomial(xs, erv_inv_imp_cn) / evaluatePolynomial(xs, erv_inv_imp_cd);
result = (y * x) + (r * x);
}
else if (x < 6)
{
// Max error found: 8.389174e-21
const float y = 0.93995571136474609375f;
double xs = x - 3;
double r = evaluatePolynomial(xs, erv_inv_imp_dn) / evaluatePolynomial(xs, erv_inv_imp_dd);
result = (y * x) + (r * x);
}
else if (x < 18)
{
// Max error found: 1.481312e-19
const float y = 0.98362827301025390625f;
double xs = x - 6;
double r = evaluatePolynomial(xs, erv_inv_imp_en) / evaluatePolynomial(xs, erv_inv_imp_ed);
result = (y * x) + (r * x);
}
else if (x < 44)
{
// Max error found: 5.697761e-20
const float y = 0.99714565277099609375f;
double xs = x - 18;
double r = evaluatePolynomial(xs, erv_inv_imp_fn) / evaluatePolynomial(xs, erv_inv_imp_fd);
result = (y * x) + (r * x);
}
else
{
// Max error found: 1.279746e-20
const float y = 0.99941349029541015625f;
double xs = x - 44;
double r = evaluatePolynomial(xs, erv_inv_imp_gn) / evaluatePolynomial(xs, erv_inv_imp_gd);
result = (y * x) + (r * x);
}
}
return s * result;
}
/// <summary>
/// Evaluate a polynomial at point x.
/// Coefficients are ordered ascending by power with power k at index k.
/// Example: coefficients [3,-1,2] represent y=2x^2-x+3.
/// </summary>
/// <param name="z">The location where to evaluate the polynomial at.</param>
/// <param name="coefficients">The coefficients of the polynomial, coefficient for power k at index k.</param>
/// <exception cref="ArgumentNullException">
/// <paramref name="coefficients"/> is a null reference.
/// </exception>
private static double evaluatePolynomial(double z, params double[] coefficients)
{
// 2020-10-07 jbialogrodzki #730 Since this is public API we should probably
// handle null arguments? It doesn't seem to have been done consistently in this class though.
if (coefficients == null)
{
throw new ArgumentNullException(nameof(coefficients));
}
// 2020-10-07 jbialogrodzki #730 Zero polynomials need explicit handling.
// Without this check, we attempted to peek coefficients at negative indices!
int n = coefficients.Length;
if (n == 0)
{
return 0;
}
double sum = coefficients[n - 1];
for (int i = n - 2; i >= 0; --i)
{
sum *= z;
sum += coefficients[i];
}
return sum;
}
}
}