consider changes for sincosd

[stonylohr]

Some changes to consider for geomath::sincosd, starting with the ones that seem most solid.

Calculate "s" and "c" using sin_cos rather than sin and cos: let (s, c) = r.sin_cos();

It should be possible to declare "q" non-mut, and skip the later changes to it, because the initial calculation should already result in the desired value range.

Use native "%" operator instead of fmod.

Consider using "match" rather than "if/else if" when checking "q":

let (s, c) = match (q as usize) & 3usize {
    0usize => { (s, c) },
    1usize => { (c, -s) },
    2usize => { (-s, -c) },
    _ => { (-c, s) } // 3
};

Contrary to what I said in #3, this is a place where I failed to find a clean way around shadowing or an equivalent. It should be possible to remove it if and when support for destructuring assignment is added. I also can't remember whether I confirmed that all my casts to usize are needed, or whether that's just paranoia on my part.

Remove the second round of shadowing for the "s" and "c" variables.

Consider replacing "90.0 q as f64" with "(90 q) as f64" to favor integer math where practical.

Review any possible difference between "(r / 90.0 + (0.5)).floor()" and "(r / 90.0).round()". I think that both are available in C++, and Charles Karney chose to use floor over round. Review should make sure to include both positive and negative boundary cases (like 0.5 and -0.5). Behavior would be my main concern, but also performance.

Review handling of NaN and infinite inputs, to make sure behavior is consistent with other versions of geographiclib.

Consider adding a comment about future use of remquo. The C++ version uses this when available, so I expect it allows some performance improvement. The nagisa crate currently has this on a todo list.

Putting it all together, it might look something like this...

pub fn sincosd(x: f64) -> (f64, f64) {
    // In order to minimize round-off errors, this function exactly reduces
    // the argument to the range [-45, 45] before converting it to radians.
    // Consider using remquo here, if and when it becomes available. See C++ version for details. --SL, 2020/07/17
    let mut r = if x.is_finite() {
        x % 360.0
    } else {
        f64::NAN
    };

    let q = if r.is_nan() {
        0
    } else {
        (r / 90.0).round() as isize
    };

    r -= (90 * q) as f64;
    r = r.to_radians();

    let (s, c) = r.sin_cos();
    // Modify to avoid shadowing, if and when destructuring assignment becomes available. --SL, 2020/07/17
    let (mut s, mut c) = match (q as usize) & 3usize {
        0usize => { (s, c) },
        1usize => { (c, -s) },
        2usize => { (-s, -c) },
        _ => { (-c, s) } // 3
    };

    // # Remove the minus sign on -0.0 except for sin(-0.0).
    // # See also Math::sincosd in the C++ library.
    // # AngNormalize has a similar fix.
    //     s, c = (x, c) if x == 0 else (0.0+s, 0.0+c)
    // return s, c
    if x != 0.0 {
        s += 0.0;
        c += 0.0;
    }
    (s, c)
}

[michaelkirk]

All of these changes sound like reasonable improvements to me. :+1:

[valarauca]

Having done this change locally.

This will break 2 assertions within geodesics.rs #[test] fn test_std_geodesic_geodsolve84

I believe they are

https://github.com/georust/geographiclib-rs/blob/43a1febc531d39f87f1bcbc28579338737b63552/src/geodesic.rs#L2566

https://github.com/georust/geographiclib-rs/blob/43a1febc531d39f87f1bcbc28579338737b63552/src/geodesic.rs#L2570

Instead of the 0.0 they'll return a something like 0.5e-15, which in practice doesn't change/break the results of the inverse/direct equations, but as it is a broken assertion I guess it should be noted.

[michaelkirk]

I believe these are all covered by #61, but let me know if I'm wrong @stonylohr.