/*	fastmath.h

This file is part of a program that implements a Software-Defined Radio.

Copyright (C) 2013, 2024, 2025 Warren Pratt, NR0V

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

The author can be reached by email at

warren@wpratt.com

*/

#ifndef _fastmath_h
#define _fastmath_h

#include <stdint.h>
#include <string.h>

/* log10 for strictly positive, finite, normal arguments; about 2.4x the
   throughput of libm's.

   Accurate to within 2.0e-13 absolute over (0, 1], which is far tighter than
   mlog10() in meterlog10.c (2.1e-4 absolute -- a 2048-entry table with no
   interpolation).  That table is fine for driving a dB meter display, but too
   coarse where the result feeds arithmetic rather than a readout.

   Caller must guarantee x > 0 and normal.  Zero, negatives, denormals,
   infinities and NaN are not handled.  */
static inline double wdsp_log10 (double x)
{
	uint64_t bits;
	double m, s, s2, p;
	int e;

	memcpy (&bits, &x, sizeof (bits));
	e = (int)((bits >> 52) & 0x7FF) - 1023;
	/* clear the exponent field, leaving the mantissa in [1, 2) */
	bits = (bits & 0x000FFFFFFFFFFFFFULL) | 0x3FF0000000000000ULL;
	memcpy (&m, &bits, sizeof (m));

	/* Recentre onto [sqrt(1/2), sqrt(2)) so the series stays in its
	   fast-converging range; |s| <= 0.1716 afterwards. */
	if (m > 1.4142135623730951)
	{
		m *= 0.5;
		e += 1;
	}

	/* log(m) = 2 * atanh(s), s = (m-1)/(m+1) */
	s = (m - 1.0) / (m + 1.0);
	s2 = s * s;
	p = 2.0 * (s + s * s2 * (3.3333333333333331e-01 + s2 * (2.0000000000000001e-01
		+ s2 * (1.4285714285714285e-01 + s2 * (1.1111111111111110e-01
		+ s2 * (9.0909090909090912e-02 + s2 * 7.6923076923076927e-02))))));

	/* log10(x) = (log(m) + e * ln2) / ln10 */
	return (p + (double)e * 6.9314718055994531e-01) * 4.3429448190325182e-01;
}

#endif