Saturday, 18 January 2025

Calculating day of the week from Lewis Carrol's algorithm

This article will present the Lewis Carrol's algorithm to calculate the day of the week.

Lewis Carrol was a mathematician and logician at Oxford University. He is a wellknown name today for his children novels such as "Alice in Wonderland", published in 1871.

The name Lewis Carrol is a pseudonym, his real name was Charles Ludwidge Dodgson. The algorithm presented here was released in 1883.

Let's look at the algorithm next for calculating the day of week. Please note that this algorithm got inaccuracies for leap years and february, it is about 90-95% correct for all dates. In many cases it will give correct results for other dates. The algorithm presented here is for dates after September 12, 1752 due to changes in our calendar system at that date.

Construct a Number n consisting of four parts where

n = CCYYMMDD | CC = century, YY = year, MM = month , DD = day

For example, the date 18th of January, 2025 is written as

n = 202501818 , CC = 20, YY = 25, MM = 01, DD = 18

The calculation will sum up the contributions from the numbers into a sum S as :

S = Sum_cc + Sum_yy + Sum_mm + Sum_dd

The sum S is then taken the modulo of 7, where the number 0 means Sunday and number 6 means Saturday. The day and week is therefore the modulo 7 of S where starting from Sunday, the remainder from S mod 7 is the weekday, starting with index 0 for Monday.

Let's look at the different part sums next.

Sum_cc is calculated as the 3 subtracted by the century modulo 4 and this is multiplied by two. Sum_yy is the year divided by 12 plus the year modulo 12 and the year divided by 4. Sum_mm is defined as the month looked up in a lookup table of values shown in the source code. Sum_dd is just the date of the month.

The summmation looks like this:


private int DayOfWeekLevisCarrol(int date){
 int century = date / 1_000_000;
 int year = (date / 10_000) % 100;
 int month = (date / 100) % 100;
 int dayValue = (date % 100); 
 int[] monthValues = [0, 3, 3, 6, 1, 4, 6, 2, 5, 0, 3, 12]; 
 int centuryValue = (3-(century % 4))*2;
 int yearValue = (year/12) + (year % 12) + ((year % 12)/4);
 int monthValue = monthValues[month-1]; 
 var sumValues= (centuryValue+yearValue+monthValue+dayValue) % 7;
 return sumValues;
}


This helper method shows the weekday presented as text:


public string WeekDayLevisCarrol(int date) => DayOfWeekLevisCarrol(date) switch {
 0 => "Sunday",
 1 => "Monday",
 2 => "Tuesday",
 3 => "Wednesday",
 4 => "Thursday",
 5 => "Friday",
 6 => "Saturday",
 _ => "Impossible"	
};


A more compressed way of expressing this is the following:


/// <summary>
/// Calculates the day of the week using Lewis Carroll's method.
/// </summary>
/// <param name="date">The date in yyyymmdd format.</param>
/// <returns>The day of the week as an integer (0 for Sunday, 1 for Monday, etc.).</returns>
private int DayOfWeekLevisCarrolV2(int date) =>
    ((3 - (date / 1_000_000 % 4)) * 2 +
    (date / 10_000 % 100 / 12) + (date / 10_000 % 100 % 12) + ((date / 10_000 % 100 % 12) / 4) +
    new int[] { 0, 3, 3, 6, 1, 4, 6, 2, 5, 0, 3, 12 }[(date / 100 % 100) - 1] +
    (date % 100)) % 7;



Running the metod for today (I used 18th of January 2025):

void Main()
{
	int todayDate = int.Parse(DateTime.Today.ToString("yyyyMMdd"));
	string todayWeekdate = WeekDayLevisCarrol(todayDate);
	Console.WriteLine($"Levis-Carrol's method calculated weekday for {DateTime.Today.ToShortDateString()} written as number {todayDate} is: {todayWeekdate}");
}


Output is:

Levis-Carrol's method calculated weekday for 18.01.2025 written as number 20250118 is: Saturday As noted, this algorithm is about 90-95% correct since it is not precise in leap years and in February. There are more precise calculation algorithms today, the algorithm is to be considered a crude approach and more modern algorithms exists today. In .NET, the calculation of the weekday is done inside the DateTime struct calculated from the start of Epoch at year 1 and there are precise calculations considering shifts in calendars and leap years, so of course you should use DayOfWeek of a DateTime in .NET. But it is a fun trivia that the author of "Alice in Wonderland" was a also a noted mathematician and according to history, he calculated weekdays for dates in few seconds using this algorithm, often being correct.

Monday, 13 January 2025

Diophantine linear equation solver in C#

Diophantine linear equations are equations of two unkowns on the format

ax + by = c

The code below shows how we can solve Diophantine equations with a sample equation :

7x - 9y = 3

DiophantineSolver.cs



class Program
{
	static void Main()
	{
		int a = 7,
		b = -9, 
		c = 3;
		DiophantineSolver.DiophantineSolution solution = DiophantineSolver.ExtendedGcd(a, b, c);

		Console.WriteLine(solution.ToString());
	}
}

/// <summary>
/// Class containing methods to solve Diophantine equations using the Extended Euclidean Algorithm.
/// </summary>
class DiophantineSolver
{
	/// <summary>
	/// Struct to hold the solution of a Diophantine equation.
	/// </summary>
	public struct DiophantineSolution
	{
		/// <summary>
		/// Indicates whether the equation has a solution.
		/// </summary>
		public bool IsSolvable { get; set; }

		/// <summary>
		/// X0, a articular solution for variable x, or null if no solution exists. Other solutions of X can be derived from this.
		/// </summary>
		public int? X0 { get; set; }

		/// <summary>
		/// Y0, a particular solution, or null if no solution exists. Other solutions of Y can be derived from this.
		/// </summary>
		public int? Y0 { get; set; }

		/// <summary>
		/// The coefficient of x in the equation.
		/// </summary>
		public readonly int A;

		/// <summary>
		/// The coefficient of y in the equation.
		/// </summary>
		public readonly int B;

		/// <summary>
		/// The constant term in the equation.
		/// </summary>
		public readonly int C;
		
		/// <summary>D: The greatest common divisor gcd(a,b) = d</summary>
		/// <remarks>Solutions to the Diophantine linear equation exists only if C|D: D is divisible by C.</remarks>
		public readonly int? D;

		/// <summary>
		/// Initializes a new instance of the DiophantineSolution struct.
		/// </summary>
		/// <param name="isSolvable">Indicates whether the equation has a solution.</param>
		/// <param name="x0">The x-coordinate of the particular solution.</param>
		/// <param name="y0">The y-coordinate of the particular solution.</param>
		/// <param name="a">The coefficient of x in the equation.</param>
		/// <param name="b">The coefficient of y in the equation.</param>
		/// <param name="c">The constant term in the equation.</param>
		/// <param name="d">The d= gcd(a,b) greatest common divisor found for a and b</param>
		public DiophantineSolution(bool isSolvable, int? x0, int? y0, int a, int b, int c, int d)
		{
			IsSolvable = isSolvable;
			X0 = x0;
			Y0 = y0;
			A = a;
			B = b;
			C = c;
			D = d;
		}

		/// <summary>
		/// Returns a string representation of the solution.
		/// </summary>
		/// <returns>A string representation of the solution.</returns>
		public override string ToString()
		{
			if (IsSolvable)
			{
				string result = $"""
				Diophantine Linear Equation Solver result for equation 
				AX + BY = C
				{A}x + {B}y = {C} 
				A={A}, B={B}, C={C}
				Greatest common divisor :
				d = gcd(a,b) = {D?.ToString() ?? "No gcd found"}
				Solutions exists? {IsSolvable}, since d|c = {D}|{C}: {C} is divisible by {D}
				Particular solution A(x0) + B(y0) = C
				x0 = {X0}, y0 = {Y0}
				({A}*{X0}) + ({B}*{Y0}) = {C}
				({A * X0}) + ({B * Y0}) = {C}
				Generalized solution:				
				{GetParameterizedSolution()}
				Additional solutions:
				""";
				if (D > 0)
				{
					for (int t = 1; t <= 3; t++)
					{
						int xn = X0!.Value + (B / D!.Value) * t;
						int yn = Y0!.Value - (A / D!.Value) * t;
						result += $"\nSolution {t}: x = {xn}, y = {yn} (t = {t}) , since : ({xn}*{X0}) + ({yn}*{Y0}) = ({xn * X0}) + ({yn * Y0}) = {C}";
					}
				}
				else {
					result += "No additional solutions exists"; 
				}
				return result;
			}
			else
			{
				return $"Diophantine Linear Equation Solver result\nNo solution exists. d = gcd(a,b) = {D}. Solution exists if: {C}|{D}. This gives a remainder of {C % D} != 0. {C} is not divisible by {D}. C = {C} = n*D + r = ({C / D}*{D}) + {C % D}.";
			}
		}

		/// <summary>
		/// Returns the parameterized solution expression as a string.
		/// </summary>
		/// <returns>A string representation of the parameterized solution.</returns>
		public string GetParameterizedSolution()
		{
			if (D == 0){
				return "No solution exists, so there is not general parameterized solution";
			}
			string generalEquation = "x' = x_0 + (B/D), y' = y_0 - (A/D)";
			return $"{generalEquation}{Environment.NewLine}x' = {X0} + ({B}/{D})t, y' = {Y0} - ({A}/{D})t, where t is any integer";
		}
	}

	/// <summary>
	/// Solves the Diophantine equation ax + by = c using the Extended Euclidean Algorithm.
	/// </summary>
	/// <param name="a">The coefficient of x.</param>
	/// <param name="b">The coefficient of y.</param>
	/// <param name="c">The constant term.</param>
	/// <returns>A DiophantineSolution struct containing the solution.</returns>
	public static DiophantineSolution ExtendedGcd(int a, int b, int c)
	{
		int X0, Y0;
		int gcd = ExtendedGcd(a, b, out X0, out Y0);

		if (c % gcd != 0)
		{
			return new DiophantineSolution(false, null, null, a, b, c, gcd);
		}
		else
		{
			int x = X0 * (c / gcd);
			int y = Y0 * (c / gcd);
			return new DiophantineSolution(true, x, y, a, b, c, gcd);
		}
	}

	/// <summary>
	/// Computes the greatest common divisor of a and b, and finds coefficients X0 and Y0 such that a*X0 + b*Y0 = gcd(a, b).
	/// </summary>
	/// <param name="a">The first integer.</param>
	/// <param name="b">The second integer.</param>
	/// <param name="X0">The coefficient of a in the linear combination.</param>
	/// <param="Y0">The coefficient of b in the linear combination.</param>
	/// <returns>The greatest common divisor of a and b.</returns>
	private static int ExtendedGcd(int a, int b, out int X0, out int Y0)
	{
		if (b == 0)
		{
			X0 = 1;
			Y0 = 0;
			return a; // Base case: gcd(a, 0) = a
		}

		int X1, Y1;
		int gcd = ExtendedGcd(b, a % b, out X1, out Y1); // Recursive call

		X0 = Y1; // Back substitution: X0 = Y1
		Y0 = X1 - (a / b) * Y1; // Back substitution: Y0 = X1 - (a / b) * Y1

		return gcd; // Return the gcd
	}
}



The extended Euclidean algorithm calculates the greatest common divisor and in the recursion also applies back substitution. The coefficients X0 and Y0 founds satisfy the equation :

(a · X0) + (b · Y0) = gcd(a, b)



That is, a linear combination of the original integers a and b equaling the gcd(a,b). The output of running the code shown above is:


Diophantine Linear Equation Solver result for equation
AX + BY = C
7x + -9y = 3
A=7, B=-9, C=3
Greatest common divisor :
d = gcd(a,b) = 1
Solutions exists? True, since d|c = 1|3: 3 is divisible by 1
Particular solution A(x0) + B(y0) = C
x0 = 12, y0 = 9
(7*12) + (-9*9) = 3
(84) + (-81) = 3
Generalized solution:               
x' = x_0 + (B/D), y' = y_0 - (A/D)
x' = 12 + (-9/1)t, y' = 9 - (7/1)t, where t is any integer
Additional solutions:
Solution 1: x = 3, y = 2 (t = 1) , since : (3*12) + (2*9) = (36) + (18) = 3
Solution 2: x = -6, y = -5 (t = 2) , since : (-6*12) + (-5*9) = (-72) + (-45) = 3
Solution 3: x = -15, y = -12 (t = 3) , since : (-15*12) + (-12*9) = (-180) + (-108) = 3

So an example of a solution of the equation : 7x -9y = 3 Is: x = 12, y = 9, since (7*12) - (9*9) = 84 - 81 = 3. Note that once you have found a solution to a linear Diophantine equation, you can found infinite other solutions. More formally:

Diophantine Equation (Linear) Definition and Solvability

A linear Diophantine equation is an equation of the form:

(a · x) + (b · y) = c

where a, b, and c are integers, and x and y are unknown integers.

Solvability

The linear Diophantine equation (a · x) + (b · y) = c has integer solutions if and only if the greatest common divisor (gcd) of a and b divides c. In mathematical notation:

d is gcd. d|c => gcd(a, b) | c

If gcd(a, b) divides c (it is common to use d to notify gcd as a short term notation), then there exist integer solutions x and y such that:

(a · x) + (b · y) = c

Otherwise, there are no integer solutions.

Sunday, 5 January 2025

Euclidean algorithm in C# to find GCD and LCM

I am reading my Elementary Number Theory book by David M. Burton for a course i took at University, MNFMA104 Tallteori at NTNU in Trondheim. I like Number Theory in Math a lot and will look into writing some C# code, starting with this article. This article will look at how we can find the greatest common divisor (gcd). Mathematically speaking, the definition is the following:

Greatest Common Divisor (GCD)

Given two integers a and b, their greatest common divisor (GCD), denoted as d, is the largest positive integer that divides both a and b without leaving a remainder. Note that divides here means there is no remainer.

Formally, if d = gcd(a, b), then:

  1. d divides both a and b. This means d | a and d | b.
  2. If there is any other integer c that divides both a and b, then c ≤ d. This ensures that d is the greatest such integer.

For example the numbers 12 and 8 can be divided by 4. And we will see that d = gcd(12,8) = 4 next. But first, let's look at the Euclidean algorithm mathematical definition.

Euclidean Algorithm for GCD

Given two integers a and b (where a ≥ b > 0), their greatest common divisor (GCD), denoted as d, can be found using the Euclidean algorithm. The process is as follows:

  1. Let a and b be two integers such that a ≥ b and b > 0.
  2. Define a sequence of equations r0, r1, r2, …, where r0 = a and r1 = b.
  3. For i ≥ 0, compute ri+2 using the division algorithm:
    ri+2 = ri mod ri+1
  4. Continue this process until rn+1 = 0 for some integer n. At this point:
    d = rn = gcd(a, b)

In summary, the GCD of a and b is the last non-zero remainder obtained through this iterative process:

gcd(a, b) = d


Let's look at C# code to calculate the GDC, we will use the Euclidean algorithm to calculate it.

int Gcd(int a, int b)
{
	int q_n = Math.Abs(a);
	int r_n = Math.Abs(b);
	while (r_n != 0)
	{
		int remainder = q_n % r_n;
		q_n = r_n;
		r_n = remainder;
	}
	return q_n; // returning here the second-last quotient that was non-zero, which is the gcd
}

Calculating the gcd of 12 and 8 gives:

void Main()
{
	int a = 8;
	int b = 12;
	int gcd = Gcd(a, b);
	Console.WriteLine($"Demo - The greatest common divisor - GCD - using the Euclidean algorithm of the numbers : ");
	Console.WriteLine($"d = gcd({a}, {b}) = {gcd}");
    
}

//output 
// Demo - The greatest common divisor - GCD - using the Euclidean algorithm of the numbers : 
// d= gcd(8, 12) = 4


We can also calculate gcd of some bigger numbers. For example, the gcd of the numbers a= 5040 and b = 3780 are 1260. Both numbers 5040 and 3780 are divisible by 1260. We can also use C# patterns in case we want to apply a more functional approach. Together with tuples and C# patterns, the code becomes very compact, especially when we use recursion and tuples to compose the current state of the currend dividend and divisor, the two numbers we consider and get the remainder from.



<summary>
Calculate the greatest common divisor. Recursive C# pattern of GCD (using 'state approach' for tuple)
</summary>
int Gcd(int a, int b) => (a, b) switch
{
    (0, _) => Math.Abs(b),
    (_, 0) => Math.Abs(a),
    _ => Gcd(b, a % b)
};



Calculating if two number are coprimes

Let's turn attention to coprime numbers a bit. It is a consequence that if the gcd of two numbers a and b, they are said to be relatively prime, or coprime.

In mathematical terms, two integers a and b are said to be coprime (or relatively prime) if their greatest common divisor (gcd) is 1.

Formally, given two integers a and b, a and b are coprime if:

gcd(a, b) = 1

This definition means that the only positive integer that divides both a and b is 1. In other words, they have no common prime factors.

For example, the numbers a = 1234, b = 3417 can be tested if they are coprimes, by calculating the gcd and if the gcd is 1, they are coprime. They have no other common factors than 1, i.e. they are both coprimes or relatively prime. Let's first define a method in C# :

bool AreCoprime(int a, int b) => Gcd(a,b) == 1; 

We can then verify that a = 1234, b = 3417 are relatively prime and the two numbers are relatively prime.

int a1 = 1234, b1 = 3417;
Console.WriteLine($"\nDemo - If two numbers - GCD - are equal to 1, they are co-primes (relative primes)  : ");
Console.WriteLine($"Are {a1} and {b1} coprime? {AreCoprime(a1, b1)}");

The output is:

Demo - If two numbers - GCD - are equal to 1, they are co-primes (relative primes)  : 
Are 1234 and 3417 coprime? True

Calculating the least common multiple (lcm)

The least common multiple - lcm - can be calculated from the gcd for two numbers a and b. The mathematical definition is this:

In mathematics, the least common multiple (LCM) of two integers a and b is the smallest positive integer m such that both a and b divide m without leaving a remainder. Formally, it can be defined as:

LCM(a, b) = |a · b| / GCD(a, b)

where GCD(a, b) is the greatest common divisor of a and b.

Consider the gcd of a= 5040 and b = 3780. We found that the gcd for the two numbers is : gcd(5040, 3780) = 1260 Here is the code for calculating the lcm:

int Lcm(int a, int b) => (a*b)/Gcd(a,b);


int lcm = Lcm(a,b);
Console.WriteLine($"Demo - Least common multiple - lcm: ");	
Console.WriteLine($"lcmd({a}, {b}) = {lcm}");

This outputs:

Demo - Least common multiple - lcm: 
lcm(5040, 3780) = 15120

The number 15120 above is divisible by both a and b, giving 3 and 4 and no remainder, and it is the lowest number which is both divisible, that is - no integer remainder from the division. We write this if we use m for lcm as : m | a and m | b. To calculate the lcm for multiple numbers, we can use:

int Lcm(int[] numbers) { 
	return numbers.Aggregate((a,b) => Lcm(a,b)); 
}

Note that Aggregate method here (from Linq) does not sum, but fold In functional programming, the term "fold" is often used to describe the process of reducing a collection of elements to a single value by iteratively applying a function. In C#, the Aggregate method is essentially a fold operation, as it reduces the array of numbers to a single value (in this case, the GCD or LCM) by applying the specified function. So, using Aggregate to calculate the GCD or LCM of multiple numbers is indeed an example of a fold operation. Example:

int a = 5040;
int b = 3780;
int c = 2520;
Console.WriteLine($"Demo - Least common multiple - lcm: ");
Console.WriteLine($"lcm({a}, {b}, {c}) = {lcm}");

Output is: Demo - Least common multiple - lcm: lcm(5040, 3780, 2520) = 15120