16
\$\begingroup\$

Description

Let a permutation of the integers {1, 2, ..., n} be called minimally interpolable if no set of k+2 points (together with their indices) fall on a polynomial of degree k. That is,

  1. No two points fall on a horizontal line (0-degree polynomial)
  2. No three points fall on a line (1-degree polynomial)
  3. No four points fall on a parabola (2-degree polynomial)
  4. Et cetera.

Challenge

Write a program that computes OEIS sequence A301802(n), the number of minimally interpolable permutations of {1, 2, ..., n} for n as a large as possible.


Scoring

I will time your code on my computer (2.3 GHz Intel Core i5, 8 GB RAM) with increasing inputs. Your score will be the greatest input that takes less than 1 minute to output the correct value.


Example

For example, the permutation [1, 2, 4, 3] is minimally interpolable because

the terms together with their indices 
[(1, 1), (2, 2), (3, 4), (4, 3)] 
have the property that
  (0) No two points have the same y-value.
  (1) No three points lie on a line.
  (2) No four points lie on a parabola.

Example illustrating that [1,2,4,3] is minimally interpolable. In the illustration, you can see that the horizontal lines (red) have at most one point on them, the lines (blue) have at most two points on them, and the parabolas (green) have three points on them.


Data

Here are the minimally interpolable permutations for n=3, n=4, and n=5:

n = 3: [1,3,2],[2,1,3],[2,3,1],[3,1,2]
n = 4: [1,2,4,3],[1,3,2,4],[1,3,4,2],[1,4,2,3],[2,1,3,4],[2,1,4,3],[2,3,1,4],[2,4,1,3],[2,4,3,1],[3,1,2,4],[3,1,4,2],[3,2,4,1],[3,4,1,2],[3,4,2,1],[4,1,3,2],[4,2,1,3],[4,2,3,1],[4,3,1,2]
n = 5: [1,2,5,3,4],[1,3,2,5,4],[1,3,4,2,5],[1,4,2,3,5],[1,4,3,5,2],[1,4,5,2,3],[1,4,5,3,2],[1,5,3,2,4],[2,1,4,3,5],[2,3,1,4,5],[2,3,5,1,4],[2,3,5,4,1],[2,4,1,5,3],[2,4,3,1,5],[2,4,5,1,3],[2,5,1,3,4],[2,5,1,4,3],[2,5,3,4,1],[2,5,4,1,3],[3,1,4,5,2],[3,1,5,2,4],[3,1,5,4,2],[3,2,5,1,4],[3,2,5,4,1],[3,4,1,2,5],[3,4,1,5,2],[3,5,1,2,4],[3,5,1,4,2],[3,5,2,1,4],[4,1,2,5,3],[4,1,3,2,5],[4,1,5,2,3],[4,1,5,3,2],[4,2,1,5,3],[4,2,3,5,1],[4,2,5,1,3],[4,3,1,2,5],[4,3,1,5,2],[4,3,5,2,1],[4,5,2,3,1],[5,1,3,4,2],[5,2,1,3,4],[5,2,1,4,3],[5,2,3,1,4],[5,2,4,3,1],[5,3,2,4,1],[5,3,4,1,2],[5,4,1,3,2]

If my program is correct, the first few values of a(n), the number of minimally interpolable permutations of {1, 2, ..., n}:

a(1) = 1
a(2) = 2
a(3) = 4
a(4) = 18
a(5) = 48
a(6) = 216
a(7) = 584
a(8) = 2870
\$\endgroup\$
7
  • \$\begingroup\$ Nice sequence number! | Although you specified fastest-code, you didn't specify which machine is it fastest on. What is exactly the winning criteria? \$\endgroup\$
    – DELETE_ME
    Mar 27, 2018 at 0:36
  • 3
    \$\begingroup\$ To add to user202729's comment, I suggest some tags that you can use to determine the winning criteria: fastest-code requires that submissions be tested on the same machine to compare the runtime (usually the OP of the challenge does this). fastest-algorithm would ask answerers to come up with code with the lowest time complexity as possible. code-golf would ask users to come up with code with the shortest source code (or equivalent) as possible. Other than that, this is indeed a nice challenge. \$\endgroup\$ Mar 27, 2018 at 0:55
  • \$\begingroup\$ Your example text uses zero-indexing though the image uses one-indexing. \$\endgroup\$ Mar 27, 2018 at 20:54
  • \$\begingroup\$ Since all points are defined by permuations of the first natural numbers, is it not impossible for any two points to occupy the same height? \$\endgroup\$ Mar 27, 2018 at 20:57
  • \$\begingroup\$ @JonathanFrech, indeed, it should be 1-indexed because these are permutations. And you're correct! Because we're dealing with permutations, the 0-degree polynomial condition comes for free. \$\endgroup\$ Mar 27, 2018 at 22:04

2 Answers 2

5
\$\begingroup\$

C#

using System;
using System.Diagnostics;
using BigInteger = System.Int32;

namespace Sandbox
{
    class PPCG160382
    {
        public static void Main(params string[] args)
        {
            if (args.Length != 0)
            {
                foreach (var arg in args) Console.WriteLine(CountValidPerms(int.Parse(arg)));
            }
            else
            {
                int[] smallValues = new int[] { 1, 1, 2, 4, 18, 48 };
                for (int n = 0; n < smallValues.Length; n++)
                {
                    var observed = CountValidPerms(n);
                    var expected = smallValues[n];
                    Console.WriteLine(observed == expected ? $"{n}: Ok" : $"{n}: expected {expected}, observed {observed}, error {observed - expected}");
                }
                for (int n = smallValues.Length; n < 13; n++)
                {
                    Stopwatch sw = new Stopwatch();
                    sw.Start();
                    Console.WriteLine($"{n}: {CountValidPerms(n)} in {sw.ElapsedMilliseconds}ms");
                }
            }
        }

        private static long CountValidPerms(int n)
        {
            // We work on the basis of exclusion by extrapolation.
            var unused = (1 << n) - 1;
            var excluded = new int[n];
            int[] perm = new int[n];

            // Symmetry exclusion: perm[0] < (n+1) / 2
            if (n > 1) excluded[0] = (1 << n) - (1 << ((n + 1) / 2));

            long count = 0;
            CountValidPerms(ref count, perm, 0, unused, excluded);
            return count;
        }

        private static void CountValidPerms(ref long count, int[] perm, int off, int unused, int[] excluded)
        {
            int n = perm.Length;
            if (off == n)
            {
                count += CountSymmetries(perm);
                return;
            }

            // Quick-aborts
            var completelyExcluded = excluded[off];
            for (int i = off + 1; i < n; i++)
            {
                if ((unused & ~excluded[i]) == 0) return;
                completelyExcluded &= excluded[i];
            }
            if ((unused & completelyExcluded) != 0) return;

            // Consider each unused non-excluded value as a candidate for perm[off]
            var candidates = unused & ~excluded[off];
            for (int val = 0; candidates > 0; val++, candidates >>= 1)
            {
                if ((candidates & 1) == 0) continue;

                perm[off] = val;

                var nextUnused = unused & ~(1 << val);

                var nextExcluded = (int[])excluded.Clone();
                // For each (non-trivial) subset of smaller indices, combine with off and extrapolate to off+1 ... excluded.Length-1
                if (off < n - 1 && off > 0)
                {
                    var points = new Point[off + 1];
                    var denoms = new BigInteger[off + 1];
                    points[0] = new Point { X = off, Y = perm[off] };
                    denoms[0] = 1;
                    ExtendExclusions(perm, off, 0, points, 1, denoms, nextExcluded);
                }

                // Symmetry exclusion: perm[0] < perm[-1] < n - 1 - perm[0]
                if (off == 0 && n > 1)
                {
                    nextExcluded[n - 1] |= (1 << n) - (2 << (n - 1 - val));
                    nextExcluded[n - 1] |= (2 << val) - 1;
                }

                CountValidPerms(ref count, perm, off + 1, nextUnused, nextExcluded);
            }
        }

        private static void ExtendExclusions(int[] perm, int off, int idx, Point[] points, int numPoints, BigInteger[] denoms, int[] excluded)
        {
            if (idx == off) return;

            // Subsets without
            ExtendExclusions(perm, off, idx + 1, points, numPoints, denoms, excluded);

            // Just add this to the subset
            points[numPoints] = new Point { X = idx, Y = perm[idx] };
            denoms = (BigInteger[])denoms.Clone();
            // Update invariant: denoms[s] = prod_{t != s} points[s].X - points[t].X
            denoms[numPoints] = 1;
            for (int s = 0; s < numPoints; s++)
            {
                denoms[s] *= points[s].X - points[numPoints].X;
                denoms[numPoints] *= points[numPoints].X - points[s].X;
            }
            numPoints++;

            for (int target = off + 1; target < excluded.Length; target++)
            {
                BigInteger prod = 1;
                for (int t = 0; t < numPoints; t++) prod *= target - points[t].X;

                Rational sum = new Rational(0, 1);
                for (int s = 0; s < numPoints; s++) sum += new Rational(prod / (target - points[s].X) * points[s].Y, denoms[s]);

                if (sum.Denom == 1 && sum.Num >= 0 && sum.Num < excluded.Length) excluded[target] |= 1 << (int)sum.Num;
            }

            // Subsets with
            ExtendExclusions(perm, off, idx + 1, points, numPoints, denoms, excluded);
        }

        private static int CountSymmetries(int[] perm)
        {
            if (perm.Length < 2) return 1;

            int cmp = 0;
            for (int i = 0, j = perm.Length - 1; i <= j; i++, j--)
            {
                cmp = perm.Length - 1 - perm[i] - perm[j];
                if (cmp != 0) break;
            }

            return cmp > 0 ? 4 : cmp == 0 ? 2 : 0;
        }

        public struct Point
        {
            public int X;
            public int Y;
        }

        public struct Rational
        {
            public Rational(BigInteger num, BigInteger denom)
            {
                if (denom == 0) throw new ArgumentOutOfRangeException(nameof(denom));

                if (denom < 0) { num = -num; denom = -denom; }

                var g = _Gcd(num, denom);
                Num = num / g;
                Denom = denom / g;
            }

            private static BigInteger _Gcd(BigInteger a, BigInteger b)
            {
                if (a < 0) a = -a;
                if (b < 0) b = -b;
                while (a != 0)
                {
                    var tmp = b % a;
                    b = a;
                    a = tmp;
                }
                return b;
            }

            public BigInteger Num;
            public BigInteger Denom;

            public static Rational operator +(Rational a, Rational b) => new Rational(a.Num * b.Denom + a.Denom * b.Num, a.Denom * b.Denom);
        }
    }
}

Takes values of n as command-line arguments, or if run without arguments times itself up to n=10. Compiling as "Release" in VS 2017 and running on an Intel Core i7-6700 I calculate n=9 in 1.2 seconds, and n=10 in 13.6 seconds. n=11 is just over 2 minutes.

FWIW:

n    a(n)
9    10408
10   45244
11   160248
12   762554
\$\endgroup\$
0
0
\$\begingroup\$

Rust

\$\color{red}{\text{I was trying to rewrite @Peter Taylor's C# answer in Rust. But the Rust code here has some mistakes.}}\$

Any help would be appreciated.

src/main.rs

use std::time::Instant;
use num_bigint::BigInt;
use num_traits::Zero;
use num_traits::One;
use num_traits::Signed;
use num_traits::ToPrimitive;


#[derive(Clone, Copy)]
struct Point {
    x: i32,
    y: i32,
}

#[derive(Clone)]
struct Rational {
    num: BigInt,
    denom: BigInt,
}


impl Rational {
    fn new(num: BigInt, denom: BigInt) -> Self {
        let g = gcd(num.clone(), denom.clone());
        Rational { num: num / g.clone(), denom: denom / g }
    }
}

impl std::ops::Add for Rational {
    type Output = Rational;

    fn add(self, other: Rational) -> Rational {
        let self_denom_clone = self.denom.clone();
        Rational::new(self.num.clone() * other.denom.clone() + self_denom_clone * other.num.clone(),
            self.denom * other.denom.clone())
    }
    
}

fn main() {
    let args: Vec<String> = std::env::args().collect();
    if args.len() > 1 {
        for arg in &args[1..] {
            let n = arg.parse::<i32>().unwrap();
            println!("{}", count_valid_perms(n));
        }
    } else {
        let small_values = vec![1, 1, 2, 4, 18, 48];
        for n in 0..small_values.len() {
            let observed = count_valid_perms(n as i32);
            let expected = small_values[n];
            if observed == expected {
                println!("{}: Ok", n);
            } else {
                println!("{}: expected {}, observed {}, error {}", n, expected, observed, observed - expected);
            }
        }
        for n in small_values.len()..13 {
            let start = Instant::now();
            println!("{}: {} in {}ms", n, count_valid_perms(n as i32), start.elapsed().as_millis());
        }
    }
}



fn gcd(a: BigInt, b: BigInt) -> BigInt {
    let mut a = a;
    let mut b = b;
    while !a.is_zero() {
        let tmp = b.clone() % a.clone();
        b = a;
        a = tmp;
    }
    b.abs()
}




fn count_valid_perms(n: i32) -> i64 {
    let unused = (1 << n) - 1;
    let mut excluded = vec![0; n as usize];
    let mut perm = vec![0; n as usize];
    if n > 1 { excluded[0] = (1 << n) - (1 << ((n + 1) / 2)); }
    let mut count: i64 = 0;
    count_valid_perms_recursive(&mut count, &mut perm, 0, unused, &mut excluded);
    count
}

// fn count_valid_perms_recursive(count: &mut i64, perm: &mut Vec<i32>, off: i32, mut unused: i32, excluded: &mut Vec<i32>) {
fn count_valid_perms_recursive(count: &mut i64, perm: &mut Vec<i32>, off: i32, unused: i32, excluded: &mut Vec<i32>) {

    let n = perm.len();
    if off == n as i32 {
        *count += count_symmetries(perm);
        return;
    }
    let mut completely_excluded = excluded[off as usize];
    for i in (off + 1)..n as i32 {
        if (unused & !excluded[i as usize]) == 0 { return; }
        completely_excluded &= excluded[i as usize];
    }
    if (unused & completely_excluded) != 0 { return; }
    let mut candidates = unused & !excluded[off as usize];
    for val in 0.. {
        if candidates > 0 { 
            candidates >>= 1;
            if (candidates & 1) == 0 { continue; }
            perm[off as usize] = val;
            let next_unused = unused & !(1 << val);
            let mut next_excluded = excluded.clone();
            if off < n as i32 - 1 && off > 0 {
                let mut points = vec![Point { x: 0, y: 0 }; (off + 1) as usize];
                let mut denoms = vec![BigInt::one(); (off + 1) as usize];
                points[0] = Point { x: off, y: perm[off as usize] };
                extend_exclusions(perm, off, 0, &mut points, 1, &mut denoms, &mut next_excluded);
            }
            if off == 0 && n > 1 {
                next_excluded[n - 1] |= (1 << n) - (2 << (n - 1 - (val as usize)));
                next_excluded[n - 1] |= (2 << val) - 1;
            }
            count_valid_perms_recursive(count, perm, off + 1, next_unused, &mut next_excluded);
        } else {
            break;
        }
    }
}

fn extend_exclusions(perm: &mut Vec<i32>, off: i32, idx: i32, points: &mut Vec<Point>, num_points: usize, denoms: &mut Vec<BigInt>, excluded: &mut Vec<i32>) {
    if idx == off { return; }
    extend_exclusions(perm, off, idx + 1, points, num_points, denoms, excluded);
    points[num_points] = Point { x: idx, y: perm[idx as usize] };
    let denoms_clone = denoms.clone();
    denoms[num_points] = BigInt::one();
    for s in 0..num_points {
        denoms[s] *= BigInt::from(points[s].x - points[num_points].x);
        denoms[num_points] *= BigInt::from(points[num_points].x - points[s].x);
    }
    let num_points = num_points + 1;
    for target in ((off + 1) as usize)..excluded.len() {
        let mut prod = BigInt::one();
        for t in 0..num_points {
            prod *= BigInt::from(target as i32 - points[t].x);
        }
        let mut sum = Rational { num: BigInt::zero(), denom: BigInt::one() };
        for s in 0..num_points {
            let prod_s = prod.clone() / BigInt::from(target as i32 - points[s].x) * BigInt::from(points[s].y);
            sum = sum + Rational::new(prod_s, denoms_clone[s].clone());
        }
        if sum.denom == BigInt::one() && sum.num >= BigInt::zero() && sum.num.to_usize().unwrap() < excluded.len() {
            excluded[target as usize] |= 1 << sum.num.to_usize().unwrap();
        }
    }
    extend_exclusions(perm, off, idx + 1, points, num_points, denoms, excluded);
}

fn count_symmetries(perm: &Vec<i32>) -> i64 {
    if perm.len() < 2 { return 1; }
    let mut cmp = 0;
    let (mut i, mut j) = (0, perm.len() - 1);
    while i <= j {
        cmp = perm.len() as i32 - 1 - perm[i] - perm[j];
        if cmp != 0 { break; }
        i += 1; j -= 1;
    }
    match cmp {
        cmp if cmp > 0 => 4,
        cmp if cmp == 0 => 2,
        _ => 0,
    }
}

Cargo.toml

[package]
name = "rust_hello"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
num-bigint = "0.4.2"
num-traits = "0.2.14"

Build and running

# win10 environment

$ cargo build
$ target\debug\rust_hello.exe
0: Ok
1: expected 1, observed 0, error -1
2: expected 2, observed 0, error -2
3: expected 4, observed 8, error 4
4: expected 18, observed 36, error 18
5: expected 48, observed 124, error 76
6: 1008 in 100ms
7: 5024 in 1864ms
8: 64224 in 57079ms
```
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.