From e9d7842caf2f1cc698efc8ebea1672db5a1fb666 Mon Sep 17 00:00:00 2001
From: Jean-Marc Valin <jeanmarcv@google.com>
Date: Sun, 30 Jun 2024 01:42:06 -0400
Subject: [PATCH] More accurate fixed-point sqrt() approximation

---
 celt/mathops.c | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/celt/mathops.c b/celt/mathops.c
index accb4700..64c9f457 100644
--- a/celt/mathops.c
+++ b/celt/mathops.c
@@ -126,7 +126,10 @@ opus_val32 celt_sqrt(opus_val32 x)
    int k;
    opus_val16 n;
    opus_val32 rt;
-   static const opus_val16 C[5] = {23175, 11561, -3011, 1699, -664};
+   /* These coeffs are optimized in fixed-point to minimize both RMS and max error
+      of sqrt(x) over .25<x<1 without exceeding 32767.
+      The RMS error is 3.4e-5 and the max is 8.2e-5. */
+   static const opus_val16 C[6] = {23171, 11574, -2901, 1592, -1002, 336};
    if (x==0)
       return 0;
    else if (x>=1073741824)
@@ -134,8 +137,8 @@ opus_val32 celt_sqrt(opus_val32 x)
    k = (celt_ilog2(x)>>1)-7;
    x = VSHR32(x, 2*k);
    n = x-32768;
-   rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
-              MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
+   rt = ADD32(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
+              MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, ADD16(C[4], MULT16_16_Q15(n, (C[5])))))))))));
    rt = VSHR32(rt,7-k);
    return rt;
 }