divisors.cpp

Go to the documentation of this file.

/* -*- coding: latin-1 -*- */
/** \file common/sigmaodd/divisors.cpp (December 20, 2017)
 *
 * GPLv3 --- Copyright (C) 2017 Olivier Pirson
 * http://www.opimedia.be/
 */

// \cond
#include <cassert>
#include <cfenv>
#include <cstdlib>

#include <utility>
// \endcond

#include "../helper/helper.hpp"
#include "divisors.hpp"
#include "harmonic.hpp"
#include "primes.hpp"


namespace sigmaodd {

  /* *******************
   * Private prototype *
   *********************/

  std::vector<nat_type>
  coprime_factors_(nat_type a, nat_type b);



  /* ******************
   * Private function *
   ********************/

  /** \brief
   * Return a list of factors such that all factors are coprimes
   * and divide a*b,
   * and these factors are (possibly) sufficient for coprime_factor_exps().
   *
   * @param a > 1
   * @param b > 1
   */
  std::vector<nat_type>
  coprime_factors_(nat_type a, nat_type b) {
    assert(a > 1);
    assert(b > 1);

    // Preliminary special cases
    if (a == b) {
      return std::vector<nat_type>{a};
    }

    if (a > b) {
      std::swap(a, b);
    }


    // Preliminary simplifications
    do {
      while (is_divide(a, b)) {
        b /= a;
      }

      if (b == 1) {
        assert(a != 1);

        return std::vector<nat_type>{a};
      }

      if (a > b) {
        std::swap(a, b);
      }
      else {
        break;
      }
    } while (is_divide(a, b));

    assert(a > 1);
    assert(b > 1);
    assert(a < b);


    // Identify {a, b} coprime or divide a and b by their GCD
    const nat_type d = gcd(a, b);

    if (d == 1) {  // {a, b} are coprimes
      return std::vector<nat_type>{a, b};
    }
    else {         // {a, b} are NOT coprimes
      do {
        a /= d;
      } while (is_divide(d, a));

      do {
        b /= d;
      } while (is_divide(d, b));

      if (a > b) {
        std::swap(a, b);
      }

      if (a == 1) {
        return std::vector<nat_type>{d, b};
      }

      assert(a > 1);
      assert(b > 1);
      assert(a < b);

      // These {a, b} are coprimes

      bool is_coprime_d_a = is_coprime(d, a);
      bool is_coprime_d_b = is_coprime(d, b);

      if (is_coprime_d_a && is_coprime_d_b) {  // {d, a, b} are coprime
        return std::vector<nat_type>{d, a, b};
      }
      else {    // {d, a} are NOT coprime and/or {d, b} are NOT coprime
        if (is_coprime_d_b) {
          std::swap(is_coprime_d_a, is_coprime_d_b);
          std::swap(a, b);
        }

        if (is_coprime_d_a) {  // {d, a} are coprime and {d, b} are NOT coprime
          std::vector<nat_type> factors = (d == b
                                           ? std::vector<nat_type>{b}
                                           : coprime_factors_(d, b));  // d*b < initial product

          factors.push_back(a);

          return factors;
        }
        else {                 // {d, a} and {d, b} are NOT coprime
          return coprime_factors_(d, a*b);  // d*a*b < initial product
        }
      }
    }
  }



  /* ***********
   * Functions *
   *************/
  std::ostream&
  operator<<(std::ostream &out, const FactorExp &factor_exp) {
    assert(factor_exp.exp != 0);

    out << factor_exp.factor;
    if (factor_exp.exp > 1) {
      out << '^' << factor_exp.exp;
    }

    return out;
  }



  std::vector<FactorExp>
  coprime_factor_exps(nat_type a, nat_type b) {
    assert(a > 1);
    assert(b > 1);

    nat_type n = a*b;

    std::vector<FactorExp> coprimes;

    for (nat_type factor : coprime_factors_(a, b)) {
      // Use collected coprime factors to rebuild product
      assert(factor != 1);
      assert(factor < a*b);

      unsigned int alpha = 0;

      while (is_divide(factor, n)) {
        n /= factor;
        ++alpha;
      }

      // factor^alpha is a (possibly) correct coprime factor
      coprimes.push_back({factor, alpha});
    }

    if (n != 1) {  // the collected factors don't complete divide a*b in coprime factors
      for (auto it = coprimes.begin(); it < coprimes.end(); ++it) {
        if (!is_coprime(it->factor, n)) {
          // the remain factor is NOT coprime with this collected coprime factor, heuristic failed
          return std::vector<FactorExp>();
        }
      }

      // the remain factor is coprime with collected coprime factors
      coprimes.push_back({n, 1});
    }

    assert(!coprimes.empty());

    return (coprimes.size() > 1
            ? coprimes  // success
            : std::vector<FactorExp>());  // heuristic failed to find at least two coprime factors
  }


  std::pair<nat_type, unsigned int>
  divide_until_nb(nat_type n, nat_type d) {
    assert(n != 0);
    assert(d > 1);

    unsigned int alpha = 0;

    while (is_divide(d, n)) {
      n /= d;
      ++alpha;
    }

    return std::make_pair(n, alpha);
  }


  nat_type
  first_divisor(nat_type n) {
    if (is_even(n)) {
      return (n != 0
              ? 2
              : 1);
    }
    else {
      const nat_type sqrt_n = floor_square_root(n);

      // Search prime divisor present in the precalculated table
      for (unsigned int i = 0; i < odd_primes_table_nb(); ++i) {
        const nat_type prime = odd_primes_table_by_index(i);

        if (is_divide(prime, n)) {
          return prime;
        }
        else if (prime > sqrt_n) {  // n is prime
          return n;
        }
      }

#ifdef PRIME16
      // Search other prime divisor by iteration on potential primes
      assert(odd_primes_table_last() % potential_prime_offsets_table_modulo() == 1);

      nat_type p = odd_primes_table_last();

      while (p < sqrt_n) {
        for (unsigned int i = 0; i < potential_prime_offsets_table_nb(); ++i) {
          p += potential_prime_offsets_table_by_index(i);

          if (is_divide(p, n)) {
            return p;
          }
        }
      }
#endif

      return n;  // n is prime
    }
  }


  bool
  is_square(nat_type n) {
    if (n >= 18446744065119617025ul) {  // >= (2^(64/2) - 1)^2 = 18446744065119617025
      return (n == 18446744065119617025ul);
    }
    else {
      nat_type sqrt_n = static_cast<nat_type>(floor(sqrt(static_cast<double>(n))));

      /* Correct possible rounding error due to floating point computation */
      if (square(sqrt_n) < n) {
        do {
          ++sqrt_n;
        } while (square(sqrt_n) < n);
      }
      else {
        while (square(sqrt_n) > n) {
          --sqrt_n;
        }
      }

      return (square(sqrt_n) == n);
    }
  }


  nat_type
  pollard_rho(nat_type n) {
    assert(n != 0);

    return pollard_rho(n, static_cast<unsigned int>(std::rand()) % n, floor_square_root(n));
  }


  nat_type
  pollard_rho(nat_type n, nat_type random, nat_type max_iteration) {
    assert(n != 0);
    assert(random < n);

    /*
      Adapted from:
      Thomas H. Cormen, Charles E. Leiserson, Ron Rivest, and Clifford Stein.
      Introduction to Algorithms. 3rd 2009
    */
    nat_type x = random;
    nat_type y = x;
    nat_type k = 2;

    for (nat_type i = 1; i < max_iteration; ++i) {
      x = (x != 0
           ? static_cast<nat_type>((static_cast<tmp_uint128_type>(x)*x - 1) % n)
           : n - 1);

      if (is_even(i)) {
        const nat_type d = gcd((y >= x
                                ? y - x
                                : x - y), n);

        if ((d != 1) && (d != n)) {
          assert(is_divide(d, n));

          return d;
        }
      }

      if (i == k) {
        y = x;
        k <<= 1;
      }
    }

    return 0;
  }


  nat_type
  pollard_rho_repeat(nat_type n, nat_type nb_tries) {
    assert(n != 0);

    while (nb_tries-- > 0) {
      const nat_type divisor = pollard_rho(n);

      if (divisor != 0) {
        return divisor;
      }
    }

    return 0;
  }


  nat_type
  sum_divisors__factorize(nat_type n) {
    assert(n != 0);

    const std::pair<nat_type, unsigned int> result_alpha = divide_until_odd_nb(n);

    n = result_alpha.first;

    nat_type sum = (static_cast<nat_type>(1) << (result_alpha.second + 1)) - 1;
    nat_type sqrt_n = floor_square_root(n);

    // Search prime factors present in the precalculated table
    for (unsigned int i = 0; (n > 1) && (i < odd_primes_table_nb()); ++i) {
      const nat_type prime = odd_primes_table_by_index(i);

      if (prime > sqrt_n) {  // remain n is prime
        return sum*(n + 1);
      }

      unsigned int alpha = 0;

      while (is_divide(prime, n)) {
        n /= prime;
        sqrt_n = floor_square_root(n);
        ++alpha;
      }

      if (alpha > 0) {
        sum *= sum_geometric_progression_strict(prime, alpha);
      }
    }

#ifdef PRIME16
    // Search other prime factors by iteration on potential primes
    assert(odd_primes_table_last() % potential_prime_offsets_table_modulo() == 1);

    nat_type p = odd_primes_table_last();

    while (n > 1) {
      for (unsigned int i = 0; i < potential_prime_offsets_table_nb(); ++i) {
        p += potential_prime_offsets_table_by_index(i);

        if (p > sqrt_n) {  // remain n is prime
          return sum*(n + 1);
        }

        unsigned int alpha = 0;

        while (is_divide(p, n)) {
          n /= p;
          sqrt_n = floor_square_root(n);
          ++alpha;
        }

        if (alpha > 0) {
          sum *= sum_geometric_progression_strict(p, alpha);
        }
      }
    }
#endif

    return sum;
  }


  nat_type
  sum_divisors__naive(nat_type n) {
    assert(n != 0);

    const nat_type sqrt_n = floor_square_root(n);

    nat_type sum = 0;

    for (nat_type i = 1; i <= sqrt_n; ++i) {
      if (is_divide(i, n)) {
        sum += i + n/i;
      }
    }

    if (square(sqrt_n) == n) {
      sum -= sqrt_n;
    }

    return sum;
  }


  nat_type
  sum_odd_divisors__factorize(nat_type n) {
    assert(n != 0);

    n = divide_until_odd(n);

    nat_type sum = 1;
    nat_type sqrt_n = floor_square_root(n);

    // Search prime factors present in the precalculated table
    for (unsigned int i = 0; (n > 1) && (i < odd_primes_table_nb()); ++i) {
      const nat_type prime = odd_primes_table_by_index(i);

      if (prime > sqrt_n) {  // remain n is prime
        return sum*(n + 1);
      }

      unsigned int alpha = 0;

      while (is_divide(prime, n)) {
        n /= prime;
        sqrt_n = floor_square_root(n);
        ++alpha;
      }

      if (alpha > 0) {
        sum *= sum_geometric_progression_strict(prime, alpha);
      }
    }

#ifdef PRIME16
    // Search other prime factors by iteration on potential primes
    assert(odd_primes_table_last() % potential_prime_offsets_table_modulo() == 1);

    nat_type p = odd_primes_table_last();

    while (n > 1) {
      for (unsigned int i = 0; i < potential_prime_offsets_table_nb(); ++i) {
        p += potential_prime_offsets_table_by_index(i);

        if (p > sqrt_n) {  // remain n is prime
          return sum*(n + 1);
        }

        unsigned int alpha = 0;

        while (is_divide(p, n)) {
          n /= p;
          sqrt_n = floor_square_root(n);
          ++alpha;
        }

        if (alpha > 0) {
          sum *= sum_geometric_progression_strict(p, alpha);
        }
      }
    }
#endif

    return sum;
  }


  nat_type
  sum_odd_divisors__naive(nat_type n) {
    assert(n != 0);

    n = divide_until_odd(n);

    const nat_type sqrt_n = floor_square_root(n);

    nat_type sum = 0;

    for (nat_type i = 1; i <= sqrt_n; i += 2) {
      if (is_divide(i, n)) {
        sum += i + n/i;
      }
    }

    if (square(sqrt_n) == n) {
      sum -= sqrt_n;
    }

    return sum;
  }


  nat_type
  sum_odd_divisors_upper_bound(nat_type n) {
    assert(n != 0);

    n = divide_until_odd(n);

    nat_type bound = sum_odd_divisors_upper_bound__DeTemple(n);

    if (n <= 10) {
      return bound;
    }

    const nat_type sqrt_n = floor_square_root(n);
    double diff = 0;

    const int round = fegetround();

    fesetround(FE_DOWNWARD);

    unsigned int i = 0;

    if (true) {
      // Remove until first prime divisor
      for ( ; i < 10; ++i) {
        const nat_type prime = odd_primes_table_by_index(i);

        if (prime > sqrt_n) {
          fesetround(round);

          return n + 1;
        }

        if (is_divide(prime, n)) {
          const nat_type k = prime - 2;

          // ???
          // bound = (k + 1)*k/2 - 1 + n - harmonic(k);
          for (nat_type d = 3; d <= k; d += 2) {
            bound -= d;
            diff += static_cast<double>(n)/static_cast<double>(d);
          }

          break;
        }
      }
    }

    if (true) {
      // Remove if not divide
      for ( ; i < 10; ++i) {
        const nat_type prime = odd_primes_table_by_index(i);

        if (prime > sqrt_n) {
          break;
        }

        if (!is_divide(prime, n)) {
          bound -= prime;
          diff += static_cast<double>(n)/static_cast<double>(prime);
        }
      }
    }

    fesetround(round);

    assert(diff >= 0);

    return bound - static_cast<nat_type>(std::floor(diff));
  }


  nat_type
  sum_odd_divisors_upper_bound__DeTemple(nat_type n) {
    assert(n != 0);

    // ????????????????? n = divide_until_odd(n);  // ???? virer ???

    if (n > 1) {
      const nat_type sqrt_n = floor_square_root(n);

      // floor(sqrt(n - 1))
      const nat_type sqrt_n1 = (square(sqrt_n) == n
                                ? sqrt_n - 1
                                : sqrt_n);

      // ((S + 1)/2)^2 + n (H_S1 - 1/2 H_{S1/2})   with S = floor(sqrt(n)) and S1 = floor(sqrt(n - 1))
      const int round = fegetround();

      fesetround(FE_UPWARD);

      const nat_type harmonic_part
        = static_cast<nat_type>(std::floor((diff_half_harmonic_upper_bound(sqrt_n1))
                                           * static_cast<double>(n)));

      fesetround(round);

      return square((sqrt_n + 1)/2) + harmonic_part;
    }
    else {
      return 1;
    }
  }


  nat_type
  sum_odd_divisors_upper_bound__DeTemple_ceil(nat_type n) {
    assert(n != 0);
    assert(is_odd(n));

    if (n > 1) {
      const nat_type sqrt_n = floor_square_root(n);

      // floor(sqrt(n - 1))
      const nat_type sqrt_n1 = (square(sqrt_n) == n
                                ? sqrt_n - 1
                                : sqrt_n);

      // ((S + 1)/2)^2 + n (H_S1 - 1/2 H_{S1/2})   with S = floor(sqrt(n)) and S1 = floor(sqrt(n - 1))
      const int round = fegetround();

      fesetround(FE_UPWARD);

      const nat_type harmonic_part
        = static_cast<nat_type>(std::ceil(diff_half_harmonic_upper_bound(sqrt_n1))) * n;

      fesetround(round);

      return square((sqrt_n + 1)/2) + harmonic_part;
    }
    else {
      return 1;
    }
  }


  nat_type
  sum_odd_divisors_upper_bound__DeTemple(nat_type n, nat_type k) {
    assert(n != 0);
    assert(is_odd(n));

    assert(is_odd(k));
    assert(k <= floor_square_root(n));

#ifndef NDEBUG
    for (nat_type i = 3; i <= k; i += 2) {
      assert(!is_divide(i, n));
    }
#endif

    nat_type bound = sum_odd_divisors_upper_bound__DeTemple(n);
    double diff = diff_half_harmonic_upper_bound(k);  // ???

    bound += n + 1;
    bound -= square(k + 1)/4;
    bound -= static_cast<nat_type>((std::floor(static_cast<double>(n)*diff)));

    return bound;
  }

}  // namespace sigmaodd