sigmaodd.cpp

Go to the documentation of this file.

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

// \cond
#include <cmath>
#include <cstdlib>

#include <algorithm>
#include <iostream>
#include <vector>
// \endcond

#include "../helper/helper.hpp"
#include "pentagonal.hpp"
#include "primes.hpp"

#include "sigmaodd.hpp"


namespace sigmaodd {

  /* ********************
   * Private prototypes *
   **********************/

  void
  exit_if_found_exception_(nat_type start_n, unsigned int length, nat_type n);


  void
  exit_if_impossible_to_check_(nat_type start_n, unsigned int length, nat_type n);


  double
  percent_(nat_type n, nat_type d);


  void
  print_data_(nat_type n,
              nat_type result, unsigned int length,
              bool print_only_longer, bool print_path);


  void
  print_lengths_(nat_type lengths[]);


  void
  print_path_(nat_type n);



  /* *******************
   * Private functions *
   *********************/

  void
  exit_if_found_exception_(nat_type start_n, unsigned int length, nat_type n) {
    if (n == start_n) {
      std::cout << std::endl
                << "! Found not trivial cycle "
                << start_n << '\t' << length << '\t' << n << std::endl;

      exit(EXIT_FAILURE);
    }
  }


  void
  exit_if_impossible_to_check_(nat_type start_n, unsigned int length, nat_type n) {
    if (n > MAX_POSSIBLE_N) {
      std::cout << std::endl
                << "! Impossible to check "
                << start_n << '\t' << length << '\t' << n << std::endl;

      exit(EXIT_FAILURE);
    }
  }


  double
  percent_(nat_type n, nat_type d) {
    assert(d != 0);

    return static_cast<double>(n*100) / static_cast<double>(d);
  }


  void
  print_data_(nat_type n,
              nat_type result, unsigned int length,
              bool print_only_longer, bool print_path) {
    if (!print_only_longer || (length > 1)) {
      std::cout << n << '\t'
                << result << '\t'
                << length;

      if (length > 1) {
        std::cout << ' ' << (length == 2
                             ? '*'
                             : '#');
      }

      if (print_path && (n > 1)) {
        print_path_(n);
      }

      std::cout << std::endl;
    }
  }


  void
  print_lengths_(nat_type lengths[]) {
    nat_type nb = 0;

    for (unsigned int i = 0; i < 4; ++i) {
      nb += lengths[i];
    }

    assert(nb != 0);

    std::cout << std::endl
              << "# length 1:   " << lengths[0]
              << " \t= " << percent_(lengths[0], nb) << '%' << std::endl
              << "# length 2:   " << lengths[1]
              << " \t= " << percent_(lengths[1], nb) << '%' << std::endl
              << "# length 3:   " << lengths[2]
              << " \t= " << percent_(lengths[2], nb) << '%' << std::endl
              << "# length > 3: " << lengths[3]
              << " \t= " << percent_(lengths[3], nb) << '%' << std::endl;
  }


  void
  print_path_(nat_type n) {
    std::cout << '\t' << n;

    nat_type value = n;

    while (value > 1) {
      const nat_type new_value = sum_odd_divisors_divided_until_odd__factorize(value);

      std::cout << ' ' << (value > new_value
                           ? '>'
                           : (value < new_value
                              ? '<'
                              : '=')) << ' ' << new_value;
      value = new_value;
    }
  }



  /* ***********
   * Functions *
   *************/

  void
  check_varsigma_odd(nat_type first_n, nat_type last_n) {
    assert(first_n >= 3);
    assert(is_odd(first_n));

    for (nat_type start_n = first_n; start_n <= last_n; start_n += 2) {
      if (is_little_mersenne_prime_unitary_divide(start_n)) {
        continue;
      }

      nat_type n = start_n;
      unsigned int length = 0;

      do {
        ++length;
        n = varsum_odd(n);
        exit_if_impossible_to_check_(start_n, length, n);
      } while (n > start_n);

      exit_if_found_exception_(start_n, length, n);

      if (length > 1) {
        std::cout << start_n << '\t' << length << '\t' << n << std::endl;
        std::cout.flush();
      }
    }
  }


  std::string
  path_to_string(const std::vector<nat_type> &path) {
    std::string s;

    if (!path.empty()) {
      s += std::to_string(path[0]);
      for (unsigned int i = 1; i < path.size(); ++i) {
        s += (path[i -1] > path[i]
              ? " > "
              : (path[i -1] < path[i]
                 ? " < "
                 : " = "));
        s += std::to_string(path[i]);
      }
    }

    return s;
  }


  void
  print_long(nat_type first, nat_type last, bool print_path) {
    assert(first >= 3);
    assert(is_odd(first));

    for (nat_type start_n = first; start_n <= last; start_n += 2) {
      if (is_little_mersenne_prime_unitary_divide(start_n) || is_square(start_n)) {
        continue;
      }

      const nat_type n = varsum_odd(start_n);

      if (n == start_n) {
        std::cout << "! Found not trivial cycle "
                  << start_n << '\t' << n << std::endl;

        exit(EXIT_FAILURE);
      }
      else if (n > start_n) {
        std::cout << start_n;
        for (unsigned int d = 3; d <= 9; d += 2) {
          if (!is_divide(d, start_n)) {
            std::cout << " !" << d;
          }
        }
        std::cout << '\t' << n << (is_square(n)
                                   ? " #"
                                   : "");
        if (print_path) {
          print_path_(start_n);
        }
        std::cout << std::endl;
        std::cout.flush();
      }
    }
  }


  void
  print_path(const std::vector<nat_type> &path, std::ostream &out) {
    assert(!path.empty());
    assert(path.back() == 1);

    print_path(path.cbegin(), path.cend(), out);
  }


  void
  print_path(const std::vector<nat_type>::const_iterator &path_begin,
             const std::vector<nat_type>::const_iterator &path_end,
             std::ostream &out) {
    if (path_begin != path_end) {
      std::vector<nat_type>::const_iterator it = path_begin;

      out << *it;
      for ( ; it + 1 != path_end; ++it) {
        out << ' ' << (*it > *(it + 1)
                       ? '>'
                       : (*it < *(it + 1)
                          ? '<'
                          : '=')) << ' ' << *(it + 1);
      }
    }
  }


  void
  print_path_infos(const std::vector<nat_type> &path,
                   bool print_bad,
                   bool print_all,
                   bool print_category,
                   bool print_lower,
                   bool print_length,
                   bool print_path,
                   std::ostream &out) {
    assert(!path.empty());

    const nat_type n = path.front();
    const unsigned int length = static_cast<unsigned int>(path.size() - 1);
    const std::vector<nat_type>::const_iterator lower_cit
      = std::find_if(path.cbegin(), path.cend(),
                     [n] (const nat_type &x) { return x < n; });
    const unsigned int partial_length
      = static_cast<unsigned int>(lower_cit - path.cbegin());
    const bool is_square = sigmaodd::is_square(n);
    const bool is_bad = is_odd(n) && !is_square && (partial_length > 1);

    if (print_all
        || (print_bad && is_bad)) {
      assert(lower_cit != path.cend());

      std::cout << n;
      if (is_bad && !is_divide(9, n)) {
        // Exception for specific conjecture about the divisibility of bad numbers
        std::cout << " *";
      }

      if (print_category) {
        std::cout << '\t' << (is_even(n)
                              ? 'E'
                              : (n == path.back()
                                 ? 'F'
                                 : (is_square
                                    ? 'S'
                                    : (is_bad
                                       ? 'B'
                                       : 'G'))));
      }


      // Partial path
      if (print_lower) {
        std::cout << '\t' << *lower_cit;
      }
      if (print_length) {
        std::cout << '\t' << partial_length;
        if (length > n) {  // exception for specific conjecture about the length
          std::cout << " *";
        }
      }
      if (print_path) {
        std::cout << '\t';
        sigmaodd::print_path(path.cbegin(), lower_cit + 1, out);
      }


      // Complete path
      if (partial_length != length) {
        assert(path.back() == 1);

        if (print_length) {
          std::cout << '\t' << length;
        }
        if (print_path) {
          std::cout << '\t';
          sigmaodd::print_path(path, out);
        }
      }

      std::cout << std::endl;
    }
  }


  void
  print_sigmaodd__factorize(nat_type first, nat_type last,
                            bool print_only_longer, bool print_path) {
    assert(first >= 3);
    assert(is_odd(first));

    nat_type lengths[4] = {0, 0, 0, 0};

    for (nat_type n = first; n <= last; n += 2) {
      const std::pair<nat_type, unsigned int> result_length =
        sum_odd_divisors_divided_until_odd_iterate_until_lower__factorize(n);

      assert(result_length.second > 0);

      ++lengths[std::min(4u, result_length.second) - 1];

      print_data_(n,
                  result_length.first, result_length.second,
                  print_only_longer, print_path);
    }

    print_lengths_(lengths);
  }


  void
  print_sigmaodd__factorize_bound(nat_type first, nat_type last,
                                  bool print_only_longer, bool print_path) {
    assert(first >= 3);
    assert(is_odd(first));

    nat_type lengths[4] = {0, 0, 0, 0};

    for (nat_type n = first; n <= last; n += 2) {
      const std::pair<nat_type, unsigned int> result_length =
        sum_odd_divisors_divided_until_odd_iterate_until_lower__factorize_bound(n);

      assert(result_length.second > 0);

      ++lengths[std::min(4u, result_length.second) - 1];

      print_data_(n,
                  result_length.first, result_length.second,
                  print_only_longer, print_path);
    }

    print_lengths_(lengths);
  }


  void
  print_sigmaodd__naive(nat_type first, nat_type last,
                        bool print_only_longer, bool print_path) {
    assert(first >= 3);
    assert(is_odd(first));

    nat_type lengths[4] = {0, 0, 0, 0};

    for (nat_type n = first; n <= last; n += 2) {
      const std::pair<nat_type, unsigned int> result_length =
        sum_odd_divisors_divided_until_odd_iterate_until_lower__naive(n);

      assert(result_length.second > 0);

      ++lengths[std::min(4u, result_length.second) - 1];

      print_data_(n,
                  result_length.first, result_length.second,
                  print_only_longer, print_path);
    }

    print_lengths_(lengths);
  }


  void
  print_square(nat_type first, nat_type last, bool print_path) {
    assert(first >= 3);
    assert(is_odd(first));

    nat_type i_first = floor_square_root(first);

    while (square(i_first) < first) {
      i_first += 2;
    }

    const nat_type i_last = floor_square_root(last);

    for (nat_type i = i_first; i <= i_last; i += 2) {
      const nat_type start_n = square(i);
      const nat_type n = varsum_odd(start_n);

      if (n == start_n) {
        std::cout << "! Found not trivial cycle "
                  << start_n << '\t' << n << std::endl;

        exit(EXIT_FAILURE);
      }

      std::cout << start_n << '\t' << n << (is_square(n)
                                            ? " #"
                                            : "");
      if (print_path) {
        print_path_(start_n);
      }
      std::cout << std::endl;
      std::cout.flush();
    }
  }


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

    return divide_until_odd(sum_divisors__euler(divide_until_odd(n)));
  }


  nat_type
  sum_odd_divisors_divided_until_odd__factorize(nat_type n, bool skip_primes_table) {
    assert(n != 0);

    n = divide_until_odd(n);

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

    if (!skip_primes_table) {
      // 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*divide_until_odd(n + 1);
        }

        unsigned int alpha = 0;

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

        if (alpha > 0) {
          sqrt_n = floor_square_root(n);
          sum *= divide_until_odd(sum_geometric_progression_strict(prime, alpha));
        }
      }
    }

    // Try to find a divisor by Pollard's rho heuristic
    // and use it to try to find coprime factors.
    {
      const nat_type divisor = pollard_rho(n,
                                           static_cast<unsigned int>(std::rand()) % n,
                                           floor_square_root(sqrt_n));

      if (divisor != 0) {
        assert(1 < divisor);
        assert(divisor < n);
        assert(is_divide(divisor, n));

        const std::vector<FactorExp> coprimes = coprime_factor_exps(divisor, n/divisor);

        if (!coprimes.empty()) {
          for (const FactorExp factor_exp : coprimes) {
            sum *= sum_odd_divisors_divided_until_odd__factorize(factor_exp.value(), true);
          }

          return sum;
        }
      }
    }

#ifdef PRIME16
    // Search other prime factors by iteration on potential primes
    nat_type p = odd_primes_table_last();

    assert(p % potential_prime_offsets_table_modulo() == 1);

    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*divide_until_odd(n + 1);
        }

        unsigned int alpha = 0;

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

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

    return sum;
  }


  nat_type
  sum_odd_divisors_divided_until_odd__factorize_bound(nat_type n, nat_type start_n) {
    assert(n != 0);

    n = divide_until_odd(n);

    unsigned int nb = 0;
    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*divide_until_odd(n + 1);
      }

      unsigned int alpha = 0;

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

      if (alpha > 0) {
        sqrt_n = floor_square_root(n);
        sum *= divide_until_odd(sum_geometric_progression_strict(prime, alpha));
        if (prime > sqrt_n) {  // remain n is prime
          return sum*divide_until_odd(n + 1);
        }
        nb = 1000;
      }

      if ((n != start_n) && (++nb > 1000)) {
        assert(n < start_n);

        nb = 0;

        const nat_type bound = sum*sum_odd_divisors_upper_bound__DeTemple(n, prime);

        if (bound < start_n) {
          return bound;  // ???
        }
      }
    }

    // Try to find a divisor by Pollard's rho heuristic
    // and use it to try to find coprime factors.
    {
      const nat_type divisor = pollard_rho(n,
                                           static_cast<unsigned int>(std::rand()) % n,
                                           floor_square_root(sqrt_n));

      if (divisor != 0) {
        assert(1 < divisor);
        assert(divisor < n);
        assert(is_divide(divisor, n));

        const std::vector<FactorExp> coprimes = coprime_factor_exps(divisor, n/divisor);

        if (!coprimes.empty()) {
          for (const FactorExp factor_exp : coprimes) {
            sum *= sum_odd_divisors_divided_until_odd__factorize(factor_exp.value(), true);
          }

          return sum;
        }
      }
    }

#ifdef PRIME16
    // Search other prime factors by iteration on potential primes
    nat_type p = odd_primes_table_last();

    assert(p % potential_prime_offsets_table_modulo() == 1);

    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*divide_until_odd(n + 1);
        }

        unsigned int alpha = 0;

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

        if (alpha > 0) {
          sqrt_n = floor_square_root(n);
          sum *= divide_until_odd(sum_geometric_progression_strict(p, alpha));
          if (p > sqrt_n) {  // remain n is prime ???
            return sum*divide_until_odd(n + 1);
          }
          nb = 1000;
        }

        if ((n != start_n) && (++nb > 1000)) {
          assert(n < start_n);

          nb = 0;

          const nat_type bound = sum*sum_odd_divisors_upper_bound__DeTemple(n, p);

          if (bound < start_n) {
            return bound;  // ???
          }
        }
      }
    }
#endif

    return sum;
  }


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

    return divide_until_odd(sum_odd_divisors__naive(n));
  }


  std::pair<nat_type, unsigned int>
  sum_odd_divisors_divided_until_odd_iterate_until_lower__euler(nat_type n) {
    assert(n > 1);

    const nat_type start_n = n;
    unsigned int length = 0;

    do {
      ++length;
      n = sum_odd_divisors_divided_until_odd__euler(n);
      exit_if_impossible_to_check_(start_n, length, n);
    } while (n > start_n);

    exit_if_found_exception_(start_n, length, n);

    return std::make_pair(n, length);  // ??? swap params?
  }


  std::pair<nat_type, unsigned int>
  sum_odd_divisors_divided_until_odd_iterate_until_lower__factorize(nat_type n) {
    assert(n > 1);

    const nat_type start_n = n;
    unsigned int length = 0;

    do {
      ++length;
      n = sum_odd_divisors_divided_until_odd__factorize(n);
      exit_if_impossible_to_check_(start_n, length, n);
    } while (n > start_n);

    exit_if_found_exception_(start_n, length, n);

    return std::make_pair(n, length);  // ??? swap params?
  }


  std::pair<nat_type, unsigned int>
  sum_odd_divisors_divided_until_odd_iterate_until_lower__factorize_bound(nat_type n) {
    assert(n > 1);

    const nat_type start_n = n;
    unsigned int length = 0;

    do {
      ++length;
      n = sum_odd_divisors_divided_until_odd__factorize_bound(n, start_n);
      exit_if_impossible_to_check_(start_n, length, n);
    } while (n > start_n);

    exit_if_found_exception_(start_n, length, n);

    return std::make_pair(n, length);  // ??? swap params?
  }


  std::pair<nat_type, unsigned int>
  sum_odd_divisors_divided_until_odd_iterate_until_lower__naive(nat_type n) {
    assert(n > 1);

    const nat_type start_n = n;
    unsigned int length = 0;

    do {
      ++length;
      n = sum_odd_divisors_divided_until_odd__naive(n);
      exit_if_impossible_to_check_(start_n, length, n);
    } while (n > start_n);

    exit_if_found_exception_(start_n, length, n);

    return std::make_pair(n, length);  // ??? swap params?
  }


  nat_type
  varsum_odd(nat_type n) {
    assert(is_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; 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*divide_until_odd(n + 1);
      }

      unsigned int alpha = 0;

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

      if (alpha > 0) {
        sqrt_n = floor_square_root(n);
        sum *= divide_until_odd(sum_geometric_progression_strict(prime, alpha));
      }
    }

    return sum*varsum_odd_big(n, sqrt_n);
  }


  nat_type
  varsum_odd_big(nat_type n, nat_type sqrt_n) {
    assert(is_odd(n));
#ifdef PRIME16
    assert(n >= potential_prime_offsets_table_modulo());
#endif

    nat_type sum = 1;

    // Try to find a divisor by Pollard's rho heuristic
    // and use it to try to find coprime factors.
    {
      const nat_type divisor = pollard_rho(n,
                                           static_cast<unsigned int>(std::rand()) % n,
                                           floor_square_root(sqrt_n));

      if (divisor != 0) {
        assert(1 < divisor);
        assert(divisor < n);
        assert(is_divide(divisor, n));

        const std::vector<FactorExp> coprimes = coprime_factor_exps(divisor, n/divisor);

        if (!coprimes.empty()) {
          for (const FactorExp factor_exp : coprimes) {
            sum *= varsum_odd_big(factor_exp.value(), floor_square_root(factor_exp.value()));
          }

          return sum;
        }
      }
    }

#ifdef PRIME16
    // Search other prime factors by iteration on potential primes
    nat_type p = odd_primes_table_last();

    assert(p % potential_prime_offsets_table_modulo() == 1);

    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*divide_until_odd(n + 1);
        }

        unsigned int alpha = 0;

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

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

    return sum;
  }

}  // namespace sigmaodd