Repository: scipy/scipy
Pull request: BUG: special: Fix formula for division of dual numbers in xsf
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #21665 

#### What does this implement/fix?
<!--Please explain your changes.-->

This PR fixes the definition of division of dual numbers in xsf. There were some typos in this definition in #21483 that weren't caught because division of dual numbers is not used in that PR.

@izaid took the definition from the 2023 master's thesis of Songchen Tan, [Higher-Order Automatic Differentiation and its Applications](https://dspace.mit.edu/bitstream/handle/1721.1/151501/tan-songchen-csesm-ccse-2023-thesis%20Songchen%20Tan.pdf?sequence=1&isAllowed=y). The relevant formula is in item (2.9) in section 2.3.2 and is pasted below:

![image](https://github.com/user-attachments/assets/4a9ee721-4470-407b-89d2-b001a12e4b7c)

In #21483, this implemented as

```C++
dual &operator/=(const dual &other) {
        for (size_t i = 0; i <= Order; ++i) {
            T coef = 1; // binomial coefficient
            for (size_t j = 1; j <= i; ++j) {
                data[i] -= coef * other.data[j] * data[i - j];
                coef *= T(i - j) / T(j + 1);
            }

            data[i] /= other.data[0];
        }
```

One can see that: 

* In iteration $(i, j)$:`coef` contains the value ${i - 1 \choose j}$
* The outer loop is over $1 \leq i \leq \mathrm{Order}$
* For each `i` the inner loop is over $1 \leq j \leq i$.
*  The case `i = 0` is skipped because `j` starts at 1, so the loop condition `j <= i` never holds. 
* `other.data[j]` contains the value of $g^{(j)}$ for each $j$
* `data[i]` initially contains the value of $f^{(i)}$ for each `i`,and the intention is that `data[i]` is updated to contain the value of $\left(\frac{f}{g}\right)^{(i)}$ as we step through the loop.

Putting this all together, we see that the formula which has actually been implemented is

$$\left(\frac{f}{g}\right)^{(n)} = \frac{1}{g}\left[f^{(n)} - \sum_{k=1}^{n}{n - 1 \choose k} g^{(k)}\left(\frac{f}{g}\right)^{\left(n - k\right)}\right]$$

Comparing to the formula above, we see that `other.data` and `data` have been switched, and there's an off by one in the indices. 

In #21568, I replaced the `coef` which is computed on the fly directly with `binom(i, j)`, which itself was incorrect, due to not fixing the off by one error in the indices. I also didn't notice that `data` and `other.data` were swapped.

I've updated things here so that the formula matches the reference.


#### Additional information
<!--Any additional information you think is important.-->
I haven't added a test case because it would require creating a new ufunc just for testing purposes. I think it's fine to test this when we get around to adding derivatives for a special function which uses division internally.

cc @mdhaber


Problem statement:
BUG: special: autodifferentiation of division operator
### Describe your issue.

In gh-21568, we found that the implementation of `dual &operator/=(const dual &other)` may not be quite right. I thought we should open this so it's not forgotten. 
@izaid can you fix this, and @steppi can review?
In any case, we should add a test and a reference so it doesn't come up again.

Hints from issue discussion:

Commit message: Fix dual division formula

Patch to review:
diff --git a/scipy/special/xsf/dual.h b/scipy/special/xsf/dual.h
index 18aa01cd7a87..c2a60ad73c01 100644
--- a/scipy/special/xsf/dual.h
+++ b/scipy/special/xsf/dual.h
@@ -120,9 +120,8 @@ class dual<T, Order> {
 
     dual &operator/=(const dual &other) {
         for (size_t i = 0; i <= Order; ++i) {
-            // General Leibniz Rule
-            for (size_t j = 1; j <= i; ++j) {
-                data[i] -= detail::fast_binom<T>(i, j) * other.data[j] * data[i - j];
+            for (size_t j = 0; j < i; ++j) {
+                data[i] -= detail::fast_binom<T>(i, j) * data[j] * other.data[i - j];
             }
 
             data[i] /= other.data[0];
@@ -250,9 +249,8 @@ class dual<T, Order0, Order1, Orders...> {
 
     dual &operator/=(const dual &other) {
         for (size_t i = 0; i <= Order0; ++i) {
-            // General Leibniz Rule
-            for (size_t j = 1; j <= i; ++j) {
-                data[i] -= detail::fast_binom<T>(i, j) * other.data[j] * data[i - j];
+            for (size_t j = 0; j < i; ++j) {
+                data[i] -= detail::fast_binom<T>(i, j) * data[j] * other.data[i - j];
             }
 
             data[i] /= other.data[0];


Code review comments:
diff_hunk: @@ -120,9 +120,8 @@ class dual<T, Order> {
 
     dual &operator/=(const dual &other) {
         for (size_t i = 0; i <= Order; ++i) {
-            // General Leibniz Rule
-            for (size_t j = 1; j <= i; ++j) {
-                data[i] -= detail::fast_binom<T>(i, j) * other.data[j] * data[i - j];
+            for (size_t j = 0; j < i; ++j) {
+                data[i] -= detail::fast_binom<T>(i, j) * data[j] * other.data[i - j];
path: scipy/special/xsf/dual.h
text: @user1:
```suggestion
            for (size_t j = 1; j <= i; ++j) {
                data[i] -= detail::fast_binom<T>(i - 1, j) * other.data[j] * data[i - j];
```

diff_hunk: @@ -250,9 +249,8 @@ class dual<T, Order0, Order1, Orders...> {
 
     dual &operator/=(const dual &other) {
         for (size_t i = 0; i <= Order0; ++i) {
-            // General Leibniz Rule
-            for (size_t j = 1; j <= i; ++j) {
-                data[i] -= detail::fast_binom<T>(i, j) * other.data[j] * data[i - j];
+            for (size_t j = 0; j < i; ++j) {
+                data[i] -= detail::fast_binom<T>(i, j) * data[j] * other.data[i - j];
path: scipy/special/xsf/dual.h
text: @user1:
```suggestion
            for (size_t j = 1; j <= i; ++j) {
                data[i] -= detail::fast_binom<T>(i - 1, j) * other.data[j] * data[i - j];
```