pxfm 0.1.28

/*
 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::polyeval::{f_estrin_polyeval7, f_polyeval3, f_polyeval4, f_polyeval6};
use crate::sin_cosf::argument_reduction::ArgumentReducer;
use crate::sin_cosf::argument_reduction_pi::ArgumentReducerPi;

/**
Generated in SageMath:

```text
print("[")
for k in range(64):
    k = RealField(150)(k) * RealField(150).pi() / RealField(150)(32)
    print(double_to_hex(k.sin()) + ",")
print("];")
```
**/
static SIN_K_PI_OVER32: [u64; 64] = [
    0x0000000000000000,
    0x3fb917a6bc29b42c,
    0x3fc8f8b83c69a60b,
    0x3fd294062ed59f06,
    0x3fd87de2a6aea963,
    0x3fde2b5d3806f63b,
    0x3fe1c73b39ae68c8,
    0x3fe44cf325091dd6,
    0x3fe6a09e667f3bcd,
    0x3fe8bc806b151741,
    0x3fea9b66290ea1a3,
    0x3fec38b2f180bdb1,
    0x3fed906bcf328d46,
    0x3fee9f4156c62dda,
    0x3fef6297cff75cb0,
    0x3fefd88da3d12526,
    0x3ff0000000000000,
    0x3fefd88da3d12526,
    0x3fef6297cff75cb0,
    0x3fee9f4156c62dda,
    0x3fed906bcf328d46,
    0x3fec38b2f180bdb1,
    0x3fea9b66290ea1a3,
    0x3fe8bc806b151741,
    0x3fe6a09e667f3bcd,
    0x3fe44cf325091dd6,
    0x3fe1c73b39ae68c8,
    0x3fde2b5d3806f63b,
    0x3fd87de2a6aea963,
    0x3fd294062ed59f06,
    0x3fc8f8b83c69a60b,
    0x3fb917a6bc29b42c,
    0xb69f77598338bfdf,
    0xbfb917a6bc29b42c,
    0xbfc8f8b83c69a60b,
    0xbfd294062ed59f06,
    0xbfd87de2a6aea963,
    0xbfde2b5d3806f63b,
    0xbfe1c73b39ae68c8,
    0xbfe44cf325091dd6,
    0xbfe6a09e667f3bcd,
    0xbfe8bc806b151741,
    0xbfea9b66290ea1a3,
    0xbfec38b2f180bdb1,
    0xbfed906bcf328d46,
    0xbfee9f4156c62dda,
    0xbfef6297cff75cb0,
    0xbfefd88da3d12526,
    0xbff0000000000000,
    0xbfefd88da3d12526,
    0xbfef6297cff75cb0,
    0xbfee9f4156c62dda,
    0xbfed906bcf328d46,
    0xbfec38b2f180bdb1,
    0xbfea9b66290ea1a3,
    0xbfe8bc806b151741,
    0xbfe6a09e667f3bcd,
    0xbfe44cf325091dd6,
    0xbfe1c73b39ae68c8,
    0xbfde2b5d3806f63b,
    0xbfd87de2a6aea963,
    0xbfd294062ed59f06,
    0xbfc8f8b83c69a60b,
    0xbfb917a6bc29b42c,
];

pub(crate) struct SinCosf {
    pub(crate) sin_k: f64,
    pub(crate) cos_k: f64,
    pub(crate) sin_y: f64,
    pub(crate) cosm1_y: f64,
}

#[inline]
pub(crate) fn sincosf_eval(x: f64, x_abs: u32) -> SinCosf {
    let argument_reduction = ArgumentReducer { x, x_abs };

    let (y, k) = argument_reduction.reduce();

    let y_sqr = y * y;

    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
    // So k is an integer and -0.5 <= y <= 0.5.
    // Then sin(x) = sin((k + y)*pi/32)
    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)

    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
    // cos_k = cos(k * pi/32)
    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);

    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
    // with:
    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
    //
    // See ./notes/sincosf.sollya
    let sin_y = y * f_polyeval4(
        y_sqr,
        f64::from_bits(0x3fb921fb54442d18),
        f64::from_bits(0xbf24abbce625abb1),
        f64::from_bits(0x3e7466bc624f2776),
        f64::from_bits(0xbdb32c3a619d4a7e),
    );
    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
    // Note that cosm1_y = cos(y*pi/32) - 1.
    //
    // See ./notes/sincosf.sollya
    let cosm1_y = y_sqr
        * f_polyeval3(
            y_sqr,
            f64::from_bits(0xbf73bd3cc9be430b),
            f64::from_bits(0x3ed03c1f070c2e27),
            f64::from_bits(0xbe155cc84bd94200),
        );

    SinCosf {
        sin_k,
        cos_k,
        sin_y,
        cosm1_y,
    }
}

#[inline(always)]
#[allow(unused)]
pub(crate) fn sincosf_eval_fma(x: f64, x_abs: u32) -> SinCosf {
    let argument_reduction = ArgumentReducer { x, x_abs };

    let (y, k) = argument_reduction.reduce_fma();

    let y_sqr = y * y;

    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
    // So k is an integer and -0.5 <= y <= 0.5.
    // Then sin(x) = sin((k + y)*pi/32)
    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)

    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
    // cos_k = cos(k * pi/32)
    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);

    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
    // with:
    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
    //
    // See ./notes/sincosf.sollya
    use crate::polyeval::d_polyeval4;
    let sin_y = y * d_polyeval4(
        y_sqr,
        f64::from_bits(0x3fb921fb54442d18),
        f64::from_bits(0xbf24abbce625abb1),
        f64::from_bits(0x3e7466bc624f2776),
        f64::from_bits(0xbdb32c3a619d4a7e),
    );
    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
    // Note that cosm1_y = cos(y*pi/32) - 1.
    //
    // See ./notes/sincosf.sollya
    use crate::polyeval::d_polyeval3;
    let cosm1_y = y_sqr
        * d_polyeval3(
            y_sqr,
            f64::from_bits(0xbf73bd3cc9be430b),
            f64::from_bits(0x3ed03c1f070c2e27),
            f64::from_bits(0xbe155cc84bd94200),
        );

    SinCosf {
        sin_k,
        cos_k,
        sin_y,
        cosm1_y,
    }
}

#[inline(always)]
pub(crate) fn sincospif_eval(x: f64) -> SinCosf {
    let argument_reduction = ArgumentReducerPi { x };

    let (y, k) = argument_reduction.reduce();

    let y_sqr = y * y;

    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
    // So k is an integer and -0.5 <= y <= 0.5.
    // Then sin(x) = sin((k + y)*pi/32)
    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)

    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
    // cos_k = cos(k * pi/32)
    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);

    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
    // with:
    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
    //
    // See ./notes/sincosf.sollya
    let sin_y = y * f_polyeval4(
        y_sqr,
        f64::from_bits(0x3fb921fb54442d18),
        f64::from_bits(0xbf24abbce625abb1),
        f64::from_bits(0x3e7466bc624f2776),
        f64::from_bits(0xbdb32c3a619d4a7e),
    );
    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
    // Note that cosm1_y = cos(y*pi/32) - 1.
    //
    // See ./notes/sincosf.sollya
    let cosm1_y = y_sqr
        * f_polyeval3(
            y_sqr,
            f64::from_bits(0xbf73bd3cc9be430b),
            f64::from_bits(0x3ed03c1f070c2e27),
            f64::from_bits(0xbe155cc84bd94200),
        );

    SinCosf {
        sin_k,
        cos_k,
        sin_y,
        cosm1_y,
    }
}

#[inline(always)]
#[allow(unused)]
pub(crate) fn sincospif_eval_fma(x: f64) -> SinCosf {
    let argument_reduction = ArgumentReducerPi { x };

    let (y, k) = argument_reduction.reduce_fma();

    let y_sqr = y * y;

    // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
    // So k is an integer and -0.5 <= y <= 0.5.
    // Then sin(x) = sin((k + y)*pi/32)
    //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)

    let sin_k = f64::from_bits(SIN_K_PI_OVER32[((k as u64) & 63) as usize]);
    // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
    // cos_k = cos(k * pi/32)
    let cos_k = f64::from_bits(SIN_K_PI_OVER32[((k.wrapping_add(16) as u64) & 63) as usize]);

    use crate::polyeval::{d_polyeval3, d_polyeval4};

    // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
    // with:
    // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
    //
    // See ./notes/sincosf.sollya
    let sin_y = y * d_polyeval4(
        y_sqr,
        f64::from_bits(0x3fb921fb54442d18),
        f64::from_bits(0xbf24abbce625abb1),
        f64::from_bits(0x3e7466bc624f2776),
        f64::from_bits(0xbdb32c3a619d4a7e),
    );
    // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
    // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
    // Note that cosm1_y = cos(y*pi/32) - 1.
    //
    // See ./notes/sincosf.sollya
    let cosm1_y = y_sqr
        * d_polyeval3(
            y_sqr,
            f64::from_bits(0xbf73bd3cc9be430b),
            f64::from_bits(0x3ed03c1f070c2e27),
            f64::from_bits(0xbe155cc84bd94200),
        );

    SinCosf {
        sin_k,
        cos_k,
        sin_y,
        cosm1_y,
    }
}

/**
Poly for sin(pi*y) on [-0.25; 0.25].
Generated by Sollya:
```text
d = [0,0.25];
f_sin = sin(y*pi)/y;
Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|D...|], d, relative, floating);
```
See ./notes/sinpif_0p25.sollya
**/
#[inline(always)]
pub(crate) fn sinpif_eval(y: f64) -> f64 {
    let y2 = y * y;
    f_polyeval6(
        y2,
        f64::from_bits(0x400921fb54442cf8),
        f64::from_bits(0xc014abbce6257778),
        f64::from_bits(0x400466bc670fa464),
        f64::from_bits(0xbfe32d2c627f47da),
        f64::from_bits(0x3fb5071be0a2d3da),
        f64::from_bits(0xbf7dd4e0ae5b9fcd),
    ) * y
}

/**
Poly for sin(pi*y) on [-0.25; 0.25].
Generated by Sollya:
```text
d = [0,0.25];
f_sin = sin(y*pi)/y;
Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10|], [|D...|], d, relative, floating);
```
See ./notes/sinpif_0p25.sollya
**/
#[inline(always)]
#[allow(unused)]
pub(crate) fn sinpif_eval_fma(y: f64) -> f64 {
    let y2 = y * y;
    use crate::polyeval::d_polyeval6;
    d_polyeval6(
        y2,
        f64::from_bits(0x400921fb54442cf8),
        f64::from_bits(0xc014abbce6257778),
        f64::from_bits(0x400466bc670fa464),
        f64::from_bits(0xbfe32d2c627f47da),
        f64::from_bits(0x3fb5071be0a2d3da),
        f64::from_bits(0xbf7dd4e0ae5b9fcd),
    ) * y
}

/**
Poly for sin(pi*y) on [-0.25; 0.25].
Generated by Sollya:
```text
d = [0,0.25];
f_sin = sin(y*pi)/y;
Q = fpminimax(f_sin, [|0, 2, 4, 6, 8, 10, 12|], [|D...|], d, relative, floating);
```
See ./notes/sinpif_0p25_2.sollya
**/
#[inline]
pub(crate) fn sinpif_eval2(y: f64) -> f64 {
    let y2 = y * y;
    f_estrin_polyeval7(
        y2,
        f64::from_bits(0x400921fb54442d18),
        f64::from_bits(0xc014abbce625be09),
        f64::from_bits(0x400466bc67754fff),
        f64::from_bits(0xbfe32d2ccdfe9424),
        f64::from_bits(0x3fb50782d5f14825),
        f64::from_bits(0xbf7e2fe76fdffd2b),
        f64::from_bits(0x3f3e357ef99eb0bb),
    ) * y
}

/**
Poly for cos(pi*y) on [-0.25; 0.25].
Generated by Sollya:
```text
d = [0, 0.25];
f_cos = cos(y*pi);
Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|1, D...|], d, relative, floating);
```
See ./notes/cospif_0p25.sollya
**/
#[inline(always)]
pub(crate) fn cospif_eval(y: f64) -> f64 {
    let y2 = y * y;
    f_estrin_polyeval7(
        y2,
        f64::from_bits(0x3ff0000000000000),
        f64::from_bits(0xc013bd3cc9be459b),
        f64::from_bits(0x40103c1f081b0c98),
        f64::from_bits(0xbff55d3c7dc25308),
        f64::from_bits(0x3fce1f4fb6e12563),
        f64::from_bits(0xbf9a6c9c101c1182),
        f64::from_bits(0x3f5f3d9faffda250),
    )
}

/**
Poly for cos(pi*y) on [-0.25; 0.25].
Generated by Sollya:
```text
d = [0, 0.25];
f_cos = cos(y*pi);
Q = fpminimax(f_cos, [|0, 2, 4, 6, 8, 10|], [|1, D...|], d, relative, floating);
```
See ./notes/cospif_0p25.sollya
**/
#[inline(always)]
#[allow(unused)]
pub(crate) fn cospif_eval_fma(y: f64) -> f64 {
    let y2 = y * y;
    use crate::polyeval::d_estrin_polyeval7;
    d_estrin_polyeval7(
        y2,
        f64::from_bits(0x3ff0000000000000),
        f64::from_bits(0xc013bd3cc9be459b),
        f64::from_bits(0x40103c1f081b0c98),
        f64::from_bits(0xbff55d3c7dc25308),
        f64::from_bits(0x3fce1f4fb6e12563),
        f64::from_bits(0xbf9a6c9c101c1182),
        f64::from_bits(0x3f5f3d9faffda250),
    )
}