Counting Integers and their Multiplicities
 

Let [I]S[/I][sub][I]n[/I][/sub] = {1 ≤ [I]i[/I] ≤ [I]n[/I]} be the set of positive integers not exceeding [I]n[/I].
Let
[tex]M_{m,n}=\bigcup_{1\leq i \leq m}i S_n[/tex]
be the union of all integers in S[sub][I]n[/I][/sub] times 1, 2, ..., [I]m[/I].

Is there a simple expression for |[I]M[/I][sub][I]m[/I],[I]n[/I][/sub]|?
Or a simple expression for [I]r[/I]([I]m[/I]) = lim[sub][I]n[/I]→∞[/sub] |[I]M[/I][sub][I]m[/I],[I]n[/I][/sub]|/n ? For [I]m[/I]=1, 2, 3, 4, 5, 6, 7 they are [I]r[/I]([I]m[/I]) = 1, 3/2, 2, 29/12, 44/15, 33/10, 134/35 respectively. What about larger [I]m[/I]? Does the ratio [I]r[/I]([I]m[/I])/[I]m[/I] converge as [I]m[/I] increases?

Alex

Did you try inputting the first few terms into

[url]http://www.research.att.com/~njas/sequences/[/url]

If there's an easy formula, this site probably has it.

good luck!

[quote=masser;101063]Did you try inputting the first few terms into

[URL]http://www.research.att.com/~njas/sequences/[/URL]
[/quote]... after transforming all terms into integers, that is ...

I multiplied by 420 and entered:

420,630,840,1015,1232,1386,1608

Got

"I am sorry, but the terms do not match anything in the table."

Does r(m)/m converge? You can show that r(m)/m is bounded by 0 and 1; it also looks like r(m)/m is decreasing - if you can prove this; the sequence converges. Is this really a homework question? It's a nifty question. :geek:

No, it's not a homework question. I stumbled upon it during some coding work gor GMP-ECM. At some point, we'll have to scale some residues by a^1, a^2, ..., a^n, a^2, a^4, a^6, ..., etc, until a^m, a^(2m), ..., a^(mn). I had wondered how much computation I could save by avoiding recomputing of powers that were encountered before. When I asked the guys here at Loria, we got the first couple of terms very quickly, but a closed form formula for anything proved much harder to get than I had anticipated - I thought this sort of problem probably was classical with a famous solution.

If anyone has any ideas, or pointers to the literature, please let me know!

Also, this isn't a homework exercise (except maybe for the most diabolical of teachers!) so I'll move it back to Math.

Alex

masser: r(m)/m is clearly bounded by 0 and 1, however, it is not monotonically decreasing. Emmanuel Thomé wrote Magma code that estimates the ratio reasonably quickly for 1 ≤ m ≤ 200. Graph is attached.

Alex

Should have guessed it would not be that easy...

Is that limiting value e^(-1)?

Looks like 1/log(log n) to me. Can you make the data of Thomé available, Alex? Have you tried any curve-fitting to it?

[QUOTE=philmoore;101188]Looks like 1/log(log n) to me. Can you make the data of Thomé available, Alex? Have you tried any curve-fitting to it?[/QUOTE]

It is attached.

I haven't worked very hard on this problem yet as it is not strictly necessary for writing GMP-ECM, which is what keeps me busy atm. However, I think it's an interesting problem as it looks so innocent and cute at first glance, but when you try to get a grip on it, it blows up right in your face!

Alex

Here's a q-n-d manual fit with 1/(c*loglog(n)). It's close, but not quite right.

Alex

We have [tex]M_{m,n}=M_{m,n}[/tex] , right ?

You consider in particular [tex]\lim_{m\to\infty} (\lim_{n\to\infty}|M_{m,n}|/n)/m[/tex]

Which might be [tex]\lim_{n\to\infty}|M_{n,n}|/n^2[/tex]

Which is the probability that a random integer has no prime factor greater than its square root.

(I think you meant [I]M[/I][sub][I]m[/I],[I]n[/I][/sub] = [I]M[/I][sub][I]n[/I],[I]m[/I][/sub])
Yes, [I]M[/I][sub][I]m[/I],[I]n[/I][/sub] is symmetric in [I]m[/I] and [I]n[/I]. Another way to define [I]M[/I][sub][I]m[/I],[I]n[/I][/sub] is
[I]M[/I][sub][I]m[/I],[I]n[/I][/sub] = {i*j, 1 ≤ i ≤ m, 1 ≤ j ≤ n}
which makes the symmetry really obvious.

> Which is the probability that a random integer has no prime factor greater than its square root.

That's a good point! That probability is ρ(2) = 1-ln(2) = 0.30685...., where ρ() is the Dickman function.
The argument sounds plausible, and the y = 1-ln(2) doesn't look all wrong as an asymptote in the graph, either.

Alex

[quote=akruppa;101034]
Or a simple expression for [I]r[/I]([I]m[/I]) = lim[sub][I]n[/I]→∞[/sub] |[I]M[/I][sub][I]m[/I],[I]n[/I][/sub]|/n ? For [I]m[/I]=1, 2, 3, 4, 5, 6, 7 they are [I]r[/I]([I]m[/I]) = 1, 3/2, 2, 29/12, 44/15, 33/10, 134/35 respectively. What about larger [I]m[/I]? [/quote]
Did you compute these limits analytically and/or are you sure of these values to be exact ?
By doing the calculation for the first few, one sees that it makes no sense to use the common denominator, but rather m! to make it an integer sequence.
So I looked up the sequence r(m)*m!

(14:08) gp > for(i=1,length(ak),print1(ak[i]*i!,","))
1,3,12,58,352,2376,19296,

Unfortunately there is no hit if you keep >4 terms, but for the first 4, there are 10 hits:
[URL="http://www.research.att.com/%7Enjas/sequences/?q=1%2C3%2C12%2C58"]http://www.research.att.com/~njas/sequences/?q=1%2C3%2C12%2C58[/URL]

I could imagine that one of those might be correct...

PS: you can easily plot that function with Maple:(takes < 1sec for 8*500)
ak:=(m,n)->nops({ ''i*j'$i=1..m'$j=1..n});
plot( [ '[i,ak(8,i)/i ]'$i=400..500 ]);

PPS: I confirm your 44/15 by doing the calculation, which shows (a) that it's a new sequence, and (b) how a recursive formula can be obtained : to go from r(m-1) to r(m), one has to add:
1/m for the multiples of m in interval ( (m-1)n, mn]
+(m-2)/(m-1)/m for the multiples of m in interval ( (m-2)n, (m-1)n] which are not multiples of m-1, [since gcd(m,m-1) =1 ]
x/m for the multiples of m in interval ( (m-3)n, (m-2)n] which are not multiples of m-1 nor of m-2,
etc. [<= well, this has to be spelled out...]

[quote=akruppa;101034]For [I]m[/I]=1, 2, 3, 4, 5, 6, 7 they are [I]r[/I]([I]m[/I]) = 1, 3/2, 2, 29/12, 44/15, 33/10, 134/35 respectively. What about larger [I]m[/I]?[/quote]
Not very elegant, but...
[code]
p(n)={local( L=1, s=1 ); for( i=1,n-2,
S = Set( vector( L=lcm( L,n-i ), j, n*j ));
for( k=n-i, n-1, S = setminus( S, Set( vector(n*L/k, j, j*k)) ));
if( length(S), s += length(S)/L, break );
);
s/n
}
a = vector(13, k, sum( i=1,k, p(i))
vector(13, k, k!*a[k] )[/code]gives

a=[1, 3/2, 2, 29/12, 44/15, 33/10, 134/35, 1173/280, 967/210, 625/126, 3818/693, 16097/2772, 19100/3003]

and the associated integer sequence

1, 3, 12, 58, 352, 2376, 19296, 168912, 1670976, 18000000, 219916800, 2781561600, 39605760000.

Let [tex] f(n)=|M_{n,n}|/n^2 [/tex].

[tex] \lim_{n\to\infty}f(n) [/tex] is NOT the probability 1-ln(2)=0.306852819 that
a random integer has no prime factor greater than its square root, sorry.
The code below shows that it is smaller: f(10000)=0.225043 and f(150000)=0.207736.

[code]
#include<stdio.h>
#include<stdlib.h>
#include<stdint.h>

int main(int ac, char **av){
int64_t i, j, s = 0, n, n2, p;
char *t;

n = atoll(av[1]);
n2 = n*n/8+1;
t = (char*)malloc(n2*sizeof(char));
for(i = 1; i <= n2; i++) t[i] = 0;
if(!t){ printf("malloc failed\n"); exit(0); }
for(i = 1; i <= n; i++) for(p = i*i, j = i; j <= n; p += i, j++) t[p>>3] |= 1<<(p&7);
n2 = n*n;
for(i = 1; i <= n2; i++) s += (t[i>>3]>>(i&7))&1;
printf("n=%lld s=%lld s/n/n=%f\n", n, s, (1.0)*s/n/n);
}
/*
fado:~$ ./akruppa 10000
n=10000 s=22504348 s/n/n=0.225043
fado:~$ ./akruppa 150000
n=150000 s=4674065960 s/n/n=0.207736
*/
[/code]

We noticed the problem as well, when Alin Bostan pointed out [url]http://www.research.att.com/~njas/sequences/A027424[/url] to us today. Emmanuel Thomé noted that, for example, [I]p[/I][sup]3[/sup] may be <[I]n[/I][sup]2[/sup] and [I]n[/I]-smooth, but isn't present in the multiplication table.

The problem keeps being surprisingly elusive!

Alex

[quote=Pascal Ochem;101488][code](...)
t = (char*)malloc(n2*sizeof(char));
for(i = 1; i <= n2; i++) t[i] = 0;
if(!t){ printf("malloc failed\n"); exit(0); }
(...)[/code][/quote]
Boom! too late : segmentation fault caused airplane crash - millions of innocent citizens died.:mad:

I used Excel and played around a bit, this is the best I could come up with.

Regards
Patrick

It seems to me that
r_m = lim |M_mn|/n = |M_mn|/n whenever n is a multiple of all k<=m,
i.e. in particular, r_m = |M_{m,L}|/L for L = lcm(1,2,...,m) = m# (primorial)
(of course also for L=m! ...)

Thus r_m = sum( p(k) / k , k=1..m )
where p(k) = 1 + sum( | { k . j , j=1..k# } \ (k-1) Z \ (k-2) Z \ ... \ (k-i) Z | , i=1..k-1 ) / k#
where \ p Z means : multiples of p removed. (hope I made no typo...)

As to the cited sequence, which equals |M_nn|, its equal to the # of different numbers in lines 1..n of the following list (the part of the multiplication table below and including the diagonal, or list of multiples of n up to n²) :
[code]
1
2 4
3 6 9
4 8 12 16
5 10 15 20 25
6 12 18 24 30 36
7 14 21 28 35 42 49
...
[/code]
Its easy to see that if n is a prime number, the whole line consists of new elements,
whereas for highly composite (n=6, n=12,...), it may be only the second half (or less?).
If p is the minimum percentage of new numbers in each line,
then we know that M_nn increases on each step k by an amount between kp and k
giving the estimate p (m+1)/2 <= r_m <= (m+1)/2
thus p/2+p/2m <= r_m/m <= 1/2+1/2m.

(just noticed my new avatar... does this mean 'hacker' ?)

[QUOTE=akruppa;101165]masser: r(m)/m is clearly bounded by 0 and 1, however, it is not monotonically decreasing. Emmanuel Thomé wrote Magma code that estimates the ratio reasonably quickly for 1 ≤ m ≤ 200. Graph is attached.

Alex[/QUOTE]

If I understand the original question, then:

I think I have a (partial) understanding of why this ratio will NOT
be uniformly decreasing. ---> Any asymptotic expression for it is likely
to contain the Buchstab function.

Notation: Let S_n = {1,2,3.....n}, let k_n = {k, 2k, 3k, 4k,....} = kS_n
and we consider the Union of k_1, k_2, k_3, .....

The size of this set may be obtained by considering the collection [not
set] of all of the elements { {k_1}, {k_2}, {k_3} ....}.

We want to remove duplicates from this collection. It can be done by
ordering them, and then *sieving* out k, 2k, 3k, ... for the various k in a
manner similar to the sieve of Eratosthenes. We want, essentially "what is
left" [i.e. the unique elements] after the sieve. We now look at a paper by
DeBruijn entitled "The Number of Uncancelled Elements in the Sieve of
Eratosthenes". The Buchstab function [a decaying oscillatory function]
arises naturally in this paper. It's been a long time since I looked at
this paper and my memory may be [probably is] shaky, but I have an
intuitive feeling that the Buchstab function will arise in this context as well.

I

nobody reads me but anyway...
 

Well since nobody protested upon " lcm($1..n) = n# ", nobody seems to read my posts :-( !
In spite of this, I'm nice enough to give you my latest calculation of the sequence r :
r = [1, 3/2, 2, 29/12, 44/15, 33/10, 134/35, 1173/280, 967/210, 625/126,
3818/693, 16097/2772, 19100/3003, 40295/6006, 15178/2145,
355741/48048, 406792/51051, 843925/102102, 2855098/323323,
59051781/6466460, 199307/20995, 1468825/149226, 1624834/156009,
3395027/318136, 562654436/50702925, 386907201/33801950,
138165574/11700675, 166202057/13724100, 1262914124/99499725,
16314534181/1260329850, 17616875026/1302340845,
192071069881/13891635680, 7777852141/548354040,
19218419059/1322376858, 498449878618/33426748355,
18233134502501/1203362940780 ]

computed using
[code]
p := proc(n) option remember; local s,t,i,j: s:=1; t:={}:
for i from n-1 by -1 to 1+n/(min@op@eval@numtheory[factorset])(n) do
t := t union { ilcm(n,i)/n };
t := select( x-> numtheory[divisors](x) intersect t = { x }, t ):
for j in combinat[powerset](t) do s := s+(-1)^nops(j)/ilcm(op(j)) od:
r := m->add( p(n), n=1..m ); map(r,[$1..36]);
[/code]

p(n) = r(n)-r(n-1) (of which r(n)/n is the mean value) is given by
[tex]
p(n) = \sum_{i=1}^n
\sum_{L\subset \{i...n-1\}} \frac{(-1)^{{\rm card}L}}{{\rm lcm} L\cup\{n\}}
[/tex]
The sum over L is in principle over all subsets of {i...n-1}, but one can easily see that one can remove elements of this set which are a multiple of another element;
also, the sum does not contribute anything if {i,...,n-1} contains a divisor of n, so it's sufficient to sum over i larger than the biggest proper divisor of n.

The sums are then of the form, for example
1/x + 1/y + 1/z - 1/lcm(x,y) - 1/lcm(y,z) - 1/lcm(x,z) + 1/lcm(x,y,z)

I did not yet find the formula I need to calculate such sums ; I suppose it is possible at least in special cases like the following:
1/ if n is a prime, then p(n) is bigger than any previous value of p(n) (n>1). In fact, the graph of p(n) looks much like 1/tau(n) (the number of divisors of n).

2/ if n has many divisors (e.g. multiples of 6, 12,...), p(n) is "locally minimal"
(usually globally: "Conjecture: if n>10, p(n)=min{p(1),..,p(n)}, then n=0 mod 6")
and gets close to 1/log( n log n ) which I suggest as candidate for its asymptotic behaviour.


[code]
p(n)={ local( P, L, t=[], s=1, last=0 );if(n<2,return(1));
forstep( i=n-1, n/factor(n)[1,1]+1, -1,
L = t; t = [ P = lcm(n,i)/n ];
/* combine with the old list, removing multiples */
for( i=1, #L, if( P >= L[i] && P % L[i] == 0, t=L; s += last; next(2));
if( L[i] % P, t = concat( t, [L[i]] ));
);
if( #L && (gcd(lcm(L), P)==1), last *= (1-1/P); /* speedup */,
P = lcm(t); /* to use integer arithmetics in sum */
last = sum( j=0, 2^#t -1, (-1)^#( L=vecextract( t,j ))*P/lcm(L))/P;
); s+=last;
); s/n
}
[/code]

Just for the records...
 

I don't know if anybody is still interested in that, and if my results are completely obsolete in view of the MAGMA code you (A.K.) mentioned.
However, just for the records, before this thread gets completely "archived", let me mention that I had a look at the function

f : powerset(N*) -> Q
f(S) = sum_{A subset S} (-1)^#A / lcm(A)

Example: f({})=1, f({m})=1-1/m, f({m,n})=1-1/m-1/n+1/lcm(m,n),...

This is the function that must be evaluated for S={n,n-1}, {n,n-1,n-2},..., {n,n-1,...,2} to calculate p(n) (= r(n)-r(n-1)).

f has lots of interesting properties which are immediate to show:[LIST][*]if 1\in S, f(S)=0[*]if m=gcd(S)>1 then f(S)=1+( f(S/m)-1)/m[*]f(S) = f(S') if S contains in addition to elements of S' only multiples thereof[*]if p\in S is relatively prime to all other elements of S, then f(S)=(1-1/p) f(S)[*]Corollary: if the elements of S are pairwise coprime, then f(S) = product( 1-1/p, p\in S)[*]Corollary: if S contains the first n primes p_1...p_n (and maybe multiples thereof, in particular any number < p[n+1] and >1), then f(S) = \phi(p_n#) / p_n#
where # is primorial (product of these primes) and \phi is Euler's totient function (OEIS A5876)
Example: f({2,3,4,5})=(2-1)(3-1)(5-1)/(2*3*5) = 4/15.[*]More generally, if S can be partitioned such that elements of different subsets never have a common factor, then f(S) is the product of the f(S_i)[/LIST]I suppose the latter is a known and well studied property (we could call it "multiplicative" in this arithmetical sense), but I never heard of before (sorry for my lack of culture ; I'd appreciate any reference)...
Using such decompositions of S, one can calculate f({n,n-1,n-2,...}) and thus r(n) quite far rather quickly. If someone is interested in, I can post the PARI code.
However, I did not yet have the time to search a formula for
f( {n,n-1,...,k} ) (at least for prime n) - I think it should not be too difficult to find, and with this, calculation of r(n) would be immediate for any n.

For the sake of completeness, I also recall that p(n)=r(n)-r(n-1) takes its maximal values for n prime, then p(n)=1/2+something growing slowly like sqrt(n) or log(n) (at these points r(n)/n increases); and it takes smallest values for highly composite n (something like 1/sigma(n), the sum of divisors of n).

If someone is interested in studying this further (I don't know if this could merit writing a paper...), feel free to contact me.

PS:
There are now also 2 OEIS sequences,
[URL="http://www.research.att.com/%7Enjas/sequences/A101459"]A101459[/URL]
a(k) = card { i*j , i <= k, j <= lcm(1,2,3...,k) }.and
[URL="http://www.research.att.com/%7Enjas/sequences/A126959"]A126959[/URL]
a(k) = k! * lim_{n->oo} card({ i*j; i=1..k, j=1..n })/n.
(The latter is r(n)*n!, the former r(n)*lcm(2,...,n).)

[QUOTE=m_f_h;104124]I don't know if anybody is still interested in that, and if my results are completely obsolete in view of the MAGMA code you (A.K.) mentioned. [/QUOTE]

I'm certainly following the thread, but I'm afraid I can't put a lot of time into this problem myself atm, so I don't really have anything interesting to add... :sad:

Alex