Any suggestions are much appreciated.
Cristina
a1 := 1/3*Exp[I*Pi/6]
a2 := 2/3*Exp[I*3*Pi/4]
n := 6
a1c := Conjugate[a1]
a2c := Conjugate[a2]
r1 := Abs[a1]
r2 := Abs[a2]
t1 := Arg[a1]
t2 := Arg[a2]
B[z_] := ((a1c/(r1))*(a1 - z)/(1 - a1c*z))^
n*((a2c/(r2))*(a2 - z)/(1 - a2c*z))^n
a := a1*a2 (a1c + a2c) - (a1 + a2)
b := (1 - (r1)^2*(r2)^2 - ((1 - (Abs[a1*a2])^2)^2 - (Abs[
a])^2)^(1/2))/(a1c*(1 - (r2)^2) + a2c*(1 - (r1)^2))
alpha = Arg[N[B[b]]]
u0[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*0*Pi)/n )]
u1[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*1*Pi)/n )]
u2[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*2*Pi)/n )]
u3[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*3*Pi)/n )]
u4[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*4*Pi)/n )]
u5[rho_, v_] := (rho)^(1/n)*Exp[I*((v + 2*5*Pi)/n )]
sol0 :=
Solve[(1 - u0[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u0[rho, v] - r2)*Exp[-I*t1] + (r2*u0[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u0[rho, v] == 0, z]
sol1 := Solve[(1 - u1[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u1[rho, v] - r2)*Exp[-I*t1] + (r2*u1[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u1[rho, v] == 0, z]
sol2 := Solve[(1 - u2[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u2[rho, v] - r2)*Exp[-I*t1] + (r2*u2[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u2[rho, v] == 0, z]
sol3 := Solve[(1 - u3[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u3[rho, v] - r2)*Exp[-I*t1] + (r2*u3[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u3[rho, v] == 0, z]
sol4 := Solve[(1 - u4[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u4[rho, v] - r2)*Exp[-I*t1] + (r2*u4[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u4[rho, v] == 0, z]
sol5 := Solve[(1 - u5[rho, v]*r1*r2)*Exp[-I (t1 + t2)]*
z^2 + ((r1*u5[rho, v] - r2)*Exp[-I*t1] + (r2*u5[rho, v] - r1)*
Exp[-I*t2])*z + r1*r2 - u5[rho, v] == 0, z]
plotpi0[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol0], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
plotpi1[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol1], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
plotpi2[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol2], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
plotpi3[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol3], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
plotpi4[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol4], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
plotpi5[tmin_, tmax_, c_] :=
ParametricPlot[
Evaluate[{Re[z], Im[z]} /. sol5], {v, 0, 2*Pi}, {rho, tmin, tmax},
ColorFunction -> Function[{x, y, v, rho}, Hue[c, v, rho]],
PlotPoints -> 25, PlotRange -> All]
With[{tmin = .0015, tmax = .01, c = .8},
Show[plotpi0[tmin, tmax, c], plotpi1[tmin, tmax, c],
plotpi2[tmin, tmax, c], plotpi3[tmin, tmax, c],
plotpi4[tmin, tmax, c], plotpi5[tmin, tmax, c], PlotRange -> All]]
When going around the circle, at some point you run into the branch
cut of the square root. Set ExclusionsStyle -> Red to see the
discontinuity. You can choose the branch cut differently for that
portion of the plot:
n = 6; {r1, r2} = {1/3, 2/3}; {t1, t2} = {Pi/6, 3 Pi/4};
u[i_][rho_, v_] := rho^(1/n) Exp[I (v + 2 i Pi)/n]
sol[i_][rho_, v_] := Module[{z}, z /. Solve[
(1 - u[i][rho, v] r1 r2) Exp[-I (t1 + t2)] z^2 +
((r1 u[i][rho, v] - r2) Exp[-I t1] + (r2 u[i][rho, v] - r1) Exp[-I
t2]) z +
r1 r2 - u[i][rho, v] == 0, z]]
Manipulate[ParametricPlot[{Re@ #, Im@ #}& /@ sol[#][rho, v] /.
Sqrt[z_] :> If[1 <= # <= 4, Sqrt[z],
Sqrt[Abs[z]] E^(I Mod[Arg[z], 2 Pi]/2)] // Evaluate,
{v, 0, 2 Pi}, {rho, uu, vv},
ColorFunction -> (Function[{}, Hue[#/6]]), PlotPoints -> {10, 5},
BoundaryStyle -> None, Mesh -> None, Exclusions -> None]& /@
Range[0, 5] // Show[#, PlotRange -> All] &,
{{uu, .0015}, .0001, .002}, {vv, .01, .2}]
Maxim Rytin
m...@inbox.ru
Best wishes,
Cristina
{m1, m2, m3, m4, m5, m6, m7, m8, m9} = {-1, -1, -1, -1, -1, -1, -1, -1,
-1}; Do[
If[Mod[n, 2] == 0, m1 = m1 + 1,
If[Mod[n, 3] == 0, m2 = m2 + 1,
If[Mod[n, 5] == 0, m3 = m3 + 1,
If[Mod[n, 7] == 0, m4 = m4 + 1,
If[Mod[n, 11] == 0, m5 = m5 + 1,
If[Mod[n, 13] == 0, m6 = m6 + 1,
If[Mod[n, 17] == 0, m7 = m7 + 1,
If[Mod[n, 19] == 0, m8 = m8 + 1,
If[Mod[n, 23] == 0, m9 = m9 + 1]]]]]]]]], {n, 1, 6!}];
Print[{m1, m2, m3, m4, m5, m6, m7, m8, m9}]
I want nested 9 times If[Mod[n,Prime[k]]==0,m[k]=m[k]+1],{k,1,9}]
I will be greatfull for any idea!
Best wishes
Artur
len = 9;
mlist = ConstantArray[-1,len]
Could do this procedurally with a nested loop.
Do [Do [If [Mod[n,Prime[k]]==0, mlist[[k]]=mlist[[k]]+1;Break[]],
{k,len}], {n,6!}]
Or use NestWhile to keep looking for a prime divisor,
Do[NestWhile[#+1&,1,(Mod[n,Prime[#]]!=0||(mlist[[#]]=mlist[[#]]+1;False))&,
1,len], {n,6!}]
Me, I'd do it the first way.
Daniel Lichtblau
Wolfram Research
t1 = f1[6!] // Timing
{0.003438,{359,119,47,26,14,11,7,5,3}}
This can be made more flexible; however, much slower.
f2[k_Integer?Positive, i_Integer?Positive] := Module[{f, m},
f[x_] := Catch[Fold[
If[Mod[x, Prime[#2]] == 0, Throw[m[#2] = m[#2] + 1], #1] &,
If[Mod[x, 2] == 0, Throw[m[1] = m[1] + 1]],
Range[2, i]]];
Table[m[n] = -1, {n, i}];
f /@ Range[k];
Table[m[n], {n, i}]];
t2 = f2[6!, 9] // Timing
{0.018913,{359,119,47,26,14,11,7,5,3}}
t2[[1]]/t1[[1]]
5.50116
This can be made slightly faster
f3[k_Integer?Positive, i_Integer?Positive] := Module[{f, m, r, p},
p = Transpose[{r = Range[2, i], Prime[r]}];
f[x_] := Catch[Fold[
If[Mod[x, #2[[2]]] == 0, Throw[m[#2[[1]]] = m[#2[[1]]] + 1], #1] &,
If[Mod[x, 2] == 0, Throw[m[1] = m[1] + 1]],
p]];
Table[m[n] = -1, {n, i}];
f /@ Range[k];
Table[m[n], {n, i}]];
t3 = f3[6!, 9] // Timing
{0.01468,{359,119,47,26,14,11,7,5,3}}
t3[[1]]/{t1[[1]], t2[[1]]}
{4.26992,0.776186}
Verifying that the results are the same:
t1[[2]] == t2[[2]] == t3[[2]]
True
For a larger range and more primes:
f3[7!, 15]
{2519,839,335,191,104,79,57,49,39,31,27,21,18,17,14}
Bob Hanlon
---- Artur <gra...@csl.pl> wrote:
=============
Dear Mathematica Gurus,
Who know how nested or folded multiple If procedure in following:
{m1, m2, m3, m4, m5, m6, m7, m8, m9} = {-1, -1, -1, -1, -1, -1, -1, -1,
-1}; Do[
If[Mod[n, 2] == 0, m1 = m1 + 1,
If[Mod[n, 3] == 0, m2 = m2 + 1,
If[Mod[n, 5] == 0, m3 = m3 + 1,
If[Mod[n, 7] == 0, m4 = m4 + 1,
If[Mod[n, 11] == 0, m5 = m5 + 1,
If[Mod[n, 13] == 0, m6 = m6 + 1,
If[Mod[n, 17] == 0, m7 = m7 + 1,
If[Mod[n, 19] == 0, m8 = m8 + 1,
If[Mod[n, 23] == 0, m9 = m9 + 1]]]]]]]]], {n, 1, 6!}];
Print[{m1, m2, m3, m4, m5, m6, m7, m8, m9}]
I want nested 9 times If[Mod[n,Prime[k]]==0,m[k]=m[k]+1],{k,1,9}]
I will be greatfull for any idea!
Best wishes
Artur
--
Bob Hanlon
Somewhat faster:
countFirstDivisors = Compile[{{max,_Integer},{ndivs,_Integer}},
Module[{mlist=ConstantArray[-1,ndivs],primes=Prime[Range[ndivs]]},
Do [Do [If [Mod[n,primes[[k]]]==0,mlist[[k]]++;Break[]],
{k,ndivs}], {n,max}];
mlist
]]
Example:
In[51]:= Timing[countFirstDivisors[10!,9]]
Out[51]= {2.18167, {1814399, 604799, 241919, 138239, 75402,
58003, 40941, 34478, 26982}}
Much faster is to work out the right formula for these counts. Lucky me,
I did that in a MathGroup thread around a decade ago.
frac[k_] := Product[(Prime[j]-1)/Prime[j],{j,k-1}]/Prime[k]
countFirstDivisors2[max_,ndivs_] :=
Table[-1+Floor[max*frac[k]], {k,ndivs}]
In[55]:= Timing[countFirstDivisors2[10!,9]] // InputForm
Out[55]//InputForm=
{8.895661984809067*^-15, {1814399, 604799, 241919, 138239,
75402, 58001, 40942, 34477, 26982}}
Daniel Lichtblau
Wolfram Research
His suggestion of using Table[-1, {ndivs}] instead of
ConstantArray[-1, ndivs] for Mathematica V6 works fine as seen below.
countFirstDivisors =
Compile[{{max, _Integer}, {ndivs, _Integer}},
Module[{mlist = Table[-1, {ndivs}], primes = =
Prime[Range[ndivs]]},
Do[Do[If[Mod[n, primes[[k]]] == 0, mlist[[k]]++; Break[]], {k, =
ndivs}], {n,
max}];
mlist]]
Timing[countFirstDivisors[10!, 9]]
{2.23923, {1814399, 604799, 241919, 138239, 75402, 58003, 40941, 34478,
26982}}
Cheers ... Syd
Syd Geraghty B.Sc, M.Sc.
Mathematica 6.0.3 for Mac OS X x86 (64 - bit) (21st May, 2008)
MacOS X V 10.5.4
MacBook Pro 2.33 Ghz Intel Core 2 Duo 2GB RAM
On Oct 17, 2008, at 8:33 AM, Daniel Lichtblau wrote:
> Syd Geraghty wrote:
>> Hi Daniel,
>> I have tried to use and understand your input below but I cannot
>> see where the bug is to generate the Error message below.
>> Attached is a Mathematica file that I used to experiment with.
>> Appreciate your help.
>> Cheers Syd
>> In[7]:= countFirstDivisors =
>> Compile[{{max, _Integer}, {ndivs, _Integer}},
>> Module[{mlist = ConstantArray[-1, ndivs], primes =
>> Prime[Range[ndivs]]},
>> Do[Do[If[Mod[n, primes[[k]]] == 0, mlist[[k]]++; Break[]], {k, =
>> ndivs}], {n,
>> max}];
>> mlist]]
>> During evaluation of In[7]:= Compile::"part" : "Part specification =
>> mlist=CE k=CE=A4 cannot be compiled since the \
>> argument is not a tensor of sufficient rank. Evaluation will use
>> the \
>> uncompiled function. ButtonBox["",
>> BaseStyle->"Link",
>> ButtonData:>"paclet:ref/Compile",
>> ButtonFrame->None,
>> .....SNIP
>> Syd Geraghty B.Sc, M.Sc.
>> sydge...@mac.com <mailto:sydge...@mac.com>
>> Mathematica 6.0.3 for Mac OS X x86 (64 - bit) (21st May, 2008)
>> MacOS X V 10.5.4
>> MacBook Pro 2.33 Ghz Intel Core 2 Duo 2GB RAM
>
> I use a development kernel. Apparently ConstantArray is unfamiliar
> to Compile in version 6. Try it instead with
>
> mlist = Table[-1,{ndivs}]
>
> and see if that works better.
>
> Daniel
>
Best wishes
Artur
Arthur,
Calling the above expression "expr", we have
In[5]:= FullSimplify[expr]
Out[5]= Pi
If I have understood your correctly, you wish to have something more
complicated than Pi, yet simpler than the original expr, which begs the
question: What do you expect? Perhaps *PowerExpand* is what you are
looking for (though the leaf count is not that much different)?
In[7]:= pw = PowerExpand[expr];
In[8]:= LeafCount /@ {expr, pw}
Out[8]= {590, 582}
Or you may want to tweak/build your own *ComplexityFunction*. See
http://reference.wolfram.com/mathematica/ref/ComplexityFunction.html
Regards,
-- Jean-Marc
>Dear Mathematica Gurus, Who know which another function as Simplify
>or FullSimplify to use to following formula (Simplify do nothing but
>FullSimplify simplify too much).
<expression snipped>
What is it you are trying to achieve? I assume the result
returned by FullSimplify, Pi, is correct. If so, why would you
want a more complex result than Pi?
Best wishes
Artur
Jean-Marc Gulliet pisze:
> Arthur,
>
> Calling the above expression "expr", we have
>
> In[5]:= FullSimplify[expr]
>
> Out[5]= Pi
>
> If I have understood your correctly, you wish to have something more
> complicated than Pi, yet simpler than the original expr, which begs
> the question: What do you expect? Perhaps *PowerExpand* is what you
> are looking for (though the leaf count is not that much different)?
>
> In[7]:= pw = PowerExpand[expr];
>
> In[8]:= LeafCount /@ {expr, pw}
>
> Out[8]= {590, 582}
>
> Or you may want to tweak/build your own *ComplexityFunction*. See
>
> http://reference.wolfram.com/mathematica/ref/ComplexityFunction.html
>
>
> Regards,
> -- Jean-Marc
>
> __________ Information from ESET NOD32 Antivirus, version of virus
> signature database 3535 (20081018) __________
>
> The message was checked by ESET NOD32 Antivirus.
>
> http://www.eset.com
>
>
>
Bill Rowe pisze:
> On 10/19/08 at 5:41 AM, gra...@csl.pl (Artur) wrote:
>
>
>> Dear Mathematica Gurus, Who know which another function as Simplify
>> or FullSimplify to use to following formula (Simplify do nothing but
>> FullSimplify simplify too much).
>>
>
> <expression snipped>
>
> What is it you are trying to achieve? I assume the result
> returned by FullSimplify, Pi, is correct. If so, why would you
> want a more complex result than Pi?
>
>
> __________ Information from ESET NOD32 Antivirus, version of virus signature database 3537 (20081020) __________