Cite where BN_div actually comes from

After spending a while trying to divine where all the bounds came from,
and coming up with some of the messy proofs for why it works, I found
this exact algorithm in Knuth, Volume 2, with... different messy proofs.
Sadly, this algorithm seems to just be messy. Cite it as reference
rather than trying to repeat it in code.

As part of this, update the discussion on branches. That was added in
https://boringssl-review.googlesource.com/c/boringssl/+/9105, back when
BN_div was used on secret inputs. It no longer is and, back then, the
function still wasn't constant-time anyway.

We could, in principle, restore the special cases now. But this would be
more complicated and diverge from Knuth's formulation, so let's just
keep it simple. (Although it might actually be a hair faster. We care
about this function to compute R^2 mod n, and the special case would
save an extra iteration through the loop. Though I think that
optimization could actually be restored with much, much less code than
OpenSSL originally did it. Probably not worth the fuss.)

Subsequent CLs will clean this code up in reference to Knuth's
formulation.

Bug: 358687140
Change-Id: I56da99c560b845f1736ab86edc79b8e711890fe3
Reviewed-on: https://boringssl-review.googlesource.com/c/boringssl/+/70170
Reviewed-by: Bob Beck <bbe@google.com>
Commit-Queue: David Benjamin <davidben@google.com>

crypto/fipsmodule/bn/div.c

diff --git a/crypto/fipsmodule/bn/div.c b/crypto/fipsmodule/bn/div.c
index 23e1b63..0410fb7 100644
--- a/crypto/fipsmodule/bn/div.c
+++ b/crypto/fipsmodule/bn/div.c
@@ -175,23 +175,29 @@
 #endif
 }
 
-// This was specifically designed to contain fewer branches that may leak
-// sensitive information; see "New Branch Prediction Vulnerabilities in OpenSSL
-// and Necessary Software Countermeasures" by Onur Acıçmez, Shay Gueron, and
-// Jean-Pierre Seifert.
 int BN_div(BIGNUM *quotient, BIGNUM *rem, const BIGNUM *numerator,
            const BIGNUM *divisor, BN_CTX *ctx) {
+  // This function implements long division, per Knuth, The Art of Computer
+  // Programming, Volume 2, Chapter 4.3.1, Algorithm D. This algorithm only
+  // divides non-negative integers, but we round towards zero, so we divide
+  // absolute values and adjust the signs separately.
+  //
+  // Inputs to this function are assumed public and may be leaked by timing and
+  // cache side channels. Division with secret inputs should use other
+  // implementation strategies such as Montgomery reduction.
+  //
+  // Historically, this function diverged from Knuth's algorithm with some
+  // shortcuts in some cases. Those have been removed per "New Branch Prediction
+  // Vulnerabilities in OpenSSL and Necessary Software Countermeasures" by Onur
+  // Acıçmez, Shay Gueron, and Jean-Pierre Seifert. We continue to omit them for
+  // simplicity, but this function is no longer used with secret inputs. (We
+  // implement a variation on "Smooth CRT-RSA" as described in the paper.)
   int norm_shift, loop;
   BIGNUM wnum;
   BN_ULONG *resp, *wnump;
   BN_ULONG d0, d1;
   int num_n, div_n;
 
-  // This function relies on the historical minimal-width |BIGNUM| invariant.
-  // It is already not constant-time (constant-time reductions should use
-  // Montgomery logic), so we shrink all inputs and intermediate values to
-  // retain the previous behavior.
-
   // Invalid zero-padding would have particularly bad consequences.
   int numerator_width = bn_minimal_width(numerator);
   int divisor_width = bn_minimal_width(divisor);