Sultan's Dowry Problem

Author

Eric W. Weisstein
September 22, 2007

This notebook downloaded from http://mathworld.wolfram.com/notebooks/RecreationalMath/SultansDowryProblem.nb.

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

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

Code

Daughter[n_]:=Module[{x},Ceiling[x/.FindRoot[HarmonicNumber[x]==HarmonicNumber[n]-1,{x,1},WorkingPrecision->20]]
]

SultansDowry[1]:=0
SultansDowry[n_]:=Ceiling[x/.FindRoot[HarmonicNumber[n]-HarmonicNumber[x]==1,{x,1,100}]]

Table

Daughter/@Range[20]

Table[SultansDowry[n],{n,20}]

Solving

Series of Probability

1st Order

FullSimplify[Normal[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)/.v:HarmonicNumber[q_]:>Series[v,{Variables[q][[1]],∞,1}]],n>0&&r≥0]

D[((n-r+2 n r Log[(n/r)])/(2 (n^2))),r]

Solve[((-1-2 n+2 n Log[(n/r)])/(2 (n^2)))==0,r]

FullSimplify[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)/.r->(E^(-1-(1/(2 n)))) n,n>0]

Limit[%,n->∞]

Series[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)/.r->(E^(-1-(1/(2 n)))) n,{n,∞,1}]

Floor[Table[{(E^(-1-(1/(2 n)))) n,Daughter[n]},{n,10}]]//Timing

Round[Table[{(E^(-1-(1/(2 n)))) n,Daughter[n]},{n,10}]]//Timing

Ceiling[Table[{(E^(-1-(1/(2 n)))) n,Daughter[n]},{n,10}]]//Timing

Subtract@@@Ceiling[Table[{(E^(-1-(1/(2 n)))) n,Daughter[n]},{n,100}]]//Union//Timing

2nd Order

FullSimplify[Normal[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)/.v:HarmonicNumber[q_]:>Series[v,{Variables[q][[1]],∞,2}]],n>0&&r≥0]

D[%,r]

Solve[%==0,r]

FullSimplify[%]

Series of Equality

FullSimplify[Sum[1/(i+1),{i,x,n-1}]]

Solve[HarmonicNumber[x]==HarmonicNumber[n]-1,x]

FullSimplify[Normal[Series[HarmonicNumber[x],{x,∞,1}]],x>0]

FullSimplify[Normal[Series[HarmonicNumber[x],{x,∞,2}]],x>0]

FullSimplify[Solve[EulerGamma+(1/(2 x))+Log[x]==EulerGamma+(1/(2 n))+Log[n]-1,x],n>0]

Solve[EulerGamma+((-1+6 x)/(12 (x^2)))+Log[x]==EulerGamma+((-1+6 n)/(12 (n^2)))+Log[n]-1,x]

Table[-(1/(2 ProductLog[-((E^(1-(1/(2 n))))/(2 n))])),{n,10}]//N

Subtract@@@Ceiling[Table[{-(1/(2 ProductLog[-((E^(1-(1/(2 n))))/(2 n))])),Daughter[n]},{n,4,500}]]//Union//Timing

Plots of Approximations

With[{max=30},
Plot[{-(1/(2 ProductLog[-((E^(1-(1/(2 n))))/(2 n))])),n Exp[-(1+1/(2n))]},{n,1,max},PlotStyle->{Red,Blue},Epilog->{PointSize[.02],Point[{#,Daughter[#]}]&/@Range[max]},
AxesLabel->TraditionalForm/@{n,{StyleForm[Subscript[Overscript[x,"⋀"],1],FontColor->Red],
StyleForm[Subscript[Overscript[x,"⋀"],2],FontColor->Blue]}}
]
]

With[{max=100},
Plot[{-(1/(2 ProductLog[-((E^(1-(1/(2 n))))/(2 n))])),n Exp[-(1+1/(2n))]},{n,1,max},PlotStyle->{Red,Blue},Epilog->{PointSize[.02],Point[{#,Daughter[#]}]&/@Range[max]},
AxesLabel->TraditionalForm/@{n,{StyleForm[Subscript[Overscript[x,"⋀"],1],FontColor->Red],
StyleForm[Subscript[Overscript[x,"⋀"],2],FontColor->Blue]}}
]
]

<<Graphics`Colors`

With[{max=30},
Show[Graphics[{

Join@@@Transpose[{{{Black,PointSize[.03]},{Red,PointSize[.02]},{Blue,PointSize[.015]}},Transpose@Table[Point[{n,Ceiling[#]}]&/@{Daughter[n],-(1/(2 ProductLog[-((E^(1-(1/(2 n))))/(2 n))])),n Exp[-(1+1/(2n))]},{n,4,max}]}]
}],Axes->True]
]

Series Expansions of Approximations

Fit[Daughter/@Range[100],{1,x},x]

Fit[Daughter/@Range[100],{1,x,1/x},x]

1./E

Normal[Series[n Exp[-(1+1/(2n))],{n,∞,2}]]

N[{-1/(2E),1/(8E)}]

Normal[Series[-(1/(2*ProductLog[-(E^(1-1/(2*n))/(2*n))])),{n,Infinity,2}]]

...hangs in v5...

v6

Normal[Series[-(1/(2*ProductLog[-(E^(1-1/(2*n))/(2*n))])),{n,Infinity,2}]]

N[{1/(2E)-1/2,1/(8E)-E/8}]

Probability for n==100

With[{n=100},NMaximize[{((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),0<=r<=n},r]]

((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)/.{r->37,n->100}//N

(cf=FromContinuedFraction/@Reverse[NestList[Most,ContinuedFraction[1/E,10],9]])//InputForm

{Daughter[Denominator[#]],Numerator[#]}&/@cf

100./E

N[1/E,20]

Exact

With[{n=100},NMaximize[{((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),0<=r<=n},r,WorkingPrecision->20]]

D[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),r]

FullSimplify[D[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),r]]

With[{n=100},((EulerGamma-HarmonicNumber[r]+PolyGamma[0,n]-r PolyGamma[1,1+r])/n)]//FullSimplify

With[{n=100},Reduce[(1/100) ((360968703235711654233892612988250163157207/69720375229712477164533808935312303556800)-HarmonicNumber[r]-r PolyGamma[1,1+r])==0&&0<r<50,r,Reals]]

Refine[%,HarmonicNumber[r]∈Reals]//Simplify

N[Root[{-360968703235711654233892612988250163157207+69720375229712477164533808935312303556800 HarmonicNumber[#1]+69720375229712477164533808935312303556800 PolyGamma[1,1+#1] #1&,36.6018821363300057693427494790956692824630336258673054108419`37.06165044590902}],20]

Plots

g1=Show[Block[{$DisplayFunction=Identity},
{
Table[Plot[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),{r,0,n},PlotStyle->Hue[.5+n/50/2]],{n,10,50,5}],
Table[Plot[((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n),{r,0,n},PlotStyle->Hue[n/18]],{n,9}]
}],PlotRange->All,AxesLabel->TraditionalForm/@{x,P[n,x]}]

g2=Plot3D[If[r>n,-(10^10),((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)],{n,1,10},{r,0,10},Mesh->False,PlotPoints->100,ClipFill->None,PlotRange->{0,1},AxesLabel->TraditionalForm/@{n,x,P[n,x]}]

Show[GraphicsArray[{g1,g2}],GraphicsSpacing->-.1]

ContourPlot[If[r>n,0,((1+r(HarmonicNumber[n-1]-HarmonicNumber[r]))/n)],{n,1,10},{r,0,10},ColorFunction->Hue,Contours->50,ContourLines->False,PlotPoints->100]