Re: [PATCH 1/3] fixp-arith: replace sin/cos table by a better precision one

[Dmitry Torokhov <dmitry.torokhov@gmail.com>]

On Wed, Feb 04, 2015 at 10:07:30AM +0100, Hans Verkuil wrote:
> From: Mauro Carvalho Chehab <mchehab@osg.samsung.com>
> 
> The cos table used at fixp-arith.h has only 8 bits of precision.
> That causes problems if it is reused on other drivers.
> 
> As some media drivers require a higher precision sin/cos
> implementation, replace the current implementation by one that
> will provide 32 bits precision.
> 
> The values generated by the new implementation matches the
> 32 bit precision of glibc's sin for an angle measured in
> integer degrees.
> 
> It also provides support for fractional angles via linear
> interpolation. On experimental calculus, when used a table
> with a 0.001 degree angle, the maximum error for sin is
> 0.000038, which is likely good enough for practical purposes.
> 
> There are some logic there that seems to be specific to the
> usage inside ff-memless.c. Move those logic to there, as they're
> not needed elsewhere.
> 
> Cc: Hans de Goede <hdegoede@redhat.com>
> Cc: Dmitry Torokhov <dmitry.torokhov@gmail.com>
> Cc: linux-input@vger.kernel.org <linux-input@vger.kernel.org>
> Signed-off-by: Mauro Carvalho Chehab <mchehab@osg.samsung.com>
> Signed-off-by: Prashant Laddha <prladdha@cisco.com>
> Signed-off-by: Hans Verkuil <hans.verkuil@cisco.com>

For input bit:

Acked-by: Dmitry Torokhov <dmitry.torokhov@gmail.com>

> ---
>  drivers/input/ff-memless.c      |  18 ++++-
>  drivers/media/usb/gspca/ov534.c |  11 +--
>  include/linux/fixp-arith.h      | 145 +++++++++++++++++++++++++++++-----------
>  3 files changed, 125 insertions(+), 49 deletions(-)
> 
> diff --git a/drivers/input/ff-memless.c b/drivers/input/ff-memless.c
> index 74c0d8c..fcc6c33 100644
> --- a/drivers/input/ff-memless.c
> +++ b/drivers/input/ff-memless.c
> @@ -237,6 +237,18 @@ static u16 ml_calculate_direction(u16 direction, u16 force,
>  		(force + new_force)) << 1;
>  }
>  
> +#define FRAC_N 8
> +static inline s16 fixp_new16(s16 a)
> +{
> +	return ((s32)a) >> (16 - FRAC_N);
> +}
> +
> +static inline s16 fixp_mult(s16 a, s16 b)
> +{
> +	a = ((s32)a * 0x100) / 0x7fff;
> +	return ((s32)(a * b)) >> FRAC_N;
> +}
> +
>  /*
>   * Combine two effects and apply gain.
>   */
> @@ -247,7 +259,7 @@ static void ml_combine_effects(struct ff_effect *effect,
>  	struct ff_effect *new = state->effect;
>  	unsigned int strong, weak, i;
>  	int x, y;
> -	fixp_t level;
> +	s16 level;
>  
>  	switch (new->type) {
>  	case FF_CONSTANT:
> @@ -255,8 +267,8 @@ static void ml_combine_effects(struct ff_effect *effect,
>  		level = fixp_new16(apply_envelope(state,
>  					new->u.constant.level,
>  					&new->u.constant.envelope));
> -		x = fixp_mult(fixp_sin(i), level) * gain / 0xffff;
> -		y = fixp_mult(-fixp_cos(i), level) * gain / 0xffff;
> +		x = fixp_mult(fixp_sin16(i), level) * gain / 0xffff;
> +		y = fixp_mult(-fixp_cos16(i), level) * gain / 0xffff;
>  		/*
>  		 * here we abuse ff_ramp to hold x and y of constant force
>  		 * If in future any driver wants something else than x and y
> diff --git a/drivers/media/usb/gspca/ov534.c b/drivers/media/usb/gspca/ov534.c
> index a9c866d..146071b 100644
> --- a/drivers/media/usb/gspca/ov534.c
> +++ b/drivers/media/usb/gspca/ov534.c
> @@ -816,21 +816,16 @@ static void sethue(struct gspca_dev *gspca_dev, s32 val)
>  		s16 huesin;
>  		s16 huecos;
>  
> -		/* fixp_sin and fixp_cos accept only positive values, while
> -		 * our val is between -90 and 90
> -		 */
> -		val += 360;
> -
>  		/* According to the datasheet the registers expect HUESIN and
>  		 * HUECOS to be the result of the trigonometric functions,
>  		 * scaled by 0x80.
>  		 *
> -		 * The 0x100 here represents the maximun absolute value
> +		 * The 0x7fff here represents the maximum absolute value
>  		 * returned byt fixp_sin and fixp_cos, so the scaling will
>  		 * consider the result like in the interval [-1.0, 1.0].
>  		 */
> -		huesin = fixp_sin(val) * 0x80 / 0x100;
> -		huecos = fixp_cos(val) * 0x80 / 0x100;
> +		huesin = fixp_sin16(val) * 0x80 / 0x7fff;
> +		huecos = fixp_cos16(val) * 0x80 / 0x7fff;
>  
>  		if (huesin < 0) {
>  			sccb_reg_write(gspca_dev, 0xab,
> diff --git a/include/linux/fixp-arith.h b/include/linux/fixp-arith.h
> index 3089d73..d4686fe 100644
> --- a/include/linux/fixp-arith.h
> +++ b/include/linux/fixp-arith.h
> @@ -1,6 +1,8 @@
>  #ifndef _FIXP_ARITH_H
>  #define _FIXP_ARITH_H
>  
> +#include <linux/math64.h>
> +
>  /*
>   * Simplistic fixed-point arithmetics.
>   * Hmm, I'm probably duplicating some code :(
> @@ -29,59 +31,126 @@
>  
>  #include <linux/types.h>
>  
> -/* The type representing fixed-point values */
> -typedef s16 fixp_t;
> +static const s32 sin_table[] = {
> +	0x00000000, 0x023be165, 0x04779632, 0x06b2f1d2, 0x08edc7b6, 0x0b27eb5c,
> +	0x0d61304d, 0x0f996a26, 0x11d06c96, 0x14060b67, 0x163a1a7d, 0x186c6ddd,
> +	0x1a9cd9ac, 0x1ccb3236, 0x1ef74bf2, 0x2120fb82, 0x234815ba, 0x256c6f9e,
> +	0x278dde6e, 0x29ac379f, 0x2bc750e8, 0x2ddf003f, 0x2ff31bdd, 0x32037a44,
> +	0x340ff241, 0x36185aee, 0x381c8bb5, 0x3a1c5c56, 0x3c17a4e7, 0x3e0e3ddb,
> +	0x3fffffff, 0x41ecc483, 0x43d464fa, 0x45b6bb5d, 0x4793a20f, 0x496af3e1,
> +	0x4b3c8c11, 0x4d084650, 0x4ecdfec6, 0x508d9210, 0x5246dd48, 0x53f9be04,
> +	0x55a6125a, 0x574bb8e5, 0x58ea90c2, 0x5a827999, 0x5c135399, 0x5d9cff82,
> +	0x5f1f5ea0, 0x609a52d1, 0x620dbe8a, 0x637984d3, 0x64dd894f, 0x6639b039,
> +	0x678dde6d, 0x68d9f963, 0x6a1de735, 0x6b598ea1, 0x6c8cd70a, 0x6db7a879,
> +	0x6ed9eba0, 0x6ff389de, 0x71046d3c, 0x720c8074, 0x730baeec, 0x7401e4bf,
> +	0x74ef0ebb, 0x75d31a5f, 0x76adf5e5, 0x777f903b, 0x7847d908, 0x7906c0af,
> +	0x79bc384c, 0x7a6831b8, 0x7b0a9f8c, 0x7ba3751c, 0x7c32a67c, 0x7cb82884,
> +	0x7d33f0c8, 0x7da5f5a3, 0x7e0e2e31, 0x7e6c924f, 0x7ec11aa3, 0x7f0bc095,
> +	0x7f4c7e52, 0x7f834ecf, 0x7fb02dc4, 0x7fd317b3, 0x7fec09e1, 0x7ffb025e,
> +	0x7fffffff
> +};
>  
> -#define FRAC_N 8
> -#define FRAC_MASK ((1<<FRAC_N)-1)
> +/**
> + * __fixp_sin32() returns the sin of an angle in degrees
> + *
> + * @degrees: angle, in degrees, from 0 to 360.
> + *
> + * The returned value ranges from -0x7fffffff to +0x7fffffff.
> + */
> +static inline s32 __fixp_sin32(int degrees)
> +{
> +	s32 ret;
> +	bool negative = false;
>  
> -/* Not to be used directly. Use fixp_{cos,sin} */
> -static const fixp_t cos_table[46] = {
> -	0x0100,	0x00FF,	0x00FF,	0x00FE,	0x00FD,	0x00FC,	0x00FA,	0x00F8,
> -	0x00F6,	0x00F3,	0x00F0,	0x00ED,	0x00E9,	0x00E6,	0x00E2,	0x00DD,
> -	0x00D9,	0x00D4,	0x00CF,	0x00C9,	0x00C4,	0x00BE,	0x00B8,	0x00B1,
> -	0x00AB,	0x00A4,	0x009D,	0x0096,	0x008F,	0x0087,	0x0080,	0x0078,
> -	0x0070,	0x0068,	0x005F,	0x0057,	0x004F,	0x0046,	0x003D,	0x0035,
> -	0x002C,	0x0023,	0x001A,	0x0011,	0x0008, 0x0000
> -};
> +	if (degrees > 180) {
> +		negative = true;
> +		degrees -= 180;
> +	}
> +	if (degrees > 90)
> +		degrees = 180 - degrees;
>  
> +	ret = sin_table[degrees];
>  
> -/* a: 123 -> 123.0 */
> -static inline fixp_t fixp_new(s16 a)
> -{
> -	return a<<FRAC_N;
> +	return negative ? -ret : ret;
>  }
>  
> -/* a: 0xFFFF -> -1.0
> -      0x8000 -> 1.0
> -      0x0000 -> 0.0
> -*/
> -static inline fixp_t fixp_new16(s16 a)
> +/**
> + * fixp_sin32() returns the sin of an angle in degrees
> + *
> + * @degrees: angle, in degrees. The angle can be positive or negative
> + *
> + * The returned value ranges from -0x7fffffff to +0x7fffffff.
> + */
> +static inline s32 fixp_sin32(int degrees)
>  {
> -	return ((s32)a)>>(16-FRAC_N);
> +	degrees = (degrees % 360 + 360) % 360;
> +
> +	return __fixp_sin32(degrees);
>  }
>  
> -static inline fixp_t fixp_cos(unsigned int degrees)
> +/* cos(x) = sin(x + 90 degrees) */
> +#define fixp_cos32(v) fixp_sin32((v) + 90)
> +
> +/*
> + * 16 bits variants
> + *
> + * The returned value ranges from -0x7fff to 0x7fff
> + */
> +
> +#define fixp_sin16(v) (fixp_sin32(v) >> 16)
> +#define fixp_cos16(v) (fixp_cos32(v) >> 16)
> +
> +/**
> + * fixp_sin32_rad() - calculates the sin of an angle in radians
> + *
> + * @radians: angle, in radians
> + * @twopi: value to be used for 2*pi
> + *
> + * Provides a variant for the cases where just 360
> + * values is not enough. This function uses linear
> + * interpolation to a wider range of values given by
> + * twopi var.
> + *
> + * Experimental tests gave a maximum difference of
> + * 0.000038 between the value calculated by sin() and
> + * the one produced by this function, when twopi is
> + * equal to 360000. That seems to be enough precision
> + * for practical purposes.
> + *
> + * Please notice that two high numbers for twopi could cause
> + * overflows, so the routine will not allow values of twopi
> + * bigger than 1^18.
> + */
> +static inline s32 fixp_sin32_rad(u32 radians, u32 twopi)
>  {
> -	int quadrant = (degrees / 90) & 3;
> -	unsigned int i = degrees % 90;
> +	int degrees;
> +	s32 v1, v2, dx, dy;
> +	s64 tmp;
>  
> -	if (quadrant == 1 || quadrant == 3)
> -		i = 90 - i;
> +	/*
> +	 * Avoid too large values for twopi, as we don't want overflows.
> +	 */
> +	BUG_ON(twopi > 1 << 18);
>  
> -	i >>= 1;
> +	degrees = (radians * 360) / twopi;
> +	tmp = radians - (degrees * twopi) / 360;
>  
> -	return (quadrant == 1 || quadrant == 2)? -cos_table[i] : cos_table[i];
> -}
> +	degrees = (degrees % 360 + 360) % 360;
> +	v1 = __fixp_sin32(degrees);
>  
> -static inline fixp_t fixp_sin(unsigned int degrees)
> -{
> -	return -fixp_cos(degrees + 90);
> -}
> +	v2 = fixp_sin32(degrees + 1);
>  
> -static inline fixp_t fixp_mult(fixp_t a, fixp_t b)
> -{
> -	return ((s32)(a*b))>>FRAC_N;
> +	dx = twopi / 360;
> +	dy = v2 - v1;
> +
> +	tmp *= dy;
> +
> +	return v1 +  div_s64(tmp, dx);
>  }
>  
> +/* cos(x) = sin(x + pi/2 radians) */
> +
> +#define fixp_cos32_rad(rad, twopi)	\
> +	fixp_sin32_rad(rad + twopi / 4, twopi)
> +
>  #endif
> -- 
> 2.1.4
> 

-- 
Dmitry