# Record Least Uncommon Multiple Counts

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).

• 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.
– xnor
Mar 27, 2023 at 0:00
• $\operatorname{lucm}(n)$ is sum(2*d for d in divisors(n) if d * d < n and gcd(d, n // d) == 1). Mar 27, 2023 at 3:02
• @alephalpha with a special case for n==1
– G B
Mar 27, 2023 at 6:34
• @alephalpha Consistent with that all primorial indices appear to generate new maxima. Mar 27, 2023 at 7:19
• 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. Apr 11, 2023 at 15:18

# 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, $$\\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, except the special case $$\n=1\$$.

• I've got lucm[20]=10. Is it right? Mar 28, 2023 at 5:24
• @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$. Mar 28, 2023 at 5:42
• Aha I mistook values of $n$ with $H$ in the question,thanks Mar 28, 2023 at 5:48
• 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) Mar 28, 2023 at 6:06
• @lesobrod I use the timeout command in linux: timeout 1m python a.py. Mar 28, 2023 at 6:11

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

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.

# Mathematica 13.2112/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

## 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])


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

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")

• 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 Apr 16, 2023 at 7:36