Record Least Uncommon Multiple Counts

Question

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

• xnor: 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.
• alephalpha: \$\operatorname{lucm}(n)\$ is sum(2*d for d in divisors(n) if d * d < n and gcd(d, n // d) == 1).
• G B: @alephalpha with a special case for n==1
• loopy walt: @alephalpha Consistent with that all primorial indices appear to generate new maxima.
• loopy walt: 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.

Answer 1

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\$.

• lesobrod: I've got lucm[20]=10. Is it right?
• alephalpha: @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\$.
• lesobrod: Aha I mistook values of \$n\$ with \$H\$ in the question,thanks
• lesobrod: 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)
• alephalpha: @lesobrod I use the timeout command in linux: timeout 1m python a.py.

Answer 2

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.

Answer 3

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

Answer 4

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 around (削除ここまで)as 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])

Answer 5

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

• lesobrod: 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

Answer 6

Rust, \$ 153/50 \approx 3.06 \$ on my laptop

use std::time::{Duration, Instant};
/// Computes the greatest common divisor of two numbers using the Euclidean algorithm.
fn gcd(mut a: u64, mut b: u64) -> u64 {
 while b != 0 {
 let temp = b;
 b = a % b;
 a = temp;
 }
 a
}
/// Returns a vector of all unitary divisors of `n`.
fn udivisors(n: u64) -> Vec<u64> {
 let mut divisors = Vec::new();
 let sqrt_n = (n as f64).sqrt() as u64;
 for d in 1..=sqrt_n {
 if n % d == 0 {
 let e = n / d;
 if gcd(d, e) == 1 {
 divisors.push(d);
 if e != d {
 divisors.push(e);
 }
 }
 }
 }
 divisors
}
/// Computes the sum of `2 * d` for each unitary divisor `d` of `n` where `d * d < n`.
fn lucm(n: u64) -> u64 {
 udivisors(n)
 .into_iter()
 .filter(|&d| d * d < n)
 .map(|d| 2 * d)
 .sum()
}
fn main() {
 let time_limit = Duration::from_secs(60); // 60-second time limit
 let start_time = Instant::now();
 let mut counter: u64 = 1;
 let mut prev_prog_max_val: u64 = 1;
 // Print the initial state
 println!("{} {} {}", counter, 1, 1);
 let mut ctr: u64 = 2;
 // Infinite loop with time limit
 loop {
 // Check if the time limit has been reached
 if start_time.elapsed() >= time_limit {
 println!("Time limit of 60 seconds reached. Terminating program.");
 break;
 }
 let count = lucm(ctr);
 if count > prev_prog_max_val {
 counter += 1;
 println!("{} {} {}", counter, ctr, count);
 prev_prog_max_val = count;
 }
 ctr += 1;
 }
}