diff options
author | Michael Sevakis <jethead71@rockbox.org> | 2010-06-08 04:51:00 +0000 |
---|---|---|
committer | Michael Sevakis <jethead71@rockbox.org> | 2010-06-08 04:51:00 +0000 |
commit | 67277f16888598ab211ad0f84cbfe490effbc9e1 (patch) | |
tree | 8bdd46fedb23857f534893ce3a9fd791441224a8 /apps | |
parent | 88f4a24c3a6723bd7951656e1c7827d44bbb7e6d (diff) | |
download | rockbox-67277f16888598ab211ad0f84cbfe490effbc9e1.tar.gz rockbox-67277f16888598ab211ad0f84cbfe490effbc9e1.zip |
Improve accuracy of NR-based fp_sqrt with better initial estimation and using one more bit internally. More reliable early termination. Good enough until better method is completed.
git-svn-id: svn://svn.rockbox.org/rockbox/trunk@26681 a1c6a512-1295-4272-9138-f99709370657
Diffstat (limited to 'apps')
-rw-r--r-- | apps/fixedpoint.c | 35 |
1 files changed, 21 insertions, 14 deletions
diff --git a/apps/fixedpoint.c b/apps/fixedpoint.c index fb89a8d30f..a8f606a5db 100644 --- a/apps/fixedpoint.c +++ b/apps/fixedpoint.c | |||
@@ -152,29 +152,36 @@ long fp_sincos(unsigned long phase, long *cos) | |||
152 | * @return Square root of argument in same fixed point format as input. | 152 | * @return Square root of argument in same fixed point format as input. |
153 | * | 153 | * |
154 | * This routine has been modified to run longer for greater precision, | 154 | * This routine has been modified to run longer for greater precision, |
155 | * but cuts calculation short if the answer is reached sooner. In | 155 | * but cuts calculation short if the answer is reached sooner. |
156 | * general, the closer x is to 1, the quicker the calculation. | ||
157 | */ | 156 | */ |
158 | long fp_sqrt(long x, unsigned int fracbits) | 157 | long fp_sqrt(long x, unsigned int fracbits) |
159 | { | 158 | { |
160 | long b = x/2 + BIT_N(fracbits); /* initial approximation */ | 159 | unsigned long xfp, b; |
161 | long c; | 160 | int n = 8; /* iteration limit (should terminate earlier) */ |
162 | unsigned n; | 161 | |
163 | const unsigned iterations = 8; | 162 | if (x <= 0) |
164 | 163 | return 0; /* no sqrt(neg), or just sqrt(0) = 0 */ | |
165 | for (n = 0; n < iterations; ++n) | 164 | |
165 | /* Increase working precision by one bit */ | ||
166 | xfp = x << 1; | ||
167 | fracbits++; | ||
168 | |||
169 | /* Get the midpoint between fracbits index and the highest bit index */ | ||
170 | b = ((sizeof(xfp)*8-1) - __builtin_clzl(xfp) + fracbits) >> 1; | ||
171 | b = BIT_N(b); | ||
172 | |||
173 | do | ||
166 | { | 174 | { |
167 | c = fp_div(x, b, fracbits); | 175 | unsigned long c = b; |
176 | b = (fp_div(xfp, b, fracbits) + b) >> 1; | ||
168 | if (c == b) break; | 177 | if (c == b) break; |
169 | b = (b + c)/2; | ||
170 | } | 178 | } |
179 | while (n-- > 0); | ||
171 | 180 | ||
172 | return b; | 181 | return b >> 1; |
173 | } | 182 | } |
174 | 183 | ||
175 | /* Accurate int sqrt with only elementary operations. (the above | 184 | /* Accurate int sqrt with only elementary operations. |
176 | * routine fails badly without enough iterations, more iterations | ||
177 | * than this requires -- [give that one a FIXME]). | ||
178 | * Snagged from: | 185 | * Snagged from: |
179 | * http://www.devmaster.net/articles/fixed-point-optimizations/ */ | 186 | * http://www.devmaster.net/articles/fixed-point-optimizations/ */ |
180 | unsigned long isqrt(unsigned long x) | 187 | unsigned long isqrt(unsigned long x) |