Smarandache Function

Author

Eric W. Weisstein
January 20, 2005

This notebook downloaded from http://mathworld.wolfram.com/notebooks/IntegerSequences/SmarandacheFunction.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/SmarandacheFunction.html.

©2005 Wolfram Research, Inc. except for portions noted otherwise

Code

RSolve[{a[ρ+1]==p a[ρ]+1,a[1]==1},a[ρ],ρ]

FullSimplify[Solve[a==((-1+(p^ν))/(-1+p)),ν],a>0&&ν>0]

mu[1]:=1
mu[n_]:=Max[mu@@@FactorInteger[n]]
mu[p_,1]:=p
mu[p_,alpha_]:=mu[p,alpha]=Module[
{a,k,r,i,nu,k0=alpha(p-1)},
i=nu=Floor[Log[p,1+k0]];
a[1]=1;
a[n_]:=(p^n-1)/(p-1);
k[nu]=Quotient[alpha,a[nu]];
r[nu]=alpha-k[nu]a[nu];
While[r[i]>0,
k[i-1]=Quotient[r[i],a[i-1]];
r[i-1]=r[i]-k[i-1]a[i-1];
i--
];
k0+Plus@@k/@Range[i,nu]
]

LeastInverse[1]:=1
LeastInverse[n_Integer]:=Module[{p=First[Last[FactorInteger[n]]]},
 p^(1+Plus@@Table[Floor[(n-1)/p^k], {k,Floor[Log[p,n-1]]}])
]

Unused Code

Case 2

mu[n_]:=mu[n]=Module[{f=FactorInteger[n]},
Which[
f==={},1,
SameQ[1,Sequence@@(Last/@f)],f[[-1,1]],
Length[f]==1&&f[[1,-1]]≤f[[1,1]],Times@@f[[1]],
True,Max[mu@@@f]
]
]
mu[p_,1]:=mu[p,1]=p
mu[p_,alpha_]:=mu[p,alpha]=Module[
{a,k,r,i,nu,k0=alpha(p-1)},
nu=Floor[Log[p,1+k0]];
a[1]=1;
a[n_]:=(p^n-1)/(p-1);
k[nu]=Quotient[alpha,a[nu]];
r[nu]=alpha-k[nu]a[nu];
i=nu;
While[r[i]>0,
k[i-1]=Quotient[r[i],a[i-1]];
r[i-1]=r[i]-k[i-1]a[i-1];
i--
];
k0+Plus@@k/@Range[i,nu]
]

Case 3

mu[1]:=1
mu[n_]:=mu[n]=Max[mu@@@FactorInteger[n]]
mu[p_,1]:=mu[p,1]=p
mu[p_,alpha_]:=mu[p,alpha]=Module[
{a,k,r,i,nu,k0},
If[alpha≤p,
alpha p,
(* Else *)
k0=alpha(p-1);
nu=Floor[Log[p,1+k0]];
a[1]=1;
a[n_]:=(p^n-1)/(p-1);
k[nu]=Quotient[alpha,a[nu]];
r[nu]=alpha-k[nu]a[nu];
i=nu;
While[r[i]>0,
k[i-1]=Quotient[r[i],a[i-1]];
r[i-1]=r[i]-k[i-1]a[i-1];
i--
];
k0+Plus@@k/@Range[i,nu]
]
]

Much much much slower

GreatestPrimeFactor[1]:=1
GreatestPrimeFactor[n_]:=FactorInteger[n][[-1,1]]

Smarandache[1]:=1
Smarandache[p_?PrimeQ]:=p
Smarandache[n_]:=Module[{k=2},
While[Mod[k!,n]!=0,k++];
k
]

s1=Smarandache/@Range[10^3];//Timing

Smarandache[1]:=1
Smarandache[p_?PrimeQ]:=p
Smarandache[n_]:=Module[{k=GreatestPrimeFactor[n]},
While[Mod[k!,n]!=0,k++];
k
]

s2=Smarandache/@Range[10^3];//Timing

s1==s2

Enumeration

Examples

mu[5,1092]

mu[3,1000]

mu[5,100000]

List

mu/@Range[20]

Brute Force Check

Smarandache/@Range[20]

Special Values

mu/@PerfectNumber/@Range[10]

Mersenne/@Take[MersennePrimeIndices,10]

Exceptional Cases

n does not divide gpf (n)

Position[Equal@@@Table[{mu[n],GreatestPrimeFactor[n]},{n,5000}],False]//Flatten//Timing

Tutescu's Conjecture

(10^5)

Position[Length/@Split[mu5=mu/@Range[10^5]],_?(#>1&)]//Timing

(10^6)

Position[Length/@Split[mu6=mu/@Range[10^6]],_?(#>1&)]//Timing

Case 2

Position[Length/@Split[mu6=mu/@Range[10^6]],_?(#>1&)]//Timing

Case 3

Position[Length/@Split[mu6=mu/@Range[10^6]],_?(#>1&)]//Timing

(10^7)

Position[Length/@Split[mu7=mu/@Range[10^7]],_?(#>1&)]//Timing

MaxMemoryUsed[]/10^6.

Put[mu7,"mu7.m"]//Timing

(10^8)

Try using the caching routine

mu8=mu/@Range[10^8];//Timing

Disable caching.  This will make things (much) slower, but use less memory

mu[1]:=1
mu[n_]:=Max[mu@@@FactorInteger[n]]
mu[p_,1]:=p
mu[p_,alpha_]:=Module[
{a,k,r,i,nu,k0=alpha(p-1)},
nu=Floor[Log[p,1+k0]];
a[1]=1;
a[n_]:=(p^n-1)/(p-1);
k[nu]=Quotient[alpha,a[nu]];
r[nu]=alpha-k[nu]a[nu];
i=nu;
While[r[i]>0,
k[i-1]=Quotient[r[i],a[i-1]];
r[i-1]=r[i]-k[i-1]a[i-1];
i--
];
k0+Plus@@k/@Range[i,nu]
]

mu8=mu/@Range[10^8];//Timing

mu8=mu/@Range[10^8];//Timing

Unfortunately, it still runs out of memory, even on a 5 GB 64-bit machine

So do it in a loop and without caching.  Messy, but what can you do?  It takes >3 hours this way though...

mu8=Table[0,{10^8}];
Do[mu8[[i]]=mu[i],{i,10^8}];//Timing

{ByteCount[mu8],MaxMemoryUsed[]}/10^6. Megabyte

Runs out of memory trying to Put or Save.  So DumpSave instead

DumpSave["mu8.mx",mu8];//Timing

Get["mu8.mx"];//Timing

Position[Length/@Split[mu8],_?(#>1&)]//Timing

MatchQ[mu8,{___,x_,x_,___}]//Timing

MatchQ[mu8,{__,x_,x_,__}]//Timing

Do[If[mu8[[i]]==mu8[[i+1]],Print[{"SAME: ",i,i+1}]],{i,Length[mu8]-1}]//Timing

High-Water Marks

HighWaterMarks[mu/@Range[30]]

Incrementally Smallest Values of μ (n)/n

mu[#]/#&/@Range[20]

{lwm,pos}=LowWaterMarks[mu[#]/#&/@Range[1000]]

{lwm,pos}=LowWaterMarks[mu[#]/#&/@Range[10^6]];//Timing

{lwm,pos}=LowWaterMarks[mu[#]/#&/@Range[10^7]];//Timing

Length[lwm]

{lwm,pos}=LowWaterMarks[mu[#]/#&/@Range[1.5 10^7]];//Timing

Length[lwm]

Through[{Numerator,Denominator}[lwm]]

Requires >5 GB

{lwm,pos}=LowWaterMarks[mu[#]/#&/@Range[5 10^7]];//Timing

Least Inverse

Formula

LeastInverseSlower[1]:=1
LeastInverseSlower[n_Integer]:=Min[#^(1+IntegerExponent[(n-1)!,#])&/@First/@FactorInteger[n]]

LeastInverse[1]:=1
LeastInverse[n_Integer]:=Module[{p=First[Last[FactorInteger[n]]]},
 p^(1+Plus@@Table[Floor[(n-1)/p^k], {k,Floor[Log[p,n-1]]}])
]

Table

LeastInverse/@Range[10]

LeastInverseSlower/@Range[10]

LeastInverse/@Range[5 10^3];//Timing

LeastInverseSlower/@Range[5 10^3];//Timing

Timing[pairs=Split[Sort[MapIndexed[{#1,#2[[1]]}&,mu6]],First[#1]==First[#2]&];]

{#[[1,1]],Last/@#}&/@Take[pairs,10]//ColumnForm

Timing[max=Position[Partition[First/@pairs,2,1],_?(Subtract@@#≠-1&),1,1][[1,1]]]

Table Look Up in OEIS

http://www.research.att.com/projects/OEIS?Anum=A046021

inverse={1,2,3,4,5,9,7,32,27,25,11,243,13,49,125,4096,17,2187,19,625,343,121,23,59049,3125,169,177147,2401,29,78125,31,134217728,1331,289,16807,43046721,37,361,2197,1953125,41,117649,43, 14641,9765625,529,47}

Elegant and Efficient

But eats gobs of memory :|

Timing[pairs=First/@Split[Sort[MapIndexed[{#1,#2[[1]]}&,mu7]],First[#1]==First[#2]&];]

Timing[max=Position[Partition[First/@pairs,2,1],_?(Subtract@@#≠-1&),1,1][[1,1]]]

inverse=Last/@Take[pairs,max]

High Water Marks

HighWaterMarks[a/@Range[100]]

HighWaterMarks[first]

HighWaterMarks[first]

FactorInteger/@{1,2,3,4,5,9,32,243,4096,59049,177147,134217728}//ColumnForm

Plot

<<Utilities`Plot`

IntegerListPlot[mu/@Range[100],AxesLabel->TraditionalForm/@{n,μ[n]},PlotRange->All]

Number of n with S(n)==k

a[1]:=1
a[k_]:=DivisorSigma[0,k!]-DivisorSigma[0,(k-1)!]

a/@Range[30]

IntegerListPlot[a/@Range[50],AxesLabel->TraditionalForm/@{n,HoldForm[a[n]]},PlotRange->All]