div_chunk(): Уменьшаем возможную ошибку до 1

big_int.cpp

@@ -223,17 +223,21 @@ static vector<Digit> div_by_base(const vector<Digit>& magnitude)
 #endif
 
 /// Делит кусок числителя на знаменатель.
-/// Числитель и знаменатель таковы, что в результате всегда одна цифра.
+/// Кусок и знаменатель таковы, что в результате всегда одна цифра.
 /// К тому же аргументы должны быть нормализованы (старшая цифра знаменателя >= base / 2)
 static Digit div_chunk(const vector<Digit>& chunk, const vector<Digit>& denominator)
 {
     if (first_is_less(chunk, denominator))
         return 0;
 
-    assert(chunk.size() == denominator.size() // Длина числителя равна длине знаменателя
-           || (chunk.size() == denominator.size() + 1 // Или числитель длиннее знаменателя на 1 цифру
+    assert(chunk.size() == denominator.size() // Длина куска равна длине знаменателя
+           || (chunk.size() == denominator.size() + 1 // Или кусок длиннее знаменателя на 1 цифру
                && first_is_less(div_by_base(chunk), denominator))); // И эта цифра не лишняя
 
+    // Когда в знаменателе одна цифра, используется div_by_digit()
+    assert(denominator.size() >= 2);
+
+    // Должно быть нормализовано
     assert(denominator.back() >= base / 2);
 
     // Две старшие цифры куска.
@@ -242,21 +246,42 @@ static Digit div_chunk(const vector<Digit>& chunk, const vector<Digit>& denomina
     if (chunk.size() != denominator.size())
         chunk_2_digits = chunk_2_digits * base + chunk[chunk.size() - 2];
 
+    // Третья цифра куска
+    DDigit chunk_digit_3 = (chunk.size() == denominator.size())
+        ? chunk[chunk.size() - 2] // Старшая цифра куска - 0, но его нет в массиве
+        : chunk[chunk.size() - 3];
+
+    // Старшая цифра знаменателя
+    DDigit denominator_digit_1 = denominator.back();
+
+    // Вторая цифра знаменателя
+    DDigit denominator_digit_2 = denominator[denominator.size() - 2];
+
     // Делим две старшие цифры куска на старшую цифру знаменателя
-    Digit digit = (Digit)min(chunk_2_digits / denominator.back(), DDigit(base - 1));
+    DDigit digit = chunk_2_digits / denominator_digit_1;
 
-#ifndef NDEBUG
-    size_t misses_count = 0; // Число промахов
-#endif
+    // Остаток может быть больше base
+    DDigit remainder = chunk_2_digits - digit * denominator_digit_1;
 
-    // while (chunk < digit * denominator) // Пока цифра слишком большая
-    while (first_is_less(chunk, mul_magnitudes(vector<Digit>{digit}, denominator)))
+    // Уменьшаем цифру, если она точно слишком большая
+    if (digit == base || digit * denominator_digit_2 > remainder * base + chunk_digit_3)
     {
         --digit;
-        assert(++misses_count <= 2);
+        remainder += denominator_digit_1;
+
+        if (remainder < base && (digit == base || digit * denominator_digit_2 > remainder * base + chunk_digit_3))
+            --digit;
     }
 
-    return digit;
+    assert(digit < base);
+
+    // Теперь проверяем, нет ли промаха на 1
+    if (first_is_less(chunk, mul_magnitudes(vector<Digit>{(Digit)digit}, denominator)))
+        --digit;
+
+    assert(!first_is_less(chunk, mul_magnitudes(vector<Digit>{(Digit)digit}, denominator)));
+
+    return (Digit)digit;
 }
 
 /// Возвращает неполное частное и остаток (делит столбиком).

docs/libs.md

@@ -35,3 +35,4 @@
 * http://hvks.com/Numerical/arbitrary_precision.html
 * https://mattmccutchen.net/bigint/ (public domain)
 * https://github.com/weidai11/cryptopp
+* https://github.com/llvm/llvm-project/blob/900be9013fdc3bab9fce906f8a71e59ecd8873b4/llvm/lib/Support/APInt.cpp#L1260