ClassiCube/third_party/fix16_sqrt.c
2025-02-16 09:15:13 +11:00

71 lines
2.6 KiB
C

#include <stdint.h>
/*-
* Copyright (c) 2019 Christophe Meessen
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE. */
/*-
* Computing the square root of an integer or a fixed point integer into a fixed
* point integer. A fixed point is a 32 bit value with the comma between the
* bits 15 and 16, where bit 0 is the less significant bit of the value.
*
* The algorithm uses the property that computing x² is trivial compared to the
* sqrt. It will thus search the biggest x so that x² <= v, assuming we compute
* sqrt(v).
*
* The algorithm tests each bit of x starting with the most significant toward
* the less significant. It tests if the bit must be set or not.
*
* The algorithm uses the relation (x + a)² = x² + 2ax + a² to add the bit
* efficiently. Instead of computing x² it keeps track of (x + a)² - x².
*
* When computing sqrt(v), r = v - x², q = 2ax, b = a² and t = 2ax + a2.
*
* Note that the input integers are signed and that the sign bit is used in the
* computation. To accept unsigned integer as input, unfolding the initial loop
* is required to handle this particular case. See the usenet discussion for the
* proposed solution.
*
* Algorithm and code Author Christophe Meessen 1993. */
static int32_t sqrt_fix16(int32_t value)
{
uint32_t t;
uint32_t r = value;
uint32_t b = 0x40000000;
uint32_t q = 0;
while (b > 0x40) {
t = q + b;
if (r >= t) {
r -= t;
q = t + b; /* Equivalent to q += 2 * b */
}
r <<= 1;
b >>= 1;
}
q >>= 8;
return q;
}