August 16, 2023
Some years ago I posted my QRat implementation on GH:

	https://github.com/quickfur/qrat

QRat is a small library for exact computations of the so-called
quadratic irrationals: algebraic numbers of the form (a+b√c) where c is
a square-free integer and a,b are rationals.  It's useful for certain
computations involving regular polygons (equilateral triangles,
pentagons, octagons) whose coordinates involve numbers of this form.
All computations are exact and not subject to roundoff errors that you'd
get if you used floating-point numbers instead.

Recently I implemented .sqrt for QRat: this function will return the
square root of a given number (a+b√c) if it happens to be a perfect
square within the extension field Q(√c).  For example:

	(3 - 2*surd!2).sqrt == surd!2 - 1
	(7 - 4*surd!3).sqrt == 2 - surd!3

True, this won't help you if the input number isn't a perfect square, but it's useful in some cases where the length of certain vectors are representible within the extension field (i.e. the sum of squares of components is a perfect square).

Unfortunately I couldn't implement .sqrt for BigInt coefficients yet, because currently there isn't an efficient BigInt square root function. (I'd implement one myself, but then I wouldn't be sure whether it would perform well.)

Also, the current implementation of isSquare is kind of a hack; it uses floating-point sqrt() to compute the square root and then checks after the fact whether the input was a perfect square.  I wasn't sure whether it was worth implementing integer sqrt() from scratch if we could leverage the intrinsic floating-point hardware sqrt instruction.


T

-- 
VI = Visual Irritation
August 16, 2023

On Wednesday, 16 August 2023 at 14:09:41 UTC, H. S. Teoh wrote:

>

Unfortunately I couldn't implement .sqrt for BigInt coefficients yet, because currently there isn't an efficient BigInt square root function. (I'd implement one myself, but then I wouldn't be sure whether it would perform well.)

/// biggest integer with r^2 <= n
T sqrt(T)(T n) if(isUnsigned!T)
{
   if(n<2) return n; // trivial cases: sqrt(1) = 1 and sqrt(0) = 0

   uint pos = (bitlen(n)-1)>>1; // root has half the bitlen
   T r, s, t;
   r.bit[pos] = true;
   s.bit[pos<<1] = true; // s and t are used for extra fast squaring
                         // of a number with only one additional bit set
   while(pos--)
   {
      s.bit[pos<<1] = true;
      t = r; t <<= pos+1; t += s;
      if(t > n)
      {
         s.bit[pos<<1] = false;
         continue;
      }
      r.bit[pos] = true;
      if(!pos) return r;

      s = t; t -= n;
      if(!t) return r;

      if(pos > ((bitlen(t)-1)>>1))
         pos = (bitlen(t)-1)>>1;
   }
   return r;
}

This will work also with Bigint and will be reasonable fast.