10
\$\begingroup\$

The Greatest Common Divisor, or gcd, of two positive integers \$x\$ and \$y\$ is the largest positive integer that divides both \$x\$ and \$y\$.

The Least Common Multiple, or lcm, of two positive integers \$x\$ and \$y\$ is the smallest positive integer that is a multiple of \$x\$ and \$y\$.

We define the Least Uncommon Multiple (lum) of \$x\$ and \$y\$ to be \$\frac{\text{lcm}(x,y)}{\text{gcd}(x,y)}\$.

We define the Least Uncommon Multiple Count (lucm) of \$n\$ as the number of positive integer pairs \$x,y\$ such that \$x,y\leq n\$ and \$lum(x,y)=n\$.

For example, \$lucm(2)=2\$, as \$lum(1,2)=lum(2,1)=2\$, while \$lum(1,1)=lum(2,2)=1\$.

Another example:

\$lum(x,y)\$ 1 2 3 4 5 6
1 1 2 3 4 5 6
2 2 1 6 2 10 3
3 3 6 1 12 15 2
4 4 2 12 1 20 6
5 5 10 15 20 1 30
6 6 3 2 6 30 1

\$lucm\$(6) is the number of times 6 appears in the above table, which is 6. Thus \$lucm(6)=6\$.

We define the sequence \$H\$ as follows: \$h_i=lucm(i)\$ for all positive integers \$i\$. \$H\$ begins as follows: 1, 2, 2, 2, 2, 6, 2, 2, 2, 6...

I was playing around with \$H\$, and what caught my interest was the set of numbers for which \$lucm(n)\$ was greater than that of any of the smaller integers. For example, \$lucm(6)=6\$, which is greater than any of the previous values in the sequence. We will call these numbers the Progressive Maxima of \$H\$. The first few progressive maxima of \$H\$ occur at \$n=1,\ 2,\ 6,\ 12,\ 20,\ 30\ldots\$

The challenge, should you choose to accept it, is to write code that calculates the most Progressive Maxima of \$H\$ within 1 minute without any hard-coding of the sequence.

Scoring

Using the following Python code as the reference code:

from math import gcd
counter=1
prev_prog_max=1
prev_prog_max_val=1
counts={1:1}
print(1,1,1)
ctr=2
while True:
    for j in range(1,ctr):
        q=ctr*j/(gcd(ctr,j)**2)
        if q>=ctr and q in counts.keys():
            counts[q]+=2
        else:
            counts[q]=2
    if counts[ctr]>prev_prog_max_val:
        counter+=1
        print(counter,ctr,counts[ctr])
        prev_prog_max=ctr
        prev_prog_max_val=counts[ctr]
    ctr+=1

Take the number of progressive maxima your code returns in one minute, and divide it by the number of progressive maxima the reference code returns in one minute on your computer. The resulting value is your score. Highest score wins.

For example, say my code returned 140 progressive maxima in 1 minute, and the reference code returned 45. My score would be \$\frac{140}{45}=\frac{28}{9}=3.\overline{09}\approx3.1\$.

If there is a tie, the earliest answer wins.


My current best, which I will post as a non-competing answer in a month after the last answer for this question, returned 389 progressive maxima in 1 minute, while the reference code returned 45. My score would thus be \$\sim\!\!8.64\$. I'm completely confident that I've gotten all of the elements of \$H\$ that are less than 1 billion, and pretty confident that I have all such elements under 1 quadrillion (though it took me about 3-4 days of runtime to get them).

\$\endgroup\$
8
  • 1
    \$\begingroup\$ I was wondering if there might be a nice number-theory characterization of lucm, but the restriction to counting pairs with \$x,y \leq n\$ makes this tricky. Each number would appear on the lum table infinitely many times if it was extended indefinitely. \$\endgroup\$
    – xnor
    Mar 27, 2023 at 0:00
  • 10
    \$\begingroup\$ \$\operatorname{lucm}(n)\$ is sum(2*d for d in divisors(n) if d * d < n and gcd(d, n // d) == 1). \$\endgroup\$
    – alephalpha
    Mar 27, 2023 at 3:02
  • 2
    \$\begingroup\$ @alephalpha with a special case for n==1 \$\endgroup\$
    – G B
    Mar 27, 2023 at 6:34
  • 2
    \$\begingroup\$ @alephalpha Consistent with that all primorial indices appear to generate new maxima. \$\endgroup\$
    – loopy walt
    Mar 27, 2023 at 7:19
  • 1
    \$\begingroup\$ 83. Did you also observe that the number of prime factors without multiplicity never goes down? That is the main (unproven) conjecture my algo is based on. \$\endgroup\$
    – loopy walt
    Apr 11, 2023 at 15:18

5 Answers 5

6
\$\begingroup\$

Python + SymPy, \$119/47\approx 2.532\$ on my computer

from sympy.ntheory.factor_ import udivisors

def lucm(n):
    return sum(2*d for d in udivisors(n) if d * d < n)

counter = 1
prev_prog_max = 1
prev_prog_max_val = 1
print(1, 1, 1)
ctr = 2

while True:
    count = lucm(ctr)
    if count > prev_prog_max_val:
        counter += 1
        print(counter, ctr, count)
        prev_prog_max = ctr
        prev_prog_max_val = count
    ctr += 1

The udivisors function returns a list of unitary divisors of a number. A unitary divisor of \$n\$ is a divisor \$d\$ of \$n\$ such that \$d\$ and \$n/d\$ are coprime.

In fact, for two positive integers \$a,b\$, let \$g=\gcd(a,b)\$, then we can write \$a=gc,b=gd\$ for some coprime \$c,d\$, and we have \$\operatorname{lcm}(a,b)=gcd\$, \$\operatorname{lum}(a,b)=cd\$.

Now assume that \$0<a,b\le n\$, \$\operatorname{lum}(a,b)=n\$. Then \$c,d\$ must be unitary divisors of \$n\$, and \$g\le n/\max(c,d)\$. So, for each unitary divisor \$d\$ of \$n\$ (and \$c=n/d\$), there are \$n/\max(c,d)=\min(c,d)\$ corresponding \$(a,b)\$'s such that \$\operatorname{lum}(a,b)=n\$.

So the total number of such \$(a,b)\$'s is \$\sum_{d\ \mid\ n\atop\gcd(d,n/d)=1}\min(d,n/d)\$. This equals to \$\sum_{d\ \mid\ n,\ d^2<n\atop\gcd(d,n/d)=1}2d\$, except the special case \$n=1\$.

\$\endgroup\$
5
  • \$\begingroup\$ I've got lucm[20]=10. Is it right? \$\endgroup\$
    – lesobrod
    Mar 28, 2023 at 5:24
  • \$\begingroup\$ @lesobrod That is right. \$\operatorname{lum}(1,20)=\operatorname{lum}(4,5)=\operatorname{lum}(5,4)=\operatorname{lum}(8,10)=\operatorname{lum}(10,8)=\operatorname{lum}(12,15)=\operatorname{lum}(15,12)=\operatorname{lum}(16,20)=\operatorname{lum}(20,1)=\operatorname{lum}(20,16)=20\$. \$\endgroup\$
    – alephalpha
    Mar 28, 2023 at 5:42
  • \$\begingroup\$ Aha I mistook values of \$n\$ with \$H\$ in the question,thanks \$\endgroup\$
    – lesobrod
    Mar 28, 2023 at 5:48
  • \$\begingroup\$ Sorry for another question. I rarely use Python, so is it possible configure some "Time constrained" option for Python (3.8) (may be with PyCharm) \$\endgroup\$
    – lesobrod
    Mar 28, 2023 at 6:06
  • \$\begingroup\$ @lesobrod I use the timeout command in linux: timeout 1m python a.py. \$\endgroup\$
    – alephalpha
    Mar 28, 2023 at 6:11
5
\$\begingroup\$

J, \$178/48 \approx 3.708 \$

Based on alephalpha's comment, powered by J's fast prime factorisation.

LUCM =: {{+/(<.|.)*/y^|:#:i.2^#y}}
main =: {{
'N C M' =. 1
echo 1 1 1
while. 1 do.
  N =. N+1
  L =. LUCM ^/__ q: N
  if. L>M do.
    C =. C+1
    M =. L
    echo C,x:N,M
  end.
end.
}}
main ''

Attempt This Online!

The function LUCM takes the prime factorisation of a number as the prime powers, e.g. 120 as 8 3 5 (2^3 3^1 5^1), and computes the LUCM from that.

\$\endgroup\$
4
\$\begingroup\$

Mathematica 13.2 112/48 113/48 \$\approx 2.35\$

PyCharm 2020 + Python 3.8 was used for run Python code.
Thanks to @alephalpha for the cool math idea, please vote for his answer!


Clear[lucm]; lucm = Compile[{{n, _Integer}}, If[n == 1, 1, 
2 Total@Select[Divisors[n], #^2 < n && GCD[#, Quotient[n, #]] == 1 &]]];

TimeConstrained[n = 1; count = 1; max = 0;
 While[True,
  res = lucm[n];
  If[res > max, max = res; Print[count, " ", n, " ", res]; count += 1];
  n += 1
  ], 60]; count

Some tips for Mathematica speed optimization:

  • Switch on "Launch parallel kernels" → "At startup"
  • Use Compile[] with type declaration
  • Use real numbers for arithmetics to avoid long exact integer calculation,
    but integers for such built-ins as GCD[], Divisors[] etc.
  • In old posts you can find a statement that old-school Do[],If[] etc work faster than such built-ins as Select[], Reap[]/Sow[] etc.
    That is not the case at all!
  • Run your program several times (maybe with small values) before the final trial.
    Mathematica uses cache ;)
  • Unfortunately there is no recurrent function for Divisors[], but it is known that Mathematica effectively uses the results of previous calculations
\$\endgroup\$
4
\$\begingroup\$

Pythran 206/46 ~ 4.5 on my laptop

Please note that this approach is on my laptop limited by available memory; the 206 were computed in well under 1 minute.

@alephalpha's method with minor tweaks:

  • keep lists of primes and decompositions found so far aroundas we are spending the memory anyway a simple sieve is more efficient
  • use the rather crude bound lumc(n) <= 2^#distinct prime factors x sqrt n for some branch cutting
  • use Pythran acceleration (this is he most impactful)

Code:

accel.py (this must be Pythran compiled)

from numpy import prod

# pythran export resdivsum(int32 [:],int,int,int,int)
# pythran export decomp(int32 [:,:])
# pythran export master(int32 [:,:],int,int)

def resdivsum(pp,mxdef,thmx,P,Q):
    if mxdef<0:
        return 0,mxdef
    ppp = prod(pp)
    if ppp*P < Q:
        ref = 2**len(pp) * thmx
        out = prod(pp+1)*P
        return out,mxdef - ref + out
    elif ppp*Q > P:
        out = 0
        inc,mxdef = resdivsum(pp[:-1],mxdef,thmx,P*pp[-1],Q)
        if mxdef < 0:
            return out+inc,mxdef
        out += inc
        inc,mxdef = resdivsum(pp[:-1],mxdef,thmx,P,Q*pp[-1])
        return out+inc,mxdef
    return 0,mxdef

def decomp(dec):
    m,n = dec.shape
    for i in range(2,m):
        if dec[i,-1] <= 0:
            ib = -dec[i,-1] or i
            ii = -ib
            dec[i,-1] = 0
            for j in range(i,m,i):
                ii += 1
                if ii == 0:
                    dec[j,-1] = -ib
                elif ii == ib:
                    ii = 0
                else:
                    dec[j,dec[j,-1]] = i
                    dec[j,-1] += 1

def master(dec,j,mxv):
    m,n = dec.shape
    lumc = 0
    while lumc <= mxv and j<m:
        j+=1
        d = dec[j]
        s = d[-1]
        d = d[:s]
        thmx = int(j**0.5)
        mxdef = (thmx<<s-1) - mxv
        lumc,defct = resdivsum(d,mxdef,thmx,1,1)
    return lumc,j

toplevel script (Python 3)

import time
from numpy import zeros,int32
from accel import resdivsum,decomp,master

# total time
tt = 60

r = time.perf_counter()

decompositions = zeros((100000000,10),int32)
decomp(decompositions)
print(time.perf_counter()-r)
maxidx = [0,1]
maxval = [0,0]
n = 1
rm = 1
print(1,1,1)
while time.perf_counter() < r+tt:
    lumc,n=master(decompositions,n,maxval[rm])
    if n==len(decompositions):
        break
    rm += 1
    maxidx.append(n)
    maxval.append(lumc)
    print(rm,maxidx[rm],2*maxval[rm])
\$\endgroup\$
2
\$\begingroup\$

SageMath, \$113/35 \approx 3.23 \$ on SageMathCell

Based on @alephalpha's answer, made a little change.

Thanks to @alephalpha for the cool math idea, please vote for his answer!


Competition code, run it online!

from sage.arith.misc import divisors
import timeit
import math

def lucm(n):
    return sum(2*d for d in divisors(n) if gcd(n // d, d) == 1 and d * d < n)

counter = 1
prev_prog_max = 1
prev_prog_max_val = 1
print(1, 1, 1)
ctr = 2

start_time = timeit.default_timer()

while True:
    elapsed_time = timeit.default_timer() - start_time
    if elapsed_time > 60:
        break

    count = lucm(ctr)
    if count > prev_prog_max_val:
        counter += 1
        print(counter, ctr, count)
        prev_prog_max = ctr
        prev_prog_max_val = count
    ctr += 1

print("end")

Reference code, run it online!

from math import gcd
import timeit

counter=1
prev_prog_max=1
prev_prog_max_val=1
counts={1:1}
print(1,1,1)
ctr=2

start_time = timeit.default_timer()

while True:
    elapsed_time = timeit.default_timer() - start_time
    if elapsed_time > 60:
        break
    for j in range(1,ctr):
        q=ctr*j/(gcd(ctr,j)**2)
        if q>=ctr and q in counts.keys():
            counts[q]+=2
        else:
            counts[q]=2
    if counts[ctr]>prev_prog_max_val:
        counter+=1
        print(counter,ctr,counts[ctr])
        prev_prog_max=ctr
        prev_prog_max_val=counts[ctr]
    ctr+=1
print("end")
\$\endgroup\$
1
  • \$\begingroup\$ Interesting observation! Are you run both SageMath and Python online? So SageMath online has same speed as Mathematica on PC, and Python online slower than with PyCharm \$\endgroup\$
    – lesobrod
    Apr 16, 2023 at 7:36

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.