Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

consider simplifying polyval #6

Closed
stonylohr opened this issue Jul 16, 2020 · 3 comments
Closed

consider simplifying polyval #6

stonylohr opened this issue Jul 16, 2020 · 3 comments

Comments

@stonylohr
Copy link
Contributor

stonylohr commented Jul 16, 2020

A few changes to consider for geomath::polyval.

Consider removing the "s" argument, instead have the calling function slice the "p" argument any time "s" might be non-zero. I think this is more idiomatic rust, and allows simplifying polyval's loop. That in turn allows for dropping the shadowing of the "n" argument (see #3) without needing to even declare "n" as a mut argument. polyval's definition then becomes a more streamlined...

pub fn polyval(n: i64, p: &[f64], x: f64) -> f64 {
    let mut y = if n < 0 { 0.0 } else { p[0] };
    for i in 1..=n {
        y = y * x + p[i as usize];
    }
    y
}

Calls to polyval would need to be modified to pass a slice starting at index "s", rather than starting at 0 and passing "s".

Example 1 (called with fixed "s" of 0):
Before: polyval(m, &COEFF, 0, sq(eps))
After: polyval(m, &COEFF, sq(eps))

Example 2 (called with variable "s" value):
Before: geomath::polyval(m, &COEFF_A3, o as usize, _n)
After: geomath::polyval(m, &COEFF_A3[o as usize..], _n)
(Note that calling contexts like this would be further improved by declaring "o" as usize in the first place, to avoid the need to cast. This change should be straightforward in most cases.)

One note of caution on this change is that C++ also supports similar array offset mechanics, yet Charles Karney chose not to use them in the original C++. I suspect this is mostly a matter of style, but it could hint that there's a slight performance benefit to passing the "s" index, in which case that difference is likely (but not certain) to carry over into rust.

Depending on your change tolerance, I would further consider changing polyval's "n" argument from i64 to isize. The difference is academic for 64 bit platforms, but for values that are primarily used as an index into arrays, it's generally preferable to stick with the architecture's native size. This would involve a change at each calling context. While the changes should be fairly simple, the benefit is also modest.

@stonylohr
Copy link
Contributor Author

I take back my cautionary note about the C++ version passing the index. Apparently my memory failed me on that point. It appears to use offset mechanics similar to those I proposed.

I'm also reminded by looking back at the C++ code that it might be worth adding a comment to consider using fma in polyval at some point in the future, as hardware support for it becomes common enough (apparently using it without hardware support results in a performance hit).

@michaelkirk
Copy link
Member

Interesting, sounds good to me!

The original transcription was based on the python implementation -so I assume that explains the divergence from the C++ implementation.

Note: I edited your post to add some "```" in order to format the code blocks more readably.

This was referenced Jan 23, 2021
stonylohr added a commit to stonylohr/geographiclib-rs that referenced this issue Jan 30, 2021
Uses a slightly different implementation than my original proposal.
stonylohr added a commit to stonylohr/geographiclib-rs that referenced this issue Jan 30, 2021
Slightly different from my original proposal.
bors bot added a commit that referenced this issue Feb 1, 2021
39: Address #6: polyval r=michaelkirk a=stonylohr

Slightly different from my original proposal. There are a few pieces worth highlighting.

First, I think the most idiomatic Rust approach would actually be to drop the `n` argument as well as `s`, and let the size of the `p` argument imply the order of the polynomial. The reason I didn't take this further step is that keeping `n` requires fewer changes, and leaves the code a closer parallel to its C++ and Java counterparts. I figure we can always choose to reconsider later, perhaps if and when we focus on performance. The performance implications are likely minor, and I'm not sure which approach they favor.

Second, while keeping `n`, it's somewhat tempting to change it to `usize`, since its support for the negative special case doesn't appear to be used in any of Karney's code that I've examined or logged. It's certainly not used in any of our current geodesic functionality. It may be used in geoid, magnetic model, or projection special cases that aren't hit by any of the scenarios I've run in my instrumented C++. My current guess is that it's not actually used and was only added for completeness. Anyway, I left it signed for now, but changed it from i64 to isize since we ultimately need to cast it to usize for use as an array index, and any value too large for i32 would strongly hint that we're already off the rails. The largest `n` I've seen in instrumented scenarios is 29.

Third, I went with a slightly different polyval implementation than my original proposal. I felt that my new version is slightly more clear that it never tries to cast a negative `n` to usize (though both should actually be safe), and that it might give the optimizer more of a hint that we don't care about the `i` value except for accessing an element of `p` (though I suspect it would recognize that anyway).

Fourth, I tried to keep changes to polyval's callers minimal in most cases, but there were two I couldn't resist: Both `_C3f` and `_C4f` were using floating point values to represent array indices. I switched these to usize, also tweaking the types of some variables they interact with.

Finally, these changes should be behavior-neutral. I think they'll help me with arcdirect troubleshooting because polyval is called extensively in some sections of code used to calculate Geodesic member variables that end up with suspect values, and these changes will make it easier for me to compare the C++ and Rust code for these sections. However, they don't do anything to resolve the arcdirect problems themselves.

I'd be happy to revise my approach if you'd prefer something different for any of these various items.

- [ ] I agree to follow the project's [code of conduct](https://github.com/georust/geo/blob/master/CODE_OF_CONDUCT.md).
- [ ] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---



Co-authored-by: Stony <stony.lohr@gmail.com>
@michaelkirk
Copy link
Member

Fixed by #39

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants