mirror of
https://github.com/torvalds/linux
synced 2024-09-23 21:07:52 +00:00
3d0d14f983
lindent these files:
errors lines of code errors/KLOC
arch/x86/math-emu/ 2236 9424 237.2
arch/x86/math-emu/ 128 8706 14.7
no other changes. No code changed:
text data bss dec hex filename
5589802 612739 3833856 10036397 9924ad vmlinux.before
5589802 612739 3833856 10036397 9924ad vmlinux.after
the intent of this patch is to ease the automated tracking of kernel
code quality - it's just much easier for us to maintain it if every file
in arch/x86 is supposed to be clean.
NOTE: it is a known problem of lindent that it causes some style damage
of its own, but it's a safe tool (well, except for the gcc array range
initializers extension), so we did the bulk of the changes via lindent,
and did the manual fixups in a followup patch.
the resulting math-emu code has been tested by Thomas Gleixner on a real
386 DX CPU as well, and it works fine.
Signed-off-by: Ingo Molnar <mingo@elte.hu>
Signed-off-by: Thomas Gleixner <tglx@linutronix.de>
379 lines
11 KiB
C
379 lines
11 KiB
C
/*---------------------------------------------------------------------------+
|
|
| poly_sin.c |
|
|
| |
|
|
| Computation of an approximation of the sin function and the cosine |
|
|
| function by a polynomial. |
|
|
| |
|
|
| Copyright (C) 1992,1993,1994,1997,1999 |
|
|
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
|
|
| E-mail billm@melbpc.org.au |
|
|
| |
|
|
| |
|
|
+---------------------------------------------------------------------------*/
|
|
|
|
#include "exception.h"
|
|
#include "reg_constant.h"
|
|
#include "fpu_emu.h"
|
|
#include "fpu_system.h"
|
|
#include "control_w.h"
|
|
#include "poly.h"
|
|
|
|
#define N_COEFF_P 4
|
|
#define N_COEFF_N 4
|
|
|
|
static const unsigned long long pos_terms_l[N_COEFF_P] = {
|
|
0xaaaaaaaaaaaaaaabLL,
|
|
0x00d00d00d00cf906LL,
|
|
0x000006b99159a8bbLL,
|
|
0x000000000d7392e6LL
|
|
};
|
|
|
|
static const unsigned long long neg_terms_l[N_COEFF_N] = {
|
|
0x2222222222222167LL,
|
|
0x0002e3bc74aab624LL,
|
|
0x0000000b09229062LL,
|
|
0x00000000000c7973LL
|
|
};
|
|
|
|
#define N_COEFF_PH 4
|
|
#define N_COEFF_NH 4
|
|
static const unsigned long long pos_terms_h[N_COEFF_PH] = {
|
|
0x0000000000000000LL,
|
|
0x05b05b05b05b0406LL,
|
|
0x000049f93edd91a9LL,
|
|
0x00000000c9c9ed62LL
|
|
};
|
|
|
|
static const unsigned long long neg_terms_h[N_COEFF_NH] = {
|
|
0xaaaaaaaaaaaaaa98LL,
|
|
0x001a01a01a019064LL,
|
|
0x0000008f76c68a77LL,
|
|
0x0000000000d58f5eLL
|
|
};
|
|
|
|
/*--- poly_sine() -----------------------------------------------------------+
|
|
| |
|
|
+---------------------------------------------------------------------------*/
|
|
void poly_sine(FPU_REG * st0_ptr)
|
|
{
|
|
int exponent, echange;
|
|
Xsig accumulator, argSqrd, argTo4;
|
|
unsigned long fix_up, adj;
|
|
unsigned long long fixed_arg;
|
|
FPU_REG result;
|
|
|
|
exponent = exponent(st0_ptr);
|
|
|
|
accumulator.lsw = accumulator.midw = accumulator.msw = 0;
|
|
|
|
/* Split into two ranges, for arguments below and above 1.0 */
|
|
/* The boundary between upper and lower is approx 0.88309101259 */
|
|
if ((exponent < -1)
|
|
|| ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa))) {
|
|
/* The argument is <= 0.88309101259 */
|
|
|
|
argSqrd.msw = st0_ptr->sigh;
|
|
argSqrd.midw = st0_ptr->sigl;
|
|
argSqrd.lsw = 0;
|
|
mul64_Xsig(&argSqrd, &significand(st0_ptr));
|
|
shr_Xsig(&argSqrd, 2 * (-1 - exponent));
|
|
argTo4.msw = argSqrd.msw;
|
|
argTo4.midw = argSqrd.midw;
|
|
argTo4.lsw = argSqrd.lsw;
|
|
mul_Xsig_Xsig(&argTo4, &argTo4);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
|
|
N_COEFF_N - 1);
|
|
mul_Xsig_Xsig(&accumulator, &argSqrd);
|
|
negate_Xsig(&accumulator);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
|
|
N_COEFF_P - 1);
|
|
|
|
shr_Xsig(&accumulator, 2); /* Divide by four */
|
|
accumulator.msw |= 0x80000000; /* Add 1.0 */
|
|
|
|
mul64_Xsig(&accumulator, &significand(st0_ptr));
|
|
mul64_Xsig(&accumulator, &significand(st0_ptr));
|
|
mul64_Xsig(&accumulator, &significand(st0_ptr));
|
|
|
|
/* Divide by four, FPU_REG compatible, etc */
|
|
exponent = 3 * exponent;
|
|
|
|
/* The minimum exponent difference is 3 */
|
|
shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
|
|
|
|
negate_Xsig(&accumulator);
|
|
XSIG_LL(accumulator) += significand(st0_ptr);
|
|
|
|
echange = round_Xsig(&accumulator);
|
|
|
|
setexponentpos(&result, exponent(st0_ptr) + echange);
|
|
} else {
|
|
/* The argument is > 0.88309101259 */
|
|
/* We use sin(st(0)) = cos(pi/2-st(0)) */
|
|
|
|
fixed_arg = significand(st0_ptr);
|
|
|
|
if (exponent == 0) {
|
|
/* The argument is >= 1.0 */
|
|
|
|
/* Put the binary point at the left. */
|
|
fixed_arg <<= 1;
|
|
}
|
|
/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
|
|
fixed_arg = 0x921fb54442d18469LL - fixed_arg;
|
|
/* There is a special case which arises due to rounding, to fix here. */
|
|
if (fixed_arg == 0xffffffffffffffffLL)
|
|
fixed_arg = 0;
|
|
|
|
XSIG_LL(argSqrd) = fixed_arg;
|
|
argSqrd.lsw = 0;
|
|
mul64_Xsig(&argSqrd, &fixed_arg);
|
|
|
|
XSIG_LL(argTo4) = XSIG_LL(argSqrd);
|
|
argTo4.lsw = argSqrd.lsw;
|
|
mul_Xsig_Xsig(&argTo4, &argTo4);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
|
|
N_COEFF_NH - 1);
|
|
mul_Xsig_Xsig(&accumulator, &argSqrd);
|
|
negate_Xsig(&accumulator);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
|
|
N_COEFF_PH - 1);
|
|
negate_Xsig(&accumulator);
|
|
|
|
mul64_Xsig(&accumulator, &fixed_arg);
|
|
mul64_Xsig(&accumulator, &fixed_arg);
|
|
|
|
shr_Xsig(&accumulator, 3);
|
|
negate_Xsig(&accumulator);
|
|
|
|
add_Xsig_Xsig(&accumulator, &argSqrd);
|
|
|
|
shr_Xsig(&accumulator, 1);
|
|
|
|
accumulator.lsw |= 1; /* A zero accumulator here would cause problems */
|
|
negate_Xsig(&accumulator);
|
|
|
|
/* The basic computation is complete. Now fix the answer to
|
|
compensate for the error due to the approximation used for
|
|
pi/2
|
|
*/
|
|
|
|
/* This has an exponent of -65 */
|
|
fix_up = 0x898cc517;
|
|
/* The fix-up needs to be improved for larger args */
|
|
if (argSqrd.msw & 0xffc00000) {
|
|
/* Get about 32 bit precision in these: */
|
|
fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
|
|
}
|
|
fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
|
|
|
|
adj = accumulator.lsw; /* temp save */
|
|
accumulator.lsw -= fix_up;
|
|
if (accumulator.lsw > adj)
|
|
XSIG_LL(accumulator)--;
|
|
|
|
echange = round_Xsig(&accumulator);
|
|
|
|
setexponentpos(&result, echange - 1);
|
|
}
|
|
|
|
significand(&result) = XSIG_LL(accumulator);
|
|
setsign(&result, getsign(st0_ptr));
|
|
FPU_copy_to_reg0(&result, TAG_Valid);
|
|
|
|
#ifdef PARANOID
|
|
if ((exponent(&result) >= 0)
|
|
&& (significand(&result) > 0x8000000000000000LL)) {
|
|
EXCEPTION(EX_INTERNAL | 0x150);
|
|
}
|
|
#endif /* PARANOID */
|
|
|
|
}
|
|
|
|
/*--- poly_cos() ------------------------------------------------------------+
|
|
| |
|
|
+---------------------------------------------------------------------------*/
|
|
void poly_cos(FPU_REG * st0_ptr)
|
|
{
|
|
FPU_REG result;
|
|
long int exponent, exp2, echange;
|
|
Xsig accumulator, argSqrd, fix_up, argTo4;
|
|
unsigned long long fixed_arg;
|
|
|
|
#ifdef PARANOID
|
|
if ((exponent(st0_ptr) > 0)
|
|
|| ((exponent(st0_ptr) == 0)
|
|
&& (significand(st0_ptr) > 0xc90fdaa22168c234LL))) {
|
|
EXCEPTION(EX_Invalid);
|
|
FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
|
|
return;
|
|
}
|
|
#endif /* PARANOID */
|
|
|
|
exponent = exponent(st0_ptr);
|
|
|
|
accumulator.lsw = accumulator.midw = accumulator.msw = 0;
|
|
|
|
if ((exponent < -1)
|
|
|| ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54))) {
|
|
/* arg is < 0.687705 */
|
|
|
|
argSqrd.msw = st0_ptr->sigh;
|
|
argSqrd.midw = st0_ptr->sigl;
|
|
argSqrd.lsw = 0;
|
|
mul64_Xsig(&argSqrd, &significand(st0_ptr));
|
|
|
|
if (exponent < -1) {
|
|
/* shift the argument right by the required places */
|
|
shr_Xsig(&argSqrd, 2 * (-1 - exponent));
|
|
}
|
|
|
|
argTo4.msw = argSqrd.msw;
|
|
argTo4.midw = argSqrd.midw;
|
|
argTo4.lsw = argSqrd.lsw;
|
|
mul_Xsig_Xsig(&argTo4, &argTo4);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
|
|
N_COEFF_NH - 1);
|
|
mul_Xsig_Xsig(&accumulator, &argSqrd);
|
|
negate_Xsig(&accumulator);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
|
|
N_COEFF_PH - 1);
|
|
negate_Xsig(&accumulator);
|
|
|
|
mul64_Xsig(&accumulator, &significand(st0_ptr));
|
|
mul64_Xsig(&accumulator, &significand(st0_ptr));
|
|
shr_Xsig(&accumulator, -2 * (1 + exponent));
|
|
|
|
shr_Xsig(&accumulator, 3);
|
|
negate_Xsig(&accumulator);
|
|
|
|
add_Xsig_Xsig(&accumulator, &argSqrd);
|
|
|
|
shr_Xsig(&accumulator, 1);
|
|
|
|
/* It doesn't matter if accumulator is all zero here, the
|
|
following code will work ok */
|
|
negate_Xsig(&accumulator);
|
|
|
|
if (accumulator.lsw & 0x80000000)
|
|
XSIG_LL(accumulator)++;
|
|
if (accumulator.msw == 0) {
|
|
/* The result is 1.0 */
|
|
FPU_copy_to_reg0(&CONST_1, TAG_Valid);
|
|
return;
|
|
} else {
|
|
significand(&result) = XSIG_LL(accumulator);
|
|
|
|
/* will be a valid positive nr with expon = -1 */
|
|
setexponentpos(&result, -1);
|
|
}
|
|
} else {
|
|
fixed_arg = significand(st0_ptr);
|
|
|
|
if (exponent == 0) {
|
|
/* The argument is >= 1.0 */
|
|
|
|
/* Put the binary point at the left. */
|
|
fixed_arg <<= 1;
|
|
}
|
|
/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
|
|
fixed_arg = 0x921fb54442d18469LL - fixed_arg;
|
|
/* There is a special case which arises due to rounding, to fix here. */
|
|
if (fixed_arg == 0xffffffffffffffffLL)
|
|
fixed_arg = 0;
|
|
|
|
exponent = -1;
|
|
exp2 = -1;
|
|
|
|
/* A shift is needed here only for a narrow range of arguments,
|
|
i.e. for fixed_arg approx 2^-32, but we pick up more... */
|
|
if (!(LL_MSW(fixed_arg) & 0xffff0000)) {
|
|
fixed_arg <<= 16;
|
|
exponent -= 16;
|
|
exp2 -= 16;
|
|
}
|
|
|
|
XSIG_LL(argSqrd) = fixed_arg;
|
|
argSqrd.lsw = 0;
|
|
mul64_Xsig(&argSqrd, &fixed_arg);
|
|
|
|
if (exponent < -1) {
|
|
/* shift the argument right by the required places */
|
|
shr_Xsig(&argSqrd, 2 * (-1 - exponent));
|
|
}
|
|
|
|
argTo4.msw = argSqrd.msw;
|
|
argTo4.midw = argSqrd.midw;
|
|
argTo4.lsw = argSqrd.lsw;
|
|
mul_Xsig_Xsig(&argTo4, &argTo4);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
|
|
N_COEFF_N - 1);
|
|
mul_Xsig_Xsig(&accumulator, &argSqrd);
|
|
negate_Xsig(&accumulator);
|
|
|
|
polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
|
|
N_COEFF_P - 1);
|
|
|
|
shr_Xsig(&accumulator, 2); /* Divide by four */
|
|
accumulator.msw |= 0x80000000; /* Add 1.0 */
|
|
|
|
mul64_Xsig(&accumulator, &fixed_arg);
|
|
mul64_Xsig(&accumulator, &fixed_arg);
|
|
mul64_Xsig(&accumulator, &fixed_arg);
|
|
|
|
/* Divide by four, FPU_REG compatible, etc */
|
|
exponent = 3 * exponent;
|
|
|
|
/* The minimum exponent difference is 3 */
|
|
shr_Xsig(&accumulator, exp2 - exponent);
|
|
|
|
negate_Xsig(&accumulator);
|
|
XSIG_LL(accumulator) += fixed_arg;
|
|
|
|
/* The basic computation is complete. Now fix the answer to
|
|
compensate for the error due to the approximation used for
|
|
pi/2
|
|
*/
|
|
|
|
/* This has an exponent of -65 */
|
|
XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
|
|
fix_up.lsw = 0;
|
|
|
|
/* The fix-up needs to be improved for larger args */
|
|
if (argSqrd.msw & 0xffc00000) {
|
|
/* Get about 32 bit precision in these: */
|
|
fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
|
|
fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
|
|
}
|
|
|
|
exp2 += norm_Xsig(&accumulator);
|
|
shr_Xsig(&accumulator, 1); /* Prevent overflow */
|
|
exp2++;
|
|
shr_Xsig(&fix_up, 65 + exp2);
|
|
|
|
add_Xsig_Xsig(&accumulator, &fix_up);
|
|
|
|
echange = round_Xsig(&accumulator);
|
|
|
|
setexponentpos(&result, exp2 + echange);
|
|
significand(&result) = XSIG_LL(accumulator);
|
|
}
|
|
|
|
FPU_copy_to_reg0(&result, TAG_Valid);
|
|
|
|
#ifdef PARANOID
|
|
if ((exponent(&result) >= 0)
|
|
&& (significand(&result) > 0x8000000000000000LL)) {
|
|
EXCEPTION(EX_INTERNAL | 0x151);
|
|
}
|
|
#endif /* PARANOID */
|
|
|
|
}
|