<commit_before>// @(#)root/mathcore:$Id$
// Author: <NAME> 11/2020

/*************************************************************************
 * Copyright (C) 1995-2020, <NAME> and <NAME>.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include <cstdint>

/// Multiply two 576 bit numbers, stored as 9 numbers of 64 bits each
///
/// \param[in] in1 first factor as 9 numbers of 64 bits each
/// \param[in] in2 second factor as 9 numbers of 64 bits each
/// \param[out] out result with 18 numbers of 64 bits each
void multiply9x9(const uint64_t *in1, const uint64_t *in2, uint64_t *out)
{
   uint64_t next = 0;
   unsigned nextCarry = 0;

#if defined(__clang__) || defined(__INTEL_COMPILER) || defined(__CUDACC__)
#pragma unroll
#elif defined(__GNUC__) && __GNUC__ >= 8
// This pragma was introduced in GCC version 8.
#pragma GCC unroll 18
#endif
   for (int i = 0; i < 18; i++) {
      uint64_t current = next;
      unsigned carry = nextCarry;

      next = 0;
      nextCarry = 0;

#if defined(__clang__) || defined(__INTEL_COMPILER) || defined(__CUDACC__)
#pragma unroll
#elif defined(__GNUC__) && __GNUC__ >= 8
// This pragma was introduced in GCC version 8.
#pragma GCC unroll 9
#endif
      for (int j = 0; j < 9; j++) {
         int k = i - j;
         if (k < 0 || k >= 9)
            continue;

         uint64_t upper1 = in1[j] >> 32;
         uint64_t lower1 = static_cast<uint32_t>(in1[j]);

         uint64_t upper2 = in2[k] >> 32;
         uint64_t lower2 = static_cast<uint32_t>(in2[k]);

         // Multiply 32-bit parts.
         uint64_t upper = upper1 * upper2;
         uint64_t middle1 = upper1 * lower2;
         uint64_t middle2 = lower1 * upper2;
         uint64_t middle = middle1 + middle2;
         if (middle < middle1) {
            // This can never overflow because the maximum value of upper is
            // (2 ** 32 - 1) ** 2 = 2 ** 64 - 2 * 2 ** 32 + 1. When now adding
            // another 2 ** 32, the result 2 ** 64 - 2 ** 32 + 1 is still smaller
            // than the maximum 2 ** 64 - 1 that can be stored in a uint64_t.
            upper += uint64_t(1) << 32;
         }
         uint64_t lower = lower1 * lower2;

         uint64_t middle_upper = middle >> 32;
         uint64_t middle_lower = middle << 32;

         // Add to current, remember carry.
         current += lower;
         if (current < lower)
            carry++;
         current += middle_lower;
         if (current < middle_lower)
            carry++;

         // Add to next, remember nextCarry.
         next += middle_upper;
         if (next < middle_upper)
            nextCarry++;
         next += upper;
         if (next < upper)
            nextCarry++;
      }

      next += carry;
      if (next < carry)
         nextCarry++;

      out[i] = current;
   }
}

/// Compute a value congruent to mul modulo m less than 2 ** 576
///
/// \param[in] mul product from multiply9x9 with 18 numbers of 64 bits each
/// \param[out] out result with 9 numbers of 64 bits each
///
/// \f$ m = 2^{576} - 2^{240} + 1 \f$
///
/// Note that this function does *not* return the smallest value congruent to<commit_msg>Optimize overflow handling in multiply9x9<commit_after>// @(#)root/mathcore:$Id$
// Author: <NAME> 11/2020

/*************************************************************************
 * Copyright (C) 1995-2020, <NAME> and <NAME>.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include <cstdint>

/// Multiply two 576 bit numbers, stored as 9 numbers of 64 bits each
///
/// \param[in] in1 first factor as 9 numbers of 64 bits each
/// \param[in] in2 second factor as 9 numbers of 64 bits each
/// \param[out] out result with 18 numbers of 64 bits each
void multiply9x9(const uint64_t *in1, const uint64_t *in2, uint64_t *out)
{
   uint64_t next = 0;
   unsigned nextCarry = 0;

#if defined(__clang__) || defined(__INTEL_COMPILER) || defined(__CUDACC__)
#pragma unroll
#elif defined(__GNUC__) && __GNUC__ >= 8
// This pragma was introduced in GCC version 8.
#pragma GCC unroll 18
#endif
   for (int i = 0; i < 18; i++) {
      uint64_t current = next;
      unsigned carry = nextCarry;

      next = 0;
      nextCarry = 0;

#if defined(__clang__) || defined(__INTEL_COMPILER) || defined(__CUDACC__)
#pragma unroll
#elif defined(__GNUC__) && __GNUC__ >= 8
// This pragma was introduced in GCC version 8.
#pragma GCC unroll 9
#endif
      for (int j = 0; j < 9; j++) {
         int k = i - j;
         if (k < 0 || k >= 9)
            continue;

         uint64_t upper1 = in1[j] >> 32;
         uint64_t lower1 = static_cast<uint32_t>(in1[j]);

         uint64_t upper2 = in2[k] >> 32;
         uint64_t lower2 = static_cast<uint32_t>(in2[k]);

         // Multiply 32-bit parts, each product has a maximum value of
         // (2 ** 32 - 1) ** 2 = 2 ** 64 - 2 * 2 ** 32 + 1.
         uint64_t upper = upper1 * upper2;
         uint64_t middle1 = upper1 * lower2;
         uint64_t middle2 = lower1 * upper2;
         uint64_t lower = lower1 * lower2;

         // When adding the two products, the maximum value for middle is
         // 2 * 2 ** 64 - 4 * 2 ** 32 + 2, which exceeds a uint64_t.
         uint64_t middle = middle1 + middle2;
         // Handling the overflow by a multiplication with 0 or 1 is cheaper
         // than branching with an if statement, which the compiler does not
         // optimize to this equivalent code. Note that we could do entirely
         // without this overflow handling when summing up the intermediate
         // products differently as described in the following SO answer:
         //    https://stackoverflow.com/a/51587262
         // However, this approach takes at least the same amount of thinking
         // why a) the code gives the same results without b) overflowing due
         // to the mixture of 32 bit arithmetic. Moreover, my tests show that
         // the scheme implemented here is actually slightly more performant.
         int overflow = (middle < middle1);
         // This addition can never overflow because the maximum value of upper
         // is 2 ** 64 - 2 * 2 ** 32 + 1 (see above). When now adding another
         // 2 ** 32, the result is 2 ** 64 - 2 ** 32 + 1 and still smaller than
         // the maximum 2 ** 64 - 1 that can be stored in a uint64_t.
         upper += overflow * (uint64_t(1) << 32);

         uint64_t middle_upper = middle >> 32;
         uint64_t middle_lower = middle << 32;

         lower += middle_lower;
         if (lower < middle_lower)
            upper++;

         // This still can't overflow since the maximum of middle_upper is
         //  - 2 ** 32 - 4 if there was an overflow for middle above, bringing
         //    the maximum value of upper to 2 ** 64 - 2.
         //  - otherwise upper still has the initial maximum value given above
         //    and the addition of a value smaller than 2 ** 32 brings it to
         //    a maximum value of 2 ** 64 - 2 ** 32 + 2.
         // (Both cases include the increment to handle the overflow in lower.)
         //
         // All the reasoning makes perfect sense given that the product of two
         // 64 bit numbers is smaller than or equal to
         //     (2 ** 64 - 1) ** 2 = 2 ** 128 - 2 * 2 ** 64 + 1
         // with the upper bits matching the 2 ** 64 - 2 of the first case.
         upper += middle_upper;

         // Add to current, remember carry.
         current += lower;
         if (current < lower)
            carry++;

         // Add to next, remember nextCarry.
         next += upper;
         if (next < upper)
            nextCarry++;
      }

      next += carry;
      if (next < carry)
         nextCarry++;

      out[i] = current;
   }
}

/// Compute a value congruent to mul modulo m less than 2 ** 576
///
/// \param[in] mul product from multiply9x9 with 18 numbers of 64 bits each
/// \param[out] out result with 9 numbers of 64 bits each
///
/// \f$ m = 2^{576} - 2^{240} + 1 \f$
///
/// Note that this function does *not* return the smallest value congruent to
{'fraction_non_alphanumeric': 0.13619124238593266, 'fraction_numerical': 0.048844960349385125, 'mean_word_length': 1.9201342281879195, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 33, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2070274', 'n_tokens_mistral': 2815, 'n_tokens_neox': 2421, 'n_words': 983}