/* * java-math-library is a Java library focused on number theory, but not necessarily limited to it. It is based on the PSIQS 4.0 factoring project. * Copyright (C) 2018 Tilman Neumann (www.tilman-neumann.de) * * This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. * * You should have received a copy of the GNU General Public License along with this program; * if not, see <http://www.gnu.org/licenses/>. */ package de.tilman_neumann.jml; import static de.tilman_neumann.jml.base.BigIntConstants.*; import java.math.BigInteger; import java.util.ArrayList; import java.util.HashSet; import java.util.List; import java.util.Map; import java.util.TreeMap; import java.util.TreeSet; import org.apache.log4j.Logger; import de.tilman_neumann.util.ConfigUtil; /** * Test Collatz or 3n+1 problem. * @author Tilman Neumann */ public class CollatzSequenceTest { private static final Logger LOG = Logger.getLogger(CollatzSequenceTest.class); private static final int N_COUNT = 30000000; private static final int START_PROGRESSION = 2; // 2 is best for dropping sequences; 2, 3, 9 are good for repeat sequences private static final int MAX_PROGRESSION = 1<<20; /** * Run the 3n+1 sequence until we find some n<nStart. */ private static void test3nPlus1DroppingSequence() { TreeSet<Integer> lengths = new TreeSet<>(); int maxLength = 0; ArrayList<Integer> recordLengths = new ArrayList<Integer>(); ArrayList<BigInteger> recordLengthsN = new ArrayList<BigInteger>(); // convention: length(1)=0 recordLengths.add(0); recordLengthsN.add(I_1); // analyze the length of all dropping sequences, arranged by length TreeMap<Integer, ArrayList<BigInteger>> length2NStartList = new TreeMap<Integer, ArrayList<BigInteger>>(); // Define arithmetic progressions of dropping sequences: progr(mul, add) = {mul*n+add, n=0,1,2,...} // Then we analyze the maximal length of sequences of progression progr(mul, add). // Many such progressions have a small maximal length. // Progressions that do not have a small maximal length are called "unbounded". // Finally we analyze the counts of unbounded progressions for sets {progr(mul, add), add=0..(mul-1)} // The smallest numbers of unbounded proressions ar achieved for mul being powers of 2. TreeMap<Integer, TreeMap<Integer, Integer>> progressionToMaxLength = new TreeMap<Integer, TreeMap<Integer, Integer>>(); HashSet<BigInteger> resolved = new HashSet<>(); resolved.add(I_1); for (int i=2; i<N_COUNT; i++) { ArrayList<BigInteger> sequence = new ArrayList<>(); BigInteger nStart = BigInteger.valueOf(i); BigInteger n = nStart; sequence.add(n); while (true) { if ((n.intValue()&1)==1) { n = n.multiply(I_3).add(I_1); } else { n = n.shiftRight(1); } sequence.add(n); if (resolved.contains(n)) { resolved.add(nStart); // sequence length is the number of operations, not the number of elements in the list int length = sequence.size() - 1; //LOG.info("len=" + length + ": " + sequence); lengths.add(length); if (length > maxLength) { maxLength = length; LOG.info("record length=" + length + " at n=" + nStart); // sequence of records = A217934 = 0, 1, 6, 11, 96, 132, 171, 220, 267, 269, 282, 287, 298, 365, 401, 468, 476, 486, 502, 613, 644, 649, 706, 729, 892, 897, 988, 1122, 1161, 1177, 1187, 1445, 1471, 1575, 1614, 1639 recordLengths.add(length); recordLengthsN.add(nStart); } addToLength2NMap(length2NStartList, length, nStart); addToProgressionToMaxLengthMap(progressionToMaxLength, nStart, length); break; } } } LOG.info("lengths = " + lengths); // A122437 // data below comes from N_COUNT=100M (16GB RAM hardly enough) LOG.info("record lengths = " + recordLengths); // A217934 = (0,) 1, 6, 11, 96, 132, 171, 220, 267, 269, 282, 287, 298, 365, 401, 468, 476, 486, 502, 613 LOG.info("N(record lengths)= " + recordLengthsN); // A060412 = (1,) 2, 3, 7, 27, 703, 10087, 35655, 270271, 362343, 381727, 626331, 1027431, 1126015, 8088063, 13421671, 20638335, 26716671, 56924955, 63728127 analyzeNStartSequences(length2NStartList); analyzeProgressions(progressionToMaxLength); } /** * Run the 2n+1 sequence until we find some n that was already contained in the sequence of some smaller n. */ @SuppressWarnings("unused") private static void test3nPlus1RepeatSequence() { TreeSet<Integer> lengths = new TreeSet<>(); int maxLength = 0; ArrayList<Integer> recordLengths = new ArrayList<Integer>(); ArrayList<BigInteger> recordLengthsN = new ArrayList<BigInteger>(); // convention: length(1)=0 recordLengths.add(0); recordLengthsN.add(I_1); // analyze the length of all repeat sequences, arranged by length TreeMap<Integer, ArrayList<BigInteger>> length2NStartList = new TreeMap<Integer, ArrayList<BigInteger>>(); // Define arithmetic progressions of "repeat" sequences: progr(mul, add) = {mul*n+add, n=0,1,2,...} // Then we analyze the maximal length of sequences of progression progr(mul, add). // Many such progressions have a small maximal length. // Progressions that do not have a small maximal length are called "unbounded". // Finally we analyze the counts of unbounded progressions for sets {progr(mul, add), add=0..(mul-1)} // The smallest numbers of unbounded progressions are achieved for mul of the form 2^r*3^s, s small. TreeMap<Integer, TreeMap<Integer, Integer>> progressionToMaxLength = new TreeMap<Integer, TreeMap<Integer, Integer>>(); HashSet<BigInteger> resolved = new HashSet<>(); resolved.add(I_1); for (int i=2; i<N_COUNT; i++) { ArrayList<BigInteger> sequence = new ArrayList<>(); BigInteger nStart = BigInteger.valueOf(i); BigInteger n = nStart; sequence.add(n); resolved.add(n); while (true) { if ((n.intValue()&1)==1) { n = n.multiply(I_3).add(I_1); } else { n = n.shiftRight(1); } sequence.add(n); if (resolved.contains(n)) { // sequence length is the number of operations, not the number of elements in the list int length = sequence.size() - 1; //LOG.info("len=" + length + ": " + sequence); lengths.add(length); if (length > maxLength) { maxLength = length; LOG.info("record length=" + length + " at n=" + nStart); recordLengths.add(length); recordLengthsN.add(nStart); } addToLength2NMap(length2NStartList, length, nStart); addToProgressionToMaxLengthMap(progressionToMaxLength, nStart, length); break; } resolved.add(n); } } LOG.info("lengths = " + lengths); // all naturals except for 2, 4 // the following data stems from N_COUNT=30M, which required about 12GB of RAM LOG.info("record lengths = " + recordLengths); // not on OEIS: 0, 1, 6, 10, 95, 132, 219, 262, 269, 271, 297, 305, 343, 357, 400, 468, 485 LOG.info("N(record lengths)= " + recordLengthsN); // not on OEIS: 1, 2, 3, 7, 27, 703, 35655, 270271, 362343, 401151, 1027431, 1327743, 1394431, 6206655, 8088063, 13421671, 26716671 analyzeNStartSequences(length2NStartList); analyzeProgressions(progressionToMaxLength); } private static void addToLength2NMap(TreeMap<Integer, ArrayList<BigInteger>> length2NList, int length, BigInteger nStart) { ArrayList<BigInteger> nList = length2NList.get(length); if (nList==null) nList = new ArrayList<BigInteger>(); nList.add(nStart); length2NList.put(length, nList); } private static void analyzeNStartSequences(TreeMap<Integer, ArrayList<BigInteger>> length2NList) { for (int length : length2NList.keySet()) { ArrayList<BigInteger> nList = length2NList.get(length); //LOG.debug("length=" + length + ": n0List = " + nList); int period = findPeriod(nList); if (period>0) { List<BigInteger> nPeriodList = nList.subList(0, period); BigInteger periodDiff = nList.get(period).subtract(nList.get(0)); LOG.info("length=" + length + ": #n0=" + nList.size() + ", period=" + period + ", periodDiff=" + periodDiff + ", nPeriodList = " + nPeriodList); } else { LOG.info("length=" + length + ": #n0=" + nList.size() + ", period not identified"); } } } private static int findPeriod(ArrayList<BigInteger> nList) { int nCount = nList.size(); ArrayList<BigInteger> diffs = new ArrayList<BigInteger>(); for (int i=1; i<nCount; i++) { diffs.add(nList.get(i).subtract(nList.get(i-1))); } //LOG.debug("diffs = " + diffs); int diffCount = diffs.size(); // should be nCount-1 for (int period=1; period<diffCount; period++) { int i; for (i=period; i<diffCount; i++) { if (!diffs.get(i).equals(diffs.get(i-period))) { // not the correct period break; } } if (i==diffCount && diffCount/period >= 2) { // all differences confirmed -> we found the correct period! return period; } } // nothing found return 0; } private static void addToProgressionToMaxLengthMap(TreeMap<Integer, TreeMap<Integer, Integer>> progressionToMaxLength, BigInteger nStart, int length) { for (int mul=START_PROGRESSION; mul<=MAX_PROGRESSION; mul<<=1) { int add = nStart.mod(BigInteger.valueOf(mul)).intValue(); TreeMap<Integer, Integer> add2MaxLength = progressionToMaxLength.get(mul); if (add2MaxLength==null) add2MaxLength = new TreeMap<Integer, Integer>(); Integer maxLength = add2MaxLength.get(add); if (maxLength==null || length > maxLength) { add2MaxLength.put(add, length); progressionToMaxLength.put(mul, add2MaxLength); } } } private static void analyzeProgressions(TreeMap<Integer, TreeMap<Integer, Integer>> progressionToMaxLength) { ArrayList<Integer> unboundedProgressionCounts = new ArrayList<Integer>(); for (Integer mul : progressionToMaxLength.keySet()) { // Find maximal bounded maxLength for arithmetic progressions with given multiplier TreeSet<Integer> maxLengths = new TreeSet<Integer>(); for (Map.Entry<Integer, Integer> add2MaxLength : progressionToMaxLength.get(mul).entrySet()) { maxLengths.add(add2MaxLength.getValue()); } int maxBoundedLength = 2; for (int maxLength : maxLengths) { // bottom-up if (maxLength - maxBoundedLength <= 3) { maxBoundedLength = maxLength; } } // count unbounded arithmetic progressions int unboundedProgressionsCount = 0; for (Map.Entry<Integer, Integer> add2MaxLength : progressionToMaxLength.get(mul).entrySet()) { int progrMaxLength = add2MaxLength.getValue(); //LOG.info("Progression " + mul + "*n + " + add2MaxLength.getKey() + ": maxLength = " + progrMaxLength); if (progrMaxLength > maxBoundedLength) { unboundedProgressionsCount++; } } LOG.info("mul=" + mul + ": maxBoundedLength = " + maxBoundedLength + ", unboundedProgressionsCount = " + unboundedProgressionsCount); unboundedProgressionCounts.add(unboundedProgressionsCount); } LOG.info("sequence of unbounded progression counts: " + unboundedProgressionCounts); // A076227 = Number of surviving Collatz residues mod 2^n = 1, 1, 1, 2, 3, 4, 8, 13, 19, 38, 64, 128, 226, 367, 734, 1295, 2114, 4228, 7495, 14990, 27328, 46611, 93222, 168807, 286581, 573162, 1037374, 1762293, 3524586, 6385637, 12771274, 23642078, 41347483, 82694966, 151917636, 263841377, 527682754, 967378591, 1934757182, 3611535862 } public static void main(String[] args) { ConfigUtil.initProject(); test3nPlus1DroppingSequence(); // test3nPlus1RepeatSequence(); } }