mirror of
https://github.com/MindWorkAI/AI-Studio.git
synced 2026-11-11 17:47:21 +00:00
Added the Mathematics tools as a prototype (#1162)
This commit is contained in:
Thorsten Sommer authored and 
GitHub committed
1 parent
ce846e540e
commit
9c98bdae86
49 files changed
+4373
-11
No files matched your search
@@ -0,0 +1,260 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Integrates a function numerically by double exponential quadrature, also over infinite intervals.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// A substitution moves the points of a trapezoidal sum doubly exponentially close to the bounds of
|
||||
/// the interval, or out to infinity: tanh-sinh on [a, b], exp-sinh on [a, oo) and (-oo, b], and
|
||||
/// sinh-sinh on the whole line. Integrands which are analytic inside the interval converge to about
|
||||
/// fourteen digits within a few hundred points, singularities at the bounds included, such as ln(x) or
|
||||
/// 1/sqrt(x) at 0. The points never lie on a bound, where an integrand such as x^3/(e^x - 1) has no
|
||||
/// value.<br/><br/>
|
||||
/// Each level halves the step, and the difference between two levels estimates the error. The terms
|
||||
/// left out at the ends add to it whenever the sum had to stop before they became negligible: at the
|
||||
/// largest parameter, at a bound which a double cannot tell apart from its neighbors, or where the
|
||||
/// integrand has no value. That is what keeps an integral like sin(x)/x over [0, oo), which the
|
||||
/// substitution cannot settle, from converging by chance. A singularity at a bound other than 0 is
|
||||
/// reached only to about seven digits, since a double near such a bound carries few digits of its
|
||||
/// distance to it; the estimated error says so.<br/><br/>
|
||||
/// The integrand is a double function, compiled by AngouriMath, so the result carries at most the
|
||||
/// digits of a double.
|
||||
/// </remarks>
|
||||
public static class DoubleExponentialQuadrature
|
||||
{
|
||||
/// <summary>
|
||||
/// The finest step is 2^-MAX_LEVEL.
|
||||
/// </summary>
|
||||
private const int MAX_LEVEL = 10;
|
||||
|
||||
/// <summary>
|
||||
/// The coarsest step which may end the integration: a coarser sum can match the next one by chance.
|
||||
/// </summary>
|
||||
private const int MIN_CONVERGED_LEVEL = 3;
|
||||
|
||||
/// <summary>
|
||||
/// Two levels which differ by less than this, relative to the integral, end the refinement.
|
||||
/// </summary>
|
||||
private const double RELATIVE_TOLERANCE = 1e-13;
|
||||
|
||||
/// <summary>
|
||||
/// The largest estimated relative error of a converged result: six digits at least.
|
||||
/// </summary>
|
||||
private const double ACCEPTABLE_RELATIVE_ERROR = 1e-6;
|
||||
|
||||
/// <summary>
|
||||
/// A term this much smaller than the sum so far no longer counts.
|
||||
/// </summary>
|
||||
private const double NEGLIGIBLE = 1e-18;
|
||||
|
||||
/// <summary>
|
||||
/// How many negligible terms in a row end one side of a sum.
|
||||
/// </summary>
|
||||
private const int NEGLIGIBLE_IN_A_ROW = 3;
|
||||
|
||||
/// <summary>
|
||||
/// The largest parameter on a finite interval: there, the points lie within 1e-37 of the bounds.
|
||||
/// </summary>
|
||||
private const double FINITE_MAX_PARAMETER = 4.0;
|
||||
|
||||
/// <summary>
|
||||
/// The largest parameter towards infinity: there, the points lie beyond 1e50.
|
||||
/// </summary>
|
||||
private const double INFINITE_MAX_PARAMETER = 5.0;
|
||||
|
||||
/// <summary>
|
||||
/// How much more the terms left out at an end may add up to than estimated. The last terms before
|
||||
/// a bound are also the least accurate ones: 1 - x carries few digits for an x close to 1, which a
|
||||
/// singularity such as 1/sqrt(1 - x) magnifies.
|
||||
/// </summary>
|
||||
private const double TAIL_SAFETY_FACTOR = 10;
|
||||
|
||||
private const int EVALUATIONS_PER_CANCELLATION_CHECK = 256;
|
||||
|
||||
private const double HALF_PI = Math.PI / 2;
|
||||
|
||||
/// <summary>
|
||||
/// Integrates a function over an interval.
|
||||
/// </summary>
|
||||
/// <param name="integrand">The function, which returns NaN or an infinity where it has no value.</param>
|
||||
/// <param name="lower">The lower bound, which may be negative infinity.</param>
|
||||
/// <param name="upper">The upper bound, which may be positive infinity; below the lower bound, the integral changes its sign.</param>
|
||||
/// <param name="token">The cancellation token.</param>
|
||||
/// <returns>The result, with its estimated error and how the integration ended.</returns>
|
||||
public static QuadratureResult Integrate(Func<double, double> integrand, double lower, double upper, CancellationToken token = default)
|
||||
{
|
||||
if (double.IsNaN(lower) || double.IsNaN(upper))
|
||||
throw new ArgumentException("A bound of the integral is not a number.");
|
||||
|
||||
if (lower.Equals(upper))
|
||||
return new(0, 0, QuadratureStatus.CONVERGED, null, 0);
|
||||
|
||||
if (lower > upper)
|
||||
{
|
||||
var reversed = Integrate(integrand, upper, lower, token);
|
||||
return reversed with { Value = -reversed.Value };
|
||||
}
|
||||
|
||||
var interval = new Interval(lower, upper);
|
||||
var evaluations = 0;
|
||||
var previous = double.NaN;
|
||||
LevelSum sum = default;
|
||||
for (var level = 0; level <= MAX_LEVEL; level++)
|
||||
{
|
||||
sum = Sum(integrand, interval, Math.ScaleB(1.0, -level), ref evaluations, token);
|
||||
if (sum.NoValueAt is { } point)
|
||||
return new(double.NaN, double.PositiveInfinity, QuadratureStatus.NO_VALUE, point, evaluations);
|
||||
|
||||
var difference = Math.Abs(sum.Integral - previous);
|
||||
var scale = Math.Max(Math.Abs(sum.Integral), sum.AbsoluteIntegral);
|
||||
var settled = level >= MIN_CONVERGED_LEVEL && difference <= RELATIVE_TOLERANCE * scale;
|
||||
if (settled || level == MAX_LEVEL)
|
||||
return Result(sum, difference + sum.TailError, scale, evaluations);
|
||||
|
||||
previous = sum.Integral;
|
||||
}
|
||||
|
||||
return Result(sum, double.PositiveInfinity, 0, evaluations);
|
||||
}
|
||||
|
||||
private static QuadratureResult Result(LevelSum sum, double error, double scale, int evaluations)
|
||||
{
|
||||
if (!double.IsFinite(sum.Integral) || double.IsNaN(error))
|
||||
return new(sum.Integral, double.PositiveInfinity, QuadratureStatus.NOT_CONVERGED, null, evaluations);
|
||||
|
||||
var status = error <= ACCEPTABLE_RELATIVE_ERROR * scale ? QuadratureStatus.CONVERGED : QuadratureStatus.NOT_CONVERGED;
|
||||
return new(sum.Integral, error, status, null, evaluations);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The trapezoidal sum of one level, from the middle of the parameter range outwards on both sides.
|
||||
/// </summary>
|
||||
private static LevelSum Sum(Func<double, double> integrand, Interval interval, double step, ref int evaluations, CancellationToken token)
|
||||
{
|
||||
var center = interval.NodeAt(0);
|
||||
var centerValue = integrand(center.X);
|
||||
evaluations++;
|
||||
if (!double.IsFinite(centerValue))
|
||||
return new(0, 0, 0, center.X);
|
||||
|
||||
var sum = center.Weight * centerValue;
|
||||
var absoluteSum = Math.Abs(sum);
|
||||
var tail = 0.0;
|
||||
foreach (var direction in (ReadOnlySpan<int>)[1, -1])
|
||||
{
|
||||
// Unknown before the first term of this side, so a missing value there is never negligible:
|
||||
var lastTerm = double.PositiveInfinity;
|
||||
var lastParameter = 0.0;
|
||||
var negligibleInARow = 0;
|
||||
for (var index = 1; ; index++)
|
||||
{
|
||||
var parameter = direction * index * step;
|
||||
if (Math.Abs(parameter) > interval.MaxParameter)
|
||||
{
|
||||
tail += Tail(lastTerm, lastParameter, step);
|
||||
break;
|
||||
}
|
||||
|
||||
var node = interval.NodeAt(parameter);
|
||||
if (node.AtBound)
|
||||
{
|
||||
tail += Tail(lastTerm, lastParameter, step);
|
||||
break;
|
||||
}
|
||||
|
||||
if (++evaluations % EVALUATIONS_PER_CANCELLATION_CHECK == 0)
|
||||
token.ThrowIfCancellationRequested();
|
||||
|
||||
var value = integrand(node.X);
|
||||
if (!double.IsFinite(value))
|
||||
{
|
||||
// Near a bound, an integrand such as x^3/(e^x - 1) loses its value to rounding, where its terms no longer count:
|
||||
if (Math.Abs(lastTerm) <= NEGLIGIBLE * Math.Abs(sum))
|
||||
{
|
||||
tail += Tail(lastTerm, lastParameter, step);
|
||||
break;
|
||||
}
|
||||
|
||||
return new(0, 0, 0, node.X);
|
||||
}
|
||||
|
||||
var term = node.Weight * value;
|
||||
sum += term;
|
||||
absoluteSum += Math.Abs(term);
|
||||
lastTerm = term;
|
||||
lastParameter = parameter;
|
||||
if (Math.Abs(term) <= NEGLIGIBLE * Math.Abs(sum) && sum != 0)
|
||||
{
|
||||
if (++negligibleInARow >= NEGLIGIBLE_IN_A_ROW)
|
||||
break;
|
||||
}
|
||||
else
|
||||
negligibleInARow = 0;
|
||||
}
|
||||
}
|
||||
|
||||
return new(step * sum, step * tail, step * absoluteSum, null);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// What the terms beyond the last one of a side may add up to, before the sum is multiplied by the
|
||||
/// step: as much as the last term, or the integral of terms which decay from it at the pace of the
|
||||
/// substitution, whichever is more, times TAIL_SAFETY_FACTOR.
|
||||
/// </summary>
|
||||
private static double Tail(double lastTerm, double lastParameter, double step) => TAIL_SAFETY_FACTOR * Math.Abs(lastTerm) * Math.Max(1, 1 / (HALF_PI * Math.Cosh(lastParameter) * step));
|
||||
|
||||
/// <summary>
|
||||
/// The sum of one level: the integral, the estimated error of what its sides left out, and the
|
||||
/// integral of the magnitude, which scales the tolerance of an integral close to zero.
|
||||
/// </summary>
|
||||
private readonly record struct LevelSum(double Integral, double TailError, double AbsoluteIntegral, double? NoValueAt);
|
||||
|
||||
/// <summary>
|
||||
/// A point of the sum and its weight; at a bound when a double cannot tell the point from the bound.
|
||||
/// </summary>
|
||||
private readonly record struct QuadratureNode(double X, double Weight, bool AtBound);
|
||||
|
||||
/// <summary>
|
||||
/// The substitution which fits the bounds.
|
||||
/// </summary>
|
||||
private readonly record struct Interval(double Lower, double Upper)
|
||||
{
|
||||
public double MaxParameter => double.IsFinite(this.Lower) && double.IsFinite(this.Upper) ? FINITE_MAX_PARAMETER : INFINITE_MAX_PARAMETER;
|
||||
|
||||
public QuadratureNode NodeAt(double parameter)
|
||||
{
|
||||
var inner = HALF_PI * Math.Sinh(parameter);
|
||||
var derivative = HALF_PI * Math.Cosh(parameter);
|
||||
if (double.IsFinite(this.Lower) && double.IsFinite(this.Upper))
|
||||
{
|
||||
// tanh-sinh, with the distance to the nearer bound computed directly: 1 - tanh(u) = 2 / (e^(2u) + 1).
|
||||
var halfWidth = (this.Upper - this.Lower) / 2;
|
||||
var distance = halfWidth * 2 / (Math.Exp(2 * Math.Abs(inner)) + 1);
|
||||
var x = parameter >= 0 ? this.Upper - distance : this.Lower + distance;
|
||||
var secant = 2 / (Math.Exp(inner) + Math.Exp(-inner));
|
||||
var atBound = distance == 0 || x.Equals(parameter >= 0 ? this.Upper : this.Lower);
|
||||
return new(x, halfWidth * derivative * secant * secant, atBound);
|
||||
}
|
||||
|
||||
if (double.IsFinite(this.Lower))
|
||||
{
|
||||
// exp-sinh towards positive infinity:
|
||||
var offset = Math.Exp(inner);
|
||||
var x = this.Lower + offset;
|
||||
return new(x, derivative * offset, !double.IsFinite(x) || x.Equals(this.Lower));
|
||||
}
|
||||
|
||||
if (double.IsFinite(this.Upper))
|
||||
{
|
||||
// exp-sinh towards negative infinity, mirrored:
|
||||
var offset = Math.Exp(inner);
|
||||
var x = this.Upper - offset;
|
||||
return new(x, derivative * offset, !double.IsFinite(x) || x.Equals(this.Upper));
|
||||
}
|
||||
|
||||
// sinh-sinh over the whole line:
|
||||
var point = Math.Sinh(inner);
|
||||
return new(point, derivative * Math.Cosh(inner), !double.IsFinite(point));
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,12 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// A computation which the Mathematics tools refused or could not finish.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// The message goes to the model as the result of the tool, in English, like the other messages
|
||||
/// about the arguments of a tool. It states what happened and, where the input was the cause, how
|
||||
/// AngouriMath reads the input and how it is written instead.
|
||||
/// </remarks>
|
||||
/// <param name="message">What happened.</param>
|
||||
public sealed class MathComputationException(string message) : Exception(message);
|
||||
@@ -0,0 +1,134 @@
|
||||
using System.Globalization;
|
||||
|
||||
using AIStudio.Chat;
|
||||
|
||||
using AngouriMath;
|
||||
using AngouriMath.Core.Exceptions;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Runs the computations of the Mathematics tools, one at a time and within a time budget.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath parses and simplifies recursively. About 500 nested parentheses overflow the stack of
|
||||
/// a thread pool thread, which ends the whole process, so every computation runs on a thread of its
|
||||
/// own with a stack of 64 MB. That holds the deepest nesting an expression within
|
||||
/// MathLimits.MAX_EXPRESSION_LENGTH can have.<br/><br/>
|
||||
/// Only one computation runs at a time, because AI Studio also runs on a Raspberry Pi. AngouriMath
|
||||
/// stops at the cancellation token of its thread wherever it checks it, and that token ends with the
|
||||
/// budget or with the request. Some computations never check it, such as an exact power with
|
||||
/// millions of digits; MathGuards refuses those which can be told beforehand. One which still runs
|
||||
/// past its budget keeps the slot until it ends: the call gives up and says so, and the next
|
||||
/// computations wait for the slot instead of taking another core.<br/><br/>
|
||||
/// The settings of AngouriMath belong to the flow which sets them. Each computation reads
|
||||
/// expressions explicitly and solves without Newton's method, see AngouriMathBehaviorTests.
|
||||
/// </remarks>
|
||||
/// <param name="logger">The logger.</param>
|
||||
public sealed class MathComputationService(ILogger<MathComputationService> logger)
|
||||
{
|
||||
private const int STACK_SIZE_BYTES = 64 * 1024 * 1024;
|
||||
|
||||
private const int MAX_MESSAGE_LENGTH = 400;
|
||||
|
||||
/// <summary>
|
||||
/// How long a computation may take, and how long it waits for the slot before that.
|
||||
/// </summary>
|
||||
public static readonly TimeSpan DEFAULT_TIME_BUDGET = TimeSpan.FromSeconds(10);
|
||||
|
||||
/// <summary>
|
||||
/// How long a call waits beyond the budget for a computation to stop at its cancellation token.
|
||||
/// </summary>
|
||||
private static readonly TimeSpan STOP_GRACE = TimeSpan.FromSeconds(1);
|
||||
|
||||
private readonly SemaphoreSlim slot = new(1, 1);
|
||||
|
||||
/// <summary>
|
||||
/// How long a computation may take. Tests shorten it.
|
||||
/// </summary>
|
||||
internal TimeSpan TimeBudget { get; init; } = DEFAULT_TIME_BUDGET;
|
||||
|
||||
/// <summary>
|
||||
/// Runs a computation on a thread of its own, once the slot is free.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Parse the expressions inside the computation, too: the parser needs the stack of its thread.
|
||||
/// An exception of AngouriMath comes back as a MathComputationException with its message.
|
||||
/// </remarks>
|
||||
/// <param name="computation">The computation, which gets the token of its budget.</param>
|
||||
/// <param name="token">The cancellation token of the request.</param>
|
||||
/// <typeparam name="T">The type of the result.</typeparam>
|
||||
/// <returns>The result of the computation.</returns>
|
||||
/// <exception cref="MathComputationException">The computation was refused, failed, or did not finish within its budget.</exception>
|
||||
public async Task<T> ComputeAsync<T>(Func<CancellationToken, T> computation, CancellationToken token = default)
|
||||
{
|
||||
if (!await this.slot.WaitAsync(this.TimeBudget, token))
|
||||
throw new MathComputationException($"Another computation is still running, so this one could not start within {this.BudgetInSeconds} seconds.");
|
||||
|
||||
// The thread disposes the budget and frees the slot once its computation ended. Should the thread not
|
||||
// even start, the slot is freed below, and the budget is left to the garbage collector:
|
||||
var budget = CancellationTokenSource.CreateLinkedTokenSource(token);
|
||||
var completion = new TaskCompletionSource<T>(TaskCreationOptions.RunContinuationsAsynchronously);
|
||||
var thread = new Thread(() => this.Run(computation, budget, completion, token), STACK_SIZE_BYTES)
|
||||
{
|
||||
IsBackground = true,
|
||||
Name = "Mathematics",
|
||||
};
|
||||
|
||||
try
|
||||
{
|
||||
budget.CancelAfter(this.TimeBudget);
|
||||
thread.Start();
|
||||
}
|
||||
catch
|
||||
{
|
||||
this.slot.Release();
|
||||
throw;
|
||||
}
|
||||
|
||||
try
|
||||
{
|
||||
return await completion.Task.WaitAsync(this.TimeBudget + STOP_GRACE, token);
|
||||
}
|
||||
catch (TimeoutException)
|
||||
{
|
||||
logger.LogWarning("A computation of the Mathematics tools ran past its time budget of {Seconds} seconds without stopping. It keeps the computing slot until it ends.", this.BudgetInSeconds);
|
||||
throw new MathComputationException($"The computation did not finish within {this.BudgetInSeconds} seconds and could not be stopped. It keeps running until it ends, and other computations wait for it meanwhile.");
|
||||
}
|
||||
}
|
||||
|
||||
private string BudgetInSeconds => this.TimeBudget.TotalSeconds.ToString(CultureInfo.InvariantCulture);
|
||||
|
||||
private void Run<T>(Func<CancellationToken, T> computation, CancellationTokenSource budget, TaskCompletionSource<T> completion, CancellationToken token)
|
||||
{
|
||||
try
|
||||
{
|
||||
MathS.Multithreading.SetLocalCancellationToken(budget.Token);
|
||||
using var explicitParsing = MathS.Settings.ExplicitParsingOnly.Set(true);
|
||||
using var withoutNewton = MathS.Settings.AllowNewton.Set(false);
|
||||
completion.TrySetResult(computation(budget.Token));
|
||||
}
|
||||
catch (OperationCanceledException) when (token.IsCancellationRequested)
|
||||
{
|
||||
completion.TrySetCanceled(token);
|
||||
}
|
||||
catch (OperationCanceledException) when (budget.IsCancellationRequested)
|
||||
{
|
||||
completion.TrySetException(new MathComputationException($"The computation did not finish within {this.BudgetInSeconds} seconds."));
|
||||
}
|
||||
catch (AngouriMathBaseException exception)
|
||||
{
|
||||
logger.LogInformation("AngouriMath could not compute an expression: {Message}", exception.Message);
|
||||
completion.TrySetException(new MathComputationException($"The tool could not compute this:{exception.Message.Shorten(MAX_MESSAGE_LENGTH)}"));
|
||||
}
|
||||
catch (Exception exception)
|
||||
{
|
||||
completion.TrySetException(exception);
|
||||
}
|
||||
finally
|
||||
{
|
||||
budget.Dispose();
|
||||
this.slot.Release();
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,165 @@
|
||||
using System.Globalization;
|
||||
using System.Text.RegularExpressions;
|
||||
|
||||
using AIStudio.Chat;
|
||||
|
||||
using AngouriMath;
|
||||
using AngouriMath.Core.Exceptions;
|
||||
|
||||
using PeterO.Numbers;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Reads an expression of the model into AngouriMath, strictly.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath reads some spellings as something else than they say, silently, see
|
||||
/// AngouriMathBehaviorTests: without explicit parsing, x2 is x^2 and foo(3) is foo times 3, and
|
||||
/// even with it, a != b is the factorial of a equated to b. This reader refuses such spellings
|
||||
/// instead of letting a model compute what it did not mean, and states how AngouriMath writes the
|
||||
/// thing instead. Its parse errors list over a thousand characters of expected tokens; the reader
|
||||
/// keeps the column and the cause.<br/><br/>
|
||||
/// Call it inside a computation of MathComputationService: the parser needs the stack of its thread.
|
||||
/// </remarks>
|
||||
public static partial class MathExpressionReader
|
||||
{
|
||||
private const int MAX_REASON_LENGTH = 200;
|
||||
|
||||
/// <summary>
|
||||
/// Reads an expression, or refuses it.
|
||||
/// </summary>
|
||||
/// <param name="text">The expression as the model wrote it.</param>
|
||||
/// <param name="argumentName">The argument the expression came in, named in every message.</param>
|
||||
/// <returns>The expression as AngouriMath reads it.</returns>
|
||||
/// <exception cref="MathComputationException">The expression is empty, too long, uses a spelling AngouriMath reads otherwise, or cannot be read.</exception>
|
||||
public static Entity Read(string text, string argumentName)
|
||||
{
|
||||
if (string.IsNullOrWhiteSpace(text))
|
||||
throw new MathComputationException($"Argument '{argumentName}' is empty.");
|
||||
|
||||
if (text.Length > MathLimits.MAX_EXPRESSION_LENGTH)
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"Argument '{argumentName}' has {text.Length:N0} characters, but at most {MathLimits.MAX_EXPRESSION_LENGTH:N0} are read."));
|
||||
|
||||
if (text.Contains("!=", StringComparison.Ordinal))
|
||||
throw new MathComputationException($"Argument '{argumentName}' contains '!=', which the tool reads as a factorial followed by '='. Not-equal is written '<>'.");
|
||||
|
||||
if (text.Contains('|'))
|
||||
throw new MathComputationException($"Argument '{argumentName}' contains '|', which the tool reads as divisibility. The absolute value is written abs(x), and 'or' joins two statements.");
|
||||
|
||||
var equalsSigns = CountEqualsSigns(text);
|
||||
if (equalsSigns > 1)
|
||||
throw new MathComputationException($"Argument '{argumentName}' contains {equalsSigns} equals signs, but an expression holds one equation at most. The tool would read a chain of them as several equations joined by 'and'.");
|
||||
|
||||
var expression = Parse(text, argumentName);
|
||||
if (expression.Vars.Any(variable => variable.Name == "oo"))
|
||||
throw new MathComputationException($"Argument '{argumentName}' uses 'oo', which the tool reads as a variable. Infinity is written +oo or -oo.");
|
||||
|
||||
return expression;
|
||||
}
|
||||
|
||||
private static Entity Parse(string text, string argumentName)
|
||||
{
|
||||
try
|
||||
{
|
||||
using var explicitParsing = MathS.Settings.ExplicitParsingOnly.Set(true);
|
||||
Entity parsed;
|
||||
|
||||
// Without the cache, an expression parsed elsewhere with other settings cannot come back here:
|
||||
using (MathS.Settings.DowncastingEnabled.Set(false))
|
||||
parsed = MathS.FromString(text, useCache: false);
|
||||
|
||||
return WithExactLiterals(parsed);
|
||||
}
|
||||
catch (MissingOperatorParseException exception)
|
||||
{
|
||||
throw new MathComputationException(DescribeMissingOperator(exception.Message, argumentName));
|
||||
}
|
||||
catch (FunctionArgumentCountException exception)
|
||||
{
|
||||
throw new MathComputationException($"Argument '{argumentName}' cannot be read: {exception.Message}.");
|
||||
}
|
||||
catch (ParseException exception)
|
||||
{
|
||||
throw new MathComputationException(DescribeParseError(exception.Message, argumentName));
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Turns every decimal literal into the fraction it writes, e.g., 1.5e-3 into 3/2000.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath downcasts while it parses: it reads 0.1 as 1/10, but 6.62607015e-34 as an inexact
|
||||
/// decimal and 1e-60 as 0, see AngouriMathBehaviorTests. Simplifying then drops an inexact decimal
|
||||
/// from a sum, so 1 + 1e-40 would be exactly 1. The reader therefore parses without downcasting,
|
||||
/// which leaves every number a decimal, and turns each into its fraction or integer. A literal
|
||||
/// with an exponent beyond MathLimits.MAX_EXACT_DIGITS stays a decimal.
|
||||
/// </remarks>
|
||||
private static Entity WithExactLiterals(Entity expression) => expression.Replace(node => node is Entity.Number.Real literal and not Entity.Number.Rational ? ExactValueOf(literal) : node);
|
||||
|
||||
private static Entity ExactValueOf(Entity.Number.Real literal)
|
||||
{
|
||||
var value = literal.EDecimal;
|
||||
if (value.IsNaN() || value.IsInfinity() || value.Exponent.Abs().CompareTo(EInteger.FromInt32(MathLimits.MAX_EXACT_DIGITS)) > 0)
|
||||
return literal;
|
||||
|
||||
var fraction = ERational.FromEDecimal(value);
|
||||
return fraction.Denominator.Equals(EInteger.One) ? Entity.Number.Integer.Create(fraction.Numerator) : Entity.Number.Rational.Create(fraction.Numerator, fraction.Denominator);
|
||||
}
|
||||
|
||||
private static string DescribeMissingOperator(string message, string argumentName)
|
||||
{
|
||||
var match = MissingOperatorTokens().Match(message);
|
||||
if (!match.Success)
|
||||
return $"Argument '{argumentName}' lacks an operator between two operands. Multiplication is written '*', as in 2*x.";
|
||||
|
||||
var left = match.Groups["left"].Value;
|
||||
var right = match.Groups["right"].Value;
|
||||
var column = int.Parse(match.Groups["column"].Value, CultureInfo.InvariantCulture) + 1;
|
||||
if (right == "(" && IsName(left))
|
||||
return $"Argument '{argumentName}' has no operator between '{left}' and '(' at column {column}. The tool knows no function named '{left}', so it would read the name as a variable times the parentheses. Multiplication is written '*'.";
|
||||
|
||||
return $"Argument '{argumentName}' has no operator between '{left}' and '{right}' at column {column}. Multiplication is written '*', as in 2*x, and a function takes its argument in parentheses, as in sin(x).";
|
||||
}
|
||||
|
||||
private static string DescribeParseError(string message, string argumentName)
|
||||
{
|
||||
var match = ParseErrorPosition().Match(message);
|
||||
if (!match.Success)
|
||||
return $"Argument '{argumentName}' cannot be read: {message.Shorten(MAX_REASON_LENGTH)}";
|
||||
|
||||
var column = int.Parse(match.Groups["column"].Value, CultureInfo.InvariantCulture) + 1;
|
||||
var reason = match.Groups["reason"].Value.Shorten(MAX_REASON_LENGTH);
|
||||
return $"Argument '{argumentName}' cannot be read at column {column}: {reason}.";
|
||||
}
|
||||
|
||||
private static bool IsName(string token) => token.Length > 0 && token.All(character => char.IsLetterOrDigit(character) || character == '_');
|
||||
|
||||
/// <summary>
|
||||
/// Counts the equals signs which stand for an equation, leaving out those of <= and >=.
|
||||
/// </summary>
|
||||
private static int CountEqualsSigns(string text)
|
||||
{
|
||||
var count = 0;
|
||||
for (var index = 0; index < text.Length; index++)
|
||||
{
|
||||
if (text[index] != '=')
|
||||
continue;
|
||||
|
||||
if (index > 0 && text[index - 1] is '<' or '>')
|
||||
continue;
|
||||
|
||||
count++;
|
||||
}
|
||||
|
||||
return count;
|
||||
}
|
||||
|
||||
// There should be an operator between [@0,0:0='2',<185>,1:0] and [@1,1:1='x',<189>,1:1], while parsing `2x`
|
||||
[GeneratedRegex(@"='(?<left>.*?)',<\d+>,\d+:\d+\] and \[@\d+,\d+:\d+='(?<right>.*?)',<\d+>,\d+:(?<column>\d+)\]")]
|
||||
private static partial Regex MissingOperatorTokens();
|
||||
|
||||
// line 1:4 extraneous input '*' expecting {'-', '+', ...}, while parsing `2 + * 3`
|
||||
[GeneratedRegex(@"^line \d+:(?<column>\d+) (?<reason>.*?)(?: expecting .*)?(?:, while parsing `.*`)?$", RegexOptions.Singleline)]
|
||||
private static partial Regex ParseErrorPosition();
|
||||
}
|
||||
@@ -0,0 +1,208 @@
|
||||
using System.Globalization;
|
||||
using System.Text;
|
||||
|
||||
using AngouriMath;
|
||||
|
||||
using PeterO.Numbers;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Writes the results of the Mathematics tools for the model: as text, as LaTeX, and as decimals.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Every number is written the same way whatever the language of the user, because it goes to a
|
||||
/// model and comes back in the plain syntax of the tools.
|
||||
/// </remarks>
|
||||
public static class MathFormatting
|
||||
{
|
||||
/// <summary>
|
||||
/// The smallest decimal exponent written without scientific notation.
|
||||
/// </summary>
|
||||
private const int MIN_PLAIN_EXPONENT = -6;
|
||||
|
||||
/// <summary>
|
||||
/// The largest decimal exponent written without scientific notation.
|
||||
/// </summary>
|
||||
private const int MAX_PLAIN_EXPONENT = 20;
|
||||
|
||||
/// <summary>
|
||||
/// An expression in the plain syntax which the tools also read.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <param name="what">What the expression is, named when it is too long.</param>
|
||||
/// <returns>The expression as text.</returns>
|
||||
/// <exception cref="MathComputationException">The text is longer than MathLimits.MAX_RESULT_LENGTH.</exception>
|
||||
public static string ToText(Entity expression, string what) => EnsureLength(expression.Stringize(), what);
|
||||
|
||||
/// <summary>
|
||||
/// An expression as LaTeX, which the chat renders.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <param name="what">What the expression is, named when it is too long.</param>
|
||||
/// <returns>The expression as LaTeX.</returns>
|
||||
/// <exception cref="MathComputationException">The LaTeX is longer than MathLimits.MAX_RESULT_LENGTH.</exception>
|
||||
public static string ToLatex(Entity expression, string what) => EnsureLength(expression.Latexize(), what);
|
||||
|
||||
/// <summary>
|
||||
/// Refuses a text longer than MathLimits.MAX_RESULT_LENGTH.
|
||||
/// </summary>
|
||||
/// <param name="text">The text of a result.</param>
|
||||
/// <param name="what">What the text is, named when it is too long.</param>
|
||||
/// <returns>The text.</returns>
|
||||
/// <exception cref="MathComputationException">The text is too long.</exception>
|
||||
public static string EnsureLength(string text, string what)
|
||||
{
|
||||
if (text.Length > MathLimits.MAX_RESULT_LENGTH)
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"{what} has {text.Length:N0} characters, but at most {MathLimits.MAX_RESULT_LENGTH:N0} are returned."));
|
||||
|
||||
return text;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// A decimal, rounded to the given number of significant digits.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Trailing zeros are left out, so 0.5 stays 0.5 at twenty digits; how many digits the value
|
||||
/// carries is reported next to it. Between 1e-6 and 1e21 the decimal is written out, beyond that
|
||||
/// in scientific notation such as 6.02214076e+23.
|
||||
/// </remarks>
|
||||
/// <param name="value">The value.</param>
|
||||
/// <param name="digits">How many significant digits to keep.</param>
|
||||
/// <returns>The decimal.</returns>
|
||||
public static string FormatDecimal(EDecimal value, int digits)
|
||||
{
|
||||
if (value.IsNaN())
|
||||
return "NaN";
|
||||
|
||||
if (value.IsInfinity())
|
||||
return value.IsNegative ? "-oo" : "+oo";
|
||||
|
||||
var rounded = value.RoundToPrecision(EContext.ForPrecisionAndRounding(digits, ERounding.HalfEven)).Reduce(null);
|
||||
if (rounded.IsZero)
|
||||
return "0";
|
||||
|
||||
var mantissa = rounded.Mantissa.Abs().ToString();
|
||||
var exponent = rounded.Exponent.ToInt64Checked() + mantissa.Length - 1;
|
||||
if (exponent is >= MIN_PLAIN_EXPONENT and <= MAX_PLAIN_EXPONENT)
|
||||
return rounded.ToPlainString();
|
||||
|
||||
var scientific = new StringBuilder();
|
||||
if (rounded.IsNegative)
|
||||
scientific.Append('-');
|
||||
|
||||
scientific.Append(mantissa[0]);
|
||||
if (mantissa.Length > 1)
|
||||
scientific.Append('.').Append(mantissa, 1, mantissa.Length - 1);
|
||||
|
||||
scientific.Append('e').Append(exponent < 0 ? '-' : '+').Append(Math.Abs(exponent).ToString(CultureInfo.InvariantCulture));
|
||||
return scientific.ToString();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// A double, rounded to the given number of significant digits.
|
||||
/// </summary>
|
||||
/// <param name="value">The value.</param>
|
||||
/// <param name="digits">How many significant digits to keep, 17 at most.</param>
|
||||
/// <returns>The decimal.</returns>
|
||||
public static string FormatDecimal(double value, int digits)
|
||||
{
|
||||
if (double.IsNaN(value))
|
||||
return "NaN";
|
||||
|
||||
if (double.IsInfinity(value))
|
||||
return value < 0 ? "-oo" : "+oo";
|
||||
|
||||
return FormatDecimal(EDecimal.FromDouble(value), Math.Min(digits, 17));
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// A real or complex number, rounded to the given number of significant digits, as in 1.5 - 2i.
|
||||
/// </summary>
|
||||
/// <param name="value">The number.</param>
|
||||
/// <param name="digits">How many significant digits to keep in each part.</param>
|
||||
/// <returns>The decimal.</returns>
|
||||
public static string FormatNumber(Entity.Number.Complex value, int digits)
|
||||
{
|
||||
var real = value.RealPart.EDecimal;
|
||||
var imaginary = value.ImaginaryPart.EDecimal;
|
||||
if (imaginary.IsZero)
|
||||
return FormatDecimal(real, digits);
|
||||
|
||||
var imaginaryText = FormatDecimal(imaginary.Abs(), digits) + "i";
|
||||
if (real.IsZero)
|
||||
return imaginary.IsNegative ? "-" + imaginaryText : imaginaryText;
|
||||
|
||||
return $"{FormatDecimal(real, digits)} {(imaginary.IsNegative ? '-' : '+')} {imaginaryText}";
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// How far apart two numbers lie, as the sum of the distances of their parts.
|
||||
/// </summary>
|
||||
/// <param name="first">One number.</param>
|
||||
/// <param name="second">The other number.</param>
|
||||
/// <returns>The distance, as a double.</returns>
|
||||
public static double Distance(Entity.Number.Complex first, Entity.Number.Complex second)
|
||||
{
|
||||
var context = EContext.ForPrecision(20);
|
||||
var real = first.RealPart.EDecimal.Subtract(second.RealPart.EDecimal, context).Abs();
|
||||
var imaginary = first.ImaginaryPart.EDecimal.Subtract(second.ImaginaryPart.EDecimal, context).Abs();
|
||||
return real.Add(imaginary, context).ToDouble();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// How many significant digits two real or complex values of the same quantity have in common,
|
||||
/// measured against the larger of them, so a tiny imaginary part does not count for more than it is.
|
||||
/// </summary>
|
||||
/// <param name="first">One value.</param>
|
||||
/// <param name="second">The other value.</param>
|
||||
/// <param name="maximum">The most digits to report, for values which are equal.</param>
|
||||
/// <returns>The number of digits, between 0 and maximum.</returns>
|
||||
public static int CountMatchingDigits(Entity.Number.Complex first, Entity.Number.Complex second, int maximum)
|
||||
{
|
||||
if (first.ImaginaryPart.EDecimal.IsZero && second.ImaginaryPart.EDecimal.IsZero)
|
||||
return CountMatchingDigits(first.RealPart.EDecimal, second.RealPart.EDecimal, maximum);
|
||||
|
||||
var context = EContext.ForPrecision(20);
|
||||
var firstSize = first.RealPart.EDecimal.Abs().Add(first.ImaginaryPart.EDecimal.Abs(), context);
|
||||
var secondSize = second.RealPart.EDecimal.Abs().Add(second.ImaginaryPart.EDecimal.Abs(), context);
|
||||
var scale = EDecimal.Max(firstSize, secondSize).ToDouble();
|
||||
var distance = Distance(first, second);
|
||||
if (!double.IsFinite(distance) || !double.IsFinite(scale))
|
||||
return 0;
|
||||
|
||||
if (distance <= 0)
|
||||
return maximum;
|
||||
|
||||
return (int)Math.Clamp(Math.Floor(-Math.Log10(distance / scale)), 0, maximum);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// How many significant digits two values of the same quantity have in common.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Counted from their relative difference, as is usual for correct digits: 0.99999 and 1.00001
|
||||
/// agree to four digits although no digit is the same. Evaluated at two working precisions, the
|
||||
/// digits they share are the ones the computation carries.
|
||||
/// </remarks>
|
||||
/// <param name="first">One value.</param>
|
||||
/// <param name="second">The other value.</param>
|
||||
/// <param name="maximum">The most digits to report, for values which are equal.</param>
|
||||
/// <returns>The number of digits, between 0 and maximum.</returns>
|
||||
public static int CountMatchingDigits(EDecimal first, EDecimal second, int maximum)
|
||||
{
|
||||
if (first.IsNaN() || second.IsNaN() || first.IsInfinity() || second.IsInfinity())
|
||||
return 0;
|
||||
|
||||
if (first.CompareTo(second) == 0)
|
||||
return maximum;
|
||||
|
||||
var context = EContext.ForPrecision(20);
|
||||
var scale = first.Abs().CompareTo(second.Abs()) >= 0 ? first.Abs() : second.Abs();
|
||||
var relative = first.Subtract(second, context).Abs().Divide(scale, context).ToDouble();
|
||||
if (relative <= 0)
|
||||
return maximum;
|
||||
|
||||
return (int)Math.Clamp(Math.Floor(-Math.Log10(relative)), 0, maximum);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,160 @@
|
||||
using System.Globalization;
|
||||
|
||||
using AIStudio.Chat;
|
||||
|
||||
using AngouriMath;
|
||||
|
||||
using PeterO.Numbers;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Refuses the computations which AngouriMath cannot stop once they started.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath computes an exact power or factorial in full, and an expansion term by term, without
|
||||
/// checking its cancellation token meanwhile: 9^9^9, (10^9)!, and the expansion of (x + y)^500 run
|
||||
/// on long past any budget, see AngouriMathBehaviorTests. The guards estimate the size of such a
|
||||
/// result from the expression alone, before anything is computed, and refuse it when it would be
|
||||
/// larger than any result the tools return anyway.<br/><br/>
|
||||
/// The estimate works from the inside out, so the inner power of 9^9^9^9 is known before the outer
|
||||
/// one is judged. Only numbers of known value count: a power of a variable is never computed as a
|
||||
/// number. Call the guards inside a computation of MathComputationService, since they evaluate the
|
||||
/// exponents and the arguments of factorials, which have passed the guards before.
|
||||
/// </remarks>
|
||||
public static class MathGuards
|
||||
{
|
||||
private const int MAX_SHOWN_LENGTH = 120;
|
||||
|
||||
/// <summary>
|
||||
/// Refuses an expression whose exact evaluation would produce a number of more than
|
||||
/// MathLimits.MAX_EXACT_DIGITS digits.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression to be evaluated or simplified.</param>
|
||||
/// <exception cref="MathComputationException">A part of the expression would grow too large.</exception>
|
||||
public static void CheckExactSize(Entity expression) => EstimateDigits(expression);
|
||||
|
||||
/// <summary>
|
||||
/// Refuses an expression whose expansion would take a power of a sum beyond
|
||||
/// MathLimits.MAX_EXPANSION_EXPONENT.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression to be expanded.</param>
|
||||
/// <exception cref="MathComputationException">A power of a sum has too large an exponent.</exception>
|
||||
public static void CheckExpansion(Entity expression)
|
||||
{
|
||||
foreach (var node in expression.Nodes)
|
||||
{
|
||||
if (node is not Entity.Powf { Base: Entity.Sumf or Entity.Minusf } power)
|
||||
continue;
|
||||
|
||||
var exponent = ValueOf(power.Exponent);
|
||||
if (exponent is null || Math.Abs(exponent.Value) <= MathLimits.MAX_EXPANSION_EXPONENT)
|
||||
continue;
|
||||
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"'{Show(power)}' raises a sum to the power of {exponent.Value:0.###}. Where the tool may expand, a sum is raised to exponents up to {MathLimits.MAX_EXPANSION_EXPONENT} only, since it cannot stop a larger expansion midway."));
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Estimates how many digits the exact value of a node needs, and refuses it beyond the limit.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// An upper bound for sums, products, and quotients: the digits of numerator and denominator
|
||||
/// add up. A number which is not exact, such as a decimal of the current precision, needs none.
|
||||
/// </remarks>
|
||||
private static double EstimateDigits(Entity node)
|
||||
{
|
||||
var digits = node switch
|
||||
{
|
||||
Entity.Number.Rational rational => DigitsOf(rational.ERational.Numerator) + DigitsOf(rational.ERational.Denominator),
|
||||
Entity.Powf power => EstimatePower(power),
|
||||
Entity.Factorialf factorial => EstimateFactorial(factorial),
|
||||
_ => node.DirectChildren.Sum(EstimateDigits),
|
||||
};
|
||||
|
||||
if (digits > MathLimits.MAX_EXACT_DIGITS)
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"'{Show(node)}' would have about {digits:N0} digits as an exact number. The tool computes such a number in full and cannot stop midway, so exact numbers are limited to {MathLimits.MAX_EXACT_DIGITS:N0} digits."));
|
||||
|
||||
return digits;
|
||||
}
|
||||
|
||||
private static double EstimatePower(Entity.Powf power)
|
||||
{
|
||||
var baseDigits = EstimateDigits(power.Base);
|
||||
var exponentDigits = EstimateDigits(power.Exponent);
|
||||
if (MathVariables.HasAny(power.Base))
|
||||
return baseDigits + exponentDigits;
|
||||
|
||||
var exponent = ValueOf(power.Exponent);
|
||||
if (exponent is null)
|
||||
return baseDigits + exponentDigits;
|
||||
|
||||
var magnitude = Math.Abs(exponent.Value);
|
||||
if (power.Base.Evaled is Entity.Number.Rational rational)
|
||||
{
|
||||
var size = Log10(rational.ERational.Numerator) + Log10(rational.ERational.Denominator);
|
||||
var digits = magnitude * size;
|
||||
if (digits > MathLimits.MAX_EXACT_DIGITS && rational.ERational.Sign > 0)
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"'{Show(power)}' would have about {digits:N0} digits as an exact number; its decimal logarithm is about {exponent.Value * (Log10(rational.ERational.Numerator) - Log10(rational.ERational.Denominator)):0.###}. The tool computes such a number in full and cannot stop midway, so exact numbers are limited to {MathLimits.MAX_EXACT_DIGITS:N0} digits."));
|
||||
|
||||
return digits;
|
||||
}
|
||||
|
||||
// A base which is exact without being rational, such as sqrt(2), can still become one, as 2^500000:
|
||||
return magnitude * Math.Max(baseDigits, 0);
|
||||
}
|
||||
|
||||
private static double EstimateFactorial(Entity.Factorialf factorial)
|
||||
{
|
||||
var argumentDigits = EstimateDigits(factorial.Argument);
|
||||
var argument = ValueOf(factorial.Argument);
|
||||
if (argument is null)
|
||||
return argumentDigits;
|
||||
|
||||
// Stirling: log10(n!) is about n * log10(n / e) + log10(2 * pi * n) / 2.
|
||||
var n = Math.Abs(argument.Value);
|
||||
if (n < 2)
|
||||
return argumentDigits + 1;
|
||||
|
||||
return n * Math.Log10(n / Math.E) + Math.Log10(2 * Math.PI * n) / 2 + 1;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The value of a part without variables, as a double, or null when it has none.
|
||||
/// </summary>
|
||||
private static double? ValueOf(Entity part)
|
||||
{
|
||||
if (MathVariables.HasAny(part))
|
||||
return null;
|
||||
|
||||
var value = part.Evaled switch
|
||||
{
|
||||
Entity.Number.Real real => real.EDecimal.ToDouble(),
|
||||
Entity.Number.Complex complex => Math.Abs(complex.RealPart.EDecimal.ToDouble()) + Math.Abs(complex.ImaginaryPart.EDecimal.ToDouble()),
|
||||
_ => double.NaN,
|
||||
};
|
||||
|
||||
return double.IsNaN(value) ? null : value;
|
||||
}
|
||||
|
||||
private static double DigitsOf(EInteger value) => value.IsZero ? 1 : value.Abs().GetDigitCountAsInt64();
|
||||
|
||||
/// <summary>
|
||||
/// The decimal logarithm of an integer's magnitude, also beyond the range of a double.
|
||||
/// </summary>
|
||||
private static double Log10(EInteger value)
|
||||
{
|
||||
var magnitude = value.Abs();
|
||||
if (magnitude.IsZero)
|
||||
return 0;
|
||||
|
||||
var digits = magnitude.GetDigitCountAsInt64();
|
||||
if (digits <= 15)
|
||||
return Math.Log10(magnitude.ToInt64Checked());
|
||||
|
||||
var leading = double.Parse(magnitude.ToString()[..15], CultureInfo.InvariantCulture);
|
||||
return Math.Log10(leading) + (digits - 15);
|
||||
}
|
||||
|
||||
private static string Show(Entity node) => node.ToString().Shorten(MAX_SHOWN_LENGTH);
|
||||
}
|
||||
@@ -0,0 +1,49 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// The limits of the Mathematics tools, in one place.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AI Studio also runs on a Raspberry Pi, and AngouriMath cannot stop every computation midway. The
|
||||
/// limits keep a single call from taking the device: they bound what is read, what is computed
|
||||
/// exactly, and what goes back to the model. "Mathematics" in documentation/Tools.md lists them.
|
||||
/// </remarks>
|
||||
public static class MathLimits
|
||||
{
|
||||
/// <summary>
|
||||
/// The longest expression read, in characters.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// It also bounds how deeply an expression can nest: 1,000 levels of parentheses at most, which
|
||||
/// the stack of a computation holds, see MathComputationService.
|
||||
/// </remarks>
|
||||
public const int MAX_EXPRESSION_LENGTH = 2_000;
|
||||
|
||||
/// <summary>
|
||||
/// The longest result returned, in characters, for each form of it.
|
||||
/// </summary>
|
||||
public const int MAX_RESULT_LENGTH = 20_000;
|
||||
|
||||
/// <summary>
|
||||
/// The most digits an exact number may have while it is computed.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath computes an exact power or factorial in full and checks no cancellation meanwhile,
|
||||
/// so MathGuards refuses one which would grow beyond this, before it starts.
|
||||
/// </remarks>
|
||||
public const int MAX_EXACT_DIGITS = 20_000;
|
||||
|
||||
/// <summary>
|
||||
/// The largest exponent of a sum which may be expanded.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// (x + y)^200 expands in a few milliseconds into about 13,000 characters, while (x + y)^500 runs
|
||||
/// on without checking for cancellation.
|
||||
/// </remarks>
|
||||
public const int MAX_EXPANSION_EXPONENT = 200;
|
||||
|
||||
/// <summary>
|
||||
/// The most significant digits a decimal result may have.
|
||||
/// </summary>
|
||||
public const int MAX_DIGITS = 100;
|
||||
}
|
||||
@@ -0,0 +1,32 @@
|
||||
using AngouriMath;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// The variables of an expression, as the Mathematics tools count them.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath counts the constants pi and e among the free variables, so ln(0) has the free
|
||||
/// variable e, the base of its logarithm. Asked for the variables of an expression, the tools mean
|
||||
/// those which have no value, see AngouriMathBehaviorTests.
|
||||
/// </remarks>
|
||||
public static class MathVariables
|
||||
{
|
||||
/// <summary>
|
||||
/// The free variables of an expression without the constants, in the order of their names.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <returns>Each variable once.</returns>
|
||||
public static IReadOnlyList<Entity.Variable> FreeIn(Entity expression) => expression.FreeVariables
|
||||
.Where(variable => variable is not Entity.Constant)
|
||||
.Distinct()
|
||||
.OrderBy(variable => variable.Name, StringComparer.Ordinal)
|
||||
.ToList();
|
||||
|
||||
/// <summary>
|
||||
/// Whether an expression has a free variable other than the constants.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <returns>True when a variable without a value occurs freely.</returns>
|
||||
public static bool HasAny(Entity expression) => expression.FreeVariables.Any(variable => variable is not Entity.Constant);
|
||||
}
|
||||
@@ -0,0 +1,212 @@
|
||||
using AngouriMath;
|
||||
|
||||
using PeterO.Numbers;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Checks numerically whether two expressions are equal, at fixed sample points.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// A symbolic answer is only as good as the library which gave it, and AngouriMath changes answers
|
||||
/// between releases. The tools therefore check what they can: a simplification against its input, an
|
||||
/// antiderivative by its derivative, a solution by inserting it.<br/><br/>
|
||||
/// The sample points are exact rationals between 0.6 and 2.3, so a polynomial is evaluated exactly,
|
||||
/// however much its terms cancel, and everything else at 50 digits. They are fixed, so the same
|
||||
/// question gets the same answer every time; the providers cache a conversation from its beginning.
|
||||
/// A point where either side has no value, such as a pole, is skipped. Agreement at the points
|
||||
/// supports the equality; it does not prove it.<br/><br/>
|
||||
/// A derivative is compared with the central difference of its function, over a step of 1e-20 at 60
|
||||
/// digits. The integer parameter of periodic solutions, which AngouriMath names n_1, takes the values
|
||||
/// -1, 0, and 1 instead, since such a solution holds at integers only.<br/><br/>
|
||||
/// Call it inside a computation of MathComputationService.
|
||||
/// </remarks>
|
||||
public static class NumericCheck
|
||||
{
|
||||
private const int SAMPLE_POINTS = 6;
|
||||
|
||||
private const int MIN_VALID_POINTS = 3;
|
||||
|
||||
private const int PRECISION = 50;
|
||||
|
||||
/// <summary>
|
||||
/// The values the variables take, in hundredths: irregular, so that no point lies on a pattern
|
||||
/// such as a multiple of pi, and eleven of them, a prime, so that the variables of a point differ.
|
||||
/// </summary>
|
||||
private static readonly int[] SAMPLE_HUNDREDTHS = [61, 73, 89, 107, 119, 131, 147, 163, 179, 197, 223];
|
||||
|
||||
private static readonly EDecimal TOLERANCE = EDecimal.FromString("1e-30");
|
||||
|
||||
/// <summary>
|
||||
/// The precision of a derivative check: the central difference divides by its tiny step.
|
||||
/// </summary>
|
||||
private const int DERIVATIVE_PRECISION = 60;
|
||||
|
||||
/// <summary>
|
||||
/// The step of the central difference: rounding at 60 digits costs about 1e-40 of it, and its own
|
||||
/// error, about the square of the step, 1e-40 as well.
|
||||
/// </summary>
|
||||
private static readonly Entity.Number.Rational DERIVATIVE_STEP = Entity.Number.Rational.Create(EInteger.One, EInteger.FromInt32(10).Pow(20));
|
||||
|
||||
private static readonly Entity.Number.Rational TWO_DERIVATIVE_STEPS = Entity.Number.Rational.Create(EInteger.FromInt32(2), EInteger.FromInt32(10).Pow(20));
|
||||
|
||||
private static readonly EDecimal DERIVATIVE_TOLERANCE = EDecimal.FromString("1e-20");
|
||||
|
||||
/// <summary>
|
||||
/// Checks whether two expressions agree at every sample point where both have a value.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// A condition of AngouriMath, as in x + 1 provided not x - 1 = 0, is left aside on both sides:
|
||||
/// the expressions are compared where they have a value.
|
||||
/// </remarks>
|
||||
/// <param name="left">One expression.</param>
|
||||
/// <param name="right">The other expression.</param>
|
||||
/// <returns>What the check found.</returns>
|
||||
public static NumericCheckOutcome AreEqual(Entity left, Entity right)
|
||||
{
|
||||
left = WithoutCondition(left);
|
||||
right = WithoutCondition(right);
|
||||
var variables = MathVariables.FreeIn(left).Concat(MathVariables.FreeIn(right)).Distinct().OrderBy(variable => variable.Name, StringComparer.Ordinal).ToArray();
|
||||
|
||||
using var precision = MathS.Settings.DecimalPrecisionContext.Set(EContext.ForPrecision(PRECISION));
|
||||
var validPoints = 0;
|
||||
for (var point = 0; point < SAMPLE_POINTS; point++)
|
||||
{
|
||||
var assignment = SampleAssignment(point, variables);
|
||||
var leftValue = ValueAt(left, assignment);
|
||||
var rightValue = ValueAt(right, assignment);
|
||||
if (leftValue is null || rightValue is null)
|
||||
continue;
|
||||
|
||||
validPoints++;
|
||||
if (!AreClose(leftValue, rightValue, TOLERANCE))
|
||||
return NumericCheckOutcome.FAILED;
|
||||
}
|
||||
|
||||
return validPoints >= MIN_VALID_POINTS ? NumericCheckOutcome.PASSED : NumericCheckOutcome.NOT_POSSIBLE;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Checks whether one expression is the derivative of another, by a central difference at each sample point.
|
||||
/// </summary>
|
||||
/// <param name="derivative">The claimed derivative.</param>
|
||||
/// <param name="function">The function.</param>
|
||||
/// <param name="variable">The variable of the derivative.</param>
|
||||
/// <returns>What the check found.</returns>
|
||||
public static NumericCheckOutcome IsDerivativeOf(Entity derivative, Entity function, Entity.Variable variable)
|
||||
{
|
||||
derivative = WithoutCondition(derivative);
|
||||
function = WithoutCondition(function);
|
||||
var variables = MathVariables.FreeIn(derivative).Concat(MathVariables.FreeIn(function)).Append(variable).Distinct().OrderBy(free => free.Name, StringComparer.Ordinal).ToArray();
|
||||
|
||||
using var precision = MathS.Settings.DecimalPrecisionContext.Set(EContext.ForPrecision(DERIVATIVE_PRECISION));
|
||||
var validPoints = 0;
|
||||
for (var point = 0; point < SAMPLE_POINTS; point++)
|
||||
{
|
||||
var assignment = SampleAssignment(point, variables);
|
||||
var at = assignment[variable];
|
||||
var claimed = ValueAt(derivative, assignment);
|
||||
var above = ValueAt(function, With(assignment, variable, (at + DERIVATIVE_STEP).InnerSimplified));
|
||||
var below = ValueAt(function, With(assignment, variable, (at - DERIVATIVE_STEP).InnerSimplified));
|
||||
if (claimed is null || above is null || below is null)
|
||||
continue;
|
||||
|
||||
var difference = NumericEvaluator.EvaluateAt(((Entity)above - below) / TWO_DERIVATIVE_STEPS, DERIVATIVE_PRECISION);
|
||||
validPoints++;
|
||||
if (!AreClose(claimed, difference, DERIVATIVE_TOLERANCE))
|
||||
return NumericCheckOutcome.FAILED;
|
||||
}
|
||||
|
||||
return validPoints >= MIN_VALID_POINTS ? NumericCheckOutcome.PASSED : NumericCheckOutcome.NOT_POSSIBLE;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Whether AngouriMath named a variable as the integer parameter of periodic solutions, as in pi / 6 + 2 * pi * n_1.
|
||||
/// </summary>
|
||||
/// <param name="variable">The variable.</param>
|
||||
/// <returns>True for n_1, n_2, and so on.</returns>
|
||||
public static bool IsIntegerParameter(Entity.Variable variable) => variable.Name.StartsWith("n_", StringComparison.Ordinal) && variable.Name.Length > 2 && variable.Name[2..].All(char.IsAsciiDigit);
|
||||
|
||||
/// <summary>
|
||||
/// The values of the variables at one sample point.
|
||||
/// </summary>
|
||||
/// <param name="point">The number of the point, from zero.</param>
|
||||
/// <param name="variables">The variables, in a fixed order.</param>
|
||||
/// <returns>Each variable with its exact value.</returns>
|
||||
public static IReadOnlyDictionary<Entity.Variable, Entity> SampleAssignment(int point, IReadOnlyList<Entity.Variable> variables)
|
||||
{
|
||||
var assignment = new Dictionary<Entity.Variable, Entity>(variables.Count);
|
||||
for (var index = 0; index < variables.Count; index++)
|
||||
{
|
||||
// An integer parameter takes -1, 0, and 1, since a periodic solution holds only at integers:
|
||||
if (IsIntegerParameter(variables[index]))
|
||||
{
|
||||
assignment[variables[index]] = Entity.Number.Integer.Create((point + index) % 3 - 1);
|
||||
continue;
|
||||
}
|
||||
|
||||
var hundredths = SAMPLE_HUNDREDTHS[(point * 3 + index * 5) % SAMPLE_HUNDREDTHS.Length];
|
||||
assignment[variables[index]] = Entity.Number.Rational.Create(hundredths, 100);
|
||||
}
|
||||
|
||||
return assignment;
|
||||
}
|
||||
|
||||
private static Dictionary<Entity.Variable, Entity> With(IReadOnlyDictionary<Entity.Variable, Entity> assignment, Entity.Variable variable, Entity value) => new(assignment) { [variable] = value };
|
||||
|
||||
/// <summary>
|
||||
/// The value of an expression with its variables replaced, or null when it has none there.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <param name="assignment">The value of each variable.</param>
|
||||
/// <returns>The value, or null for no value, an infinite one, or one AngouriMath cannot evaluate.</returns>
|
||||
public static Entity.Number.Complex? ValueAt(Entity expression, IReadOnlyDictionary<Entity.Variable, Entity> assignment)
|
||||
{
|
||||
try
|
||||
{
|
||||
var substituted = expression;
|
||||
foreach (var (variable, value) in assignment)
|
||||
substituted = substituted.Substitute(variable, value);
|
||||
|
||||
if (!substituted.EvaluableNumerical)
|
||||
return null;
|
||||
|
||||
//
|
||||
// A polynomial at a sample point stays exact. Anything else is evaluated without the downcasting which turns
|
||||
// tiny values into 0, and so is an exact 0: simplifying takes a tiny value such as 6.62607015e-34 * sin(0.61)
|
||||
// for 0 at this precision, see AngouriMathBehaviorTests:
|
||||
//
|
||||
var exact = NumericEvaluator.ExactValueOf(substituted);
|
||||
var number = exact is not null && !IsZero(exact) ? exact : NumericEvaluator.EvaluateAt(substituted, MathS.Settings.DecimalPrecisionContext.Value.Precision.ToInt32Checked());
|
||||
return IsFinite(number) ? number : null;
|
||||
}
|
||||
catch (Exception exception) when (exception is not OperationCanceledException)
|
||||
{
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
private static Entity WithoutCondition(Entity expression) => expression is Entity.Providedf provided ? provided.Expression : expression;
|
||||
|
||||
private static bool IsFinite(Entity.Number.Complex number) => IsFinite(number.RealPart.EDecimal) && IsFinite(number.ImaginaryPart.EDecimal);
|
||||
|
||||
private static bool IsZero(Entity.Number.Complex number) => number.RealPart.EDecimal.IsZero && number.ImaginaryPart.EDecimal.IsZero;
|
||||
|
||||
private static bool IsFinite(EDecimal value) => !value.IsNaN() && !value.IsInfinity();
|
||||
|
||||
private static bool AreClose(Entity.Number.Complex left, Entity.Number.Complex right, EDecimal tolerance)
|
||||
{
|
||||
var context = EContext.ForPrecision(PRECISION);
|
||||
var difference = left.RealPart.EDecimal.Subtract(right.RealPart.EDecimal, context).Abs().Add(left.ImaginaryPart.EDecimal.Subtract(right.ImaginaryPart.EDecimal, context).Abs(), context);
|
||||
var leftSize = left.RealPart.EDecimal.Abs().Add(left.ImaginaryPart.EDecimal.Abs(), context);
|
||||
var rightSize = right.RealPart.EDecimal.Abs().Add(right.ImaginaryPart.EDecimal.Abs(), context);
|
||||
|
||||
//
|
||||
// Against 0, the tolerance is absolute: an expression whose terms cancel, such as sin(2*x) - 2*sin(x)*cos(x),
|
||||
// comes only as close as its rounding noise. Between two other values, it is relative, so that tiny values,
|
||||
// such as those of a formula with the Planck constant, are still told apart:
|
||||
//
|
||||
var scale = leftSize.IsZero || rightSize.IsZero ? EDecimal.One : EDecimal.Max(leftSize, rightSize);
|
||||
return difference.CompareTo(tolerance.Multiply(scale, context)) <= 0;
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,22 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// What a numeric check found.
|
||||
/// </summary>
|
||||
public enum NumericCheckOutcome
|
||||
{
|
||||
/// <summary>
|
||||
/// Too few sample points had a value on both sides, or one side cannot be evaluated as a number.
|
||||
/// </summary>
|
||||
NOT_POSSIBLE = 0,
|
||||
|
||||
/// <summary>
|
||||
/// Both sides agree at every sample point where both have a value.
|
||||
/// </summary>
|
||||
PASSED,
|
||||
|
||||
/// <summary>
|
||||
/// The sides differ at a sample point.
|
||||
/// </summary>
|
||||
FAILED,
|
||||
}
|
||||
@@ -0,0 +1,327 @@
|
||||
using System.Globalization;
|
||||
|
||||
using AIStudio.Chat;
|
||||
|
||||
using AngouriMath;
|
||||
|
||||
using PeterO.Numbers;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// Evaluates an expression without variables to a number, exactly where it can, and states how many
|
||||
/// of its digits are certain.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// A rational value stays exact, so 547 / 57576 comes back as that fraction. Every other value is
|
||||
/// evaluated twice, at ten and at twenty digits beyond those asked for, and only the digits both agree
|
||||
/// on are reported: AngouriMath rounds at every step, and a compound expression such as pi^4/15 loses
|
||||
/// its last digits at the precision asked for, see AngouriMathBehaviorTests.<br/><br/>
|
||||
/// A definite integral is integrated numerically, never by AngouriMath: its antiderivative may jump
|
||||
/// between the bounds, and Evaled turns the constants inside an integral it cannot take into decimals.
|
||||
/// Its value then carries the digits of a double at most. The two evaluations shift every integral
|
||||
/// by its estimated error, down in one and up in the other, so the digits they agree on are those the
|
||||
/// integrals vouch for, too.<br/><br/>
|
||||
/// A tiny value, such as one computed with the Planck constant, gets as many more digits as it has
|
||||
/// leading zeros, see EvaluateAt.<br/><br/>
|
||||
/// Call it inside a computation of MathComputationService.
|
||||
/// </remarks>
|
||||
public static class NumericEvaluator
|
||||
{
|
||||
private const int GUARD_DIGITS = 10;
|
||||
|
||||
private const int MAX_SHOWN_LENGTH = 120;
|
||||
|
||||
/// <summary>
|
||||
/// The most digits an evaluation computes with. A value which needs more, because it is smaller
|
||||
/// than about 1e-900, loses its digits.
|
||||
/// </summary>
|
||||
private const int MAX_PRECISION = 1_000;
|
||||
|
||||
/// <summary>
|
||||
/// The precision of a bound of an integral, which ends up as a double.
|
||||
/// </summary>
|
||||
private const int BOUND_PRECISION = 30;
|
||||
|
||||
/// <summary>
|
||||
/// Evaluates an expression to a number.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression, read by MathExpressionReader.</param>
|
||||
/// <param name="digits">How many significant digits are asked for.</param>
|
||||
/// <param name="token">The cancellation token of the computation.</param>
|
||||
/// <returns>The value with its certain digits.</returns>
|
||||
/// <exception cref="MathComputationException">The expression has no numeric value, or cannot be integrated.</exception>
|
||||
public static NumericValue Evaluate(Entity expression, int digits, CancellationToken token)
|
||||
{
|
||||
RefuseWhatIsNoNumber(expression);
|
||||
MathGuards.CheckExactSize(expression);
|
||||
|
||||
var integrals = new List<NumericIntegral>();
|
||||
var issues = new List<string>();
|
||||
var withPlaceholders = expression.Replace(node => node is Entity.Integralf integral ? Integrate(integral, integrals, issues, token) : node);
|
||||
|
||||
var exact = integrals.Count == 0 ? ExactValueOf(expression) : null;
|
||||
if (exact is not null && IsZero(exact) && MayBeTiny(expression, digits))
|
||||
exact = null;
|
||||
|
||||
if (exact is not null)
|
||||
return new(exact, DecimalOf(exact, digits + GUARD_DIGITS), digits, false, null, issues);
|
||||
|
||||
var substituted = Substituted(withPlaceholders, integrals, 0);
|
||||
if (!substituted.EvaluableNumerical)
|
||||
throw new MathComputationException($"'{Show(expression)}' has no value the tool can compute as a single number, as for a sum without a closed form, an infinite sum, or a matrix.");
|
||||
|
||||
var central = EvaluateAt(substituted, digits + GUARD_DIGITS);
|
||||
if (IsNaN(central))
|
||||
throw new MathComputationException($"'{Show(expression)}' has no value: it evaluates to NaN, as for a division by zero, a pole, or the factorial of a negative integer.");
|
||||
|
||||
if (IsInfinite(central))
|
||||
{
|
||||
issues.Add("The value is infinite.");
|
||||
return new(null, central, 0, integrals.Count > 0, null, issues);
|
||||
}
|
||||
|
||||
//
|
||||
// Where AngouriMath rounds to decimal places, a value far below 1 keeps only some of its digits, and a
|
||||
// tiny one none at all. The precision therefore grows by the leading zeros of the value, and a value of
|
||||
// zero is evaluated again with twice the digits, until MAX_PRECISION. More digits cannot change an
|
||||
// integral's value of zero, which comes from the quadrature:
|
||||
//
|
||||
var shift = 0;
|
||||
var maxShift = MAX_PRECISION - digits - 2 * GUARD_DIGITS;
|
||||
while (true)
|
||||
{
|
||||
int wantedShift;
|
||||
if (!IsZero(central))
|
||||
wantedShift = LeadingZerosOf(central);
|
||||
else if (integrals.Count == 0)
|
||||
wantedShift = 2 * shift + digits + GUARD_DIGITS;
|
||||
else
|
||||
break;
|
||||
|
||||
wantedShift = Math.Min(wantedShift, maxShift);
|
||||
if (wantedShift <= shift + GUARD_DIGITS)
|
||||
break;
|
||||
|
||||
token.ThrowIfCancellationRequested();
|
||||
shift = wantedShift;
|
||||
central = EvaluateAt(substituted, digits + GUARD_DIGITS + shift);
|
||||
}
|
||||
|
||||
var low = integrals.Count == 0 ? central : EvaluateAt(Substituted(withPlaceholders, integrals, -1), digits + GUARD_DIGITS + shift);
|
||||
var high = EvaluateAt(Substituted(withPlaceholders, integrals, 1), digits + 2 * GUARD_DIGITS + shift);
|
||||
var certainDigits = MathFormatting.CountMatchingDigits(low, high, digits);
|
||||
double? estimatedError = integrals.Count == 0 ? null : MathFormatting.Distance(low, high) / 2;
|
||||
if (IsZero(central) && integrals.Count == 0)
|
||||
{
|
||||
certainDigits = 0;
|
||||
issues.Add(string.Create(CultureInfo.InvariantCulture, $"The value evaluates to 0 with {digits + GUARD_DIGITS + shift} digits, but the tool cannot show that it is exactly 0: a value too small for that many digits comes out as 0, too."));
|
||||
}
|
||||
else if (certainDigits == 0)
|
||||
issues.Add(string.Create(CultureInfo.InvariantCulture, $"No digit of the value is certain: it lies within about {MathFormatting.FormatDecimal(MathFormatting.Distance(low, high), 2)} of the decimal."));
|
||||
else if (certainDigits < digits)
|
||||
issues.Add(string.Create(CultureInfo.InvariantCulture, $"Only {certainDigits} of the {digits} significant digits asked for are certain, so the decimal has {certainDigits}."));
|
||||
|
||||
return new(null, central, certainDigits, integrals.Count > 0, estimatedError, issues);
|
||||
}
|
||||
|
||||
private static void RefuseWhatIsNoNumber(Entity expression)
|
||||
{
|
||||
if (expression is Entity.Statement)
|
||||
throw new MathComputationException($"'{Show(expression)}' is a statement, such as a comparison, and has no numeric value. An equation with one '=' is checked by comparing its sides.");
|
||||
|
||||
if (expression is Entity.Set)
|
||||
throw new MathComputationException($"'{Show(expression)}' is a set, and has no numeric value.");
|
||||
|
||||
var variables = MathVariables.FreeIn(expression).Select(variable => variable.Name).ToList();
|
||||
if (variables.Count > 0)
|
||||
throw new MathComputationException($"'{Show(expression)}' contains the variable{(variables.Count == 1 ? string.Empty : "s")} {string.Join(", ", variables)}, so it has no numeric value. compute_symbolically works with variables.");
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The exact value of an expression when it is rational, or complex with rational parts.
|
||||
/// </summary>
|
||||
/// <param name="expression">The expression.</param>
|
||||
/// <returns>The exact value, or null for an expression without one.</returns>
|
||||
public static Entity.Number.Complex? ExactValueOf(Entity expression) => expression.InnerSimplified switch
|
||||
{
|
||||
Entity.Number.Rational rational => rational,
|
||||
Entity.Number.Complex { RealPart: Entity.Number.Rational, ImaginaryPart: Entity.Number.Rational } complex => complex,
|
||||
_ => null,
|
||||
};
|
||||
|
||||
/// <summary>
|
||||
/// Integrates a definite integral numerically, and puts a variable of its own in its place.
|
||||
/// </summary>
|
||||
private static Entity Integrate(Entity.Integralf integral, List<NumericIntegral> integrals, List<string> issues, CancellationToken token)
|
||||
{
|
||||
if (integral.Range is not { } range)
|
||||
throw new MathComputationException($"'{Show(integral)}' has no bounds, so it has no numeric value. compute_symbolically integrates without bounds.");
|
||||
|
||||
var (from, to) = range;
|
||||
|
||||
if (integral.Var is not Entity.Variable variable)
|
||||
throw new MathComputationException($"'{Show(integral)}' does not integrate over a variable.");
|
||||
|
||||
if (MathVariables.FreeIn(integral.Expression).Any(free => free != variable) || MathVariables.HasAny(from) || MathVariables.HasAny(to))
|
||||
throw new MathComputationException($"'{Show(integral)}' depends on the variable of an outer integral. Nested integrals are not integrated numerically.");
|
||||
|
||||
var lower = BoundOf(from, integral);
|
||||
var upper = BoundOf(to, integral);
|
||||
Func<double, double> integrand;
|
||||
try
|
||||
{
|
||||
integrand = integral.Expression.Compile<double, double>(variable);
|
||||
}
|
||||
catch (Exception exception) when (exception is not OperationCanceledException)
|
||||
{
|
||||
throw new MathComputationException($"The integrand of '{Show(integral)}' cannot be evaluated numerically: {exception.Message.Shorten(MAX_SHOWN_LENGTH)}");
|
||||
}
|
||||
|
||||
var result = DoubleExponentialQuadrature.Integrate(integrand, lower, upper, token);
|
||||
switch (result.Status)
|
||||
{
|
||||
case QuadratureStatus.NO_VALUE:
|
||||
throw new MathComputationException(string.Create(CultureInfo.InvariantCulture, $"The integrand of '{Show(integral)}' has no finite value at {variable.Name} = {MathFormatting.FormatDecimal(result.FailurePoint ?? double.NaN, 6)}, inside the interval of integration."));
|
||||
|
||||
case QuadratureStatus.NOT_CONVERGED:
|
||||
throw new MathComputationException($"The numerical integration of '{Show(integral)}' did not settle. The integral may diverge, converge only conditionally such as sin(x)/x over [0, +oo), or have a singularity inside the interval, which a split at the singularity moves to a bound.");
|
||||
}
|
||||
|
||||
// The expression has no free variables at this point, so no name can clash with one of them:
|
||||
var placeholder = MathS.Var($"Q_{integrals.Count + 1}");
|
||||
integrals.Add(new(placeholder, result.Value, Math.Max(result.EstimatedError, Math.Abs(result.Value) * Math.ScaleB(1.0, -52))));
|
||||
issues.Add(string.Create(CultureInfo.InvariantCulture, $"'{Show(integral)}' was integrated numerically, with an estimated error of {MathFormatting.FormatDecimal(result.EstimatedError, 2)}."));
|
||||
return placeholder;
|
||||
}
|
||||
|
||||
private static double BoundOf(Entity bound, Entity.Integralf integral)
|
||||
{
|
||||
var value = EvaluateAt(bound, BOUND_PRECISION);
|
||||
if (!value.ImaginaryPart.EDecimal.IsZero)
|
||||
throw new MathComputationException($"A bound of '{Show(integral)}' is not a real number.");
|
||||
|
||||
var number = value.RealPart.EDecimal;
|
||||
if (number.IsNaN())
|
||||
throw new MathComputationException($"A bound of '{Show(integral)}' has no value.");
|
||||
|
||||
if (number.IsInfinity())
|
||||
return number.IsNegative ? double.NegativeInfinity : double.PositiveInfinity;
|
||||
|
||||
return number.ToDouble();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// The expression with every integral replaced by its value, shifted by its estimated error in the given direction.
|
||||
/// </summary>
|
||||
private static Entity Substituted(Entity expression, IReadOnlyList<NumericIntegral> integrals, int direction)
|
||||
{
|
||||
// Downcasting would turn the value of a tiny integral into 0, see EvaluateAt:
|
||||
using var withoutDowncasting = MathS.Settings.DowncastingEnabled.Set(false);
|
||||
foreach (var integral in integrals)
|
||||
expression = expression.Substitute(integral.Placeholder, Entity.Number.Real.Create(EDecimal.FromDouble(integral.Value + direction * integral.Error)));
|
||||
|
||||
return expression;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Evaluates an expression without variables to a number, with at least the given number of significant digits.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// By default, AngouriMath downcasts: it turns a number which looks like an integer at the precision into that
|
||||
/// integer, and so a tiny value such as 1.380649e-23 * 300 into 0, see AngouriMathBehaviorTests. This evaluation
|
||||
/// does without. Then AngouriMath rounds the quotient of two integers to as many decimal places as the precision
|
||||
/// has digits, so 21/10^50 vanishes at 30 digits. The precision therefore grows by the leading zeros of the
|
||||
/// smallest fraction in the expression, up to MAX_PRECISION.
|
||||
/// </remarks>
|
||||
/// <param name="expression">The expression, without variables.</param>
|
||||
/// <param name="precision">How many significant digits to compute with, at least.</param>
|
||||
/// <returns>The value, which may be NaN or infinite.</returns>
|
||||
public static Entity.Number.Complex EvaluateAt(Entity expression, int precision)
|
||||
{
|
||||
precision = Math.Min(precision + LeadingZerosOfFractions(expression), MAX_PRECISION);
|
||||
using var context = MathS.Settings.DecimalPrecisionContext.Set(EContext.ForPrecision(precision));
|
||||
using var withoutDowncasting = MathS.Settings.DowncastingEnabled.Set(false);
|
||||
return expression.EvalNumerical();
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Whether an exact 0 may stand for a tiny value instead.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath simplifies a value below about 1e-1000 to 0, e.g., pi/10^1500, see AngouriMathBehaviorTests.
|
||||
/// A fraction that small, such as 1e-2000, leaves the 0 in doubt. Otherwise, a true 0, as for sin(pi),
|
||||
/// evaluates to rounding noise which changes with the precision, but a tiny value, as for exp(-3000), keeps
|
||||
/// its digits. A quotient such as pi/10^1500 evaluates to 0 as well, since AngouriMath rounds it to decimal
|
||||
/// places; its 0 stands.
|
||||
/// </remarks>
|
||||
private static bool MayBeTiny(Entity expression, int digits)
|
||||
{
|
||||
if (LeadingZerosOfFractions(expression) > MAX_PRECISION - digits - 2 * GUARD_DIGITS)
|
||||
return true;
|
||||
|
||||
var value = EvaluateAt(expression, digits + GUARD_DIGITS);
|
||||
var check = EvaluateAt(expression, digits + 2 * GUARD_DIGITS);
|
||||
return !IsZero(value) && MathFormatting.CountMatchingDigits(value, check, GUARD_DIGITS) >= GUARD_DIGITS / 2;
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// How many zeros the smallest fraction of an expression has after its decimal point, e.g., 49 for 21/10^50.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// The fractions are those the expression writes, such as 1e-200, and those its simplification computes, such as
|
||||
/// 21/10^50 from 3/10^50 * 7.
|
||||
/// </remarks>
|
||||
private static int LeadingZerosOfFractions(Entity expression) => expression.Nodes
|
||||
.Concat(expression.InnerSimplified.Nodes)
|
||||
.OfType<Entity.Number.Complex>()
|
||||
.SelectMany(number => new[] { number.RealPart, number.ImaginaryPart })
|
||||
.OfType<Entity.Number.Rational>()
|
||||
.Select(fraction => fraction.ERational.Denominator.ToString().Length - fraction.ERational.Numerator.Abs().ToString().Length)
|
||||
.Append(0)
|
||||
.Max();
|
||||
|
||||
/// <summary>
|
||||
/// How many zeros a number has after its decimal point, e.g., 18 for 3.97e-19, and 0 for a number of at least 0.1.
|
||||
/// </summary>
|
||||
private static int LeadingZerosOf(Entity.Number.Complex number)
|
||||
{
|
||||
var real = number.RealPart.EDecimal;
|
||||
var imaginary = number.ImaginaryPart.EDecimal;
|
||||
var largest = real.Abs().CompareTo(imaginary.Abs()) >= 0 ? real : imaginary;
|
||||
if (largest.IsZero || largest.IsNaN() || largest.IsInfinity())
|
||||
return 0;
|
||||
|
||||
var exponentOfFirstDigit = largest.Exponent.ToInt64Checked() + largest.UnsignedMantissa.ToString().Length - 1;
|
||||
return (int)Math.Clamp(-exponentOfFirstDigit - 1, 0, MAX_PRECISION);
|
||||
}
|
||||
|
||||
private static bool IsZero(Entity.Number.Complex number) => number.RealPart.EDecimal.IsZero && number.ImaginaryPart.EDecimal.IsZero;
|
||||
|
||||
/// <summary>
|
||||
/// An exact value as a decimal with the given number of significant digits.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// AngouriMath converts a fraction to a fixed number of decimal places, so 1/10^200 would become 0,
|
||||
/// see AngouriMathBehaviorTests.
|
||||
/// </remarks>
|
||||
private static Entity.Number.Complex DecimalOf(Entity.Number.Complex exact, int precision)
|
||||
{
|
||||
var context = EContext.ForPrecision(precision);
|
||||
using var withoutDowncasting = MathS.Settings.DowncastingEnabled.Set(false);
|
||||
return Entity.Number.Complex.Create(Entity.Number.Real.Create(DecimalOf(exact.RealPart, context)), Entity.Number.Real.Create(DecimalOf(exact.ImaginaryPart, context)));
|
||||
}
|
||||
|
||||
private static EDecimal DecimalOf(Entity.Number.Real part, EContext context) => part is Entity.Number.Rational fraction ? fraction.ERational.ToEDecimal(context) : part.EDecimal;
|
||||
|
||||
private static bool IsNaN(Entity.Number.Complex number) => number.RealPart.EDecimal.IsNaN() || number.ImaginaryPart.EDecimal.IsNaN();
|
||||
|
||||
private static bool IsInfinite(Entity.Number.Complex number) => number.RealPart.EDecimal.IsInfinity() || number.ImaginaryPart.EDecimal.IsInfinity();
|
||||
|
||||
private static string Show(Entity expression) => expression.ToString().Shorten(MAX_SHOWN_LENGTH);
|
||||
|
||||
/// <summary>
|
||||
/// A definite integral, integrated numerically, and the variable which stands for it.
|
||||
/// </summary>
|
||||
private sealed record NumericIntegral(Entity.Variable Placeholder, double Value, double Error);
|
||||
}
|
||||
@@ -0,0 +1,14 @@
|
||||
using AngouriMath;
|
||||
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// The value of an expression without variables, with as many digits as the computation vouches for.
|
||||
/// </summary>
|
||||
/// <param name="Exact">The exact value when it is rational, or complex with rational parts; otherwise null.</param>
|
||||
/// <param name="Value">The value at the working precision, before it is rounded to its digits.</param>
|
||||
/// <param name="Digits">How many significant digits of the value are certain; zero when none is.</param>
|
||||
/// <param name="Approximate">Whether a part of the value was integrated numerically.</param>
|
||||
/// <param name="EstimatedError">The estimated absolute error of an approximate value; null for any other.</param>
|
||||
/// <param name="Issues">What the model has to know about the value.</param>
|
||||
public sealed record NumericValue(Entity.Number? Exact, Entity.Number.Complex Value, int Digits, bool Approximate, double? EstimatedError, IReadOnlyList<string> Issues);
|
||||
@@ -0,0 +1,11 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// The result of a numerical integration.
|
||||
/// </summary>
|
||||
/// <param name="Value">The estimate of the integral; only an estimate of it when Status is CONVERGED.</param>
|
||||
/// <param name="EstimatedError">The estimated absolute error of the value.</param>
|
||||
/// <param name="Status">How the integration ended.</param>
|
||||
/// <param name="FailurePoint">Where the integrand had no finite value, when Status is NO_VALUE.</param>
|
||||
/// <param name="Evaluations">How often the integrand was evaluated.</param>
|
||||
public sealed record QuadratureResult(double Value, double EstimatedError, QuadratureStatus Status, double? FailurePoint, int Evaluations);
|
||||
@@ -0,0 +1,23 @@
|
||||
namespace AIStudio.Tools.Mathematics;
|
||||
|
||||
/// <summary>
|
||||
/// How a numerical integration ended.
|
||||
/// </summary>
|
||||
public enum QuadratureStatus
|
||||
{
|
||||
/// <summary>
|
||||
/// The estimate settled within the tolerance; its estimated error says how far.
|
||||
/// </summary>
|
||||
CONVERGED = 0,
|
||||
|
||||
/// <summary>
|
||||
/// The estimate did not settle: the integral may diverge or oscillate, or the integrand varies too
|
||||
/// much for the finest step. The value is the last estimate and must not be taken for the integral.
|
||||
/// </summary>
|
||||
NOT_CONVERGED,
|
||||
|
||||
/// <summary>
|
||||
/// The integrand has no finite value at a point inside the interval, such as a pole.
|
||||
/// </summary>
|
||||
NO_VALUE,
|
||||
}
|
||||
Reference in new issue
Block a user