diff --git a/Algorithms/Knapsack/BranchAndBoundKnapsackSolver.cs b/Algorithms/Knapsack/BranchAndBoundKnapsackSolver.cs index a3c9c37e..c3f5b27d 100644 --- a/Algorithms/Knapsack/BranchAndBoundKnapsackSolver.cs +++ b/Algorithms/Knapsack/BranchAndBoundKnapsackSolver.cs @@ -2,164 +2,163 @@ using System.Collections.Generic; using System.Linq; -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +/// +/// Branch and bound Knapsack solver. +/// +/// Type of items in knapsack. +public class BranchAndBoundKnapsackSolver { /// - /// Branch and bound Knapsack solver. + /// Returns the knapsack containing the items that maximize value while not exceeding weight capacity. + /// Construct a tree structure with total number of items + 1 levels, each node have two child nodes, + /// starting with a dummy item root, each following levels are associated with 1 items, construct the + /// tree in breadth first order to identify the optimal item set. /// - /// Type of items in knapsack. - public class BranchAndBoundKnapsackSolver + /// All items to choose from. + /// The maximum weight capacity of the knapsack to be filled. + /// + /// A function that returns the value of the specified item + /// from the items list. + /// + /// + /// A function that returns the weight of the specified item + /// from the items list. + /// + /// + /// The array of items that provides the maximum value of the + /// knapsack without exceeding the specified weight capacity. + /// + public T[] Solve(T[] items, int capacity, Func weightSelector, Func valueSelector) { - /// - /// Returns the knapsack containing the items that maximize value while not exceeding weight capacity. - /// Construct a tree structure with total number of items + 1 levels, each node have two child nodes, - /// starting with a dummy item root, each following levels are associated with 1 items, construct the - /// tree in breadth first order to identify the optimal item set. - /// - /// All items to choose from. - /// The maximum weight capacity of the knapsack to be filled. - /// - /// A function that returns the value of the specified item - /// from the items list. - /// - /// - /// A function that returns the weight of the specified item - /// from the items list. - /// - /// - /// The array of items that provides the maximum value of the - /// knapsack without exceeding the specified weight capacity. - /// - public T[] Solve(T[] items, int capacity, Func weightSelector, Func valueSelector) - { - // This is required for greedy approach in upper bound calculation to work. - items = items.OrderBy(i => valueSelector(i) / weightSelector(i)).ToArray(); + // This is required for greedy approach in upper bound calculation to work. + items = items.OrderBy(i => valueSelector(i) / weightSelector(i)).ToArray(); - // nodesQueue --> used to construct tree in breadth first order - Queue nodesQueue = new(); + // nodesQueue --> used to construct tree in breadth first order + Queue nodesQueue = new(); - // maxCumulativeValue --> maximum value while not exceeding weight capacity. - var maxCumulativeValue = 0.0; + // maxCumulativeValue --> maximum value while not exceeding weight capacity. + var maxCumulativeValue = 0.0; - // starting node, associated with a temporary created dummy item - BranchAndBoundNode root = new(level: -1, taken: false); + // starting node, associated with a temporary created dummy item + BranchAndBoundNode root = new(level: -1, taken: false); - // lastNodeOfOptimalPat --> last item in the optimal item sets identified by this algorithm - BranchAndBoundNode lastNodeOfOptimalPath = root; + // lastNodeOfOptimalPat --> last item in the optimal item sets identified by this algorithm + BranchAndBoundNode lastNodeOfOptimalPath = root; - nodesQueue.Enqueue(root); + nodesQueue.Enqueue(root); + + while (nodesQueue.Count != 0) + { + // parent --> parent node which represents the previous item, may or may not be taken into the knapsack + BranchAndBoundNode parent = nodesQueue.Dequeue(); - while (nodesQueue.Count != 0) + // IF it is the last level, branching cannot be performed + if (parent.Level == items.Length - 1) { - // parent --> parent node which represents the previous item, may or may not be taken into the knapsack - BranchAndBoundNode parent = nodesQueue.Dequeue(); - - // IF it is the last level, branching cannot be performed - if (parent.Level == items.Length - 1) - { - continue; - } - - // create a child node where the associated item is taken into the knapsack - var left = new BranchAndBoundNode(parent.Level + 1, true, parent); - - // create a child node where the associated item is not taken into the knapsack - var right = new BranchAndBoundNode(parent.Level + 1, false, parent); - - // Since the associated item on current level is taken for the first node, - // set the cumulative weight of first node to cumulative weight of parent node + weight of the associated item, - // set the cumulative value of first node to cumulative value of parent node + value of current level's item. - left.CumulativeWeight = parent.CumulativeWeight + weightSelector(items[left.Level]); - left.CumulativeValue = parent.CumulativeValue + valueSelector(items[left.Level]); - right.CumulativeWeight = parent.CumulativeWeight; - right.CumulativeValue = parent.CumulativeValue; - - // IF cumulative weight is smaller than the weight capacity of the knapsack AND - // current cumulative value is larger then the current maxCumulativeValue, update the maxCumulativeValue - if (left.CumulativeWeight <= capacity && left.CumulativeValue > maxCumulativeValue) - { - maxCumulativeValue = left.CumulativeValue; - lastNodeOfOptimalPath = left; - } - - left.UpperBound = ComputeUpperBound(left, items, capacity, weightSelector, valueSelector); - right.UpperBound = ComputeUpperBound(right, items, capacity, weightSelector, valueSelector); - - // IF upperBound of this node is larger than maxCumulativeValue, - // the current path is still possible to reach or surpass the maximum value, - // add current node to nodesQueue so that nodes below it can be further explored - if (left.UpperBound > maxCumulativeValue && left.CumulativeWeight < capacity) - { - nodesQueue.Enqueue(left); - } - - // Cumulative weight is the same as for parent node and < capacity - if (right.UpperBound > maxCumulativeValue) - { - nodesQueue.Enqueue(right); - } + continue; } - return GetItemsFromPath(items, lastNodeOfOptimalPath); - } + // create a child node where the associated item is taken into the knapsack + var left = new BranchAndBoundNode(parent.Level + 1, true, parent); - // determine items taken based on the path - private static T[] GetItemsFromPath(T[] items, BranchAndBoundNode lastNodeOfPath) - { - List takenItems = new(); + // create a child node where the associated item is not taken into the knapsack + var right = new BranchAndBoundNode(parent.Level + 1, false, parent); - // only bogus initial node has no parent - for (var current = lastNodeOfPath; current.Parent is not null; current = current.Parent) + // Since the associated item on current level is taken for the first node, + // set the cumulative weight of first node to cumulative weight of parent node + weight of the associated item, + // set the cumulative value of first node to cumulative value of parent node + value of current level's item. + left.CumulativeWeight = parent.CumulativeWeight + weightSelector(items[left.Level]); + left.CumulativeValue = parent.CumulativeValue + valueSelector(items[left.Level]); + right.CumulativeWeight = parent.CumulativeWeight; + right.CumulativeValue = parent.CumulativeValue; + + // IF cumulative weight is smaller than the weight capacity of the knapsack AND + // current cumulative value is larger then the current maxCumulativeValue, update the maxCumulativeValue + if (left.CumulativeWeight <= capacity && left.CumulativeValue > maxCumulativeValue) { - if(current.IsTaken) - { - takenItems.Add(items[current.Level]); - } + maxCumulativeValue = left.CumulativeValue; + lastNodeOfOptimalPath = left; } - return takenItems.ToArray(); + left.UpperBound = ComputeUpperBound(left, items, capacity, weightSelector, valueSelector); + right.UpperBound = ComputeUpperBound(right, items, capacity, weightSelector, valueSelector); + + // IF upperBound of this node is larger than maxCumulativeValue, + // the current path is still possible to reach or surpass the maximum value, + // add current node to nodesQueue so that nodes below it can be further explored + if (left.UpperBound > maxCumulativeValue && left.CumulativeWeight < capacity) + { + nodesQueue.Enqueue(left); + } + + // Cumulative weight is the same as for parent node and < capacity + if (right.UpperBound > maxCumulativeValue) + { + nodesQueue.Enqueue(right); + } } - /// - /// Returns the upper bound value of a given node. - /// - /// The given node. - /// All items to choose from. - /// The maximum weight capacity of the knapsack to be filled. - /// - /// A function that returns the value of the specified item - /// from the items list. - /// - /// - /// A function that returns the weight of the specified item - /// from the items list. - /// - /// - /// upper bound value of the given node. - /// - private static double ComputeUpperBound(BranchAndBoundNode aNode, T[] items, int capacity, Func weightSelector, Func valueSelector) + return GetItemsFromPath(items, lastNodeOfOptimalPath); + } + + // determine items taken based on the path + private static T[] GetItemsFromPath(T[] items, BranchAndBoundNode lastNodeOfPath) + { + List takenItems = new(); + + // only bogus initial node has no parent + for (var current = lastNodeOfPath; current.Parent is not null; current = current.Parent) { - var upperBound = aNode.CumulativeValue; - var availableWeight = capacity - aNode.CumulativeWeight; - var nextLevel = aNode.Level + 1; + if(current.IsTaken) + { + takenItems.Add(items[current.Level]); + } + } + + return takenItems.ToArray(); + } - while (availableWeight > 0 && nextLevel < items.Length) + /// + /// Returns the upper bound value of a given node. + /// + /// The given node. + /// All items to choose from. + /// The maximum weight capacity of the knapsack to be filled. + /// + /// A function that returns the value of the specified item + /// from the items list. + /// + /// + /// A function that returns the weight of the specified item + /// from the items list. + /// + /// + /// upper bound value of the given node. + /// + private static double ComputeUpperBound(BranchAndBoundNode aNode, T[] items, int capacity, Func weightSelector, Func valueSelector) + { + var upperBound = aNode.CumulativeValue; + var availableWeight = capacity - aNode.CumulativeWeight; + var nextLevel = aNode.Level + 1; + + while (availableWeight > 0 && nextLevel < items.Length) + { + if (weightSelector(items[nextLevel]) <= availableWeight) { - if (weightSelector(items[nextLevel]) <= availableWeight) - { - upperBound += valueSelector(items[nextLevel]); - availableWeight -= weightSelector(items[nextLevel]); - } - else - { - upperBound += valueSelector(items[nextLevel]) / weightSelector(items[nextLevel]) * availableWeight; - availableWeight = 0; - } - - nextLevel++; + upperBound += valueSelector(items[nextLevel]); + availableWeight -= weightSelector(items[nextLevel]); + } + else + { + upperBound += valueSelector(items[nextLevel]) / weightSelector(items[nextLevel]) * availableWeight; + availableWeight = 0; } - return upperBound; + nextLevel++; } + + return upperBound; } } diff --git a/Algorithms/Knapsack/BranchAndBoundNode.cs b/Algorithms/Knapsack/BranchAndBoundNode.cs index 3d873739..0a65f16a 100644 --- a/Algorithms/Knapsack/BranchAndBoundNode.cs +++ b/Algorithms/Knapsack/BranchAndBoundNode.cs @@ -1,30 +1,29 @@ -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +public class BranchAndBoundNode { - public class BranchAndBoundNode - { - // isTaken --> true = the item where index = level is taken, vice versa - public bool IsTaken { get; } + // isTaken --> true = the item where index = level is taken, vice versa + public bool IsTaken { get; } - // cumulativeWeight --> um of weight of item associated in each nodes starting from root to this node (only item that is taken) - public int CumulativeWeight { get; set; } + // cumulativeWeight --> um of weight of item associated in each nodes starting from root to this node (only item that is taken) + public int CumulativeWeight { get; set; } - // cumulativeValue --> sum of value of item associated in each nodes starting from root to this node (only item that is taken) - public double CumulativeValue { get; set; } + // cumulativeValue --> sum of value of item associated in each nodes starting from root to this node (only item that is taken) + public double CumulativeValue { get; set; } - // upperBound --> largest possible value after taking/not taking the item associated to this node (fractional) - public double UpperBound { get; set; } + // upperBound --> largest possible value after taking/not taking the item associated to this node (fractional) + public double UpperBound { get; set; } - // level --> level of the node in the tree structure - public int Level { get; } + // level --> level of the node in the tree structure + public int Level { get; } - // parent node - public BranchAndBoundNode? Parent { get; } + // parent node + public BranchAndBoundNode? Parent { get; } - public BranchAndBoundNode(int level, bool taken, BranchAndBoundNode? parent = null) - { - Level = level; - IsTaken = taken; - Parent = parent; - } + public BranchAndBoundNode(int level, bool taken, BranchAndBoundNode? parent = null) + { + Level = level; + IsTaken = taken; + Parent = parent; } } diff --git a/Algorithms/Knapsack/DynamicProgrammingKnapsackSolver.cs b/Algorithms/Knapsack/DynamicProgrammingKnapsackSolver.cs index 2b5abc18..db60b55d 100644 --- a/Algorithms/Knapsack/DynamicProgrammingKnapsackSolver.cs +++ b/Algorithms/Knapsack/DynamicProgrammingKnapsackSolver.cs @@ -1,101 +1,100 @@ -using System; +using System; using System.Collections.Generic; -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +/// +/// Dynamic Programming Knapsack solver. +/// +/// Type of items in knapsack. +public class DynamicProgrammingKnapsackSolver { /// - /// Dynamic Programming Knapsack solver. + /// Returns the knapsack containing the items that + /// maximize value while not exceeding weight capacity. /// - /// Type of items in knapsack. - public class DynamicProgrammingKnapsackSolver + /// The list of items from which we select ones to be in the knapsack. + /// + /// The maximum weight capacity of the knapsack + /// to be filled. Only integer values of this capacity are tried. If + /// a greater resolution is needed, multiply the + /// weights/capacity by a factor of 10. + /// + /// + /// A function that returns the value of the specified item + /// from the items list. + /// + /// + /// A function that returns the weight of the specified item + /// from the items list. + /// + /// + /// The array of items that provides the maximum value of the + /// knapsack without exceeding the specified weight capacity. + /// + public T[] Solve(T[] items, int capacity, Func weightSelector, Func valueSelector) { - /// - /// Returns the knapsack containing the items that - /// maximize value while not exceeding weight capacity. - /// - /// The list of items from which we select ones to be in the knapsack. - /// - /// The maximum weight capacity of the knapsack - /// to be filled. Only integer values of this capacity are tried. If - /// a greater resolution is needed, multiply the - /// weights/capacity by a factor of 10. - /// - /// - /// A function that returns the value of the specified item - /// from the items list. - /// - /// - /// A function that returns the weight of the specified item - /// from the items list. - /// - /// - /// The array of items that provides the maximum value of the - /// knapsack without exceeding the specified weight capacity. - /// - public T[] Solve(T[] items, int capacity, Func weightSelector, Func valueSelector) - { - var cache = Tabulate(items, weightSelector, valueSelector, capacity); - return GetOptimalItems(items, weightSelector, cache, capacity); - } + var cache = Tabulate(items, weightSelector, valueSelector, capacity); + return GetOptimalItems(items, weightSelector, cache, capacity); + } - private static T[] GetOptimalItems(T[] items, Func weightSelector, double[,] cache, int capacity) - { - var currentCapacity = capacity; + private static T[] GetOptimalItems(T[] items, Func weightSelector, double[,] cache, int capacity) + { + var currentCapacity = capacity; - var result = new List(); - for (var i = items.Length - 1; i >= 0; i--) + var result = new List(); + for (var i = items.Length - 1; i >= 0; i--) + { + if (cache[i + 1, currentCapacity] > cache[i, currentCapacity]) { - if (cache[i + 1, currentCapacity] > cache[i, currentCapacity]) - { - var item = items[i]; - result.Add(item); - currentCapacity -= weightSelector(item); - } + var item = items[i]; + result.Add(item); + currentCapacity -= weightSelector(item); } - - result.Reverse(); // we added items back to front - return result.ToArray(); } - private static double[,] Tabulate( - T[] items, - Func weightSelector, - Func valueSelector, - int maxCapacity) + result.Reverse(); // we added items back to front + return result.ToArray(); + } + + private static double[,] Tabulate( + T[] items, + Func weightSelector, + Func valueSelector, + int maxCapacity) + { + // Store the incremental results in a bottom up manner + var n = items.Length; + var results = new double[n + 1, maxCapacity + 1]; + for (var i = 0; i <= n; i++) { - // Store the incremental results in a bottom up manner - var n = items.Length; - var results = new double[n + 1, maxCapacity + 1]; - for (var i = 0; i <= n; i++) + for (var w = 0; w <= maxCapacity; w++) { - for (var w = 0; w <= maxCapacity; w++) + if (i == 0 || w == 0) + { + // If we have no items to take, or + // if we have no capacity in our knapsack + // we cannot possibly have any value + results[i, w] = 0; + } + else if (weightSelector(items[i - 1]) <= w) { - if (i == 0 || w == 0) - { - // If we have no items to take, or - // if we have no capacity in our knapsack - // we cannot possibly have any value - results[i, w] = 0; - } - else if (weightSelector(items[i - 1]) <= w) - { - // Decide if it is better to take or not take this item - var iut = items[i - 1]; // iut = Item under test - var vut = valueSelector(iut); // vut = Value of item under test - var wut = weightSelector(iut); // wut = Weight of item under test - var valueIfTaken = vut + results[i - 1, w - wut]; - var valueIfNotTaken = results[i - 1, w]; - results[i, w] = Math.Max(valueIfTaken, valueIfNotTaken); - } - else - { - // There is not enough room to take this item - results[i, w] = results[i - 1, w]; - } + // Decide if it is better to take or not take this item + var iut = items[i - 1]; // iut = Item under test + var vut = valueSelector(iut); // vut = Value of item under test + var wut = weightSelector(iut); // wut = Weight of item under test + var valueIfTaken = vut + results[i - 1, w - wut]; + var valueIfNotTaken = results[i - 1, w]; + results[i, w] = Math.Max(valueIfTaken, valueIfNotTaken); + } + else + { + // There is not enough room to take this item + results[i, w] = results[i - 1, w]; } } - - return results; } + + return results; } } diff --git a/Algorithms/Knapsack/IHeuristicKnapsackSolver.cs b/Algorithms/Knapsack/IHeuristicKnapsackSolver.cs index b5403d2b..e73bcdd1 100644 --- a/Algorithms/Knapsack/IHeuristicKnapsackSolver.cs +++ b/Algorithms/Knapsack/IHeuristicKnapsackSolver.cs @@ -1,25 +1,24 @@ -using System; +using System; -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +/// +/// Solves knapsack problem using some heuristics +/// Sum of values of taken items -> max +/// Sum of weights of taken items. <= capacity. +/// +/// Type of items in knapsack. +public interface IHeuristicKnapsackSolver { /// /// Solves knapsack problem using some heuristics /// Sum of values of taken items -> max /// Sum of weights of taken items. <= capacity. /// - /// Type of items in knapsack. - public interface IHeuristicKnapsackSolver - { - /// - /// Solves knapsack problem using some heuristics - /// Sum of values of taken items -> max - /// Sum of weights of taken items. <= capacity. - /// - /// All items to choose from. - /// How much weight we can take. - /// Maps item to its weight. - /// Maps item to its value. - /// Items that were chosen. - T[] Solve(T[] items, double capacity, Func weightSelector, Func valueSelector); - } + /// All items to choose from. + /// How much weight we can take. + /// Maps item to its weight. + /// Maps item to its value. + /// Items that were chosen. + T[] Solve(T[] items, double capacity, Func weightSelector, Func valueSelector); } diff --git a/Algorithms/Knapsack/IKnapsackSolver.cs b/Algorithms/Knapsack/IKnapsackSolver.cs index fb1cdaa4..3eca055a 100644 --- a/Algorithms/Knapsack/IKnapsackSolver.cs +++ b/Algorithms/Knapsack/IKnapsackSolver.cs @@ -1,12 +1,11 @@ -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +/// +/// Solves knapsack problem: +/// to maximize sum of values of taken items, +/// while sum of weights of taken items is less than capacity. +/// +/// Type of items in knapsack. +public interface IKnapsackSolver : IHeuristicKnapsackSolver { - /// - /// Solves knapsack problem: - /// to maximize sum of values of taken items, - /// while sum of weights of taken items is less than capacity. - /// - /// Type of items in knapsack. - public interface IKnapsackSolver : IHeuristicKnapsackSolver - { - } } diff --git a/Algorithms/Knapsack/NaiveKnapsackSolver.cs b/Algorithms/Knapsack/NaiveKnapsackSolver.cs index f0bd1dd5..841b8750 100644 --- a/Algorithms/Knapsack/NaiveKnapsackSolver.cs +++ b/Algorithms/Knapsack/NaiveKnapsackSolver.cs @@ -1,38 +1,37 @@ -using System; +using System; using System.Collections.Generic; -namespace Algorithms.Knapsack +namespace Algorithms.Knapsack; + +/// +/// Greedy heurictic solver. +/// +/// Type of items in knapsack. +public class NaiveKnapsackSolver : IHeuristicKnapsackSolver { /// - /// Greedy heurictic solver. + /// TODO. /// - /// Type of items in knapsack. - public class NaiveKnapsackSolver : IHeuristicKnapsackSolver + /// TODO. 2. + /// TODO. 3. + /// TODO. 4. + /// TODO. 5. + /// TODO. 6. + public T[] Solve(T[] items, double capacity, Func weightSelector, Func valueSelector) { - /// - /// TODO. - /// - /// TODO. 2. - /// TODO. 3. - /// TODO. 4. - /// TODO. 5. - /// TODO. 6. - public T[] Solve(T[] items, double capacity, Func weightSelector, Func valueSelector) - { - var weight = 0d; - var left = new List(); + var weight = 0d; + var left = new List(); - foreach (var item in items) + foreach (var item in items) + { + var weightDelta = weightSelector(item); + if (weight + weightDelta <= capacity) { - var weightDelta = weightSelector(item); - if (weight + weightDelta <= capacity) - { - weight += weightDelta; - left.Add(item); - } + weight += weightDelta; + left.Add(item); } - - return left.ToArray(); } + + return left.ToArray(); } } diff --git a/Algorithms/LinearAlgebra/Distances/Euclidean.cs b/Algorithms/LinearAlgebra/Distances/Euclidean.cs index fc893036..67fd4129 100644 --- a/Algorithms/LinearAlgebra/Distances/Euclidean.cs +++ b/Algorithms/LinearAlgebra/Distances/Euclidean.cs @@ -1,31 +1,27 @@ using System; -using System.Collections.Generic; using System.Linq; -using System.Text; -using System.Threading.Tasks; -namespace Algorithms.LinearAlgebra.Distances +namespace Algorithms.LinearAlgebra.Distances; + +/// +/// Implementation for Euclidean distance. +/// +public static class Euclidean { /// - /// Implementation for Euclidean distance. + /// Calculate Euclidean distance for two N-Dimensional points. /// - public static class Euclidean + /// First N-Dimensional point. + /// Second N-Dimensional point. + /// Calculated Euclidean distance. + public static double Distance(double[] point1, double[] point2) { - /// - /// Calculate Euclidean distance for two N-Dimensional points. - /// - /// First N-Dimensional point. - /// Second N-Dimensional point. - /// Calculated Euclidean distance. - public static double Distance(double[] point1, double[] point2) + if (point1.Length != point2.Length) { - if (point1.Length != point2.Length) - { - throw new ArgumentException("Both points should have the same dimensionality"); - } - - // distance = sqrt((x1-y1)^2 + (x2-y2)^2 + ... + (xn-yn)^2) - return Math.Sqrt(point1.Zip(point2, (x1, x2) => (x1 - x2) * (x1 - x2)).Sum()); + throw new ArgumentException("Both points should have the same dimensionality"); } + + // distance = sqrt((x1-y1)^2 + (x2-y2)^2 + ... + (xn-yn)^2) + return Math.Sqrt(point1.Zip(point2, (x1, x2) => (x1 - x2) * (x1 - x2)).Sum()); } } diff --git a/Algorithms/LinearAlgebra/Distances/Manhattan.cs b/Algorithms/LinearAlgebra/Distances/Manhattan.cs index 1b47563f..ea0dd96f 100644 --- a/Algorithms/LinearAlgebra/Distances/Manhattan.cs +++ b/Algorithms/LinearAlgebra/Distances/Manhattan.cs @@ -1,35 +1,31 @@ using System; -using System.Collections.Generic; using System.Linq; -using System.Text; -using System.Threading.Tasks; -namespace Algorithms.LinearAlgebra.Distances +namespace Algorithms.LinearAlgebra.Distances; + +/// +/// Implementation fo Manhattan distance. +/// It is the sum of the lengths of the projections of the line segment between the points onto the coordinate axes. +/// In other words, it is the sum of absolute difference between the measures in all dimensions of two points. +/// +/// Its commonly used in regression analysis. +/// +public static class Manhattan { /// - /// Implementation fo Manhattan distance. - /// It is the sum of the lengths of the projections of the line segment between the points onto the coordinate axes. - /// In other words, it is the sum of absolute difference between the measures in all dimensions of two points. - /// - /// Its commonly used in regression analysis. + /// Calculate Manhattan distance for two N-Dimensional points. /// - public static class Manhattan + /// First N-Dimensional point. + /// Second N-Dimensional point. + /// Calculated Manhattan distance. + public static double Distance(double[] point1, double[] point2) { - /// - /// Calculate Manhattan distance for two N-Dimensional points. - /// - /// First N-Dimensional point. - /// Second N-Dimensional point. - /// Calculated Manhattan distance. - public static double Distance(double[] point1, double[] point2) + if (point1.Length != point2.Length) { - if (point1.Length != point2.Length) - { - throw new ArgumentException("Both points should have the same dimensionality"); - } - - // distance = |x1-y1| + |x2-y2| + ... + |xn-yn| - return point1.Zip(point2, (x1, x2) => Math.Abs(x1 - x2)).Sum(); + throw new ArgumentException("Both points should have the same dimensionality"); } + + // distance = |x1-y1| + |x2-y2| + ... + |xn-yn| + return point1.Zip(point2, (x1, x2) => Math.Abs(x1 - x2)).Sum(); } } diff --git a/Algorithms/LinearAlgebra/Eigenvalue/PowerIteration.cs b/Algorithms/LinearAlgebra/Eigenvalue/PowerIteration.cs index 51e3cf85..8e503992 100644 --- a/Algorithms/LinearAlgebra/Eigenvalue/PowerIteration.cs +++ b/Algorithms/LinearAlgebra/Eigenvalue/PowerIteration.cs @@ -2,86 +2,85 @@ using System.Linq; using Utilities.Extensions; -namespace Algorithms.LinearAlgebra.Eigenvalue +namespace Algorithms.LinearAlgebra.Eigenvalue; + +/// +/// Power iteration method - eigenvalue numeric algorithm, based on recurrent relation: +/// Li+1 = (A * Li) / || A * Li ||, where Li - eigenvector approximation. +/// +public static class PowerIteration { /// - /// Power iteration method - eigenvalue numeric algorithm, based on recurrent relation: - /// Li+1 = (A * Li) / || A * Li ||, where Li - eigenvector approximation. + /// Returns approximation of the dominant eigenvalue and eigenvector of matrix. /// - public static class PowerIteration + /// + /// + /// The algorithm will not converge if the start vector is orthogonal to the eigenvector. + /// + /// + /// The matrix must be square-shaped. + /// + /// + /// Source square-shaped matrix. + /// Start vector. + /// Accuracy of the result. + /// Dominant eigenvalue and eigenvector pair. + /// The matrix is not square-shaped. + /// The length of the start vector doesn't equal the size of the source matrix. + public static (double eigenvalue, double[] eigenvector) Dominant( + double[,] source, + double[] startVector, + double error = 0.00001) { - /// - /// Returns approximation of the dominant eigenvalue and eigenvector of matrix. - /// - /// - /// - /// The algorithm will not converge if the start vector is orthogonal to the eigenvector. - /// - /// - /// The matrix must be square-shaped. - /// - /// - /// Source square-shaped matrix. - /// Start vector. - /// Accuracy of the result. - /// Dominant eigenvalue and eigenvector pair. - /// The matrix is not square-shaped. - /// The length of the start vector doesn't equal the size of the source matrix. - public static (double eigenvalue, double[] eigenvector) Dominant( - double[,] source, - double[] startVector, - double error = 0.00001) + if (source.GetLength(0) != source.GetLength(1)) { - if (source.GetLength(0) != source.GetLength(1)) - { - throw new ArgumentException("The source matrix is not square-shaped."); - } - - if (source.GetLength(0) != startVector.Length) - { - throw new ArgumentException( - "The length of the start vector doesn't equal the size of the source matrix."); - } - - double eigenNorm; - double[] previousEigenVector; - double[] currentEigenVector = startVector; + throw new ArgumentException("The source matrix is not square-shaped."); + } - do - { - previousEigenVector = currentEigenVector; - currentEigenVector = source.Multiply( - previousEigenVector.ToColumnVector()) - .ToRowVector(); + if (source.GetLength(0) != startVector.Length) + { + throw new ArgumentException( + "The length of the start vector doesn't equal the size of the source matrix."); + } - eigenNorm = currentEigenVector.Magnitude(); - currentEigenVector = currentEigenVector.Select(x => x / eigenNorm).ToArray(); - } - while (Math.Abs(currentEigenVector.Dot(previousEigenVector)) < 1.0 - error); + double eigenNorm; + double[] previousEigenVector; + double[] currentEigenVector = startVector; - var eigenvalue = source.Multiply(currentEigenVector.ToColumnVector()).ToRowVector().Magnitude(); + do + { + previousEigenVector = currentEigenVector; + currentEigenVector = source.Multiply( + previousEigenVector.ToColumnVector()) + .ToRowVector(); - return (eigenvalue, eigenvector: currentEigenVector); + eigenNorm = currentEigenVector.Magnitude(); + currentEigenVector = currentEigenVector.Select(x => x / eigenNorm).ToArray(); } + while (Math.Abs(currentEigenVector.Dot(previousEigenVector)) < 1.0 - error); + + var eigenvalue = source.Multiply(currentEigenVector.ToColumnVector()).ToRowVector().Magnitude(); - /// - /// Returns approximation of the dominant eigenvalue and eigenvector of matrix. - /// Random normalized vector is used as the start vector to decrease chance of orthogonality to the eigenvector. - /// - /// - /// - /// The algorithm will not converge if the start vector is orthogonal to the eigenvector. - /// - /// - /// The matrix should be square-shaped. - /// - /// - /// Source square-shaped matrix. - /// Accuracy of the result. - /// Dominant eigenvalue and eigenvector pair. - /// The matrix is not square-shaped. - /// The length of the start vector doesn't equal the size of the source matrix. - public static (double eigenvalue, double[] eigenvector) Dominant(double[,] source, double error = 0.00001) => - Dominant(source, new Random().NextVector(source.GetLength(1)), error); + return (eigenvalue, eigenvector: currentEigenVector); } + + /// + /// Returns approximation of the dominant eigenvalue and eigenvector of matrix. + /// Random normalized vector is used as the start vector to decrease chance of orthogonality to the eigenvector. + /// + /// + /// + /// The algorithm will not converge if the start vector is orthogonal to the eigenvector. + /// + /// + /// The matrix should be square-shaped. + /// + /// + /// Source square-shaped matrix. + /// Accuracy of the result. + /// Dominant eigenvalue and eigenvector pair. + /// The matrix is not square-shaped. + /// The length of the start vector doesn't equal the size of the source matrix. + public static (double eigenvalue, double[] eigenvector) Dominant(double[,] source, double error = 0.00001) => + Dominant(source, new Random().NextVector(source.GetLength(1)), error); } diff --git a/Algorithms/ModularArithmetic/ChineseRemainderTheorem.cs b/Algorithms/ModularArithmetic/ChineseRemainderTheorem.cs index a8b702a6..1eeaef76 100644 --- a/Algorithms/ModularArithmetic/ChineseRemainderTheorem.cs +++ b/Algorithms/ModularArithmetic/ChineseRemainderTheorem.cs @@ -3,189 +3,188 @@ using System.Linq; using System.Numerics; -namespace Algorithms.ModularArithmetic +namespace Algorithms.ModularArithmetic; + +/// +/// Chinese Remainder Theorem: https://en.wikipedia.org/wiki/Chinese_remainder_theorem. +/// +public static class ChineseRemainderTheorem { /// - /// Chinese Remainder Theorem: https://en.wikipedia.org/wiki/Chinese_remainder_theorem. + /// The Chinese Remainder Theorem is used to compute x for given set of pairs of integers (a_i, n_i) with: + /// + /// x = a_0 mod n_0 + /// x = a_1 mod n_1 + /// ... + /// x = a_k mod n_k + /// + /// for 0 <= i < k, for some positive integer k. Additional requirements are: + /// + /// n_i > 1 for 0 <= i < k + /// n_i and n_j are coprime, for all 0 <= i < j < k + /// 0 <= a_i < n_i, for all 0 <= i < k + /// 0 <= x < n_0 * n_1 * ... * n_(k-1) + /// /// - public static class ChineseRemainderTheorem + /// An ordered list of a_0, a_1, ..., a_k. + /// An ordered list of n_0, n_1, ..., n_k. + /// The value x. + /// If any of the requirements is not fulfilled. + public static long Compute(List listOfAs, List listOfNs) { - /// - /// The Chinese Remainder Theorem is used to compute x for given set of pairs of integers (a_i, n_i) with: - /// - /// x = a_0 mod n_0 - /// x = a_1 mod n_1 - /// ... - /// x = a_k mod n_k - /// - /// for 0 <= i < k, for some positive integer k. Additional requirements are: - /// - /// n_i > 1 for 0 <= i < k - /// n_i and n_j are coprime, for all 0 <= i < j < k - /// 0 <= a_i < n_i, for all 0 <= i < k - /// 0 <= x < n_0 * n_1 * ... * n_(k-1) - /// - /// - /// An ordered list of a_0, a_1, ..., a_k. - /// An ordered list of n_0, n_1, ..., n_k. - /// The value x. - /// If any of the requirements is not fulfilled. - public static long Compute(List listOfAs, List listOfNs) + // Check the requirements for this algorithm: + CheckRequirements(listOfAs, listOfNs); + + // For performance-reasons compute the product of all n_i as prodN, because we're going to need access to (prodN / n_i) for all i: + var prodN = 1L; + foreach (var n in listOfNs) { - // Check the requirements for this algorithm: - CheckRequirements(listOfAs, listOfNs); + prodN *= n; + } - // For performance-reasons compute the product of all n_i as prodN, because we're going to need access to (prodN / n_i) for all i: - var prodN = 1L; - foreach (var n in listOfNs) - { - prodN *= n; - } + var result = 0L; + for (var i = 0; i < listOfNs.Count; i++) + { + var a_i = listOfAs[i]; + var n_i = listOfNs[i]; + var modulus_i = prodN / n_i; - var result = 0L; - for (var i = 0; i < listOfNs.Count; i++) - { - var a_i = listOfAs[i]; - var n_i = listOfNs[i]; - var modulus_i = prodN / n_i; + var bezout_modulus_i = ExtendedEuclideanAlgorithm.Compute(n_i, modulus_i).bezoutB; + result += a_i * bezout_modulus_i * modulus_i; + } - var bezout_modulus_i = ExtendedEuclideanAlgorithm.Compute(n_i, modulus_i).bezoutB; - result += a_i * bezout_modulus_i * modulus_i; - } + // Make sure, result is in [0, prodN). + result %= prodN; + if (result < 0) + { + result += prodN; + } - // Make sure, result is in [0, prodN). - result %= prodN; - if (result < 0) - { - result += prodN; - } + return result; + } - return result; - } + /// + /// The Chinese Remainder Theorem is used to compute x for given set of pairs of integers (a_i, n_i) with: + /// + /// x = a_0 mod n_0 + /// x = a_1 mod n_1 + /// ... + /// x = a_k mod n_k + /// + /// for 0 <= i < k, for some positive integer k. Additional requirements are: + /// + /// n_i > 1 for 0 <= i < k + /// n_i and n_j are coprime, for all 0 <= i < j < k + /// 0 <= a_i < n_i, for all 0 <= i < k + /// 0 <= x < n_0 * n_1 * ... * n_(k-1) + /// + /// + /// An ordered list of a_0, a_1, ..., a_k. + /// An ordered list of n_0, n_1, ..., n_k. + /// The value x. + /// If any of the requirements is not fulfilled. + public static BigInteger Compute(List listOfAs, List listOfNs) + { + // Check the requirements for this algorithm: + CheckRequirements(listOfAs, listOfNs); - /// - /// The Chinese Remainder Theorem is used to compute x for given set of pairs of integers (a_i, n_i) with: - /// - /// x = a_0 mod n_0 - /// x = a_1 mod n_1 - /// ... - /// x = a_k mod n_k - /// - /// for 0 <= i < k, for some positive integer k. Additional requirements are: - /// - /// n_i > 1 for 0 <= i < k - /// n_i and n_j are coprime, for all 0 <= i < j < k - /// 0 <= a_i < n_i, for all 0 <= i < k - /// 0 <= x < n_0 * n_1 * ... * n_(k-1) - /// - /// - /// An ordered list of a_0, a_1, ..., a_k. - /// An ordered list of n_0, n_1, ..., n_k. - /// The value x. - /// If any of the requirements is not fulfilled. - public static BigInteger Compute(List listOfAs, List listOfNs) + // For performance-reasons compute the product of all n_i as prodN, because we're going to need access to (prodN / n_i) for all i: + var prodN = BigInteger.One; + foreach (var n in listOfNs) { - // Check the requirements for this algorithm: - CheckRequirements(listOfAs, listOfNs); + prodN *= n; + } - // For performance-reasons compute the product of all n_i as prodN, because we're going to need access to (prodN / n_i) for all i: - var prodN = BigInteger.One; - foreach (var n in listOfNs) - { - prodN *= n; - } + var result = BigInteger.Zero; + for (var i = 0; i < listOfNs.Count; i++) + { + var a_i = listOfAs[i]; + var n_i = listOfNs[i]; + var modulus_i = prodN / n_i; - var result = BigInteger.Zero; - for (var i = 0; i < listOfNs.Count; i++) - { - var a_i = listOfAs[i]; - var n_i = listOfNs[i]; - var modulus_i = prodN / n_i; + var bezout_modulus_i = ExtendedEuclideanAlgorithm.Compute(n_i, modulus_i).bezoutB; + result += a_i * bezout_modulus_i * modulus_i; + } - var bezout_modulus_i = ExtendedEuclideanAlgorithm.Compute(n_i, modulus_i).bezoutB; - result += a_i * bezout_modulus_i * modulus_i; - } + // Make sure, result is in [0, prodN). + result %= prodN; + if (result < 0) + { + result += prodN; + } - // Make sure, result is in [0, prodN). - result %= prodN; - if (result < 0) - { - result += prodN; - } + return result; + } - return result; + /// + /// Checks the requirements for the algorithm and throws an ArgumentException if they are not being met. + /// + /// An ordered list of a_0, a_1, ..., a_k. + /// An ordered list of n_0, n_1, ..., n_k. + /// If any of the requirements is not fulfilled. + private static void CheckRequirements(List listOfAs, List listOfNs) + { + if (listOfAs == null || listOfNs == null || listOfAs.Count != listOfNs.Count) + { + throw new ArgumentException("The parameters 'listOfAs' and 'listOfNs' must not be null and have to be of equal length!"); } - /// - /// Checks the requirements for the algorithm and throws an ArgumentException if they are not being met. - /// - /// An ordered list of a_0, a_1, ..., a_k. - /// An ordered list of n_0, n_1, ..., n_k. - /// If any of the requirements is not fulfilled. - private static void CheckRequirements(List listOfAs, List listOfNs) + if (listOfNs.Any(x => x <= 1)) { - if (listOfAs == null || listOfNs == null || listOfAs.Count != listOfNs.Count) - { - throw new ArgumentException("The parameters 'listOfAs' and 'listOfNs' must not be null and have to be of equal length!"); - } - - if (listOfNs.Any(x => x <= 1)) - { - throw new ArgumentException($"The value {listOfNs.First(x => x <= 1)} for some n_i is smaller than or equal to 1."); - } + throw new ArgumentException($"The value {listOfNs.First(x => x <= 1)} for some n_i is smaller than or equal to 1."); + } - if (listOfAs.Any(x => x < 0)) - { - throw new ArgumentException($"The value {listOfAs.First(x => x < 0)} for some a_i is smaller than 0."); - } + if (listOfAs.Any(x => x < 0)) + { + throw new ArgumentException($"The value {listOfAs.First(x => x < 0)} for some a_i is smaller than 0."); + } - // Check if all pairs of (n_i, n_j) are coprime: - for (var i = 0; i < listOfNs.Count; i++) + // Check if all pairs of (n_i, n_j) are coprime: + for (var i = 0; i < listOfNs.Count; i++) + { + for (var j = i + 1; j < listOfNs.Count; j++) { - for (var j = i + 1; j < listOfNs.Count; j++) + long gcd; + if ((gcd = ExtendedEuclideanAlgorithm.Compute(listOfNs[i], listOfNs[j]).gcd) != 1L) { - long gcd; - if ((gcd = ExtendedEuclideanAlgorithm.Compute(listOfNs[i], listOfNs[j]).gcd) != 1L) - { - throw new ArgumentException($"The GCD of n_{i} = {listOfNs[i]} and n_{j} = {listOfNs[j]} equals {gcd} and thus these values aren't coprime."); - } + throw new ArgumentException($"The GCD of n_{i} = {listOfNs[i]} and n_{j} = {listOfNs[j]} equals {gcd} and thus these values aren't coprime."); } } } + } - /// - /// Checks the requirements for the algorithm and throws an ArgumentException if they are not being met. - /// - /// An ordered list of a_0, a_1, ..., a_k. - /// An ordered list of n_0, n_1, ..., n_k. - /// If any of the requirements is not fulfilled. - private static void CheckRequirements(List listOfAs, List listOfNs) + /// + /// Checks the requirements for the algorithm and throws an ArgumentException if they are not being met. + /// + /// An ordered list of a_0, a_1, ..., a_k. + /// An ordered list of n_0, n_1, ..., n_k. + /// If any of the requirements is not fulfilled. + private static void CheckRequirements(List listOfAs, List listOfNs) + { + if (listOfAs == null || listOfNs == null || listOfAs.Count != listOfNs.Count) { - if (listOfAs == null || listOfNs == null || listOfAs.Count != listOfNs.Count) - { - throw new ArgumentException("The parameters 'listOfAs' and 'listOfNs' must not be null and have to be of equal length!"); - } + throw new ArgumentException("The parameters 'listOfAs' and 'listOfNs' must not be null and have to be of equal length!"); + } - if (listOfNs.Any(x => x <= 1)) - { - throw new ArgumentException($"The value {listOfNs.First(x => x <= 1)} for some n_i is smaller than or equal to 1."); - } + if (listOfNs.Any(x => x <= 1)) + { + throw new ArgumentException($"The value {listOfNs.First(x => x <= 1)} for some n_i is smaller than or equal to 1."); + } - if (listOfAs.Any(x => x < 0)) - { - throw new ArgumentException($"The value {listOfAs.First(x => x < 0)} for some a_i is smaller than 0."); - } + if (listOfAs.Any(x => x < 0)) + { + throw new ArgumentException($"The value {listOfAs.First(x => x < 0)} for some a_i is smaller than 0."); + } - // Check if all pairs of (n_i, n_j) are coprime: - for (var i = 0; i < listOfNs.Count; i++) + // Check if all pairs of (n_i, n_j) are coprime: + for (var i = 0; i < listOfNs.Count; i++) + { + for (var j = i + 1; j < listOfNs.Count; j++) { - for (var j = i + 1; j < listOfNs.Count; j++) + BigInteger gcd; + if ((gcd = ExtendedEuclideanAlgorithm.Compute(listOfNs[i], listOfNs[j]).gcd) != BigInteger.One) { - BigInteger gcd; - if ((gcd = ExtendedEuclideanAlgorithm.Compute(listOfNs[i], listOfNs[j]).gcd) != BigInteger.One) - { - throw new ArgumentException($"The GCD of n_{i} = {listOfNs[i]} and n_{j} = {listOfNs[j]} equals {gcd} and thus these values aren't coprime."); - } + throw new ArgumentException($"The GCD of n_{i} = {listOfNs[i]} and n_{j} = {listOfNs[j]} equals {gcd} and thus these values aren't coprime."); } } } diff --git a/Algorithms/ModularArithmetic/ExtendedEuclideanAlgorithm.cs b/Algorithms/ModularArithmetic/ExtendedEuclideanAlgorithm.cs index 19481fcd..aed5b219 100644 --- a/Algorithms/ModularArithmetic/ExtendedEuclideanAlgorithm.cs +++ b/Algorithms/ModularArithmetic/ExtendedEuclideanAlgorithm.cs @@ -1,95 +1,94 @@ using System.Numerics; -namespace Algorithms.ModularArithmetic +namespace Algorithms.ModularArithmetic; + +/// +/// Extended Euclidean algorithm: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm. +/// +public static class ExtendedEuclideanAlgorithm { /// - /// Extended Euclidean algorithm: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm. + /// Computes the greatest common divisor (gcd) of integers a and b, also the coefficients of Bézout's identity, + /// which are integers x and y such that a*bezoutCoefficientOfA + b*bezoutCoefficientOfB = gcd(a, b). /// - public static class ExtendedEuclideanAlgorithm + /// Input number. + /// Second input number. + /// A record of ExtendedEuclideanAlgorithmResult containing the bezout coefficients of a and b as well as the gcd(a,b). + public static ExtendedEuclideanAlgorithmResult Compute(long a, long b) { - /// - /// Computes the greatest common divisor (gcd) of integers a and b, also the coefficients of Bézout's identity, - /// which are integers x and y such that a*bezoutCoefficientOfA + b*bezoutCoefficientOfB = gcd(a, b). - /// - /// Input number. - /// Second input number. - /// A record of ExtendedEuclideanAlgorithmResult containing the bezout coefficients of a and b as well as the gcd(a,b). - public static ExtendedEuclideanAlgorithmResult Compute(long a, long b) - { - long quotient; - long tmp; - var s = 0L; - var bezoutOfA = 1L; - var r = b; - var gcd = a; - var bezoutOfB = 0L; - - while (r != 0) - { - quotient = gcd / r; // integer division + long quotient; + long tmp; + var s = 0L; + var bezoutOfA = 1L; + var r = b; + var gcd = a; + var bezoutOfB = 0L; - tmp = gcd; - gcd = r; - r = tmp - quotient * r; - - tmp = bezoutOfA; - bezoutOfA = s; - s = tmp - quotient * s; - } + while (r != 0) + { + quotient = gcd / r; // integer division - if (b != 0) - { - bezoutOfB = (gcd - bezoutOfA * a) / b; // integer division - } + tmp = gcd; + gcd = r; + r = tmp - quotient * r; - return new ExtendedEuclideanAlgorithmResult(bezoutOfA, bezoutOfB, gcd); + tmp = bezoutOfA; + bezoutOfA = s; + s = tmp - quotient * s; } - /// - /// Computes the greatest common divisor (gcd) of integers a and b, also the coefficients of Bézout's identity, - /// which are integers x and y such that a*bezoutCoefficientOfA + b*bezoutCoefficientOfB = gcd(a, b). - /// - /// Input number. - /// Second input number. - /// A record of ExtendedEuclideanAlgorithmResult containing the bezout coefficients of a and b as well as the gcd(a,b). - public static ExtendedEuclideanAlgorithmResult Compute(BigInteger a, BigInteger b) + if (b != 0) { - BigInteger quotient; - BigInteger tmp; - var s = BigInteger.Zero; - var bezoutOfA = BigInteger.One; - var r = b; - var gcd = a; - var bezoutOfB = BigInteger.Zero; + bezoutOfB = (gcd - bezoutOfA * a) / b; // integer division + } - while (r != 0) - { - quotient = gcd / r; // integer division + return new ExtendedEuclideanAlgorithmResult(bezoutOfA, bezoutOfB, gcd); + } - tmp = gcd; - gcd = r; - r = tmp - quotient * r; + /// + /// Computes the greatest common divisor (gcd) of integers a and b, also the coefficients of Bézout's identity, + /// which are integers x and y such that a*bezoutCoefficientOfA + b*bezoutCoefficientOfB = gcd(a, b). + /// + /// Input number. + /// Second input number. + /// A record of ExtendedEuclideanAlgorithmResult containing the bezout coefficients of a and b as well as the gcd(a,b). + public static ExtendedEuclideanAlgorithmResult Compute(BigInteger a, BigInteger b) + { + BigInteger quotient; + BigInteger tmp; + var s = BigInteger.Zero; + var bezoutOfA = BigInteger.One; + var r = b; + var gcd = a; + var bezoutOfB = BigInteger.Zero; + + while (r != 0) + { + quotient = gcd / r; // integer division - tmp = bezoutOfA; - bezoutOfA = s; - s = tmp - quotient * s; - } + tmp = gcd; + gcd = r; + r = tmp - quotient * r; - if (b != 0) - { - bezoutOfB = (gcd - bezoutOfA * a) / b; // integer division - } + tmp = bezoutOfA; + bezoutOfA = s; + s = tmp - quotient * s; + } - return new ExtendedEuclideanAlgorithmResult(bezoutOfA, bezoutOfB, gcd); + if (b != 0) + { + bezoutOfB = (gcd - bezoutOfA * a) / b; // integer division } - /// - /// The result type for the computation of the Extended Euclidean Algorithm. - /// - /// The data type of the computation (i.e. long or BigInteger). - /// The bezout coefficient of the parameter a to the computation. - /// The bezout coefficient of the parameter b to the computation. - /// The greatest common divisor of the parameters a and b to the computation. - public record ExtendedEuclideanAlgorithmResult(T bezoutA, T bezoutB, T gcd); + return new ExtendedEuclideanAlgorithmResult(bezoutOfA, bezoutOfB, gcd); } + + /// + /// The result type for the computation of the Extended Euclidean Algorithm. + /// + /// The data type of the computation (i.e. long or BigInteger). + /// The bezout coefficient of the parameter a to the computation. + /// The bezout coefficient of the parameter b to the computation. + /// The greatest common divisor of the parameters a and b to the computation. + public record ExtendedEuclideanAlgorithmResult(T bezoutA, T bezoutB, T gcd); } diff --git a/Algorithms/ModularArithmetic/ModularMultiplicativeInverse.cs b/Algorithms/ModularArithmetic/ModularMultiplicativeInverse.cs index e1ba724c..0e42cd68 100644 --- a/Algorithms/ModularArithmetic/ModularMultiplicativeInverse.cs +++ b/Algorithms/ModularArithmetic/ModularMultiplicativeInverse.cs @@ -1,65 +1,64 @@ using System; using System.Numerics; -namespace Algorithms.ModularArithmetic +namespace Algorithms.ModularArithmetic; + +/// +/// Modular multiplicative inverse: https://en.wikipedia.org/wiki/Modular_multiplicative_inverse. +/// +public static class ModularMultiplicativeInverse { /// - /// Modular multiplicative inverse: https://en.wikipedia.org/wiki/Modular_multiplicative_inverse. + /// Computes the modular multiplicative inverse of a in Z/nZ, if there is any (i.e. if a and n are coprime). /// - public static class ModularMultiplicativeInverse + /// The number a, of which to compute the multiplicative inverse. + /// The modulus n. + /// The multiplicative inverse of a in Z/nZ, a value in the interval [0, n). + /// If there exists no multiplicative inverse of a in Z/nZ. + public static long Compute(long a, long n) { - /// - /// Computes the modular multiplicative inverse of a in Z/nZ, if there is any (i.e. if a and n are coprime). - /// - /// The number a, of which to compute the multiplicative inverse. - /// The modulus n. - /// The multiplicative inverse of a in Z/nZ, a value in the interval [0, n). - /// If there exists no multiplicative inverse of a in Z/nZ. - public static long Compute(long a, long n) - { - var eeaResult = ExtendedEuclideanAlgorithm.Compute(a, n); - - // Check if there is an inverse: - if (eeaResult.gcd != 1) - { - throw new ArithmeticException($"{a} is not invertible in Z/{n}Z."); - } - - // Make sure, inverseOfA (i.e. the bezout coefficient of a) is in the interval [0, n). - var inverseOfA = eeaResult.bezoutA; - if (inverseOfA < 0) - { - inverseOfA += n; - } + var eeaResult = ExtendedEuclideanAlgorithm.Compute(a, n); - return inverseOfA; + // Check if there is an inverse: + if (eeaResult.gcd != 1) + { + throw new ArithmeticException($"{a} is not invertible in Z/{n}Z."); } - /// - /// Computes the modular multiplicative inverse of a in Z/nZ, if there is any (i.e. if a and n are coprime). - /// - /// The number a, of which to compute the multiplicative inverse. - /// The modulus n. - /// The multiplicative inverse of a in Z/nZ, a value in the interval [0, n). - /// If there exists no multiplicative inverse of a in Z/nZ. - public static BigInteger Compute(BigInteger a, BigInteger n) + // Make sure, inverseOfA (i.e. the bezout coefficient of a) is in the interval [0, n). + var inverseOfA = eeaResult.bezoutA; + if (inverseOfA < 0) { - var eeaResult = ExtendedEuclideanAlgorithm.Compute(a, n); + inverseOfA += n; + } - // Check if there is an inverse: - if (eeaResult.gcd != 1) - { - throw new ArithmeticException($"{a} is not invertible in Z/{n}Z."); - } + return inverseOfA; + } - // Make sure, inverseOfA (i.e. the bezout coefficient of a) is in the interval [0, n). - var inverseOfA = eeaResult.bezoutA; - if (inverseOfA < 0) - { - inverseOfA += n; - } + /// + /// Computes the modular multiplicative inverse of a in Z/nZ, if there is any (i.e. if a and n are coprime). + /// + /// The number a, of which to compute the multiplicative inverse. + /// The modulus n. + /// The multiplicative inverse of a in Z/nZ, a value in the interval [0, n). + /// If there exists no multiplicative inverse of a in Z/nZ. + public static BigInteger Compute(BigInteger a, BigInteger n) + { + var eeaResult = ExtendedEuclideanAlgorithm.Compute(a, n); - return inverseOfA; + // Check if there is an inverse: + if (eeaResult.gcd != 1) + { + throw new ArithmeticException($"{a} is not invertible in Z/{n}Z."); } + + // Make sure, inverseOfA (i.e. the bezout coefficient of a) is in the interval [0, n). + var inverseOfA = eeaResult.bezoutA; + if (inverseOfA < 0) + { + inverseOfA += n; + } + + return inverseOfA; } } diff --git a/Algorithms/Numeric/AliquotSumCalculator.cs b/Algorithms/Numeric/AliquotSumCalculator.cs index b2703a3a..b74b550f 100644 --- a/Algorithms/Numeric/AliquotSumCalculator.cs +++ b/Algorithms/Numeric/AliquotSumCalculator.cs @@ -1,38 +1,37 @@ -using System; +using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// In number theory, the aliquot sum s(n) of a positive integer n is the sum of all proper divisors +/// of n, that is, all divisors of n other than n itself. For example, the proper divisors of 15 +/// (that is, the positive divisors of 15 that are not equal to 15) are 1, 3 and 5, so the aliquot +/// sum of 15 is 9 i.e. (1 + 3 + 5). Wikipedia: https://en.wikipedia.org/wiki/Aliquot_sum. +/// +public static class AliquotSumCalculator { /// - /// In number theory, the aliquot sum s(n) of a positive integer n is the sum of all proper divisors - /// of n, that is, all divisors of n other than n itself. For example, the proper divisors of 15 - /// (that is, the positive divisors of 15 that are not equal to 15) are 1, 3 and 5, so the aliquot - /// sum of 15 is 9 i.e. (1 + 3 + 5). Wikipedia: https://en.wikipedia.org/wiki/Aliquot_sum. + /// Finds the aliquot sum of an integer number. /// - public static class AliquotSumCalculator + /// Positive number. + /// The Aliquot Sum. + /// Error number is not on interval (0.0; int.MaxValue). + public static int CalculateAliquotSum(int number) { - /// - /// Finds the aliquot sum of an integer number. - /// - /// Positive number. - /// The Aliquot Sum. - /// Error number is not on interval (0.0; int.MaxValue). - public static int CalculateAliquotSum(int number) + if (number < 0) { - if (number < 0) - { - throw new ArgumentException($"{nameof(number)} cannot be negative"); - } + throw new ArgumentException($"{nameof(number)} cannot be negative"); + } - var sum = 0; - for (int i = 1, limit = number / 2; i <= limit; ++i) + var sum = 0; + for (int i = 1, limit = number / 2; i <= limit; ++i) + { + if (number % i == 0) { - if (number % i == 0) - { - sum += i; - } + sum += i; } - - return sum; } + + return sum; } } diff --git a/Algorithms/Numeric/AmicableNumbersChecker.cs b/Algorithms/Numeric/AmicableNumbersChecker.cs index ad18cd8c..3b80037b 100644 --- a/Algorithms/Numeric/AmicableNumbersChecker.cs +++ b/Algorithms/Numeric/AmicableNumbersChecker.cs @@ -1,41 +1,34 @@ -using System; -using System.Collections.Generic; -using System.Linq; -using System.Text; -using System.Threading.Tasks; +namespace Algorithms.Numeric; -namespace Algorithms.Numeric +/// +/// Amicable numbers are two different natural numbers related in such a way that the sum of the proper divisors of +/// each is equal to the other number. That is, σ(a)=b+a and σ(b)=a+b, where σ(n) is equal to the sum of positive divisors of n (see also divisor function). +/// See here for more info. +/// +public static class AmicableNumbersChecker { /// - /// Amicable numbers are two different natural numbers related in such a way that the sum of the proper divisors of - /// each is equal to the other number. That is, σ(a)=b+a and σ(b)=a+b, where σ(n) is equal to the sum of positive divisors of n (see also divisor function). - /// See here for more info. + /// Checks if two numbers are amicable or not. /// - public static class AmicableNumbersChecker + /// First number to check. + /// Second number to check. + /// True if they are amicable numbers. False if not. + public static bool AreAmicableNumbers(int x, int y) { - /// - /// Checks if two numbers are amicable or not. - /// - /// First number to check. - /// Second number to check. - /// True if they are amicable numbers. False if not. - public static bool AreAmicableNumbers(int x, int y) - { - return SumOfDivisors(x) == y && SumOfDivisors(y) == x; - } + return SumOfDivisors(x) == y && SumOfDivisors(y) == x; + } - private static int SumOfDivisors(int number) + private static int SumOfDivisors(int number) + { + var sum = 0; /* sum of its positive divisors */ + for (var i = 1; i < number; ++i) { - var sum = 0; /* sum of its positive divisors */ - for (var i = 1; i < number; ++i) + if (number % i == 0) { - if (number % i == 0) - { - sum += i; - } + sum += i; } - - return sum; } + + return sum; } } diff --git a/Algorithms/Numeric/AutomorphicNumber.cs b/Algorithms/Numeric/AutomorphicNumber.cs index c3bd9991..159f4988 100644 --- a/Algorithms/Numeric/AutomorphicNumber.cs +++ b/Algorithms/Numeric/AutomorphicNumber.cs @@ -3,72 +3,71 @@ using System.Linq; using System.Numerics; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// Calculates Automorphic numbers. A number is said to be an Automorphic number +/// if the square of the number will contain the number itself at the end. +/// +public static class AutomorphicNumber { /// - /// Calculates Automorphic numbers. A number is said to be an Automorphic number - /// if the square of the number will contain the number itself at the end. + /// Generates a list of automorphic numbers that are between and + /// inclusive. /// - public static class AutomorphicNumber + /// The lower bound of the list. + /// The upper bound of the list. + /// A list that contains all of the automorphic numbers between and inclusive. + /// If the + /// or is not greater than zero + /// or is lower than the . + public static IEnumerable GetAutomorphicNumbers(int lowerBound, int upperBound) { - /// - /// Generates a list of automorphic numbers that are between and - /// inclusive. - /// - /// The lower bound of the list. - /// The upper bound of the list. - /// A list that contains all of the automorphic numbers between and inclusive. - /// If the - /// or is not greater than zero - /// or is lower than the . - public static IEnumerable GetAutomorphicNumbers(int lowerBound, int upperBound) + if (lowerBound < 1) { - if (lowerBound < 1) - { - throw new ArgumentException($"Lower bound must be greater than 0."); - } - - if (upperBound < 1) - { - throw new ArgumentException($"Upper bound must be greater than 0."); - } + throw new ArgumentException($"Lower bound must be greater than 0."); + } - if (lowerBound > upperBound) - { - throw new ArgumentException($"The lower bound must be less than or equal to the upper bound."); - } + if (upperBound < 1) + { + throw new ArgumentException($"Upper bound must be greater than 0."); + } - return Enumerable.Range(lowerBound, upperBound).Where(IsAutomorphic); + if (lowerBound > upperBound) + { + throw new ArgumentException($"The lower bound must be less than or equal to the upper bound."); } - /// - /// Checks if a given natural number is automorphic or not. - /// - /// The number to check. - /// True if the number is automorphic, false otherwise. - /// If the number is non-positive. - public static bool IsAutomorphic(int number) + return Enumerable.Range(lowerBound, upperBound).Where(IsAutomorphic); + } + + /// + /// Checks if a given natural number is automorphic or not. + /// + /// The number to check. + /// True if the number is automorphic, false otherwise. + /// If the number is non-positive. + public static bool IsAutomorphic(int number) + { + if (number < 1) { - if (number < 1) - { - throw new ArgumentException($"An automorphic number must always be positive."); - } + throw new ArgumentException($"An automorphic number must always be positive."); + } - BigInteger square = BigInteger.Pow(number, 2); + BigInteger square = BigInteger.Pow(number, 2); - // Extract the last digits of both numbers - while (number > 0) + // Extract the last digits of both numbers + while (number > 0) + { + if (number % 10 != square % 10) { - if (number % 10 != square % 10) - { - return false; - } - - number /= 10; - square /= 10; + return false; } - return true; + number /= 10; + square /= 10; } + + return true; } } diff --git a/Algorithms/Numeric/BinomialCoefficient.cs b/Algorithms/Numeric/BinomialCoefficient.cs index cb61239c..e2fe6fed 100644 --- a/Algorithms/Numeric/BinomialCoefficient.cs +++ b/Algorithms/Numeric/BinomialCoefficient.cs @@ -1,49 +1,48 @@ using System; using System.Numerics; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// The binomial coefficients are the positive integers +/// that occur as coefficients in the binomial theorem. +/// +public static class BinomialCoefficient { /// - /// The binomial coefficients are the positive integers - /// that occur as coefficients in the binomial theorem. + /// Calculates Binomial coefficients for given input. /// - public static class BinomialCoefficient + /// First number. + /// Second number. + /// Binimial Coefficients. + public static BigInteger Calculate(BigInteger num, BigInteger k) { - /// - /// Calculates Binomial coefficients for given input. - /// - /// First number. - /// Second number. - /// Binimial Coefficients. - public static BigInteger Calculate(BigInteger num, BigInteger k) + if (num < k || k < 0) { - if (num < k || k < 0) - { - throw new ArgumentException("num ≥ k ≥ 0"); - } - - // Tricks to gain performance: - // 1. Because (num over k) equals (num over (num-k)), we can save multiplications and divisions - // by replacing k with the minimum of k and (num - k). - k = BigInteger.Min(k, num - k); + throw new ArgumentException("num ≥ k ≥ 0"); + } - // 2. We can simplify the computation of (num! / (k! * (num - k)!)) to ((num * (num - 1) * ... * (num - k + 1) / (k!)) - // and thus save some multiplications and divisions. - var numerator = BigInteger.One; - for (var val = num - k + 1; val <= num; val++) - { - numerator *= val; - } + // Tricks to gain performance: + // 1. Because (num over k) equals (num over (num-k)), we can save multiplications and divisions + // by replacing k with the minimum of k and (num - k). + k = BigInteger.Min(k, num - k); - // 3. Typically multiplication is a lot faster than division, therefore compute the value of k! first (i.e. k - 1 multiplications) - // and then divide the numerator by the denominator (i.e. 1 division); instead of performing k - 1 divisions (1 for each factor in k!). - var denominator = BigInteger.One; - for (var val = k; val > BigInteger.One; val--) - { - denominator *= val; - } + // 2. We can simplify the computation of (num! / (k! * (num - k)!)) to ((num * (num - 1) * ... * (num - k + 1) / (k!)) + // and thus save some multiplications and divisions. + var numerator = BigInteger.One; + for (var val = num - k + 1; val <= num; val++) + { + numerator *= val; + } - return numerator / denominator; + // 3. Typically multiplication is a lot faster than division, therefore compute the value of k! first (i.e. k - 1 multiplications) + // and then divide the numerator by the denominator (i.e. 1 division); instead of performing k - 1 divisions (1 for each factor in k!). + var denominator = BigInteger.One; + for (var val = k; val > BigInteger.One; val--) + { + denominator *= val; } + + return numerator / denominator; } } diff --git a/Algorithms/Numeric/Decomposition/LU.cs b/Algorithms/Numeric/Decomposition/LU.cs index eb369a85..8b4ccc6f 100644 --- a/Algorithms/Numeric/Decomposition/LU.cs +++ b/Algorithms/Numeric/Decomposition/LU.cs @@ -1,114 +1,113 @@ using System; -namespace Algorithms.Numeric.Decomposition +namespace Algorithms.Numeric.Decomposition; + +/// +/// LU-decomposition factors the "source" matrix as the product of lower triangular matrix +/// and upper triangular matrix. +/// +public static class Lu { /// - /// LU-decomposition factors the "source" matrix as the product of lower triangular matrix - /// and upper triangular matrix. + /// Performs LU-decomposition on "source" matrix. + /// Lower and upper matrices have same shapes as source matrix. + /// Note: Decomposition can be applied only to square matrices. /// - public static class Lu + /// Square matrix to decompose. + /// Tuple of lower and upper matrix. + /// Source matrix is not square shaped. + public static (double[,] L, double[,] U) Decompose(double[,] source) { - /// - /// Performs LU-decomposition on "source" matrix. - /// Lower and upper matrices have same shapes as source matrix. - /// Note: Decomposition can be applied only to square matrices. - /// - /// Square matrix to decompose. - /// Tuple of lower and upper matrix. - /// Source matrix is not square shaped. - public static (double[,] L, double[,] U) Decompose(double[,] source) + if (source.GetLength(0) != source.GetLength(1)) { - if (source.GetLength(0) != source.GetLength(1)) - { - throw new ArgumentException("Source matrix is not square shaped."); - } + throw new ArgumentException("Source matrix is not square shaped."); + } - var pivot = source.GetLength(0); - var lower = new double[pivot, pivot]; - var upper = new double[pivot, pivot]; + var pivot = source.GetLength(0); + var lower = new double[pivot, pivot]; + var upper = new double[pivot, pivot]; - for (var i = 0; i < pivot; i++) + for (var i = 0; i < pivot; i++) + { + for (var k = i; k < pivot; k++) { - for (var k = i; k < pivot; k++) + double sum = 0; + + for (var j = 0; j < i; j++) { - double sum = 0; + sum += lower[i, j] * upper[j, k]; + } - for (var j = 0; j < i; j++) - { - sum += lower[i, j] * upper[j, k]; - } + upper[i, k] = source[i, k] - sum; + } - upper[i, k] = source[i, k] - sum; + for (var k = i; k < pivot; k++) + { + if (i == k) + { + lower[i, i] = 1; } - - for (var k = i; k < pivot; k++) + else { - if (i == k) + double sum = 0; + + for (var j = 0; j < i; j++) { - lower[i, i] = 1; + sum += lower[k, j] * upper[j, i]; } - else - { - double sum = 0; - - for (var j = 0; j < i; j++) - { - sum += lower[k, j] * upper[j, i]; - } - lower[k, i] = (source[k, i] - sum) / upper[i, i]; - } + lower[k, i] = (source[k, i] - sum) / upper[i, i]; } } - - return (L: lower, U: upper); } - /// - /// Eliminates linear equations system represented as A*x=b, using LU-decomposition, - /// where A - matrix of equation coefficients, b - vector of absolute terms of equations. - /// - /// Matrix of equation coefficients. - /// Vector of absolute terms of equations. - /// Vector-solution for linear equations system. - /// Matrix of equation coefficients is not square shaped. - public static double[] Eliminate(double[,] matrix, double[] coefficients) - { - if (matrix.GetLength(0) != matrix.GetLength(1)) - { - throw new ArgumentException("Matrix of equation coefficients is not square shaped."); - } + return (L: lower, U: upper); + } - var pivot = matrix.GetLength(0); - var upperTransform = new double[pivot, 1]; // U * upperTransform = coefficients - var solution = new double[pivot]; // L * solution = upperTransform - (double[,] l, double[,] u) = Decompose(matrix); + /// + /// Eliminates linear equations system represented as A*x=b, using LU-decomposition, + /// where A - matrix of equation coefficients, b - vector of absolute terms of equations. + /// + /// Matrix of equation coefficients. + /// Vector of absolute terms of equations. + /// Vector-solution for linear equations system. + /// Matrix of equation coefficients is not square shaped. + public static double[] Eliminate(double[,] matrix, double[] coefficients) + { + if (matrix.GetLength(0) != matrix.GetLength(1)) + { + throw new ArgumentException("Matrix of equation coefficients is not square shaped."); + } - for (var i = 0; i < pivot; i++) - { - double pivotPointSum = 0; + var pivot = matrix.GetLength(0); + var upperTransform = new double[pivot, 1]; // U * upperTransform = coefficients + var solution = new double[pivot]; // L * solution = upperTransform + (double[,] l, double[,] u) = Decompose(matrix); - for (var j = 0; j < i; j++) - { - pivotPointSum += upperTransform[j, 0] * l[i, j]; - } + for (var i = 0; i < pivot; i++) + { + double pivotPointSum = 0; - upperTransform[i, 0] = (coefficients[i] - pivotPointSum) / l[i, i]; + for (var j = 0; j < i; j++) + { + pivotPointSum += upperTransform[j, 0] * l[i, j]; } - for (var i = pivot - 1; i >= 0; i--) - { - double pivotPointSum = 0; + upperTransform[i, 0] = (coefficients[i] - pivotPointSum) / l[i, i]; + } - for (var j = i; j < pivot; j++) - { - pivotPointSum += solution[j] * u[i, j]; - } + for (var i = pivot - 1; i >= 0; i--) + { + double pivotPointSum = 0; - solution[i] = (upperTransform[i, 0] - pivotPointSum) / u[i, i]; + for (var j = i; j < pivot; j++) + { + pivotPointSum += solution[j] * u[i, j]; } - return solution; + solution[i] = (upperTransform[i, 0] - pivotPointSum) / u[i, i]; } + + return solution; } } diff --git a/Algorithms/Numeric/Decomposition/ThinSVD.cs b/Algorithms/Numeric/Decomposition/ThinSVD.cs index 4af41275..8c07f2e6 100644 --- a/Algorithms/Numeric/Decomposition/ThinSVD.cs +++ b/Algorithms/Numeric/Decomposition/ThinSVD.cs @@ -1,145 +1,142 @@ using System; using Utilities.Extensions; -using M = Utilities.Extensions.MatrixExtensions; -using V = Utilities.Extensions.VectorExtensions; -namespace Algorithms.Numeric.Decomposition +namespace Algorithms.Numeric.Decomposition; + +/// +/// Singular Vector Decomposition decomposes any general matrix into its +/// singular values and a set of orthonormal bases. +/// +public static class ThinSvd { /// - /// Singular Vector Decomposition decomposes any general matrix into its - /// singular values and a set of orthonormal bases. + /// Computes a random unit vector. /// - public static class ThinSvd + /// The dimensions of the required vector. + /// The unit vector. + public static double[] RandomUnitVector(int dimensions) { - /// - /// Computes a random unit vector. - /// - /// The dimensions of the required vector. - /// The unit vector. - public static double[] RandomUnitVector(int dimensions) + Random random = new(); + double[] result = new double[dimensions]; + for (var i = 0; i < dimensions; i++) { - Random random = new(); - double[] result = new double[dimensions]; - for (var i = 0; i < dimensions; i++) - { - result[i] = 2 * random.NextDouble() - 1; - } - - var magnitude = result.Magnitude(); - result = result.Scale(1 / magnitude); - return result; + result[i] = 2 * random.NextDouble() - 1; } - /// - /// Computes a single singular vector for the given matrix, corresponding to the largest singular value. - /// - /// The matrix. - /// A singular vector, with dimension equal to number of columns of the matrix. - public static double[] Decompose1D(double[,] matrix) => - Decompose1D(matrix, 1E-5, 100); - - /// - /// Computes a single singular vector for the given matrix, corresponding to the largest singular value. - /// - /// The matrix. - /// The error margin. - /// The maximum number of iterations. - /// A singular vector, with dimension equal to number of columns of the matrix. - public static double[] Decompose1D(double[,] matrix, double epsilon, int maxIterations) + var magnitude = result.Magnitude(); + result = result.Scale(1 / magnitude); + return result; + } + + /// + /// Computes a single singular vector for the given matrix, corresponding to the largest singular value. + /// + /// The matrix. + /// A singular vector, with dimension equal to number of columns of the matrix. + public static double[] Decompose1D(double[,] matrix) => + Decompose1D(matrix, 1E-5, 100); + + /// + /// Computes a single singular vector for the given matrix, corresponding to the largest singular value. + /// + /// The matrix. + /// The error margin. + /// The maximum number of iterations. + /// A singular vector, with dimension equal to number of columns of the matrix. + public static double[] Decompose1D(double[,] matrix, double epsilon, int maxIterations) + { + var n = matrix.GetLength(1); + var iterations = 0; + double mag; + double[] lastIteration; + double[] currIteration = RandomUnitVector(n); + double[,] b = matrix.Transpose().Multiply(matrix); + do { - var n = matrix.GetLength(1); - var iterations = 0; - double mag; - double[] lastIteration; - double[] currIteration = RandomUnitVector(n); - double[,] b = matrix.Transpose().Multiply(matrix); - do + lastIteration = currIteration.Copy(); + currIteration = b.MultiplyVector(lastIteration); + currIteration = currIteration.Scale(100); + mag = currIteration.Magnitude(); + if (mag > epsilon) { - lastIteration = currIteration.Copy(); - currIteration = b.MultiplyVector(lastIteration); - currIteration = currIteration.Scale(100); - mag = currIteration.Magnitude(); - if (mag > epsilon) - { - currIteration = currIteration.Scale(1 / mag); - } - - iterations++; + currIteration = currIteration.Scale(1 / mag); } - while (lastIteration.Dot(currIteration) < 1 - epsilon && iterations < maxIterations); - return currIteration; + iterations++; } + while (lastIteration.Dot(currIteration) < 1 - epsilon && iterations < maxIterations); - public static (double[,] U, double[] S, double[,] V) Decompose(double[,] matrix) => - Decompose(matrix, 1E-5, 100); - - /// - /// Computes the SVD for the given matrix, with singular values arranged from greatest to least. - /// - /// The matrix. - /// The error margin. - /// The maximum number of iterations. - /// The SVD. - public static (double[,] U, double[] S, double[,] V) Decompose( - double[,] matrix, - double epsilon, - int maxIterations) + return currIteration; + } + + public static (double[,] U, double[] S, double[,] V) Decompose(double[,] matrix) => + Decompose(matrix, 1E-5, 100); + + /// + /// Computes the SVD for the given matrix, with singular values arranged from greatest to least. + /// + /// The matrix. + /// The error margin. + /// The maximum number of iterations. + /// The SVD. + public static (double[,] U, double[] S, double[,] V) Decompose( + double[,] matrix, + double epsilon, + int maxIterations) + { + var m = matrix.GetLength(0); + var n = matrix.GetLength(1); + var numValues = Math.Min(m, n); + + // sigmas is be a diagonal matrix, hence only a vector is needed + double[] sigmas = new double[numValues]; + double[,] us = new double[m, numValues]; + double[,] vs = new double[n, numValues]; + + // keep track of progress + double[,] remaining = matrix.Copy(); + + // for each singular value + for (var i = 0; i < numValues; i++) { - var m = matrix.GetLength(0); - var n = matrix.GetLength(1); - var numValues = Math.Min(m, n); + // compute the v singular vector + double[] v = Decompose1D(remaining, epsilon, maxIterations); + double[] u = matrix.MultiplyVector(v); - // sigmas is be a diagonal matrix, hence only a vector is needed - double[] sigmas = new double[numValues]; - double[,] us = new double[m, numValues]; - double[,] vs = new double[n, numValues]; + // compute the contribution of this pair of singular vectors + double[,] contrib = u.OuterProduct(v); - // keep track of progress - double[,] remaining = matrix.Copy(); + // extract the singular value + var s = u.Magnitude(); - // for each singular value - for (var i = 0; i < numValues; i++) + // v and u should be unit vectors + if (s < epsilon) + { + u = new double[m]; + v = new double[n]; + } + else { - // compute the v singular vector - double[] v = Decompose1D(remaining, epsilon, maxIterations); - double[] u = matrix.MultiplyVector(v); - - // compute the contribution of this pair of singular vectors - double[,] contrib = u.OuterProduct(v); - - // extract the singular value - var s = u.Magnitude(); - - // v and u should be unit vectors - if (s < epsilon) - { - u = new double[m]; - v = new double[n]; - } - else - { - u = u.Scale(1 / s); - } - - // save u, v and s into the result - for (var j = 0; j < u.Length; j++) - { - us[j, i] = u[j]; - } - - for (var j = 0; j < v.Length; j++) - { - vs[j, i] = v[j]; - } - - sigmas[i] = s; - - // remove the contribution of this pair and compute the rest - remaining = remaining.Subtract(contrib); + u = u.Scale(1 / s); } - return (U: us, S: sigmas, V: vs); + // save u, v and s into the result + for (var j = 0; j < u.Length; j++) + { + us[j, i] = u[j]; + } + + for (var j = 0; j < v.Length; j++) + { + vs[j, i] = v[j]; + } + + sigmas[i] = s; + + // remove the contribution of this pair and compute the rest + remaining = remaining.Subtract(contrib); } + + return (U: us, S: sigmas, V: vs); } } diff --git a/Algorithms/Numeric/EulerMethod.cs b/Algorithms/Numeric/EulerMethod.cs index 043fd26f..20112c2b 100644 --- a/Algorithms/Numeric/EulerMethod.cs +++ b/Algorithms/Numeric/EulerMethod.cs @@ -1,85 +1,84 @@ using System; using System.Collections.Generic; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// In mathematics and computational science, the Euler method (also called forward Euler method) +/// is a first-order numerical procedure for solving ordinary differential equations (ODEs) +/// with a given initial value (aka. Cauchy problem). It is the most basic explicit method for numerical integration +/// of ordinary differential equations. The method proceeds in a series of steps. At each step +/// the y-value is calculated by evaluating the differential equation at the previous step, +/// multiplying the result with the step-size and adding it to the last y-value: +/// y_n+1 = y_n + stepSize * f(x_n, y_n). +/// (description adapted from https://en.wikipedia.org/wiki/Euler_method ) +/// (see also: https://www.geeksforgeeks.org/euler-method-solving-differential-equation/ ). +/// +public static class EulerMethod { /// - /// In mathematics and computational science, the Euler method (also called forward Euler method) - /// is a first-order numerical procedure for solving ordinary differential equations (ODEs) - /// with a given initial value (aka. Cauchy problem). It is the most basic explicit method for numerical integration - /// of ordinary differential equations. The method proceeds in a series of steps. At each step - /// the y-value is calculated by evaluating the differential equation at the previous step, - /// multiplying the result with the step-size and adding it to the last y-value: - /// y_n+1 = y_n + stepSize * f(x_n, y_n). - /// (description adapted from https://en.wikipedia.org/wiki/Euler_method ) - /// (see also: https://www.geeksforgeeks.org/euler-method-solving-differential-equation/ ). + /// Loops through all the steps until xEnd is reached, adds a point for each step and then + /// returns all the points. /// - public static class EulerMethod + /// Initial conditions x-value. + /// Last x-value. + /// Step-size on the x-axis. + /// Initial conditions y-value. + /// The right hand side of the differential equation. + /// The solution of the Cauchy problem. + public static List EulerFull( + double xStart, + double xEnd, + double stepSize, + double yStart, + Func yDerivative) { - /// - /// Loops through all the steps until xEnd is reached, adds a point for each step and then - /// returns all the points. - /// - /// Initial conditions x-value. - /// Last x-value. - /// Step-size on the x-axis. - /// Initial conditions y-value. - /// The right hand side of the differential equation. - /// The solution of the Cauchy problem. - public static List EulerFull( - double xStart, - double xEnd, - double stepSize, - double yStart, - Func yDerivative) + if (xStart >= xEnd) { - if (xStart >= xEnd) - { - throw new ArgumentOutOfRangeException( - nameof(xEnd), - $"{nameof(xEnd)} should be greater than {nameof(xStart)}"); - } - - if (stepSize <= 0) - { - throw new ArgumentOutOfRangeException( - nameof(stepSize), - $"{nameof(stepSize)} should be greater than zero"); - } - - List points = new(); - double[] firstPoint = { xStart, yStart }; - points.Add(firstPoint); - var yCurrent = yStart; - var xCurrent = xStart; - - while (xCurrent < xEnd) - { - yCurrent = EulerStep(xCurrent, stepSize, yCurrent, yDerivative); - xCurrent += stepSize; - double[] point = { xCurrent, yCurrent }; - points.Add(point); - } + throw new ArgumentOutOfRangeException( + nameof(xEnd), + $"{nameof(xEnd)} should be greater than {nameof(xStart)}"); + } - return points; + if (stepSize <= 0) + { + throw new ArgumentOutOfRangeException( + nameof(stepSize), + $"{nameof(stepSize)} should be greater than zero"); } - /// - /// Calculates the next y-value based on the current value of x, y and the stepSize. - /// - /// Current x-value. - /// Step-size on the x-axis. - /// Current y-value. - /// The right hand side of the differential equation. - /// The next y-value. - private static double EulerStep( - double xCurrent, - double stepSize, - double yCurrent, - Func yDerivative) + List points = new(); + double[] firstPoint = { xStart, yStart }; + points.Add(firstPoint); + var yCurrent = yStart; + var xCurrent = xStart; + + while (xCurrent < xEnd) { - var yNext = yCurrent + stepSize * yDerivative(xCurrent, yCurrent); - return yNext; + yCurrent = EulerStep(xCurrent, stepSize, yCurrent, yDerivative); + xCurrent += stepSize; + double[] point = { xCurrent, yCurrent }; + points.Add(point); } + + return points; + } + + /// + /// Calculates the next y-value based on the current value of x, y and the stepSize. + /// + /// Current x-value. + /// Step-size on the x-axis. + /// Current y-value. + /// The right hand side of the differential equation. + /// The next y-value. + private static double EulerStep( + double xCurrent, + double stepSize, + double yCurrent, + Func yDerivative) + { + var yNext = yCurrent + stepSize * yDerivative(xCurrent, yCurrent); + return yNext; } } diff --git a/Algorithms/Numeric/Factorial.cs b/Algorithms/Numeric/Factorial.cs index 845a92ea..6c30f50e 100644 --- a/Algorithms/Numeric/Factorial.cs +++ b/Algorithms/Numeric/Factorial.cs @@ -1,39 +1,38 @@ -using System; +using System; using System.Numerics; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// The factorial of a positive integer n, denoted by n!, +/// is the product of all positive integers less than or equal to n. +/// +public static class Factorial { /// - /// The factorial of a positive integer n, denoted by n!, - /// is the product of all positive integers less than or equal to n. + /// Calculates factorial of a integer number. /// - public static class Factorial + /// Integer Input number. + /// Factorial of integer input number. + public static BigInteger Calculate(int inputNum) { - /// - /// Calculates factorial of a integer number. - /// - /// Integer Input number. - /// Factorial of integer input number. - public static BigInteger Calculate(int inputNum) - { - // Convert integer input to BigInteger - BigInteger num = new BigInteger(inputNum); + // Convert integer input to BigInteger + BigInteger num = new BigInteger(inputNum); - // Don't calculate factorial if input is a negative number. - if (BigInteger.Compare(num, BigInteger.Zero) < 0) - { - throw new ArgumentException("Only for num >= 0"); - } - - // Factorial of numbers greater than 0. - BigInteger result = BigInteger.One; + // Don't calculate factorial if input is a negative number. + if (BigInteger.Compare(num, BigInteger.Zero) < 0) + { + throw new ArgumentException("Only for num >= 0"); + } - for (BigInteger i = BigInteger.One; BigInteger.Compare(i, num) <= 0; i = BigInteger.Add(i, BigInteger.One)) - { - result = BigInteger.Multiply(result, i); - } + // Factorial of numbers greater than 0. + BigInteger result = BigInteger.One; - return result; + for (BigInteger i = BigInteger.One; BigInteger.Compare(i, num) <= 0; i = BigInteger.Add(i, BigInteger.One)) + { + result = BigInteger.Multiply(result, i); } + + return result; } } diff --git a/Algorithms/Numeric/Factorization/IFactorizer.cs b/Algorithms/Numeric/Factorization/IFactorizer.cs index dcce13d6..4a45f86e 100755 --- a/Algorithms/Numeric/Factorization/IFactorizer.cs +++ b/Algorithms/Numeric/Factorization/IFactorizer.cs @@ -1,16 +1,15 @@ -namespace Algorithms.Numeric.Factorization +namespace Algorithms.Numeric.Factorization; + +/// +/// Finds a factor of a given number or returns false if it's prime. +/// +public interface IFactorizer { /// /// Finds a factor of a given number or returns false if it's prime. /// - public interface IFactorizer - { - /// - /// Finds a factor of a given number or returns false if it's prime. - /// - /// Integer to factor. - /// Found factor. - /// if factor is found, if is prime. - bool TryFactor(int n, out int factor); - } + /// Integer to factor. + /// Found factor. + /// if factor is found, if is prime. + bool TryFactor(int n, out int factor); } diff --git a/Algorithms/Numeric/Factorization/TrialDivisionFactorizer.cs b/Algorithms/Numeric/Factorization/TrialDivisionFactorizer.cs index 43820e24..f105fd6d 100755 --- a/Algorithms/Numeric/Factorization/TrialDivisionFactorizer.cs +++ b/Algorithms/Numeric/Factorization/TrialDivisionFactorizer.cs @@ -1,24 +1,23 @@ using System; using System.Linq; -namespace Algorithms.Numeric.Factorization +namespace Algorithms.Numeric.Factorization; + +/// +/// Factors number using trial division algorithm. +/// +public class TrialDivisionFactorizer : IFactorizer { /// - /// Factors number using trial division algorithm. + /// Finds the smallest non trivial factor (i.e.: 1 < factor <= sqrt()) of a given number or returns false if it's prime. /// - public class TrialDivisionFactorizer : IFactorizer + /// Integer to factor. + /// Found factor. + /// if factor is found, if is prime. + public bool TryFactor(int n, out int factor) { - /// - /// Finds the smallest non trivial factor (i.e.: 1 < factor <= sqrt()) of a given number or returns false if it's prime. - /// - /// Integer to factor. - /// Found factor. - /// if factor is found, if is prime. - public bool TryFactor(int n, out int factor) - { - n = Math.Abs(n); - factor = Enumerable.Range(2, (int)Math.Sqrt(n) - 1).FirstOrDefault(i => n % i == 0); - return factor != 0; - } + n = Math.Abs(n); + factor = Enumerable.Range(2, (int)Math.Sqrt(n) - 1).FirstOrDefault(i => n % i == 0); + return factor != 0; } } diff --git a/Algorithms/Numeric/GaussJordanElimination.cs b/Algorithms/Numeric/GaussJordanElimination.cs index 27a54dd8..3ef7969b 100644 --- a/Algorithms/Numeric/GaussJordanElimination.cs +++ b/Algorithms/Numeric/GaussJordanElimination.cs @@ -1,153 +1,152 @@ -using System; +using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// Algorithm used to find the inverse of any matrix that can be inverted. +/// +public class GaussJordanElimination { + private int RowCount { get; set; } + /// - /// Algorithm used to find the inverse of any matrix that can be inverted. + /// Method to find a linear equation system using gaussian elimination. /// - public class GaussJordanElimination + /// The key matrix to solve via algorithm. + /// + /// whether the input matrix has a unique solution or not. + /// and solves on the given matrix. + /// + public bool Solve(double[,] matrix) { - private int RowCount { get; set; } - - /// - /// Method to find a linear equation system using gaussian elimination. - /// - /// The key matrix to solve via algorithm. - /// - /// whether the input matrix has a unique solution or not. - /// and solves on the given matrix. - /// - public bool Solve(double[,] matrix) + RowCount = matrix.GetUpperBound(0) + 1; + + if (!CanMatrixBeUsed(matrix)) { - RowCount = matrix.GetUpperBound(0) + 1; + throw new ArgumentException("Please use a n*(n+1) matrix with Length > 0."); + } - if (!CanMatrixBeUsed(matrix)) - { - throw new ArgumentException("Please use a n*(n+1) matrix with Length > 0."); - } + var pivot = PivotMatrix(ref matrix); + if (!pivot) + { + return false; + } - var pivot = PivotMatrix(ref matrix); - if (!pivot) - { - return false; - } + Elimination(ref matrix); - Elimination(ref matrix); + return ElementaryReduction(ref matrix); + } - return ElementaryReduction(ref matrix); - } + /// + /// To make simple validation of the matrix to be used. + /// + /// Multidimensional array matrix. + /// + /// True: if algorithm can be use for given matrix; + /// False: Otherwise. + /// + private bool CanMatrixBeUsed(double[,] matrix) => matrix?.Length == RowCount * (RowCount + 1) && RowCount > 1; - /// - /// To make simple validation of the matrix to be used. - /// - /// Multidimensional array matrix. - /// - /// True: if algorithm can be use for given matrix; - /// False: Otherwise. - /// - private bool CanMatrixBeUsed(double[,] matrix) => matrix?.Length == RowCount * (RowCount + 1) && RowCount > 1; - - /// - /// To prepare given matrix by pivoting rows. - /// - /// Input matrix. - /// Matrix. - private bool PivotMatrix(ref double[,] matrix) + /// + /// To prepare given matrix by pivoting rows. + /// + /// Input matrix. + /// Matrix. + private bool PivotMatrix(ref double[,] matrix) + { + for (var col = 0; col + 1 < RowCount; col++) { - for (var col = 0; col + 1 < RowCount; col++) + if (matrix[col, col] == 0) { - if (matrix[col, col] == 0) - { - // To find a non-zero coefficient - var rowToSwap = FindNonZeroCoefficient(ref matrix, col); + // To find a non-zero coefficient + var rowToSwap = FindNonZeroCoefficient(ref matrix, col); - if (matrix[rowToSwap, col] != 0) - { - var tmp = new double[RowCount + 1]; - for (var i = 0; i < RowCount + 1; i++) - { - // To make the swap with the element above. - tmp[i] = matrix[rowToSwap, i]; - matrix[rowToSwap, i] = matrix[col, i]; - matrix[col, i] = tmp[i]; - } - } - else + if (matrix[rowToSwap, col] != 0) + { + var tmp = new double[RowCount + 1]; + for (var i = 0; i < RowCount + 1; i++) { - // To return that the matrix doesn't have a unique solution. - return false; + // To make the swap with the element above. + tmp[i] = matrix[rowToSwap, i]; + matrix[rowToSwap, i] = matrix[col, i]; + matrix[col, i] = tmp[i]; } } + else + { + // To return that the matrix doesn't have a unique solution. + return false; + } } - - return true; } - private int FindNonZeroCoefficient(ref double[,] matrix, int col) - { - var rowToSwap = col + 1; + return true; + } - // To find a non-zero coefficient - for (; rowToSwap < RowCount; rowToSwap++) + private int FindNonZeroCoefficient(ref double[,] matrix, int col) + { + var rowToSwap = col + 1; + + // To find a non-zero coefficient + for (; rowToSwap < RowCount; rowToSwap++) + { + if (matrix[rowToSwap, col] != 0) { - if (matrix[rowToSwap, col] != 0) - { - return rowToSwap; - } + return rowToSwap; } - - return col + 1; } - /// - /// Applies REF. - /// - /// Input matrix. - private void Elimination(ref double[,] matrix) + return col + 1; + } + + /// + /// Applies REF. + /// + /// Input matrix. + private void Elimination(ref double[,] matrix) + { + for (var srcRow = 0; srcRow + 1 < RowCount; srcRow++) { - for (var srcRow = 0; srcRow + 1 < RowCount; srcRow++) + for (var destRow = srcRow + 1; destRow < RowCount; destRow++) { - for (var destRow = srcRow + 1; destRow < RowCount; destRow++) - { - var df = matrix[srcRow, srcRow]; - var sf = matrix[destRow, srcRow]; + var df = matrix[srcRow, srcRow]; + var sf = matrix[destRow, srcRow]; - for (var i = 0; i < RowCount + 1; i++) - { - matrix[destRow, i] = matrix[destRow, i] * df - matrix[srcRow, i] * sf; - } + for (var i = 0; i < RowCount + 1; i++) + { + matrix[destRow, i] = matrix[destRow, i] * df - matrix[srcRow, i] * sf; } } } + } - /// - /// To continue reducing the matrix using RREF. - /// - /// Input matrix. - /// True if it has a unique solution; false otherwise. - private bool ElementaryReduction(ref double[,] matrix) + /// + /// To continue reducing the matrix using RREF. + /// + /// Input matrix. + /// True if it has a unique solution; false otherwise. + private bool ElementaryReduction(ref double[,] matrix) + { + for (var row = RowCount - 1; row >= 0; row--) { - for (var row = RowCount - 1; row >= 0; row--) + var element = matrix[row, row]; + if (element == 0) { - var element = matrix[row, row]; - if (element == 0) - { - return false; - } - - for (var i = 0; i < RowCount + 1; i++) - { - matrix[row, i] /= element; - } + return false; + } - for (var destRow = 0; destRow < row; destRow++) - { - matrix[destRow, RowCount] -= matrix[destRow, row] * matrix[row, RowCount]; - matrix[destRow, row] = 0; - } + for (var i = 0; i < RowCount + 1; i++) + { + matrix[row, i] /= element; } - return true; + for (var destRow = 0; destRow < row; destRow++) + { + matrix[destRow, RowCount] -= matrix[destRow, row] * matrix[row, RowCount]; + matrix[destRow, row] = 0; + } } + + return true; } } diff --git a/Algorithms/Numeric/GreatestCommonDivisor/BinaryGreatestCommonDivisorFinder.cs b/Algorithms/Numeric/GreatestCommonDivisor/BinaryGreatestCommonDivisorFinder.cs index 935488d2..691a855a 100644 --- a/Algorithms/Numeric/GreatestCommonDivisor/BinaryGreatestCommonDivisorFinder.cs +++ b/Algorithms/Numeric/GreatestCommonDivisor/BinaryGreatestCommonDivisorFinder.cs @@ -1,71 +1,70 @@ using System; -namespace Algorithms.Numeric.GreatestCommonDivisor +namespace Algorithms.Numeric.GreatestCommonDivisor; + +/// +/// Finds greatest common divisor for numbers u and v +/// using binary algorithm. +/// Wiki: https://en.wikipedia.org/wiki/Binary_GCD_algorithm. +/// +public class BinaryGreatestCommonDivisorFinder : IGreatestCommonDivisorFinder { - /// - /// Finds greatest common divisor for numbers u and v - /// using binary algorithm. - /// Wiki: https://en.wikipedia.org/wiki/Binary_GCD_algorithm. - /// - public class BinaryGreatestCommonDivisorFinder : IGreatestCommonDivisorFinder + public int FindGcd(int u, int v) { - public int FindGcd(int u, int v) + // GCD(0, 0) = 0 + if (u == 0 && v == 0) { - // GCD(0, 0) = 0 - if (u == 0 && v == 0) - { - return 0; - } + return 0; + } - // GCD(0, v) = v; GCD(u, 0) = u - if (u == 0 || v == 0) - { - return u + v; - } + // GCD(0, v) = v; GCD(u, 0) = u + if (u == 0 || v == 0) + { + return u + v; + } - // GCD(-a, -b) = GCD(-a, b) = GCD(a, -b) = GCD(a, b) - u = Math.Sign(u) * u; - v = Math.Sign(v) * v; + // GCD(-a, -b) = GCD(-a, b) = GCD(a, -b) = GCD(a, b) + u = Math.Sign(u) * u; + v = Math.Sign(v) * v; - // Let shift := lg K, where K is the greatest power of 2 dividing both u and v - var shift = 0; - while (((u | v) & 1) == 0) - { - u >>= 1; - v >>= 1; - shift++; - } + // Let shift := lg K, where K is the greatest power of 2 dividing both u and v + var shift = 0; + while (((u | v) & 1) == 0) + { + u >>= 1; + v >>= 1; + shift++; + } + + while ((u & 1) == 0) + { + u >>= 1; + } - while ((u & 1) == 0) + // From here on, u is always odd + do + { + // Remove all factors of 2 in v as they are not common + // v is not zero, so while will terminate + while ((v & 1) == 0) { - u >>= 1; + v >>= 1; } - // From here on, u is always odd - do + // Now u and v are both odd. Swap if necessary so u <= v, + if (u > v) { - // Remove all factors of 2 in v as they are not common - // v is not zero, so while will terminate - while ((v & 1) == 0) - { - v >>= 1; - } - - // Now u and v are both odd. Swap if necessary so u <= v, - if (u > v) - { - var t = v; - v = u; - u = t; - } - - // Here v >= u and v - u is even - v -= u; + var t = v; + v = u; + u = t; } - while (v != 0); - // Restore common factors of 2 - return u << shift; + // Here v >= u and v - u is even + v -= u; } + while (v != 0); + + // Restore common factors of 2 + return u << shift; } } diff --git a/Algorithms/Numeric/GreatestCommonDivisor/EuclideanGreatestCommonDivisorFinder.cs b/Algorithms/Numeric/GreatestCommonDivisor/EuclideanGreatestCommonDivisorFinder.cs index fdbdafdf..df93caf9 100644 --- a/Algorithms/Numeric/GreatestCommonDivisor/EuclideanGreatestCommonDivisorFinder.cs +++ b/Algorithms/Numeric/GreatestCommonDivisor/EuclideanGreatestCommonDivisorFinder.cs @@ -1,41 +1,40 @@ -namespace Algorithms.Numeric.GreatestCommonDivisor +namespace Algorithms.Numeric.GreatestCommonDivisor; + +/// +/// TODO. +/// +public class EuclideanGreatestCommonDivisorFinder : IGreatestCommonDivisorFinder { /// - /// TODO. + /// Finds greatest common divisor for numbers a and b + /// using euclidean algorithm. /// - public class EuclideanGreatestCommonDivisorFinder : IGreatestCommonDivisorFinder + /// TODO. + /// TODO. 2. + /// Greatest common divisor. + public int FindGcd(int a, int b) { - /// - /// Finds greatest common divisor for numbers a and b - /// using euclidean algorithm. - /// - /// TODO. - /// TODO. 2. - /// Greatest common divisor. - public int FindGcd(int a, int b) + if (a == 0 && b == 0) { - if (a == 0 && b == 0) - { - return int.MaxValue; - } - - if (a == 0 || b == 0) - { - return a + b; - } + return int.MaxValue; + } - var aa = a; - var bb = b; - var cc = aa % bb; + if (a == 0 || b == 0) + { + return a + b; + } - while (cc != 0) - { - aa = bb; - bb = cc; - cc = aa % bb; - } + var aa = a; + var bb = b; + var cc = aa % bb; - return bb; + while (cc != 0) + { + aa = bb; + bb = cc; + cc = aa % bb; } + + return bb; } } diff --git a/Algorithms/Numeric/GreatestCommonDivisor/IGreatestCommonDivisorFinder.cs b/Algorithms/Numeric/GreatestCommonDivisor/IGreatestCommonDivisorFinder.cs index f0e40e61..49fa8ed5 100644 --- a/Algorithms/Numeric/GreatestCommonDivisor/IGreatestCommonDivisorFinder.cs +++ b/Algorithms/Numeric/GreatestCommonDivisor/IGreatestCommonDivisorFinder.cs @@ -1,7 +1,6 @@ -namespace Algorithms.Numeric.GreatestCommonDivisor +namespace Algorithms.Numeric.GreatestCommonDivisor; + +public interface IGreatestCommonDivisorFinder { - public interface IGreatestCommonDivisorFinder - { - int FindGcd(int a, int b); - } + int FindGcd(int a, int b); } diff --git a/Algorithms/Numeric/KeithNumberChecker.cs b/Algorithms/Numeric/KeithNumberChecker.cs index 92be23a5..e1a91f93 100644 --- a/Algorithms/Numeric/KeithNumberChecker.cs +++ b/Algorithms/Numeric/KeithNumberChecker.cs @@ -1,57 +1,56 @@ using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// In number theory, a Keith number or repfigit number is a natural number n in a given number base b with k digits such that +/// when a sequence is created such that the first k terms are the k digits of n and each subsequent term is the sum of the +/// previous k terms, n is part of the sequence. +/// +public static class KeithNumberChecker { /// - /// In number theory, a Keith number or repfigit number is a natural number n in a given number base b with k digits such that - /// when a sequence is created such that the first k terms are the k digits of n and each subsequent term is the sum of the - /// previous k terms, n is part of the sequence. + /// Checks if a number is a Keith number or not. /// - public static class KeithNumberChecker + /// Number to check. + /// True if it is a Keith number; False otherwise. + public static bool IsKeithNumber(int number) { - /// - /// Checks if a number is a Keith number or not. - /// - /// Number to check. - /// True if it is a Keith number; False otherwise. - public static bool IsKeithNumber(int number) + if (number < 0) { - if (number < 0) - { - throw new ArgumentException($"{nameof(number)} cannot be negative"); - } - - var tempNumber = number; + throw new ArgumentException($"{nameof(number)} cannot be negative"); + } - var stringNumber = number.ToString(); + var tempNumber = number; - var digitsInNumber = stringNumber.Length; + var stringNumber = number.ToString(); - /* storing the terms of the series */ - var termsArray = new int[number]; + var digitsInNumber = stringNumber.Length; - for (var i = digitsInNumber - 1; i >= 0; i--) - { - termsArray[i] = tempNumber % 10; - tempNumber /= 10; - } + /* storing the terms of the series */ + var termsArray = new int[number]; - var sum = 0; - var k = digitsInNumber; - while (sum < number) - { - sum = 0; + for (var i = digitsInNumber - 1; i >= 0; i--) + { + termsArray[i] = tempNumber % 10; + tempNumber /= 10; + } - for (var j = 1; j <= digitsInNumber; j++) - { - sum += termsArray[k - j]; - } + var sum = 0; + var k = digitsInNumber; + while (sum < number) + { + sum = 0; - termsArray[k] = sum; - k++; + for (var j = 1; j <= digitsInNumber; j++) + { + sum += termsArray[k - j]; } - return sum == number; + termsArray[k] = sum; + k++; } + + return sum == number; } } diff --git a/Algorithms/Numeric/KrishnamurthyNumberChecker.cs b/Algorithms/Numeric/KrishnamurthyNumberChecker.cs index 95a199e0..c4d245a0 100644 --- a/Algorithms/Numeric/KrishnamurthyNumberChecker.cs +++ b/Algorithms/Numeric/KrishnamurthyNumberChecker.cs @@ -1,38 +1,37 @@ using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// A Krishnamurthy number is a number whose sum of the factorial of digits +/// is equal to the number itself. +/// +/// For example, 145 is a Krishnamurthy number since: 1! + 4! + 5! = 1 + 24 + 120 = 145. +/// +public static class KrishnamurthyNumberChecker { /// - /// A Krishnamurthy number is a number whose sum of the factorial of digits - /// is equal to the number itself. - /// - /// For example, 145 is a Krishnamurthy number since: 1! + 4! + 5! = 1 + 24 + 120 = 145. + /// Check if a number is Krishnamurthy number or not. /// - public static class KrishnamurthyNumberChecker + /// The number to check. + /// True if the number is Krishnamurthy, false otherwise. + public static bool IsKMurthyNumber(int n) { - /// - /// Check if a number is Krishnamurthy number or not. - /// - /// The number to check. - /// True if the number is Krishnamurthy, false otherwise. - public static bool IsKMurthyNumber(int n) - { - int sumOfFactorials = 0; - int tmp = n; + int sumOfFactorials = 0; + int tmp = n; - if (n <= 0) - { - return false; - } - - while (n != 0) - { - int factorial = (int)Factorial.Calculate(n % 10); - sumOfFactorials += factorial; - n = n / 10; - } + if (n <= 0) + { + return false; + } - return tmp == sumOfFactorials; + while (n != 0) + { + int factorial = (int)Factorial.Calculate(n % 10); + sumOfFactorials += factorial; + n = n / 10; } + + return tmp == sumOfFactorials; } } diff --git a/Algorithms/Numeric/MillerRabinPrimalityChecker.cs b/Algorithms/Numeric/MillerRabinPrimalityChecker.cs index b458888a..0850a43c 100644 --- a/Algorithms/Numeric/MillerRabinPrimalityChecker.cs +++ b/Algorithms/Numeric/MillerRabinPrimalityChecker.cs @@ -1,84 +1,83 @@ using System; using System.Numerics; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// https://en.wikipedia.org/wiki/Miller-Rabin_primality_test +/// The Miller–Rabin primality test or Rabin–Miller primality test is a probabilistic primality test: +/// an algorithm which determines whether a given number is likely to be prime, +/// similar to the Fermat primality test and the Solovay–Strassen primality test. +/// It is of historical significance in the search for a polynomial-time deterministic primality test. +/// Its probabilistic variant remains widely used in practice, as one of the simplest and fastest tests known. +/// +public static class MillerRabinPrimalityChecker { /// - /// https://en.wikipedia.org/wiki/Miller-Rabin_primality_test - /// The Miller–Rabin primality test or Rabin–Miller primality test is a probabilistic primality test: - /// an algorithm which determines whether a given number is likely to be prime, - /// similar to the Fermat primality test and the Solovay–Strassen primality test. - /// It is of historical significance in the search for a polynomial-time deterministic primality test. - /// Its probabilistic variant remains widely used in practice, as one of the simplest and fastest tests known. - /// - public static class MillerRabinPrimalityChecker + /// Run the probabilistic primality test. + /// + /// Number to check. + /// Number of rounds, the parameter determines the accuracy of the test, recommended value is Log2(n). + /// Seed for random number generator. + /// True if is a highly likely prime number; False otherwise. + /// Error: number should be more than 3. + public static bool IsProbablyPrimeNumber(BigInteger n, BigInteger rounds, int? seed = null) + { + Random rand = seed is null + ? new() + : new(seed.Value); + return IsProbablyPrimeNumber(n, rounds, rand); + } + + private static bool IsProbablyPrimeNumber(BigInteger n, BigInteger rounds, Random rand) { - /// - /// Run the probabilistic primality test. - /// - /// Number to check. - /// Number of rounds, the parameter determines the accuracy of the test, recommended value is Log2(n). - /// Seed for random number generator. - /// True if is a highly likely prime number; False otherwise. - /// Error: number should be more than 3. - public static bool IsProbablyPrimeNumber(BigInteger n, BigInteger rounds, int? seed = null) + if (n <= 3) { - Random rand = seed is null - ? new() - : new(seed.Value); - return IsProbablyPrimeNumber(n, rounds, rand); + throw new ArgumentException($"{nameof(n)} should be more than 3"); } - private static bool IsProbablyPrimeNumber(BigInteger n, BigInteger rounds, Random rand) + // Input #1: n > 3, an odd integer to be tested for primality + // Input #2: k, the number of rounds of testing to perform, recommended k = Log2(n) + // Output: false = “composite” + // true = “probably prime” + + // write n as 2r·d + 1 with d odd(by factoring out powers of 2 from n − 1) + BigInteger r = 0; + BigInteger d = n - 1; + while (d % 2 == 0) { - if (n <= 3) - { - throw new ArgumentException($"{nameof(n)} should be more than 3"); - } + r++; + d /= 2; + } - // Input #1: n > 3, an odd integer to be tested for primality - // Input #2: k, the number of rounds of testing to perform, recommended k = Log2(n) - // Output: false = “composite” - // true = “probably prime” + // as there is no native random function for BigInteger we suppose a random int number is sufficient + int nMaxValue = (n > int.MaxValue) ? int.MaxValue : (int)n; + BigInteger a = rand.Next(2, nMaxValue - 2); // ; pick a random integer a in the range[2, n − 2] - // write n as 2r·d + 1 with d odd(by factoring out powers of 2 from n − 1) - BigInteger r = 0; - BigInteger d = n - 1; - while (d % 2 == 0) + while (rounds > 0) + { + rounds--; + var x = BigInteger.ModPow(a, d, n); + if (x == 1 || x == (n - 1)) { - r++; - d /= 2; + continue; } - // as there is no native random function for BigInteger we suppose a random int number is sufficient - int nMaxValue = (n > int.MaxValue) ? int.MaxValue : (int)n; - BigInteger a = rand.Next(2, nMaxValue - 2); // ; pick a random integer a in the range[2, n − 2] - - while (rounds > 0) + BigInteger tempr = r - 1; + while (tempr > 0 && (x != n - 1)) { - rounds--; - var x = BigInteger.ModPow(a, d, n); - if (x == 1 || x == (n - 1)) - { - continue; - } - - BigInteger tempr = r - 1; - while (tempr > 0 && (x != n - 1)) - { - tempr--; - x = BigInteger.ModPow(x, 2, n); - } - - if (x == n - 1) - { - continue; - } + tempr--; + x = BigInteger.ModPow(x, 2, n); + } - return false; + if (x == n - 1) + { + continue; } - return true; + return false; } + + return true; } } diff --git a/Algorithms/Numeric/ModularExponentiation.cs b/Algorithms/Numeric/ModularExponentiation.cs index 30c09b28..c7f4f2e9 100644 --- a/Algorithms/Numeric/ModularExponentiation.cs +++ b/Algorithms/Numeric/ModularExponentiation.cs @@ -1,43 +1,42 @@ using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// Modular exponentiation is a type of exponentiation performed over a modulus +/// Modular exponentiation c is: c = b^e mod m where b is base, e is exponent, m is modulus +/// (Wiki: https://en.wikipedia.org/wiki/Modular_exponentiation). +/// +public class ModularExponentiation { /// - /// Modular exponentiation is a type of exponentiation performed over a modulus - /// Modular exponentiation c is: c = b^e mod m where b is base, e is exponent, m is modulus - /// (Wiki: https://en.wikipedia.org/wiki/Modular_exponentiation). + /// Performs Modular Exponentiation on b, e, m. /// - public class ModularExponentiation + /// Base. + /// Exponent. + /// Modulus. + /// Modular Exponential. + public int ModularPow(int b, int e, int m) { - /// - /// Performs Modular Exponentiation on b, e, m. - /// - /// Base. - /// Exponent. - /// Modulus. - /// Modular Exponential. - public int ModularPow(int b, int e, int m) + // initialize result in variable res + int res = 1; + if (m == 1) { - // initialize result in variable res - int res = 1; - if (m == 1) - { - // 1 divides every number - return 0; - } - - if (m <= 0) - { - // exponential not defined in this case - throw new ArgumentException(string.Format("{0} is not a positive integer", m)); - } + // 1 divides every number + return 0; + } - for (int i = 0; i < e; i++) - { - res = (res * b) % m; - } + if (m <= 0) + { + // exponential not defined in this case + throw new ArgumentException(string.Format("{0} is not a positive integer", m)); + } - return res; + for (int i = 0; i < e; i++) + { + res = (res * b) % m; } + + return res; } } diff --git a/Algorithms/Numeric/NarcissisticNumberChecker.cs b/Algorithms/Numeric/NarcissisticNumberChecker.cs index 17e13960..4e6283d8 100644 --- a/Algorithms/Numeric/NarcissisticNumberChecker.cs +++ b/Algorithms/Numeric/NarcissisticNumberChecker.cs @@ -1,40 +1,39 @@ -using System; +using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// A Narcissistic number is equal to the sum of the cubes of its digits. For example, 370 is a +/// Narcissistic number because 3*3*3 + 7*7*7 + 0*0*0 = 370. +/// +public static class NarcissisticNumberChecker { /// - /// A Narcissistic number is equal to the sum of the cubes of its digits. For example, 370 is a - /// Narcissistic number because 3*3*3 + 7*7*7 + 0*0*0 = 370. + /// Checks if a number is a Narcissistic number or not. /// - public static class NarcissisticNumberChecker + /// Number to check. + /// True if is a Narcissistic number; False otherwise. + public static bool IsNarcissistic(int number) { - /// - /// Checks if a number is a Narcissistic number or not. - /// - /// Number to check. - /// True if is a Narcissistic number; False otherwise. - public static bool IsNarcissistic(int number) + var sum = 0; + var temp = number; + var numberOfDigits = 0; + while (temp != 0) { - var sum = 0; - var temp = number; - var numberOfDigits = 0; - while (temp != 0) - { - numberOfDigits++; - temp /= 10; - } - - temp = number; - while (number > 0) - { - var remainder = number % 10; - var power = (int)Math.Pow(remainder, numberOfDigits); + numberOfDigits++; + temp /= 10; + } - sum += power; - number /= 10; - } + temp = number; + while (number > 0) + { + var remainder = number % 10; + var power = (int)Math.Pow(remainder, numberOfDigits); - return sum == temp; + sum += power; + number /= 10; } + + return sum == temp; } } diff --git a/Algorithms/Numeric/PerfectNumberChecker.cs b/Algorithms/Numeric/PerfectNumberChecker.cs index 1a91def9..4abd55f8 100644 --- a/Algorithms/Numeric/PerfectNumberChecker.cs +++ b/Algorithms/Numeric/PerfectNumberChecker.cs @@ -1,37 +1,36 @@ -using System; +using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// In number theory, a perfect number is a positive integer that is equal to the sum of its positive +/// divisors, excluding the number itself.For instance, 6 has divisors 1, 2 and 3 (excluding +/// itself), and 1 + 2 + 3 = 6, so 6 is a perfect number. +/// +public static class PerfectNumberChecker { /// - /// In number theory, a perfect number is a positive integer that is equal to the sum of its positive - /// divisors, excluding the number itself.For instance, 6 has divisors 1, 2 and 3 (excluding - /// itself), and 1 + 2 + 3 = 6, so 6 is a perfect number. + /// Checks if a number is a perfect number or not. /// - public static class PerfectNumberChecker + /// Number to check. + /// True if is a perfect number; False otherwise. + /// Error number is not on interval (0.0; int.MaxValue). + public static bool IsPerfectNumber(int number) { - /// - /// Checks if a number is a perfect number or not. - /// - /// Number to check. - /// True if is a perfect number; False otherwise. - /// Error number is not on interval (0.0; int.MaxValue). - public static bool IsPerfectNumber(int number) + if (number < 0) { - if (number < 0) - { - throw new ArgumentException($"{nameof(number)} cannot be negative"); - } + throw new ArgumentException($"{nameof(number)} cannot be negative"); + } - var sum = 0; /* sum of its positive divisors */ - for (var i = 1; i < number; ++i) + var sum = 0; /* sum of its positive divisors */ + for (var i = 1; i < number; ++i) + { + if (number % i == 0) { - if (number % i == 0) - { - sum += i; - } + sum += i; } - - return sum == number; } + + return sum == number; } } diff --git a/Algorithms/Numeric/PerfectSquareChecker.cs b/Algorithms/Numeric/PerfectSquareChecker.cs index 8a8d937a..6465c3b6 100644 --- a/Algorithms/Numeric/PerfectSquareChecker.cs +++ b/Algorithms/Numeric/PerfectSquareChecker.cs @@ -1,26 +1,25 @@ -using System; +using System; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// A perfect square is an element of algebraic structure that is equal to the square of another element. +/// +public static class PerfectSquareChecker { /// - /// A perfect square is an element of algebraic structure that is equal to the square of another element. + /// Checks if a number is a perfect square or not. /// - public static class PerfectSquareChecker + /// Number too check. + /// True if is a perfect square; False otherwise. + public static bool IsPerfectSquare(int number) { - /// - /// Checks if a number is a perfect square or not. - /// - /// Number too check. - /// True if is a perfect square; False otherwise. - public static bool IsPerfectSquare(int number) + if (number < 0) { - if (number < 0) - { - return false; - } - - var sqrt = (int)Math.Sqrt(number); - return sqrt * sqrt == number; + return false; } + + var sqrt = (int)Math.Sqrt(number); + return sqrt * sqrt == number; } } diff --git a/Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs b/Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs index 5e8d6938..0e9317aa 100644 --- a/Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs +++ b/Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs @@ -1,46 +1,45 @@ -using System; +using System; using Algorithms.Numeric.Decomposition; using Utilities.Extensions; -namespace Algorithms.Numeric.Pseudoinverse +namespace Algorithms.Numeric.Pseudoinverse; + +/// +/// The Moore–Penrose pseudo-inverse A+ of a matrix A, +/// is a general way to find the solution to the following system of linear equations: +/// ~b = A ~y. ~b e R^m; ~y e R^n; A e Rm×n. +/// There are varios methods for construction the pseudo-inverse. +/// This one is based on Singular Value Decomposition (SVD). +/// +public static class PseudoInverse { /// - /// The Moore–Penrose pseudo-inverse A+ of a matrix A, - /// is a general way to find the solution to the following system of linear equations: - /// ~b = A ~y. ~b e R^m; ~y e R^n; A e Rm×n. - /// There are varios methods for construction the pseudo-inverse. - /// This one is based on Singular Value Decomposition (SVD). + /// Return the pseudoinverse of a matrix based on the Moore-Penrose Algorithm. + /// using Singular Value Decomposition (SVD). /// - public static class PseudoInverse + /// Input matrix to find its inverse to. + /// The inverse matrix approximation of the input matrix. + public static double[,] PInv(double[,] inMat) { - /// - /// Return the pseudoinverse of a matrix based on the Moore-Penrose Algorithm. - /// using Singular Value Decomposition (SVD). - /// - /// Input matrix to find its inverse to. - /// The inverse matrix approximation of the input matrix. - public static double[,] PInv(double[,] inMat) - { - // To compute the SVD of the matrix to find Sigma. - var (u, s, v) = ThinSvd.Decompose(inMat); + // To compute the SVD of the matrix to find Sigma. + var (u, s, v) = ThinSvd.Decompose(inMat); - // To take the reciprocal of each non-zero element on the diagonal. - var len = s.Length; + // To take the reciprocal of each non-zero element on the diagonal. + var len = s.Length; - var sigma = new double[len]; - for (var i = 0; i < len; i++) - { - sigma[i] = Math.Abs(s[i]) < 0.0001 ? 0 : 1 / s[i]; - } + var sigma = new double[len]; + for (var i = 0; i < len; i++) + { + sigma[i] = Math.Abs(s[i]) < 0.0001 ? 0 : 1 / s[i]; + } - // To construct a diagonal matrix based on the vector result. - var diag = sigma.ToDiagonalMatrix(); + // To construct a diagonal matrix based on the vector result. + var diag = sigma.ToDiagonalMatrix(); - // To construct the pseudo-inverse using the computed information above. - var matinv = u.Multiply(diag).Multiply(v.Transpose()); + // To construct the pseudo-inverse using the computed information above. + var matinv = u.Multiply(diag).Multiply(v.Transpose()); - // To Transpose the result matrix. - return matinv.Transpose(); - } + // To Transpose the result matrix. + return matinv.Transpose(); } } diff --git a/Algorithms/Numeric/RungeKuttaMethod.cs b/Algorithms/Numeric/RungeKuttaMethod.cs index 6f694e6e..2ef987bd 100644 --- a/Algorithms/Numeric/RungeKuttaMethod.cs +++ b/Algorithms/Numeric/RungeKuttaMethod.cs @@ -1,69 +1,68 @@ using System; using System.Collections.Generic; -namespace Algorithms.Numeric +namespace Algorithms.Numeric; + +/// +/// In numerical analysis, the Runge–Kutta methods are a family of implicit and explicit iterative methods, +/// used in temporal discretization for the approximate solutions of simultaneous nonlinear equations. +/// The most widely known member of the Runge–Kutta family is generally referred to as +/// "RK4", the "classic Runge–Kutta method" or simply as "the Runge–Kutta method". +/// +public static class RungeKuttaMethod { /// - /// In numerical analysis, the Runge–Kutta methods are a family of implicit and explicit iterative methods, - /// used in temporal discretization for the approximate solutions of simultaneous nonlinear equations. - /// The most widely known member of the Runge–Kutta family is generally referred to as - /// "RK4", the "classic Runge–Kutta method" or simply as "the Runge–Kutta method". - /// - public static class RungeKuttaMethod + /// Loops through all the steps until xEnd is reached, adds a point for each step and then + /// returns all the points. + /// + /// Initial conditions x-value. + /// Last x-value. + /// Step-size on the x-axis. + /// Initial conditions y-value. + /// The right hand side of the differential equation. + /// The solution of the Cauchy problem. + public static List ClassicRungeKuttaMethod( + double xStart, + double xEnd, + double stepSize, + double yStart, + Func function) { - /// - /// Loops through all the steps until xEnd is reached, adds a point for each step and then - /// returns all the points. - /// - /// Initial conditions x-value. - /// Last x-value. - /// Step-size on the x-axis. - /// Initial conditions y-value. - /// The right hand side of the differential equation. - /// The solution of the Cauchy problem. - public static List ClassicRungeKuttaMethod( - double xStart, - double xEnd, - double stepSize, - double yStart, - Func function) + if (xStart >= xEnd) { - if (xStart >= xEnd) - { - throw new ArgumentOutOfRangeException( - nameof(xEnd), - $"{nameof(xEnd)} should be greater than {nameof(xStart)}"); - } - - if (stepSize <= 0) - { - throw new ArgumentOutOfRangeException( - nameof(stepSize), - $"{nameof(stepSize)} should be greater than zero"); - } + throw new ArgumentOutOfRangeException( + nameof(xEnd), + $"{nameof(xEnd)} should be greater than {nameof(xStart)}"); + } - List points = new(); - double[] firstPoint = { xStart, yStart }; - points.Add(firstPoint); + if (stepSize <= 0) + { + throw new ArgumentOutOfRangeException( + nameof(stepSize), + $"{nameof(stepSize)} should be greater than zero"); + } - var yCurrent = yStart; - var xCurrent = xStart; + List points = new(); + double[] firstPoint = { xStart, yStart }; + points.Add(firstPoint); - while (xCurrent < xEnd) - { - var k1 = function(xCurrent, yCurrent); - var k2 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k1); - var k3 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k2); - var k4 = function(xCurrent + stepSize, yCurrent + stepSize * k3); + var yCurrent = yStart; + var xCurrent = xStart; - yCurrent += (1.0 / 6.0) * stepSize * (k1 + 2 * k2 + 2 * k3 + k4); - xCurrent += stepSize; + while (xCurrent < xEnd) + { + var k1 = function(xCurrent, yCurrent); + var k2 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k1); + var k3 = function(xCurrent + 0.5 * stepSize, yCurrent + 0.5 * stepSize * k2); + var k4 = function(xCurrent + stepSize, yCurrent + stepSize * k3); - double[] newPoint = { xCurrent, yCurrent }; - points.Add(newPoint); - } + yCurrent += (1.0 / 6.0) * stepSize * (k1 + 2 * k2 + 2 * k3 + k4); + xCurrent += stepSize; - return points; + double[] newPoint = { xCurrent, yCurrent }; + points.Add(newPoint); } + + return points; } } diff --git a/Algorithms/Numeric/Series/Maclaurin.cs b/Algorithms/Numeric/Series/Maclaurin.cs index cfae31a3..75924a50 100644 --- a/Algorithms/Numeric/Series/Maclaurin.cs +++ b/Algorithms/Numeric/Series/Maclaurin.cs @@ -1,140 +1,139 @@ using System; using System.Linq; -namespace Algorithms.Numeric.Series +namespace Algorithms.Numeric.Series; + +/// +/// Maclaurin series calculates nonlinear functions approximation +/// starting from point x = 0 in a form of infinite power series: +/// f(x) = f(0) + f'(0) * x + ... + (f'n(0) * (x ^ n)) / n! + ..., +/// where n is natural number. +/// +public static class Maclaurin { /// - /// Maclaurin series calculates nonlinear functions approximation - /// starting from point x = 0 in a form of infinite power series: - /// f(x) = f(0) + f'(0) * x + ... + (f'n(0) * (x ^ n)) / n! + ..., - /// where n is natural number. + /// Calculates approximation of e^x function: + /// e^x = 1 + x + x^2 / 2! + ... + x^n / n! + ..., + /// where n is number of terms (natural number), + /// and x is given point (rational number). /// - public static class Maclaurin - { - /// - /// Calculates approximation of e^x function: - /// e^x = 1 + x + x^2 / 2! + ... + x^n / n! + ..., - /// where n is number of terms (natural number), - /// and x is given point (rational number). - /// - /// Given point. - /// The number of terms in polynomial. - /// Approximated value of the function in the given point. - public static double Exp(double x, int n) => - Enumerable.Range(0, n).Sum(i => ExpTerm(x, i)); + /// Given point. + /// The number of terms in polynomial. + /// Approximated value of the function in the given point. + public static double Exp(double x, int n) => + Enumerable.Range(0, n).Sum(i => ExpTerm(x, i)); - /// - /// Calculates approximation of sin(x) function: - /// sin(x) = x - x^3 / 3! + ... + (-1)^n * x^(2*n + 1) / (2*n + 1)! + ..., - /// where n is number of terms (natural number), - /// and x is given point (rational number). - /// - /// Given point. - /// The number of terms in polynomial. - /// Approximated value of the function in the given point. - public static double Sin(double x, int n) => - Enumerable.Range(0, n).Sum(i => SinTerm(x, i)); + /// + /// Calculates approximation of sin(x) function: + /// sin(x) = x - x^3 / 3! + ... + (-1)^n * x^(2*n + 1) / (2*n + 1)! + ..., + /// where n is number of terms (natural number), + /// and x is given point (rational number). + /// + /// Given point. + /// The number of terms in polynomial. + /// Approximated value of the function in the given point. + public static double Sin(double x, int n) => + Enumerable.Range(0, n).Sum(i => SinTerm(x, i)); - /// - /// Calculates approximation of cos(x) function: - /// cos(x) = 1 - x^2 / 2! + ... + (-1)^n * x^(2*n) / (2*n)! + ..., - /// where n is number of terms (natural number), - /// and x is given point (rational number). - /// - /// Given point. - /// The number of terms in polynomial. - /// Approximated value of the function in the given point. - public static double Cos(double x, int n) => - Enumerable.Range(0, n).Sum(i => CosTerm(x, i)); + /// + /// Calculates approximation of cos(x) function: + /// cos(x) = 1 - x^2 / 2! + ... + (-1)^n * x^(2*n) / (2*n)! + ..., + /// where n is number of terms (natural number), + /// and x is given point (rational number). + /// + /// Given point. + /// The number of terms in polynomial. + /// Approximated value of the function in the given point. + public static double Cos(double x, int n) => + Enumerable.Range(0, n).Sum(i => CosTerm(x, i)); - /// - /// Calculates approximation of e^x function: - /// e^x = 1 + x + x^2 / 2! + ... + x^n / n! + ..., - /// and x is given point (rational number). - /// - /// Given point. - /// Last term error value. - /// Approximated value of the function in the given point. - /// Error value is not on interval (0.0; 1.0). - public static double Exp(double x, double error = 0.00001) => ErrorTermWrapper(x, error, ExpTerm); + /// + /// Calculates approximation of e^x function: + /// e^x = 1 + x + x^2 / 2! + ... + x^n / n! + ..., + /// and x is given point (rational number). + /// + /// Given point. + /// Last term error value. + /// Approximated value of the function in the given point. + /// Error value is not on interval (0.0; 1.0). + public static double Exp(double x, double error = 0.00001) => ErrorTermWrapper(x, error, ExpTerm); - /// - /// Calculates approximation of sin(x) function: - /// sin(x) = x - x^3 / 3! + ... + (-1)^n * x^(2*n + 1) / (2*n + 1)! + ..., - /// and x is given point (rational number). - /// - /// Given point. - /// Last term error value. - /// Approximated value of the function in the given point. - /// Error value is not on interval (0.0; 1.0). - public static double Sin(double x, double error = 0.00001) => ErrorTermWrapper(x, error, SinTerm); + /// + /// Calculates approximation of sin(x) function: + /// sin(x) = x - x^3 / 3! + ... + (-1)^n * x^(2*n + 1) / (2*n + 1)! + ..., + /// and x is given point (rational number). + /// + /// Given point. + /// Last term error value. + /// Approximated value of the function in the given point. + /// Error value is not on interval (0.0; 1.0). + public static double Sin(double x, double error = 0.00001) => ErrorTermWrapper(x, error, SinTerm); - /// - /// Calculates approximation of cos(x) function: - /// cos(x) = 1 - x^2 / 2! + ... + (-1)^n * x^(2*n) / (2*n)! + ..., - /// and x is given point (rational number). - /// - /// Given point. - /// Last term error value. - /// Approximated value of the function in the given point. - /// Error value is not on interval (0.0; 1.0). - public static double Cos(double x, double error = 0.00001) => ErrorTermWrapper(x, error, CosTerm); + /// + /// Calculates approximation of cos(x) function: + /// cos(x) = 1 - x^2 / 2! + ... + (-1)^n * x^(2*n) / (2*n)! + ..., + /// and x is given point (rational number). + /// + /// Given point. + /// Last term error value. + /// Approximated value of the function in the given point. + /// Error value is not on interval (0.0; 1.0). + public static double Cos(double x, double error = 0.00001) => ErrorTermWrapper(x, error, CosTerm); - /// - /// Wrapper function for calculating approximation with estimated - /// count of terms, where last term value is less than given error. - /// - /// Given point. - /// Last term error value. - /// Indexed term of approximation series. - /// Approximated value of the function in the given point. - /// Error value is not on interval (0.0; 1.0). - private static double ErrorTermWrapper(double x, double error, Func term) + /// + /// Wrapper function for calculating approximation with estimated + /// count of terms, where last term value is less than given error. + /// + /// Given point. + /// Last term error value. + /// Indexed term of approximation series. + /// Approximated value of the function in the given point. + /// Error value is not on interval (0.0; 1.0). + private static double ErrorTermWrapper(double x, double error, Func term) + { + if (error <= 0.0 || error >= 1.0) { - if (error <= 0.0 || error >= 1.0) - { - throw new ArgumentException("Error value is not on interval (0.0; 1.0)."); - } - - var i = 0; - var termCoefficient = 0.0; - var result = 0.0; + throw new ArgumentException("Error value is not on interval (0.0; 1.0)."); + } - do - { - result += termCoefficient; - termCoefficient = term(x, i); - i++; - } - while (Math.Abs(termCoefficient) > error); + var i = 0; + var termCoefficient = 0.0; + var result = 0.0; - return result; + do + { + result += termCoefficient; + termCoefficient = term(x, i); + i++; } + while (Math.Abs(termCoefficient) > error); - /// - /// Single term for e^x function approximation: x^i / i!. - /// - /// Given point. - /// Term index from 0 to n. - /// Single term value. - private static double ExpTerm(double x, int i) => Math.Pow(x, i) / (long)Factorial.Calculate(i); + return result; + } - /// - /// Single term for sin(x) function approximation: (-1)^i * x^(2*i + 1) / (2*i + 1)!. - /// - /// Given point. - /// Term index from 0 to n. - /// Single term value. - private static double SinTerm(double x, int i) => - Math.Pow(-1, i) / ((long)Factorial.Calculate(2 * i + 1)) * Math.Pow(x, 2 * i + 1); + /// + /// Single term for e^x function approximation: x^i / i!. + /// + /// Given point. + /// Term index from 0 to n. + /// Single term value. + private static double ExpTerm(double x, int i) => Math.Pow(x, i) / (long)Factorial.Calculate(i); - /// - /// Single term for cos(x) function approximation: (-1)^i * x^(2*i) / (2*i)!. - /// - /// Given point. - /// Term index from 0 to n. - /// Single term value. - private static double CosTerm(double x, int i) => - Math.Pow(-1, i) / ((long)Factorial.Calculate(2 * i)) * Math.Pow(x, 2 * i); - } + /// + /// Single term for sin(x) function approximation: (-1)^i * x^(2*i + 1) / (2*i + 1)!. + /// + /// Given point. + /// Term index from 0 to n. + /// Single term value. + private static double SinTerm(double x, int i) => + Math.Pow(-1, i) / ((long)Factorial.Calculate(2 * i + 1)) * Math.Pow(x, 2 * i + 1); + + /// + /// Single term for cos(x) function approximation: (-1)^i * x^(2*i) / (2*i)!. + /// + /// Given point. + /// Term index from 0 to n. + /// Single term value. + private static double CosTerm(double x, int i) => + Math.Pow(-1, i) / ((long)Factorial.Calculate(2 * i)) * Math.Pow(x, 2 * i); }