# Question
Title: Find the largest prime whose length, sum and product is prime
The number `113` is the first prime whose length `3` is prime, digital sum `5 = 1 + 1 + 3` is prime, and digital product `3 = 1 * 1 * 3` is prime.

A prime that has these 3 properties will be called *supremely prime*. The primes `11117` and `1111151` are other examples.

# Goal

Write a program that can find the largest *supremely prime* number possible in less than an hour on a decent modern personal computer (such as the preferred spec here).

You should not simply give us a large supreme prime. You need to show us your search process with code that actually works. You can build on your or other people's solutions but be sure to give them credit. We're kind of communally trying to find the largest supreme prime realizable on a normal computer in an hour.

# Scoring

The submission that finds the largest supreme prime wins. If it turns out that there are finitely many supreme primes then the first submission that generates the highest supreme prime wins.

(If you can mathematically prove that either there are or aren't infinitely many supreme primes I'll give you 200 bounty rep just because. :) )

# Details

* You may use any source to generate your primes (e.g. internet).
* You may use probabilistic prime testing methods.
* Everything is in base 10.
* Zero and one are NOT considered prime.
* Primes that contain `0` have a digital product of `0` so obviously they can't be supreme.
* To keep the page less cluttered put large (100+ digit) supreme primes in the form:

```
  {[number of 1's before the prime digit]}[prime digit]{[number of 1's after the prime digit]}

```

So `1111151` could be expressed as `{5}5{1}`.

# Answer
> 10 votes
# Perl, 15101 digits, {83}7{15017}, 8 minutes. Max found: 72227 digits

Using my module Math::Prime::Util and its GMP back end. It has a number of compositeness tests including is\_prob\_prime() with an ES BPSW test (slightly more stringent than Pari's ispseudoprime), is\_prime() which adds one random-base M-R, and is\_provable\_prime() which will run BLS75 T5 or ECPP. At these sizes and types, doing a proof is going to take a long time. I threw in another M-R test in the pedantic verifier sub. Times on a Core2 E7500 which is definitely not my fastest computer (it takes 2.5 minutes on my i7-4770K).

As Tim S. points out, we could keep searching for larger values, up to the point where a single test takes an hour. At ~15000 digits on this E7500 it takes about 26s for an M-R test and 2 minutes for the full is\_prime (trial division plus base-2 M-R plus ES Lucas plus one random-base M-R). My i7-4770K is over 3x faster. I tried a few sizes, mainly seeing how it did on other people's results. I tried 8k, 20k, and 16k, killing each after ~5 minutes. I then tried 15k in progression for ~10m each and got lucky on the 4th one.

OpenPFGW's PRP tests are certainly faster once past 4000 or so digits, and much faster indeed in the 50k+ range. Its test has a fair amount of false positives however, which makes it a great pre-test but one would still like to verify the results with something else.

This could be parallelized with perl threads or using MCE similar to the parallel Fibonacci prime finder examples in the module.

Timing and results on an idle i7-4770K using single core:

* input 3000, 16 seconds, 3019 digits, {318}5{2700}
* input 4000, 47 seconds, 4001 digits, {393}7{3607}
* input 4100, 5 seconds, 4127 digits, {29}7{4097}
* input 6217, 5 seconds, 6217 digits, {23}5{6193}
* input 6500, 5 minutes, 6547 digits, {598}5{5948}
* input 7000, 15 minutes, 7013 digits, {2411}7{4601}
* input 9000, 11 minutes, 9001 digits, {952}7{8048}
* input 12000, 10 minutes, 12007 digits, {652}5{11354}
* input 15100, 2.5 minutes, 15101 digits, {83}7{15017}
* input 24600, 47 minutes, 24671 digits, {621}7{24049}
* input 32060, 18 minutes, 32063 digits, {83}7{31979}
* input 57000, 39 minutes, 57037 digits, {112}5{56924}
* input 72225, 42 minutes, 72227 digits, {16}3{72210}

For the 32k digit result, I started 6 scripts running at the same time each with successive arguments starting at 32000. After 26.5 minutes one finished with the 32063 digit result shown. For 57k I let successive scripts run 6 at a time for an hour in input increments of 500 until the 57k result returned in 57 minutes. The 72k digit result was found by doing successive primes from 70k up, so definitely not found in an hour (though once you know where to start, it is).

The script:

```
#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
use Math::Prime::Util::GMP;  # Just to ensure it is used.

my $l = shift || 1000;  $l--;

while (1) {
  $l = next_prime($l);
  my @D = grep { is_prime($l-1 + $_) } (3,5,7);
  next unless scalar @D > 0;
  for my $s (0 .. $l-1) {
    my $e = $l-$s-1;
    warn "   checking $l $s\n" unless $s % 100;
    for my $d (@D) {
      my $n = "1"x$s . $d . "1"x$e;
      die unless length($n) == $l;
      verify_supreme($n,$s,$d,$e) if is_prime($n);  # ES BPSW + 1 rand-base M-R
    }
  }
}
sub verify_supreme {  # Be pedantic and verify the result
  my($n,$s,$d,$e) = @_;
  die "Incorrect length" unless is_prime(length($n));
  die "Incorrect sum" unless is_prime(vecsum(split(//,$n)));
  my $prod = 1; $prod *= $_ for split(//,$n);
  die "Incorrect product" unless is_prime($prod);
  die "n is not a prime!" unless miller_rabin_random($n,1);  # One more M-R test
  die "{$s} $d {$e}\n";
}

```

# Answer
> 21 votes
# Python 2.7 on PyPy, {2404}3{1596} (~10^4000)

```
11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111113111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

```

Found this one about 50 minutes after starting from 4000. Therefore, I would estimate this is the upper limit of this code approach.

Change: I've noticed that some lengths seem to be more fruitful for generating this sort of prime than others, so I've decided to only check 50 random locations of the digit that isn't 1 instead of all possible locations, before moving on. I'm not completely sure this will improve performance, or that 50 is correct, but we'll see.

Possibilities list is generated based on the fact that for the products requirement to be fulfilled, the number must be all ones except for a prime. In addition, the prime can't be 2 because of the sum and length relationship, and the digital sum must not be divisible by three, giving the %3 requirements.

is\_prime taken from http://codepad.org/KtXsydxK , written by @primo

Note: this is\_prime function is actually a Baillie-PSW pseudoprime test, but there are no known counter-examples, so I'm not going to worry about the distinction.

```
#http://codepad.org/KtXsydxK
from my_math import is_prime
import time,random
LLIMIT=2748
time.clock()
start_time=time.time()
checked=0
while time.time()-start_time<3600:
    small_primes = [a for a in range(LLIMIT,2*LLIMIT) if is_prime(a)]
    leng,dig=(0,0)
    for a in small_primes:
        if a+2 in small_primes:
            leng,dig=(a,3)
            break
        if a+4 in small_primes:
            leng,dig=(a,5)
            break
        if a+6 in small_primes:
            leng,dig=(a,7)
            break
    start=time.clock()
    print leng,dig,time.clock(),checked
    for loc in random.sample(range(leng),50):
        checked+=1
        if is_prime(int('1'*loc+str(dig)+'1'*(leng-loc-1))):
            print leng-1,loc,dig,time.clock(),time.clock()-start, \
                  int('1'*loc+str(dig)+'1'*(leng-loc-1))
            break
    LLIMIT=leng+1

```

# Answer
> 13 votes
# PARI/GP, 4127 digits

(104127-1)/9 + 2\*10515

This is a fairly straightforward search: check only prime digit lengths, then compute the possible primes to use, then iterate through all possibilities. I special-cased the common cases where there are 0 or 1 suitable prime digits to use.

```
supreme(lim,startAt=3)={
    forprime(d=startAt,lim,
        my(N=10^d\9, P=select(p->isprime(d+p),[1,2,4,6]), D, n=1);
        if(#P==0, next);
        if(#P==1,
            for(i=0,d-1,
                if (ispseudoprime(D=N+n*P[1]), print(D));
                n*=10
            );
            next
        );
        D=vector(#P-1,i,P[i+1]-P[i]);
        for(i=0,d-1,
            forstep(k=N+n*P[1],N+n*P[#P],n*D,
                if (ispseudoprime(k), print(k))
            );
            n*=10
        )
    )
};
supreme(4200, 4100)

```

This took 36 minutes to compute on one core of a fairly old machine. It would have no trouble finding such a prime over 5000 digits in an hour, I'm sure, but I'm also impatient.

A better solution would be to use any reasonable language to do everything but the innermost loop, then construct an abc file for primeform which is optimized for that particular sort of calculation. This should be able to push the calculation up to at least 10,000 digits.

**Edit**: I implemented the hybrid solution described above, but on my old machine I can't find the first term with \>= 10,000 digits in less than an hour. Unless I run it on something faster I'll have to change to a less lofty target.

# Answer
> 7 votes
# Mathematica 3181 digits

Update: There were some serious mistakes in my first submission. I was able to devote some time to checking the results for this one. The output is formatted as a list of digits. That makes for easy checking of each of the conditions.

```
f[primeDigitLength_]:=
Module[{id=ConstantArray[1,primeDigitLength-1]},
permutations=Reverse@Sort@Flatten[Table[Insert[id,#,pos],{pos,primeDigitLength}]&/@{3,5,7},1];
Flatten[Select[permutations,PrimeQ[FromDigits[#]]\[And]PrimeQ[Plus@@#]&,1],1]]

```

**Example**

This was my first test, a search for a solution with 3181 digits. It found the first case in 26 seconds.

Let's go through the reasoning. Then we'll step into the program itself.

Let's start, as I did, "What is the 450th prime?" Can we find a solution with that many digits (3181)?

```
primeDigits = Prime[450]

```

> 3181

---

The number is found by joining the digits.

```
number = FromDigits[digits];

```

But rather than display it, we can ask instead what the digits are and where they lie.

```
DigitCount[number]

```

> {3180, 0, 0, 0, 0, 0, 1, 0, 0, 0}

This means that there were 3180 instances of the digit 1, and a single instance of the digit 7.

At what position is the digit 7?

```
Position[digits, 7][[1, 1]]

```

> 142

So the digit 7 is the 142nd digit. All the others are 1's.

---

Of course, the product of the digits must be a prime, namely 7.

```
digitProduct = Times @@ digits

```

> 7

---

And the sum of the digits is also a prime.

```
digitSum = Plus @@ digits
PrimeQ[digitSum]

```

> 3187
> True

---

And we know that the number of digits is a prime. Remember, we selected the 450th prime, namely 3118.

So all the conditions have been met.

# Answer
> 7 votes
# Python 2.7, 6217 digits: {23}5{6193} 6 mins 51 secs

I was working on my own version and was disappointed to see that @issacg had beaten me to the punch with a very similar approach, albeit with is\_(very\_probably)\_prime(). However, I see that I have some significant differences that result in a better answer in less time (when I also use is\_prime). To make this clear, when also starting from 4000, I arrive at a better 4001 digit answer ({393}7{3607}) in only 26 mins, 37 seconds using the *standard* Python interpreter (also at version 2.7), **not** the PyPy version. Also, I'm not 'spot' checking the numbers. All of the candidates are checked.

Here's the primary improvements:

1. Use a prime generator (https://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618) to create a list of primes to check against and (his version of "small primes") and for generating eligible number lengths.
2. We want to spend our time looking for the largest supreme prime of a given length, not the smallest, so, I construct the largest possible numbers first for checking, not the smallest. Then, once one is found, we can immediately move on to the next length.

# EDIT: Now with multiprocessing

This is a significant change to previous versions. Before, I noticed that my 8-core machine was hardly working, so, I decided to try my hand at multiprocessing in Python (first time). The results are very nice!

In this version, 7 children processes are spawned, which grab a 'task' off a queue of potential possibilities (num\_length + eligible digits). They churn through trying different \[7,5,3\] positions until it finds one. If it does, informs the master process of the new longest length that has been found. If children are working on a num\_length that is shorter, they just bail, and go get the next length.

I started this run with 6000, and it is still running, but so far, I'm very pleased with the results.

The program doesn't yet stop correctly, but not a huge deal to me.

Now the code:

```
#!/usr/bin/env python
from __future__ import print_function

import sys
from multiprocessing import Pool, cpu_count, Value
from datetime import datetime, timedelta

# is_prime() et al from: http://codepad.org/KtXsydxK - omitted for brevity
# gen_primes() from: https://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618 - ommitted for brevity
from external_sources import is_prime, gen_primes

def gen_tasks(start_length):
    """
    A generator that produces a stream of eligible number lengths and digits
    """
    for num_length in gen_primes():
        if num_length < start_length:
            continue

        ns = [ n for n in [7,5,3] if num_length + n - 1 in prime_list ]
        if ns:
            yield (num_length, ns)

def hunt(num_length, ns):
    """
    Given the num_length and list of eligible digits to try, build combinations
    to try, and try them.
    """

    if datetime.now() > end_time or num_length <= largest.value:
        return

    print('Tasked with: {0}, {1}'.format(num_length, ns))
    sys.stdout.flush()
    template = list('1' * num_length)
    for offset in range(num_length):
        for n in ns:
            if datetime.now() > end_time or num_length <= largest.value:
                return

            num_list = template[:]
            num_list[offset] = str(n)
            num = int(''.join(num_list))

            if is_prime(num):
                elapsed = datetime.now() - start_time
                largest.value = num_length
                print('\n{0} - "{1}"\a'.format(elapsed, num))

if __name__ == '__main__':
    start_time = datetime.now()
    end_time = start_time + timedelta(seconds=3600)

    print('Starting @ {0}, will stop @ {1}'.format(start_time, end_time))

    start_length = int(sys.argv[1])

    #
    # Just create a list of primes for checking. Up to 20006 will cover the first
    # 20,000 digits of solutions
    #
    prime_list = []
    for prime in gen_primes():
        prime_list.append(prime)
        if prime > 20006:
            break;
    print('prime_list is primed.')

    largest = Value('d', 0)

    task_generator = gen_tasks(start_length)

    cores = cpu_count()
    print('Number of cores: {0}'.format(cores))

    #
    # We reduce the number of cores by 1 because __main__ is another process
    #
    pool = Pool(processes=cores - 1)

    while datetime.now() < end_time:
        pool.apply_async(hunt, next(task_generator))

```

# Answer
> 6 votes
## C, GMP \- {7224}5{564} = 7789

Kudos to @issacg and all of you guys for the inspirations and tricks.
And also the masterful question asker @Calvin's Hobbies for this question.

Compile: `gcc -I/usr/local/include -o p_out p.c -pthread -L/usr/local/lib -lgmp`

If you feel like donating your computation power or curious of the performance, feel free copy the code and compile. ;) You will need GMP installed.

```
#include<stdio.h>
#include<stdlib.h>
#include<sys/time.h>
#include<gmp.h>
#include<pthread.h>

#define THREAD_COUNT 1
#define MAX_DIGITS   7800
#define MIN_DIGITS   1000

static void huntSupremePrime(int startIndex) {

    char digits[MAX_DIGITS + 1];

    for (int i = 0; i < MAX_DIGITS; digits[i++] = '1');

    digits[MAX_DIGITS] = '\0';
    mpz_t testPrime, digitSum, digitCount, increment;

    for (int j = 0; j < MAX_DIGITS - startIndex; digits[j++] = '0');

    int step = THREAD_COUNT * 2;

    for (int i = startIndex, l = MAX_DIGITS - startIndex; i > MIN_DIGITS - 1;
        i -= step, l += step) {
        fprintf(stderr, "Testing for %d digits.\n", i);
        mpz_init_set_ui(digitCount, i);
        if (mpz_millerrabin(digitCount, 10)) {
            for (int j = 3; j < 8; j += 2) {
                mpz_init_set_ui(digitSum, i - 1 + j);
                if (mpz_millerrabin(digitSum, 10)) {
                    gmp_printf("%Zd \n", digitSum);
                    digits[MAX_DIGITS - 1] = j + 48;
                    mpz_init_set_str(testPrime, digits, 10);
                    mpz_init_set_ui(increment, (j - 1) * 99);
                    for (int k = 0; k < i/20; ++k) {
                        if (mpz_millerrabin(testPrime, 25)) {
                            i = 0;
                            j = 9;
                            k = l;
                            gmp_printf("%Zd\n", testPrime);
                            break;
                        }
                        mpz_add(testPrime, testPrime, increment);
                        mpz_mul_ui(increment, increment, 100);
                        fprintf(stderr, "TICK %d\n", k);
                    }

                }
            }
        }
        for (int j = 0; j < step; digits[l + j++] = '0');

    }
}

static void *huntSupremePrimeThread(void *p) {
    int* startIndex = (int*) p;
    huntSupremePrime(*startIndex);
    pthread_exit(NULL);
}

int main(int argc, char *argv[]) {

    int  startIndexes[THREAD_COUNT];
    pthread_t threads[THREAD_COUNT];

    int startIndex = MAX_DIGITS;
    for (int i = 0; i < THREAD_COUNT; ++i) {
        for (;startIndex % 2 == 0; --startIndex);
        startIndexes[i] = startIndex;
        int rc = pthread_create(&threads[i], NULL, huntSupremePrimeThread, (void*)&startIndexes[i]);
        if (rc) {
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
        --startIndex;
    }

    for (int i = 0; i < THREAD_COUNT; ++i) {
        void * status;
        int rc = pthread_join(threads[i], &status);
        if (rc) {
            printf("ERROR: return code from pthread_join() is %d\n", rc);
            exit(-1);
        }
    }

    pthread_exit(NULL);
    return 0;
}

```

# Answer
> 5 votes
## PFGW, 6067 digits, {5956}7{110}

Run PFGW with the following input file and `-f100` to prefactor the numbers. In about 2-3 CPU minutes on my computer (i5 Haswell), it finds the PRP (10^(6073-6)-1)/9+6\*10^110, or {5956}7{110}. I chose 6000 digits as the starting point as a nothing-up-my-sleeve number that's a little higher than all previous submissions.

```
ABC2 $a-$b & (10^($a-$b)-1)/9+$b*10^$c
a: primes from 6000 to 6200
b: in { 2 4 6 }
c: from 0 to 5990

```

Based on how quickly I was able to find this one, I could easily bump up the number of digits and still find a PRP within an hour. With how the rules are written, I might even just find the size where my CPU, running on all 4 cores, can finish one PRP test in an hour, take a *long* time to find a PRP, and have my "search" consist solely of the one PRP.

P.S. In some senses, this isn't a "code" solution because I didn't write anything besides the input file...but then, many one-line Mathematica solutions to mathematical problems could be described in the same way, as could using a library that does the hard work for you. In reality, I think it's hard to draw a good line between the two. If you like, I could write a script that creates the PFGW input file and calls PFGW. The script could even search in parallel, to use all 4 cores and speed up the search by ~4 times (on my CPU).

P.P.S. I think LLR can do the PRP tests for these numbers, and I'd expect it to be far faster than PFGW. A dedicated sieving program could also be better at factoring these numbers than PFGW's one-at-a-time. If you combined these, I'm sure you could push the bounds much higher than current solutions.

# Answer
> 4 votes
# Python 2.7, 17-19 digits

```
11111111171111111

```

Found 5111111111111 (13 digits) in 3 seconds and this 17 digit supreme prime in 3 minutes. I'll guess that the target machine could run this and get a 19 digit supreme prime in less than an hour. This approach does not scale well because it keeps primes up to half the number of target digits in memory. 17 digit search requires storing an array of 100M booleans. 19 digits would require a 1B element array, and memory would be exhausted before getting to 23 digits. Runtime probably would be, too.

Primality test approaches that don't involve a ridiculously large array of divisor primes will fare much better.

```
#!/usr/bin/env python
import math
import numpy as np
import sys

max_digits = int(sys.argv[1])
max_num = 10**max_digits

print "largest supreme prime of " + str(max_digits) + " or fewer digits"

def sum_product_digits(num):
    add = 0
    mul = 1
    while num:
         add, mul, num = add + num % 10, mul * (num % 10), num / 10
    return add, mul

def primesfrom2to(n):
    # http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n/3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)/3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

def checkprime(n):
    for divisor in primes:
        if (divisor>math.sqrt(n)):
            break
        if n%divisor==0:
            return False
    return True

# make an array of all primes we need to check as divisors of our max_num
primes = primesfrom2to(math.sqrt(max_num))
# only consider digit counts that are prime
for num_digits in primes:
    if num_digits > max_digits:
        break
    for ones_on_right in range(0,num_digits):
        for mid_prime in ['3','5','7']:
            # assemble a number of the form /1*[357]1*/
            candidate = int('1'*(num_digits-ones_on_right-1)+mid_prime+'1'*ones_on_right)
            # check for primeness of digit sum first digit product first
            add, mul = sum_product_digits(candidate)
            if add in primes and mul in primes:
                # check for primality next
                if checkprime(candidate):
                    # supreme prime!
                    print candidate

```

# Answer
> 3 votes
# Mathematica 4211 4259 digits

With the number: {1042}7{3168} {388}3{3870}

Which was generated by the following code:

```
TimeConstrained[
 Do[
  p = Prime[n];
  curlargest = Catch[
    If[PrimeQ[p + 6],
     list = ConstantArray[1, p];
     Do[
      temp = FromDigits[ReplacePart[list, i -> 7]];
      If[PrimeQ[temp],
       Throw[temp]
       ], {i, p}]
     ];

    If[PrimeQ[p + 4],
     list = ConstantArray[1, p];
     Do[
      temp = FromDigits[ReplacePart[list, i -> 5]];
      If[PrimeQ[temp],
       Throw[temp]
       ], {i, p}]
     ];
    If[PrimeQ[p + 2],
     list = ConstantArray[1, p];
     Do[
      temp = FromDigits[ReplacePart[list, i -> 3]];
      If[PrimeQ[temp],
       Throw[temp]
       ], {i, p}]
     ];
    Throw[curlargest];
    ]

  , {n, 565, 10000}]
 , 60*60]

```

The throws cause it to stop testing for other numbers with the same digits as the one currently found. since it begins testing at the most significant digit this will mean that it always returns the largest number unless the number of digits is a member of a prime triplet.

Simply started testing just below the value of one of the preceding answers :)

Once it finishes, the number is stored in the variable curlargest

# Answer
> 3 votes
# JavaScript, 3019 digits, {2,273}5{745}

This uses the MillerRabin test included in BigInteger.js by Tom Wu.

Starting from **0 =\> 2,046 digits** = {1799}7{263} **in one hour**.

Starting from **3000 =\> 3,019 digits** = {2,273}5{745} **in one hour, less 3 seconds**.

When it started from 0, the program skipped ahead and began searching again at a length of 1.5X the length of the last s-prime found. Then when I saw how fast it was running I guessed it would find one starting at 3000 in one hour - which it did with only 3 seconds to spare.

You can try it here: http://goo.gl/t3TmTk
(Set to calculate all s-primes, or skip ahead.)

The program works by constructing strings of all "1"s, but with one "3", "5", or "7". I added a quick check in the IsStrPrime function to reject numbers ending in "5".

```
if (IsIntPrime(length)) {

    var do3s = IsIntPrime(length + 2);
    var do5s = IsIntPrime(length + 4);
    var do7s = IsIntPrime(length + 6);

    if (do3s || do5s || do7s) {

        // loop through length of number
        var before, digit, after;

        for (var after = 0; after <= length - 1; after++) {

            before = length - after - 1;
            beforeStr = Ones(before);
            afterStr = Ones(after);

            if (do3s && IsStrPrime(beforeStr + (digit = "3") + afterStr)) { RecordAnswer(); if (brk) break; }
            if (AreDone()) break;

            if (do5s && IsStrPrime(beforeStr + (digit = "5") + afterStr)) { RecordAnswer(); if (brk) break; }
            if (AreDone()) break;

            if (do7s && IsStrPrime(beforeStr + (digit = "7") + afterStr)) { RecordAnswer(); if (brk) break; }
            if (AreDone()) break;

            if (after % 10 == 0) document.title = "b=" + bestLength + ", testing=" + length + "-" + after;
        }
    }
}

```

This was pretty fun. Reminds me of a puzzle I did many years back to calculate what's called a **digit removed prime**. This is a prime number that if you remove any digit, then the remaining number is still prime. For example 1037 is a digit removed prime because 1037, 037, 137, 107, and 103 are prime. I found one 84 digits long, and the longest I know of is 332 digits long. I'm sure we could find one much longer with the techniques used for this puzzle. (But choosing the trial numbers is a little bit trickier - maybe?)

# Answer
> 1 votes
# Python 3.9.2, with sympy and gmpy2, 3019 digits

Completed the search **from 3000 digits to 3019 digits in approximately 10 minutes** on my laptop.

Port of @Charles's Pari/GP answer in Python.

```
import sympy
import gmpy2
import time

def supreme(lim, startAt=3):
    for d in sympy.primerange(startAt, lim + 1):
        N = 10**d // 9
        P = [p for p in [1,2,4,6] if sympy.isprime(d + p)]
        if len(P) == 0:
            continue
        if len(P) == 1:
            n = 1
            for i in range(0, d):
                if gmpy2.is_strong_bpsw_prp(N + n * P[0]):
                    print(N + n * P[0])
                n *= 10
            continue
        D = [P[i+1] - P[i] for i in range(len(P) - 1)]
        n = 1
        for i in range(0, d):
            for k in range(N + n * P[0], N + n * P[-1] + 1, n * D[0]):
                if gmpy2.is_strong_bpsw_prp(k):
                    print(k)
            n *= 10

start_time=time.time()
supreme(3019, 3000)
print(f"The search ends in {time.time()-start_time} seconds")

```

```
[encoded data, 6038 characters]
The search ends in 606.180704832077 seconds

```

# Answer
> 0 votes
# Axiom, 3019 digits {318}5{2700}

```
)time on

-- Return true if n is pseudo prime else return false
-- **Can Fail**
prime1(n:PI):Boolean==
     n=1=>false
     n<4=>true
     i:=5;sq:=sqrt(1.*n)
     repeat
       if i>sq or i>50000 then break
       if n rem i=0       then return false
       i:=i+2
     if i~=50001        then return true
     --output("i")
     if powmod(3,n,n)=3 then return true
     --output("e")
     false

-- input  'n': must be n>1 prime
-- output [0] if not find any number, else return
-- [n,a,b,c,z] 'n' digits of solution,
-- 'a' number of '1', 'b' central digit, 'b' number of final digit '1'
-- 'z' the number found
g(n:PI):List NNI==
    x:=b:=z:=1
    for i in 1..n-1 repeat
        z:=z*10+1
        b:=b*10
    repeat
       --output b
       k:=0    -- 3 5 7 <-> 2 4 6
       for i in [2,4,6] repeat
           ~prime?(n+i)=>0 --somma
           k:=k+1
           t:=z+b*i
           if prime1(t) then return [n,x-1,i+1,n-x,t]
       --if x=1 then output ["k=", k]
       if k=0  then break
       x:=x+1
       b:=b quo 10
       if b<=0 then break
    [0]

-- start from number of digits n
-- and return g(i) with i prime i>=n
find(n:PI):List NNI==
    i:=n
    if i rem 2=0 then i:=i+1
    repeat
        if prime?(i) then --solo le lunghezze prime sono accettate
             output i
             a:=g(i)
             if a~=[0] then return a
        i:=i+2

```

result from the start value 3000 in 529 sec

```
(4) -> find(3000)
   3001
   3011
   3019

   (4)
   [3019, 318, 5, 2700, Omissis]
                                            Type: List NonNegativeInteger
       Time: 0.02 (IN) + 525.50 (EV) + 0.02 (OT) + 3.53 (GC) = 529.07 sec

```

---
Tags: primes, fastest-code
---


# Question
Title: Extending OEIS: Counting Diamond Tilings
I promise, this will be my last challenge about diamong tilings (for a while, anyway). On the bright side, this challenge doesn't have anything to do with ASCII art, and is not a code golf either, so this is actually completely different.

So just as a reminder, every hexagon can be titled with three different diamonds:

An interesting question to ask is how many of these tilings exist for a given hexagon size. It seems these numbers have been studied fairly thoroughly and can be found in OEIS A008793.

However, the problem gets trickier if we ask how many tilings exist *up to rotation and reflection*. For instance, for side length N = 2, the following 20 tilings exist:

```
   ____     ____     ____     ____     ____     ____     ____     ____     ____     ____
  /\_\_\   /\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
 /\/\_\_\ /\/_/\_\ /\/_/_/\ /\/_/\_\ /\_\/\_\ /\_\/_/\ /\/_/_/\ /\_\/\_\ /\_\/_/\ /_/\/\_\
 \/\/_/_/ \/\_\/_/ \/\_\_\/ \/_/\/_/ \/\_\/_/ \/\_\_\/ \/_/\_\/ \/_/\/_/ \/_/\_\/ \_\/\/_/
  \/_/_/   \/_/_/   \/_/_/   \_\/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/   \_\/_/   \_\/_/
   ____     ____     ____     ____     ____     ____     ____     ____     ____     ____
  /_/_/\   /\_\_\   /_/\_\   /_/_/\   /_/\_\   /_/\_\   /_/_/\   /_/_/\   /_/_/\   /_/_/\
 /\_\_\/\ /\/_/_/\ /_/\/_/\ /\_\_\/\ /\_\/_/\ /_/\/_/\ /_/\_\/\ /\_\_\/\ /_/\_\/\ /_/_/\/\
 \/\_\_\/ \/_/_/\/ \_\/\_\/ \/_/\_\/ \/_/_/\/ \_\/_/\/ \_\/\_\/ \/_/_/\/ \_\/_/\/ \_\_\/\/
  \/_/_/   \_\_\/   \_\/_/   \_\/_/   \_\_\/   \_\_\/   \_\/_/   \_\_\/   \_\_\/   \_\_\/

```

But many of these are identical under rotation and reflection. If we take these symmetries into account, only *6* distinct tilings remain:

```
   ____     ____     ____     ____     ____     ____
  /\_\_\   /\_\_\   /\_\_\   /_/\_\   /_/\_\   /_/\_\
 /\/\_\_\ /\/_/\_\ /\/_/_/\ /\_\/_/\ /\_\/_/\ /_/\/\_\
 \/\/_/_/ \/\_\/_/ \/\_\_\/ \/\_\_\/ \/_/\_\/ \_\/\/_/
  \/_/_/   \/_/_/   \/_/_/   \/_/_/   \_\/_/   \_\/_/

   2        2        6        6        1        3

```

where the numbers indicate the multiplicity of each tiling. Note that for larger hexagons there are also tilings with multiplicities 4 and 12.

It appears that the number of tilings up to symmetry has been studied less thoroughly. The OEIS entry A066931 only lists the five terms:

```
1, 1, 6, 113, 20174

```

where the first term is for side length `N = 0` and the last term for side length `N = 4`.

I'm sure we can do better than that!

Your task is to compute the number of tilings for a given side length.

This is . Your score will be the highest side-length `N` for which your code produces the correct result within *30 minutes* on my machine. In case of a tie, I will accept the submission which produces the result for *that* `N` fastest.

As usual, you must not hardcode results you already know to win the tie-breaker. The algorithm that solve `N = 3` should be identical to the one that solves `N = 5`.

Your submission must not use more than 4GB of memory. I will give some leeway on this if you're operating close to that limit, but if you're consistently above that limit, or if you spike significantly beyond it, I will not count that `N` for your submission.

I will test all submissions on my Windows 8 machine, so make sure your language of choice is freely available on Windows. The only exception to this is Mathematica (because I happen to have a licence for it). Please include instructions for how to compile/run your code.

Of course, feel free to compute more terms in your own time (for science, and for others to check their numbers against), but your answer's *score* will be determined in those 30 minutes.

# Answer
## Algebra, graph theory, Möbius inversion, research, and Java

The symmetry group of the hexagon is the dihedral group of order 12, and is generated by a 60 degree rotation and a mirror flip across a diameter. It has 16 subgroups, but some of them are in non-trivial conjugacy groups (the ones which only have reflections have 3 choices of axis), so there are 10 fundamentally different symmetries which a tiling of the hexagon can have:

The number of diamond tilings of a subset of a triangular lattice can be calculated as a determinant, so my initial approach was to set up one determinant for each of the symmetries of the hexagon, to calculate the number of tilings which have *at least* those symmetries; and then use Möbius inversion in the incidence algebra of their poset (basically a generalisation of the inclusion-exclusion principle) to work out the number of tilings whose symmetry group is *exactly* each of the 10 cases. However, some of the symmetries have nasty edge conditions, so I was forced to sum over exponentially many determinants. Fortunately, the values obtained for `n < 10` gave me enough data to be able to identify relevant sequences in OEIS and piece together a closed form (for some value of "closed" which allows finite products). There's a bit of discussion of the sequences, and references to proofs, in the formal writeup which I prepared to justify the OEIS sequence updates.

Once the double-counting is taken care of, it turns out that four of the ten values cancel out neatly, so we only have to calculate the remaining six and then do a weighted sum.

This code takes under 30 seconds for `N=1000` on my machine.

```
import java.math.BigInteger;

public class OptimisedCounter {
    private static int[] minp = new int[2];

    public static void main(String[] args) {
        if (args.length > 0) {
            for (String arg : args) System.out.println(count(Integer.parseInt(arg)));
        }
        else {
            for (int n = 0; n < 16; n++) {
                System.out.format("%d\t%s\n", n, count(n));
            }
        }
    }

    private static BigInteger count(int n) {
        if (n == 0) return BigInteger.ONE;

        if (minp.length < 3*n) {
            int[] wider = new int[3*n];
            System.arraycopy(minp, 0, wider, 0, minp.length);
            for (int x = minp.length; x < wider.length; x++) {
                // Find the smallest prime which divides x
                for (wider[x] = 2; x % wider[x] != 0; wider[x]++) { /* Do nothing */ }
            }
            minp = wider;
        }

        BigInteger E = countE(n), R2 = countR2(n), F = countF(n), R3 = countR3(n), R = countR(n), FR = countFR(n);
        BigInteger sum = E.add(R3);
        sum = sum.add(R2.add(R).multiply(BigInteger.valueOf(2)));
        sum = sum.add(F.add(FR).multiply(BigInteger.valueOf(3)));
        return sum.divide(BigInteger.valueOf(12));
    }

    private static BigInteger countE(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j <= i + n; j++) w[j]--;
            for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countR2(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            w[3*i+2]++;
            for (int j = 3*i + 1; j <= 2*i + n + 1; j++) w[j]--;
            for (int j = 2*i + n + 1; j <= i + n + n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countF(int n) {
        int[] w = new int[3*n];
        for (int i = 0; i < n; i++) {
            for (int j = 2*i + 1; j <= 2*i + n; j++) w[j]--;
            for (int j = i + n + 1; j <= i + 2*n; j++) w[j]++;
        }
        return powerProd(w);
    }

    private static BigInteger countR3(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        return countE(n / 2).pow(2);
    }

    private static BigInteger countR(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        int m = n / 2;
        int[] w = new int[3*m-1];
        for (int i = 0; i < m; i++) {
            for (int j = 1; j <= 3*i+1; j++) w[j] += 2;
            for (int j = 1; j <= i + m; j++) w[j] -= 2;
        }
        return powerProd(w);
    }

    private static BigInteger countFR(int n) {
        if ((n & 1) == 1) return BigInteger.ZERO;
        int m = n / 2;
        int[] w = new int[3*n-2];
        for (int j = 1; j <= m; j++) w[j]--;
        for (int j = 2*m; j <= 3*m-1; j++) w[j]++;
        for (int i = 0; i <= 2*m-3; i++) {
            for (int j = i + 2*m + 1; j <= i + 4*m; j++) w[j]++;
            for (int j = 2*i + 3; j <= 2*i + 2*m + 2; j++) w[j]--;
        }
        return powerProd(w);
    }

    private static BigInteger powerProd(int[] w) {
        BigInteger result = BigInteger.ONE;
        for (int x = w.length - 1; x > 1; x--) {
            if (w[x] == 0) continue;

            int p = minp[x];
            if (p == x) result = result.multiply(BigInteger.valueOf(p).pow(w[p]));
            else {
                // Redistribute it. This should ensure we avoid negatives.
                w[p] += w[x];
                w[x / p] += w[x];
            }
        }

        return result;
    }
}

```
> 89 votes

# Answer
# C

**Introduction**

As commented by David Carraher, the simplest way of analysing the hexagon tiling seemed to be to take advantage of its isomorphism with the 3 dimensional Young Diagram, essentially an x,y square filled with integer height bars whose z heights must stay the same or increase as the z axis is approached.

I started with an algorithm for finding the totals which is more amenable to adaptation for symmetry counting than the published algorithm, which is based on a bias to one of the three cartesian axes.

**Algorithm**

I start by filling the cells of the x,y, and z planes with 1's, while the rest of the area contains zeroes. Once that is done, I build up the pattern layer by layer, with each layer containing the cells which have a common 3D manhattan distance from the origin. A cell can only contain a 1 if the three cells below it also contain a 1. if any of them contains a 0, then the cell must be a 0.

The advantage of building the pattern up in this way is that each layer is symmetrical about the x=y=z line. This means that each layer can be checked independently for symmetry.

**Symmetry checking**

The symmetries of the solid are as follows: 3 fold rotation about the x=y=z line --\> 3 fold rotation about the hexagon centre; and 3 x reflections about the 3 planes containing the x=y=z line and each of the axes x,y,z --\> reflection about the lines through the hexagon corners.

This adds up only to 6 fold symmetry. To get the full symmetry of the hexagon, another kind of symmetry must be considered. Each solid (built from 1's) has a complementary solid (built from 0's). Where N is odd, the complementary solid must be different from the original solid (because it is not possible for them to have the same number of cubes). Yet when the complementary solid is turned round, it will be found that its 2D representation as a diamond tiling is identical (except for a 2-fold symmetry operation) to the original solid. Where N is even, it is possible for the solid to be self-inverse.

This can be seen in the examples for N=2 in the question. If viewed from the left, the first hexagon looks like a solid cube with 8 small cubes, while the last hexagon looks like an empty shell with 0 small cubes. If viewed from the right, the reverse is true. The 3rd, 4th and 5th hexagons and the 16th, 17th and 18th hexagons look like they contain either 2 or 6 cubes, and thus they complement each other in 3 dimensions. They are related to each other in 2 dimensions by a 2-fold symmetry operation (2-fold rotation, or reflection about an axis through the hexagon's edges.) On the other hand the 9th,10th,11th and 12th hexagons show 3D patterns which are their own complements, and therefore have a higher symmetry (these are therefore the only patterns with odd multiplicity).

Note that having (N^3)/2 cubes is a necessary condition to be self complement, but in general it is not a sufficient condition if N\>2. The result of all this is that for odd N, tilings always occur in pairs (N^3)/2 cubes must be carefully inspected.

**Current code (generates the right total for N=1,2,3,5. Error as discussed for N=4.)**

```
int n;                     //side length

char t[11][11][11];        //grid sized for N up to 10

int q[29][192], r[29];     //tables of coordinates for up to 10*3-2=28 layers

int c[9];                  //counts arrangements found by symmetry class. c[8] contains total.

//recursive layer counting function. m= manhattan distance, e= number of cells in previous layers, s=symmetry class.
void f(int m,int e,int s){

  int u[64], v[64], w[64]; //shortlists for x,y,z coordinates of cells in this layer
  int j=0;
  int x,y,z;

  for (int i=r[m]*3; i; i-=3){
    // get a set of coordinates for a cell in the current layer.
    x=q[m][i-3]; y= q[m][i-2]; z= q[m][i-1];
    // if the three cells in the previous layer are filled, add it to the shortlist u[],v[],w[]. j indicates the length of the shortlist.
    if (t[x][y][z-1] && t[x][y-1][z] && t[x-1][y][z]) u[j]=x, v[j]=y, w[j++]=z ;
  }

  // there are 1<<j possible arrangements for this layer.
  for (int i = 1 << j; i--;) {

    int d = 0;

    // for each value of i, set the 1's bits of t[] to the 1's bits of i. Count the number of 1's into d as we go.
    for (int k = j; k--;) d+=(t[u[k]][v[k]][w[k]]=(i>>k)&1);

    // we have no interest in i=0 as it is the empty layer and therefore the same as the previous recursion step.
    // Still we loop through it to ensure t[] is properly cleared.

    if(i>0){
      int s1=s;    //local copy of symmetry class. 1's bit for 3 fold rotation, 2's bit for reflection in y axis.
      int sc=0;    //symmetry of self-complement.

      //if previous layers were symmetrical, test if the symmetry has been reduced by the current layer
      if (s1) for (int k = j; k--;) s1 &= (t[u[k]][v[k]][w[k]]==t[w[k]][u[k]][v[k]]) | (t[u[k]][v[k]][w[k]]==t[w[k]][v[k]][u[k]])<<1;

      //if exactly half the cells are filled, test for self complement
      if ((e+d)*2==n*n*n){
        sc=1;
        for(int A=1; A<=(n>>1); A++)for(int B=1; B<=n; B++)for(int C=1; C<=n; C++) sc&=t[A][B][C]^t[n+1-A][n+1-B][n+1-C];
      }

      //increment counters for total and for symmetry class.
      c[8]++; c[s1+(sc<<2)]++;

      //uncomment for graphic display of each block stacking with metadata. not recommended for n>3.
      //printf("m=%d  j=%d  i=%d c1=%d-2*%d=%d c3=%d cy=%d(cs=%d) c3v=%d ctot=%d\n",m,j,i,c[0],c[2],c[0]-2*c[2],c[1],c[2],c[2]*3,c[3],c[8]);
      //printf("m=%d  j=%d  i=%d C1=%d-2*%d=%d C3=%d CY=%d(CS=%d) C3V=%d ctot=%d\n",m,j,i,c[4],c[6],c[4]-2*c[6],c[5],c[6],c[6]*3,c[7],c[8]);
      //for (int A = 0; A<4; A++, puts(""))for (int B = 0; B<4; B++, printf(" "))for (int C = 0; C<4; C++) printf("%c",34+t[A][B][C]);

      //recurse to next level.
      if(m<n*3-2)f(m + 1,e+d,s1);

    }
  }
}

main()
{
  scanf("%d",&n);

  int x,y,z;

  // Fill x,y and z planes of t[] with 1's
  for (int a=0; a<9; a++) for (int b=0; b<9; b++) t[a][b][0]= t[0][a][b]= t[b][0][a]= 1;

  // Build table of coordinates for each manhattan layer
  for (int m=1; m < n*3-1; m++){
    printf("m=%d : ",m);
    int j=0;
    for (x = 1; x <= n; x++) for (y = 1; y <= n; y++) {
      z=m+2-x-y;
      if (z>0 && z <= n) q[m][j++] = x, q[m][j++] = y, q[m][j++]=z, printf(" %d%d%d ",x,y,z);
      r[m]=j/3;
    }
    printf(" : r=%d\n",r[m]);
  }

  // Set count to 1 representing the empty box (symmetry c3v)
  c[8]=1; c[3]=1;

  // Start searching at f=1, with 0 cells occupied and symmetry 3=c3v
  f(1,0,3);

  // c[2 and 6] only contain reflections in y axis, therefore must be multiplied by 3.
  // Similarly the reflections in x and z axis must be subtracted from c[0] and c[4].
  c[0]-=c[2]*2; c[2]*=3;
  c[4]-=c[6]*2; c[6]*=3;

  int cr[9];cr[8]=0;
  printf("non self-complement                   self-complement\n");
  printf("c1  %9d/12=%9d           C1  %9d/6=%9d\n",   c[0], cr[0]=c[0]/12,     c[4], cr[4]=c[4]/6);
  if(cr[0]*12!=c[0])puts("c1 division error");if(cr[4]*6!=c[4])puts("C1 division error");

  printf("c3  %9d/4 =%9d           C3  %9d/2=%9d\n",   c[1], cr[1]=c[1]/4,      c[5], cr[5]=c[5]/2);
  if(cr[1]*4!=c[1])puts("c3 division error");if(cr[5]*2!=c[5])puts("C3 division error");

  printf("cs  %9d/6 =%9d           CS  %9d/3=%9d\n",   c[2], cr[2]=c[2]/6,      c[6], cr[6]=c[6]/3);
  if(cr[2]*6!=c[2])puts("cs division error");if(cr[6]*3!=c[6])puts("CS division error");

  printf("c3v %9d/2 =%9d           C3V %9d/1=%9d\n",   c[3], cr[3]=c[3]/2,      c[7], cr[7]=c[7]);
  if(cr[3]*2!=c[3])puts("c3v division error");

  for(int i=8;i--;)cr[8]+=cr[i];
  printf("total =%d unique =%d",c[8],cr[8]);
}

```

**Output**

The program generates an output table of 8 entries, in accordance with the 8 symmetries of the solid. The solid can have any of 4 symmetries as follows (Schoenflies notation)

```
c1: no symmetry
c3: 3-fold axis of rotation (produces 3-fold axis of rotation in hexagon tiling)
cs: plane of reflection (produces line of reflection in hexagon tiling)
c3v both of the above (produces 3-fold axis of rotation and three lines of reflection through the hexagon corners)

```

Additionally, when the solid has exactly half the cells with 1's and half with 0's, there exists the possibility of flipping all the 1's and 0's, then inverting the coordinates through the centre of the cube space. This is what I am calling self-complement, but a more mathematical term would be "antisymmetric with respect to a centre of inversion."

This symmetry operation gives a 2-fold axis of rotation in the hexagon tiling.

Patterns which have this symmetry are listed in a separate column. They only occur where N is even.

My count seems to be slightly off for N=4. In discussion with Peter Taylor it appears that I am not detecting tilings which only have a symmetry of a line through the hexagon edges. This is presumably because I have not tested for self complement (antisymmetry) for operations other than (inversion)x(identity.) Testing for self complement for the operatons (inversion)x(reflection) and (inversion)x(3-fold rotation) may uncover the missing symmetries. I would then expect the first line of the data for N=4 to look like this (16 less in c1 and 32 more in C1):

```
c1   224064/12=18672          C1  534/6=89

```

This would bring the totals in line with Peter's answer and https://oeis.org/A066931/a066931.txt

current output is as follows.

```
N=1
non self-complement     self-complement
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs      0/6 = 0           CS  0/3= 0
c3v     2/2 = 1           C3V 0/1= 0
total =2 unique =1

non self-complement     self-complement
N=2
c1      0/12= 0           C1  0/6= 0
c3      0/4 = 0           C3  0/2= 0
cs     12/6 = 2           CS  3/3= 1
c3v     4/2 = 2           C3V 1/1= 1
total =20 unique =6

N=3
non self-complement     self-complement
c1    672/12=56           C1  0/6= 0
c3      4/4 = 1           C3  0/2= 0
cs    288/6 =48           CS  0/3= 0
c3v    16/2 = 8           C3V 0/1= 0
total =980 unique =113

N=4 (errors as discussed)
non self-complement     self-complement
c1   224256/12=18688          C1  342/6=57
c3       64/4 =16             C3  2/2= 1
cs     8064/6 =1344           CS  54/3=18
c3v      64/2 =32             C3V 2/1= 2
total =232848 unique =20158

N=5
non self-complement     self-complement
c1  266774112/12=22231176        C1  0/6= 0
c3       1100/4 =275             C3  0/2= 0
cs     451968/6 =75328           CS  0/3= 0
c3v       352/2 =176             C3V 0/1= 0
total =267227532 unique =22306955

```

**To-do list (updated)**

> Tidy up current code.

Done, more or less

> Implement symmetry checking for current layer, and pass a parameter for the symmetry of the previous layer (no point in checking if the last layer was asymmetrical.)

Done, results for odd N agree with published data

> Add an option to suppress counting asymmetrical figures (should run much faster)

This can be done by adding another condition to the recursion call: `if(s1 && m<n*3-2)f(m + 1,e+d,s1)` It reduces the run time for N=5 from 5 minutes down to about a second. As a result the first line of the output becomes total garbage (as do the overall totals) but if the total is already known from OEIS the number of asymmetric tilings can be reconstituted, at least for odd N.

But for even N, the number of asymmetric (according to c3v symmetries) solids that are self-complement would be lost. For this case, a separate program dedicated to solids with exactly (N\**3)/2 cells with a 1 may be useful. With this available (and counting correctly) it may be possible to try N=6, but it will take a long time to run.

> Implement counting of cells to reduce search to up to (N^3)/2 cubes.

Not done, savings expected to be marginal

> Implement symmetry (complementary solid) checking for patterns containing exactly (N^3)/2 cubes.

Done, but seems to have omissions, see N=4.

> Find a way to pick the lexically lowest figure from an asymetrical one.

Savings are not expected to be that great. Supressing asymmetrical figures eliminates most of this. The only reflection that is checked is the plane through the y axis (x and z are calculated later by multiplying by 3.) Figures with only rotational symmetry are counted in both their enantiomeric forms. Perhaps it would run nearly twice as fast if only one was counted.

> To facilitate this, possibly improve the way the coordinates in each layer are listed (they form degenerate groups of 6 or 3, with possibly a group of 1 in the exact centre of the layer.)

Interesting but probably there are other questions on the site to explore.
> 15 votes

# Answer
# Rust, 0.03 s for n=0..=15

Port of @Peter Taylor's Java answer in Rust.

Run it on RustPlayground!

---

### `src/main.rs`

```
use std::collections::HashMap;
use num_bigint::BigUint;
use num_traits::{Zero, One, Pow};

fn main() {
    for n in 0..16 {
        println!("{}\t{}", n, count(n));
    }
}

fn count(n: usize) -> BigUint {
    if n == 0 { return BigUint::one() }

    let mut minp = vec![0; 3 * n];

    for index in 2..3 * n {
        minp[index] = (2..=index).find(|p| index % p == 0).unwrap_or(2);
    }

    let e = count_e(n, &minp);
    // println!("n={:?}  e={:?}", n, e);
    let r2 = count_r2(n, &minp);
    // println!("n={:?}  r2={:?}", n, r2);
    let f = count_f(n, &minp);
    // println!("n={:?}  f={:?}", n, f);
    let r3 = count_r3(n, &minp);
    // println!("n={:?}  r3={:?}", n, r3);
    let r = count_r(n, &minp);
    // println!("n={:?}  r={:?}", n, r);
    let fr = count_fr(n, &minp);
    // println!("n={:?}  fr={:?}", n, fr);

    let sum = e
        + r3
        + (r2 + r.clone()) * BigUint::from(2u8)
        + (f + fr) * BigUint::from(3u8);

    sum / BigUint::from(12u8)
}

fn count_e(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        for j in (i + 1)..=(i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    // println!("w={:?}  minp={:?}",w, minp);

    power_prod(&mut w, minp)
}

fn count_r2(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        w[3 * i + 2] += 1;
        for j in (3 * i + 1)..=(2 * i + n + 1) { w[j] -= 1 }
        for j in (2 * i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    power_prod(&mut w, minp)
}

fn count_f(n: usize, minp: &[usize]) -> BigUint {
    let mut w = vec![0; 3 * n];

    for i in 0..n {
        for j in (2 * i + 1)..=(2 * i + n) { w[j] -= 1 }
        for j in (i + n + 1)..=(i + 2 * n) { w[j] += 1 }
    }

    power_prod(&mut w, minp)
}

fn count_r3(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let res = count_e(n / 2, minp);
    res.pow(2u8)
}

fn count_r(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m = n / 2;
    // let mut w = vec![0; 3 * m - 1];
    let mut w = vec![0; 3 * m - 1];

    for i in 0..m {
        for j in 1..=(3 * i + 1) { w[j] += 2 }
        for j in 1..=(i + m) { w[j] -= 2 }
    }

    // power_prod(&w, minp)
    power_prod(&mut w, minp)
}

fn count_fr(n: usize, minp: &[usize]) -> BigUint {
    if n & 1 == 1 { return BigUint::zero() }
    let m:i32 = n as i32 / 2;
    let mut w = vec![0; 3 * n - 2];
    // println!("n={:?}   m={:?}", n, m);

    for j in 1..=m { w[j as usize] -= 1 }
    for j in (2 * m)..(3 * m) { w[j as usize] += 1 }
    for i in 0..=(2 * m - 3) {
        // println!("start{:?}    end{:?}", (i + 2 * m + 1),  (i + 4 * m));
        for j in (i + 2 * m + 1)..=std::cmp::min(w.len() as i32 - 1, i + 4 * m) { w[j as usize] += 1 }
        // println!("start{:?}    end{:?}", (2 * i + 3),  (2 * i + 2 * m + 2));
        for j in (2 * i + 3)..=std::cmp::min(w.len() as i32 - 1, (2 * i + 2 * m + 2)) { w[j as usize] -= 1 }
        // println!("i={:?}",i);
    }

    power_prod(&mut w, minp)
}

fn power_prod(w: &mut [i32], minp: &[usize]) -> BigUint {
    let mut result = BigUint::one();
    // println!("left={:?}   right={:?}", 2, w.len());
    for x in (2..w.len()).rev() {
        if w[x] != 0 {
            let p = minp[x];
            if p == x {
                result *= BigUint::from(p as u8).pow(w[p] as u32);
            } else {
                w[p] += w[x];
                w[x / p] += w[x];
            }
        }
        // println!("x={:?}", x);
    }
    result
}

```

### `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.3"
num-traits = "0.2"

[profile.release]
lto = "y"
panic = "abort"

```

### Build and run

```
$ cargo build --release
$ cargo run --release

    Finished release [optimized] target(s) in 0.03s
     Running `target\release\rust_hello.exe`
0       1
1       1
2       6
3       113
4       20174
5       22306955
6       123222909271
7       3283834214485890
8       421263391026827547540
9       260028731850596651411721718
10      772086476515163830856527013278243
11      11025620741283840573496993339545350520150
12      757129347300072898736973484532998417574513923224
13      249988147151001673705019407088899436686198854480028470000
14      396838668136562149340506906750722486726178589993407032745345886384
15      3028477458789741699983305448780147795991274413836547758038383842702471625120

```
> 1 votes

# Answer
# Scala

A port of @Peter Taylor's Java answer in Scala.

TIO

```
import scala.collection.mutable.ArrayBuffer
import scala.math._
import java.math.BigInteger

object OptimisedCounter extends App {
  private var minp = ArrayBuffer.fill(2)(0)

  if (args.length > 0) {
    args.map(arg => println(count(arg.toInt)))
  } else {
    (0 until 16).foreach(n => println(s"$n\t${count(n)}"))
  }

  private def count(n: Int): BigInteger = {
    if (n == 0) return BigInteger.ONE

    if (minp.length < 3 * n) {
      val wider = ArrayBuffer.fill(3 * n)(0)
      for (index <- 0 until (minp.length) by 1) wider(index) = minp(index)
            //println(minp.mkString(","))
      (minp.length until wider.length).foreach { x =>
        // Find the smallest prime which divides x
        wider(x) = (2 to x).find(x % _ == 0).getOrElse(2)
      }
      minp = wider

            //println(minp.mkString(","))
    }

    val E = countE(n)
    val R2 = countR2(n)
    val F = countF(n)
    val R3 = countR3(n)
    val R = countR(n)
    val FR = countFR(n)
    val sum = E
      .add(R3)
      .add(R2.add(R).multiply(BigInteger.valueOf(2)))
      .add(F.add(FR).multiply(BigInteger.valueOf(3)))
    sum.divide(BigInteger.valueOf(12))
  }

  private def countE(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((i + 1) to (i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countR2(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      w(3 * i + 2) += 1
      ((3 * i + 1) to (2 * i + n + 1)).foreach(j => w(j) -= 1)
      ((2 * i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countF(n: Int): BigInteger = {
    val w = ArrayBuffer.fill(3 * n)(0)
    (0 until n).foreach { i =>
      ((2 * i + 1) to (2 * i + n)).foreach(j => w(j) -= 1)
      ((i + n + 1) to (i + 2 * n)).foreach(j => w(j) += 1)
    }
    powerProd(w.toArray)
  }

  private def countR3(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    countE(n / 2).pow(2)
  }

  private def countR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * m - 1)(0)
    (0 until m).foreach { i =>
      (1 to (3 * i + 1)).foreach(j => w(j) += 2)
      (1 to (i + m)).foreach(j => w(j) -= 2)
    }
    powerProd(w.toArray)
  }

  private def countFR(n: Int): BigInteger = {
    if ((n & 1) == 1) return BigInteger.ZERO
    val m = n / 2
    val w = ArrayBuffer.fill(3 * n - 2)(0)
    (1 to m).foreach(j => w(j) -= 1)
    (2 * m until 3 * m).foreach(j => w(j) += 1)
    (0 to 2 * m - 3).foreach { i =>
      ((i + 2 * m + 1) to (i + 4 * m)).foreach(j => w(j) += 1)
      ((2 * i + 3) to (2 * i + 2 * m + 2)).foreach(j => w(j) -= 1)
    }
    powerProd(w.toArray)
  }

  private def powerProd(w: Array[Int]): BigInteger = {
    var result = BigInteger.ONE
    (2 until w.length).reverse.foreach { x =>
      if (w(x) == 0) ()
      else {
        val p = minp(x)
        if (p == x) result = result.multiply(BigInteger.valueOf(p).pow(w(p)))
        else {
          // Redistribute it. This should ensure we avoid negatives.
          w(p) += w(x)
          w(x / p) += w(x)
        }
      }
    }
    result
  }
}

```
> 0 votes

---
Tags: fastest-code, combinatorics, tiling
---


title: Shifted auto-sum

text:
# Question
Title: Shifted auto-sum
Let’s take a positive integer such as `123`. We define the shifted auto-sum of this integer as follows:

* 123 has 3 digits. We thus consider 3 copies of 123.
* We stack each copy on top of each other, shifted by 1 digit each time:

```
123
 123
  123

```
* We pad each copy with `0`s (excluding the last one) to get a standard column addition:

```
12300
 1230
  123
-----
13653

```

### Test cases

```
Input                Shifted auto-sum

1                    1
12                   132
100                  11100
123                  13653
6789                 7542579
4815162342           5350180379464981962

```

### Scoring

This is , so the shortest answer in bytes wins.

# Answer
# Japt, 3 bytes

```
ì£U

```

Try it

Found a 3 that works with numbers.

```
ì£U
ì      Convert to base-10 array of digits, do ... and convert back:
 £     Map over the array with:
  U    Return the original input

```

---

# Japt, 4 bytes

```
*sç1

```

Try it

```
*sç1
*       Multiply the input with
 s      Transform the input as a string:
  ç1    Replace each char by 1

```

If it is allowed to take the number as a string:

# Japt, 3 bytes

```
*ç1

```

Try it

```
*ç1
*      Multiply the input with (coercing both sides to number)
 ç1    Replace each char by 1

```
> 10 votes

# Answer
# Jelly, ~~4~~ 3 bytes

-1 byte thanks to Unrelated String.

```
ṁDḌ

```

Try it online!

A port of my PARI/GP answer.

```
 D     # To digits
ṁ      # Fill with the input
  Ḍ    # From digits

```
> 8 votes

# Answer
# Haskell, 23 bytes

```
f n=n*read('1'<$show n)

```

Try it online!

Multiplies `n` by the number obtained from replacing each of its digits with a `1`. Same length pointfree:

```
(*)<*>read.('1'<$).show

```

Try it online!
> 7 votes

# Answer
# , 21 bytes

```
\(x)x*10^nchar(x)%/%9
```

Attempt This Online!

Uses an '`x` multiplied by `111`' approach, where the '`111`' is constructed as the next power-of-10 greater-or-equal to `x`, integer-divided by 9.

(now reading the other answers: this is the same approach already used in xnor's Python answer)
> 7 votes

# Answer
# Python 2, 25 bytes

```
lambda n:10**len(`n`)/9*n

```

Try it online!
> 6 votes

# Answer
# Vyxal, 8 bytes

```
₌ẏL+↲›⌊∑

```

Try it Online!

(yes I know there's a 4 byter, but I liked this approach so much I posted it separately)

Kids these days with their "ten to the power of" approaches. Hasn't anyone ever heard of literal spec interpretation? :p

Takes input as a string.

## Explained

```
₌ẏL+↲›⌊∑
₌ẏL      # Push the range [0, len(in)) and len(in)
  +      # add those together
   ↲›    # pad the input left with spaces until length for each item in that list and replace all spaces with 0s
     ⌊∑  # convert each item to int and sum

```
> 5 votes

# Answer
# JavaScript (Node.js), 27 bytes

```
n=>(n+'').replace(/./g,1)*n

```

Try it online!

It is shorter than

```
n=>n*(g=i=>i<n?g(i+9):i/9)`` // 28 bytes
n=>~~(.1**~Math.log10(n)/9)*n // 29 bytes
n=>(10**-~Math.log10(n)-1)/9*n // 30 bytes

```
> 5 votes

# Answer
# Bash, 19 bytes

```
expr ${1//?/1} * $1

```

An unusual example of `expr` being shorter than regular `$[]` arithmetic expansion.

Try it online!
> 5 votes

# Answer
# Jelly, 5 bytes

```
D1€Ḍ×

```

Try it online!

I suspect there's a 4-byter but I can't find one.

```
D  Ḍ  # To digits...
 1€   # Fill with 1
    × # Multiply by original

```
> 4 votes

# Answer
# C (gcc), ~~38~~ ~~36~~ 35 bytes

```
m;f(n){n*=m=exp10(m=log10(n)+1)/9;}

```

Try it online!

*Saved ~~2~~ 3 bytes thanks to tsh!!!*
> 4 votes

# Answer
# Retina, 16 bytes

```
.+
$.($($.&*1)**

```

Try it online! Link includes test cases. Explanation: `$($.&*1)` represents the input with its characters replaced by `1`s, which is then implicitly multiplied by the original input. (Retina's multiplication is decimal times unary, so the extra `*` is required to convert one of the decimal numbers to unary for the multiplication, and the `$.(` converts the result back to decimal.)
> 4 votes

# Answer
# APL (Dyalog Unicode), 7 bytes

Anonymous tacit prefix function.

```
10⊥⊢⊣¨⍕

```

Try it online!

stringify

`⊣¨` defer each digit in favour of:

the argument

`10⊥` evaluate as (carrying) base-10 digits
> 4 votes

# Answer
# BQN, 15 bytes

```
10⊸×⊸+˜´⊢⊣¨•Fmt

```

Try it at BQN online!

### Explanation

```
10⊸×⊸+˜´⊢⊣¨•Fmt
            •Fmt  Convert the argument number to a string (list of characters)
          ⊣¨      Replace each character with
        ⊢         The argument number
       ´          Right fold
      ˜           on this function with reversed arguments:
    ⊸+              Add the right argument to
10⊸×                The left argument times 10

```

### Other versions

Another 15-byte solution that constructs a list of powers of 10:

```
+´⊢×10⋆↕∘≠∘•Fmt

```

A 16-byte solution that constructs a list of 10s and multiplication-scans it:

```
+´÷⟜10×`·10¨•Fmt

```

Another 16-byte solution that constructs a string of 1s and evaluates it:

```
⊢×·•BQN·'1'¨•Fmt

```
> 4 votes

# Answer
# Desmos, 29 bytes

```
f(k)=∑_{n=0}^{logk-.5}10^nk

```

Try It On Desmos!

Using the formula $$n\cdot\frac{10^{length(n)}-1}9$$ is surprisingly 1 byte longer:

```
f(k)=k(10^{floor(logk)+1}-1)/9

```
> 4 votes

# Answer
# K (ngn/k), ~~17~~ 14 bytes

```
{x*10/(#$x)#1}

```

Try it online!

*-3 bytes thanks to ovs!*

Explanations:

```
{x*10/(#$x)#1}  Main function. x is input (123)
            1   The number 1              (1)
           #    Duplicate by the amount of
      ( $x)     x converted to string     ("123")
       #        Length times              (3)
                (This creates an array of integers, which is (1, 1, 1))
   10/          Convert the integers to base 10
 x*             Multiply by x

```
> 4 votes

# Answer
# 05AB1E, ~~5~~ 4 bytes

```
$g×*

```

-1 byte porting emanresuA's Jelly answer

Try it online or verify all test cases.

Alternative 4-byter (port of *@Bubbler*'s Japt answer:

```
gиTβ

```

Try it online or verify all test cases.

**Explanation:**

```
$     # Push 1 and the input
 g    # Pop the input, and push its length
  ×   # Repeat 1 that many times
   *  # Multiply it to the (implicit) input-integer
      # (after which the result is output implicitly)

g     # Push the length of the (implicit) input-integer
 и    # Repeat the (implicit) input-list that many times
  Tβ  # Convert it from a base-10 list to a base-10 integer
      # (after which the result is output implicitly)

```
> 3 votes

# Answer
# Python, 57 bytes

```
lambda n:sum(n//10for i in range(len(str(n)))if(n:=n*10))

```

Attempt This Online!

# Python, 31 bytes

```
lambda n:int('1'*len(str(n)))*n

```

Attempt This Online!

# Python, 29 bytes

Port of `@xnor's` answer.

```
lambda n:10**len(str(n))//9*n

```

Attempt This Online!
> 3 votes

# Answer
# Charcoal, 8 bytes

```
Ｉ×ＮＩ×1Ｌθ

```

Try it online! Link is to verbose version of code. Works with both integer and string input. Explanation:

```
     1      Literal string `1`
    ×       Repeated by
       θ    First input
      Ｌ     Digit length
   Ｉ        Cast to integer
 ×          Multiplied by
  Ｎ         Input as an integer
Ｉ           Cast to string
            Implicitly print

```

7 bytes but only accepts numeric input in JSON format:

```
Ｉ↨χＥＩθθ

```

Try it online! Link is to verbose version of code. Explanation: Port of @alephalpha's Jelly answer.

```
     θ  Input
    Ｉ   Cast to string
   Ｅ    Map over characters
      θ Input
 ↨χ     Convert from "base 10"
Ｉ       Cast to string
        Implicitly print

```
> 3 votes

# Answer
# Nibbles, 4 bytes (8 nibbles)

```
`@~.`p$@

```

Inspired by alephalpha's Jelly answer.

```
    `p$     # convert input to string
   .        # map function over each element (characters in string):
       @    # second argument (arg1=element, arg2=input)
            # so now we have digits-of-input copies of the input
`@          # interpret as digits in base
  ~         # 10

```
> 3 votes

# Answer
# Burlesque, 8 bytes

```
Jq./]m|*

```

Try it online!

Takes input as a string (as from standard in) containing number.

```
J   # Duplicate input
q./ # Quoted isNumeric (returns 1 for each digit)
]m  # Map and concat to string
|*  # String multiply

```
> 3 votes

# Answer
# Vyxal `r`, 4 bytes

```
Lẋ₀β

```

Try it Online!

Another approach.

## Vyxal `s`, 3 bytes

```
ẏ↵*

```

Try it Online!

To get even with these 3 byters, lol.
> 3 votes

# Answer
# Pip, 7 bytes

```
a*:1X#a

```

Try it online!

### Explanation

Port of Bubbler's string-based Japt answer:

```
a*:1X#a
     #a  Length (number of digits) of cmdline argument
   1X    String of that many 1s
 *:      (Treat as a number and) Multiply by
a        Cmdline argument

```
> 3 votes

# Answer
# dc, 10

```
?dZAr^9/*p

```

Try it online!

### Explanation

```
?             # Read input
 d            # duplicate input on stack
  Z           # push count digits of input to stack
   Ar         # push 10, reverse top 2 stack elements
     ^        # calculate power of 10
      9/      # divide by 9 to get 111....1
        *     # multiply by original input
         p    # print

```
> 3 votes

# Answer
# Brachylog, 5 bytes

```
lᵐ×↙?

```

Try it online!

```
 ᵐ       For each digit of the input,
l        get its length (i.e. 1),
 ᵐ       and concatenate the lengths back into an integer.
  ×↙?    Multiply by the input.

```
> 3 votes

# Answer
# Wolfram Language (Mathematica), ~~27~~ 26 bytes

```
#(10^IntegerLength@#-1)/9&

```

Try it online!
> 3 votes

# Answer
# Perl 5 `-Minteger -pF` , 12 bytes

```
$_*=10**@F/9

```

Try it online!
> 3 votes

# Answer
# Regex  (PCRE2 v10.35+), 36 bytes

```
^(?*(\1{10}|^x{9})*x(x*$))(?*\2x+)x+

```

Attempt This Online!
Attempt This Online! \- just the test cases (that are small enough to complete in reasonable time)

Takes its input in unary, as the length of a string of `x`s. Returns its output as the number of ways the regex can match. (The rabbit emoji indicates this output method. It can yield outputs bigger than the input, and is really good at multiplying.)

```
^                     # Anchor to start; tail = N = input number
(?*                   # Non-atomic lookahead - try all possibilities to find ones
                      # resulting in a later match. This is used to multiply the
                      # number of possible matches outside by the number of
                      # possible matches inside.
    (\1{10}|^x{9})*x  # tail -= {any power of 10}
    (x*$)             # \2 = tail = N - {the power of 10}
)
(?*                   # Another non-atomic lookahead
    \2                # tail = {the power of 10 chosen above}
    x+                # Multiply the number of possible matches by tail. This
                      # applies for each power of 10 found, and those numbers of
                      # matches will add together.
)
x+                   # Add tail possible matches, where tail==N, with the above
                     # multiplier being applied.

```

# Regex  (Perl / PCRE), ~~108~~ 76 bytes

```
^(x?)(?=(x*)\2$).*(?=\2((x+)\4{8}(?=\4$))*x{5}$)()?x+.*(?=\2$)(()?x+|\1+$)|x

```

Try it online! \- Perl v5.28.2 / Attempt This Online! \- Perl v5.36+
Try it online! \- PCRE1
Try it online! \- PCRE2 v10.33 / Attempt This Online! \- PCRE2 v10.40+
Try it online! \- PCRE2 - just the test cases

Without `(?*`...`)` molecular lookahead, it's much harder to do multiplication. So we reserve half of the number to use as a counter for one operand of the multiplication, letting \\$a=\lfloor n/2\rfloor\\$ and \\$b=n\bmod 2\\$, such that \\$2a+b=n\\$, in order to emulate operations on \\$n\\$.

We don't actually need to construct the number of the form 111...111 to multiply by. We can pretty much directly sum \\$n\\$ multiplied by each power of \\$10\\$, that is to say,

$$\large\left({\sum\_{k=1}^{\lfloor log\_{10}n\rfloor}2{{10^k}\over 2}(2a+b)}\right)+n = \sum\_{k=0}^{\lfloor log\_{10}n\rfloor}10^k n$$

using emulated arithmetic on \\$2a+b=n\\$.

```
    ^                         # Anchor to start; tail = N = input number
    (x?)                      # \1 = N % 2; tail -= \1
    (?=(x*)\2$)               # \2 = floor(N / 2)
    .*                        # Find any value of tail ≤ N satisfying the
                              # following:
    (?=
        # Assert that 2*(tail-\2) is a power of 10
        \2                      # tail -= \2
        (
            (x+)\4{8}(?=\4$)  # Assert tail is divisible by 10; tail /= 10
        )*                    # Iterate the above as many times as possible,
                              # minimum zero.
        x{5}$                 # Assert tail == 5
    )
    ()?                       # possibleMatches *= 2, for the following:
    x+                        # possibleMatches *= tail - \2, for the following:
    .*(?=\2$)                 # tail = \2
    (
        ()?                   # possibleMatches *= 2, for the following:
        x+                    # possibleMatches += tail
    |
        \1+$                  # possibleMatches += \1
    )
|
                              # Unanchored
    x                         # possibleMatches += N

```

# Regex (.NET), ~~129~~ 117 bytes

`^(?=(x)*)(x?)(?=(x*)\3$)(((?=\3((x+)\7{8}(?=\7$))*x{5}$)(?=((?=.*(?=\3)(?<1>x)*){4}(?=.*(?=\2)(?<1>x)?){2}x)*\3))?x)*`

Try it online!
Try it online! \- just the test cases

Returns its output as the capture count of `\1`.

This is now a fairly straight port of the Perl/PCRE  version.

The former 129 byte version actually constructed the number whose digits are all `1` to multiply by – indirectly – by constructing the number \\$c\\$ whose every digit is a `1` and has one less digit than the input, such that \\$10c+1\\$ is the actually desired number, then returning \\$(10c)(2a+b)+n\\$ \\$=(10c+1)n\\$ as a capture count. That version couldn't return a result for an input of \\$0\\$, but this version can.

```
^                               # Anchor to start; tail = N = input number
(?=(x)*)                        # \1.captureCount = tail = N
(x?)                            # \2 = N % 2; tail = tail-\2 == \3 * 2
(?=(x*)\3$)                     # \3 = floor(N / 2)
(
    (
        (?=                       # Lookahead conditional
            # Assert that 2*(tail-\3) is a power of 10
            \3                      # tail -= \3
            (
                (x+)\7{8}(?=\7$)    # Assert tail is divisible by 10; tail /= 10
            )*                      # Iterate the above as many times as possible,
                                    # minimum zero.
            x{5}$                   # Assert tail == 5
        )
        # If the above assertion matches, do the following:
        (?=
            (
                (?=
                    .*(?=\3)    # tail = \3
                    (?<1>x)*    # \1.captureCount += tail
                ){4}            # Iterate the above 4 times
                (?=
                    .*(?=\2)    # tail = \2
                    (?<1>x)?    # \1.captureCount += tail
                ){2}            # Iterate the above 2 times
                x
            )*                  # Iterate the above as many times as possible,
                                # which will be a number of iterations equal to
                                # exactly half the current power of 10.
            \3                  # Assert tail >= \3
        )
    )?
    x
)*                              # Iterate the above as many times as possible,
                                # applying the payload to each power of 10

```

I think I may have discovered a bug in the .NET regex engine, or at least the version on TIO, while working on this. The following variant of the same length should have worked, but didn't:

`^(?=(x)*)(x?)(?=(x*)\3$)((?(?=\3((x+)\6{8}(?=\6$))*x{5}$)(?=((?=.*(?=\3)(?<1>x)*){4}(?=.*(?=\2)(?<1>x)?){2}x)*\3))x)*`
> 3 votes

# Answer
# perl -p, 12 bytes

```
$_*=s/./1/gr

```

Multiplies itself with `1`, `11`, `111`, etc, depending on the length of the input.

Try it online!
> 3 votes

# Answer
# PHP, 29 bytes

```
fn($n)=>$n*~-10**strlen($n)/9

```

Try it online!
> 3 votes

# Answer
# Knight, 12 bytes

```
O*/^10L=xP9x

```

Try it online!

Ungolfed:

```
OUTPUT (* (/ (^ 10 (LENGTH (= x PROMPT))) 9) x)

```
> 3 votes

# Answer
# Ruby, 40 36 bytes

```
p [*0...gets.size].sum{eval$_+?0*_1}
```

Attempt This Online!

-4 bytes thanks to Steffan
> 3 votes

# Answer
# \>\<\>, 25 bytes

```
::a(?va,{:@a*}!
=?n+0>~l1

```

Try it online

**Explanation**

```
:             !   # save a copy of the input on the stack
 :a(?v            # if the current number is lower than 10 move to next row
      a,          # else divide the current number by 10
        {:@a*}    # and store a copy of the original number * 10 on the stack

     >~           # discard the current number
=?n    l1         # if there is only one item on the stack, print it
   +0             # else add the top 2 numbers, and push a 0 to be discarded

```
> 3 votes

# Answer
# , ~~12~~ 11 bytes

```
10#.]#~$@":

```

Try it online!

-1 byte thanks to rdm
> 3 votes

# Answer
# GolfScript, 9 bytes

```
.,1`*~\~*

```

Try it online!

There's clearly a better way to do this, only I haven't discovered it.
> 3 votes

# Answer
# Vyxal, 4 bytes

```
ẏ↵*∑

```

Try it Online!

```
ẏ    # Exclusive range 0-length
 ↵   # 10 to the power of each of those
  *  # Multiply by the input
   ∑ # Sum

```
> 2 votes

# Answer
# PARI/GP, 21 bytes

```
n->n*(10^#Str(n)-1)/9
```

Attempt This Online!

---

# PARI/GP, 29 bytes

```
f(n,m=n)=if(n,10*f(n\10,m)+m)
```

Attempt This Online!

---

# PARI/GP, 31 bytes

```
n->fromdigits([n|i<-digits(n)])
```

Attempt This Online!
> 2 votes

# Answer
# MathGolf, 5 bytes

```
hª*y*

```

Try it online.

**Explanation:**

```
h      # Push the length of the (implicit) input-integer (without popping)
 ª*    # Push a list of the length amount of 1s:
 ª     #  Push [1]
  *    #  Python-style multiply it to the length
   y   # Join it together and implicitly convert it to an integer
    *  # Multiply it to the input
       # (after which the entire stack is output implicitly as result)

```
> 2 votes

# Answer
# Factor, 33 bytes

```
[ dup log10 1 /i 1 + 10^ 9 /i * ]

```

Try it online!

```
       ! 123
dup    ! 123 123
log10  ! 123 2.089905111439398
1      ! 123 2.089905111439398 1
/i     ! 123 2
1      ! 123 2 1
+      ! 123 3
10^    ! 123 1000
9      ! 123 1000 9
/i     ! 123 111
*      ! 13653

```
> 2 votes

# Answer
# BQN, 27 bytes

```
{≤÷10?(×)+×10;×}⟜1

```

Try it at BQN REPL

Recursive function to try to circumvent the lack of simple base-conversion, string-conversion or logarithm functions in BQN.

```
{≤÷10?(×)+×10;×}⟜1     # recursive function with input
                         ⟜1     #  and second argument  initiall set to 1
 ≤÷10?                        # is  less than or equal to ÷10?
                                # if yes:
        (×)                   #  return ×
             +                  #  plus
              ×10            #  the result of a recursive call
                                #   with  multiplied by 10
                    ;           # otherwise:
                     ×        #  just return ×

```
> 2 votes

# Answer
# Raku, 13 bytes

```
{$_*S:g/./1/}

```

Try it online!
> 2 votes

# Answer
# C++ (gcc), 42 bytes

```
[](int&n){int i=1;for(;n/i;)i*=10;n*=i/9;}

```

Try it online!
> 2 votes

# Answer
# Prolog (SWI), 40 bytes

```
X+N:-N is(10^(floor(log10(X))+1)-1)*X/9.

```

Try it online!
> 2 votes

# Answer
# ECMAScript 2016, 73 bytes

```
n=>Array(String(n).length).fill(n).map((a,b)=>a*10**b).reduce((a,b)=>a+b)

```
> 1 votes

# Answer
# PowerShell Core, 46 bytes

```
"$(1..($s="$args").Length|%{"+$s";$s+=0})"|iex

```

Try it online!

For example for `123`:
We generate an array containing: `"+123","+1230""+12300"`
Then we join them all in a space separated string `"+123 +1230 +12300"`
Then we execute it as a PowerShell expression, returning `13653`
> 1 votes

# Answer
# Thunno 2 `S`, 3 bytes

```
Że×

```

Try it online!

#### Explanation

```
Że×  # Implicit input
Ż    # Push [0..len(input))
 e   # Take 10 ** each
  ×  # Multiply by input
     # Auto-sum and output

```
> 1 votes

---
Tags: code-golf, number, integer
---