Skip to content

Instantly share code, notes, and snippets.

@micampbell
Last active December 12, 2024 01:02
Show Gist options
  • Save micampbell/4a84f8839de902f8253a2053e32c0c53 to your computer and use it in GitHub Desktop.
Save micampbell/4a84f8839de902f8253a2053e32c0c53 to your computer and use it in GitHub Desktop.
A method to find the square root of Int128 without any fancy under the hood stuff
/// <summary>
/// There currently is no method in C# to take the square root of an Int128.
/// However, we occasionally need to do this. This method finds the integer
/// component of the square root - as if you did a Floor function on the
/// actual square root.
/// </summary>
/// <param name="num"></param>
/// <returns></returns>
internal static Int128 SquareRoot(Int128 num)
{
if (num < 0)
throw new ArgumentException("Negative numbers are not " +
"allowed in this square root function.");
if (num<=long.MaxValue)
return (Int128)Math.Sqrt((long)num);
// start out with approximation when converted to double
var fracTerm = Math.Sqrt((double)num);
// get the "floor" integer from the double
var intTerm = (Int128)fracTerm;
// squaring should give us close to the input, num
// given that digits were likely lost in the conversion to
// double, then this might be larger than the expected answer
var intTermSqd = intTerm * intTerm;
if (intTermSqd < 0)
{ // this happens for large numbers where the casting
// to double and re-squaring leads to overflow (now a negative result)
// it seems a slight fractional reduction is enough to solve the problem
fracTerm /= 1.00000000001;
intTerm = (Int128)fracTerm;
intTermSqd = intTerm * intTerm;
}
// now here is the difference, generally positive but not always
var delta = num - intTerm * intTerm;
// consider the problem as (x + y)^2 = num where x is the integer
// part and y is the fractional part. if the fractional part is
// greater than 1 then we add to the integer part.
if (Int128.Abs(delta) > (intTerm << 1))
{
fracTerm = (double)delta / (double)(intTerm << 1);
intTerm += (Int128)fracTerm;
delta = num - intTerm * intTerm;
}
// occasionally this is still larger. this happens because 2xy +y^2
// produce and extra unit, but never more than once.
if (intTerm * intTerm > num)
intTerm--;
return intTerm;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment