mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   Math (https://www.mersenneforum.org/forumdisplay.php?f=8)
-   -   Counting Integers and their Multiplicities (https://www.mersenneforum.org/showthread.php?t=7530)

akruppa 2007-03-18 14:45

1 Attachment(s)
(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

m_f_h 2007-03-19 18:08

[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...]

m_f_h 2007-03-20 06:24

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

Pascal Ochem 2007-03-20 12:52

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]

akruppa 2007-03-20 14:14

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

m_f_h 2007-03-20 17:24

[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:

Patrick123 2007-03-20 22:25

1 Attachment(s)
I used Excel and played around a bit, this is the best I could come up with.

Regards
Patrick

m_f_h 2007-03-21 14:01

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.

m_f_h 2007-03-21 14:03

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

R.D. Silverman 2007-03-21 17:39

[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

m_f_h 2007-03-24 17:26

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]


All times are UTC. The time now is 22:11.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.