Возведение в степень Mathematica и нахождение указанного коэффициента

Question

У меня есть следующий код, и он делает именно то, что я хочу, за исключением того, что это смехотворно медленно. Я не был бы так обеспокоен, за исключением того, что, когда я обрабатываю код «вручную», то есть я разбиваю его на части и делаю их по отдельности, это почти мгновенно.

Вот мой код:

Coefficient[Product[Sum[x^(j*Prime[i]), {j, 0, Floor[q/Prime[i]]}], 
                        {i, 1, PrimePi[q]}], x, q]

Изображение добавлено для наглядности:

Я думаю, что он пытается оптимизировать сумму, но я не уверен. Есть ли способ остановить это?

Кроме того, поскольку все мои коэффициенты положительны, а я хочу только x ^ qth, есть ли способ заставить Mathematica отбросить все показатели, которые больше, чем это, и не делать все умножения с этими?

Answer 1

Возможно, я неправильно понимаю, что вы хотите, но, поскольку коэффициент будет зависеть от q, я предполагаю, что вы хотите, чтобы он был оценен для конкретной q. Поскольку я подозревал (как и вы), что для оптимизации продукта и суммы времени требуется время, я переписал его. У вас было что-то вроде:

With[{q = 80}, Coefficient[\!\(
\*UnderoverscriptBox[\(\[Product]\), \(i = 1\), \(PrimePi[q]\)]\((
\*UnderoverscriptBox[\(\[Sum]\), \(j = 0\), \(\[LeftFloor]
\*FractionBox[\(q\), \(Prime[i]\)]\[RightFloor]\)]
\*SuperscriptBox[\(x\), \(j*Prime[i]\)])\)\), x, q]] // Timing
(*
-> {8.36181, 10003}
*)

, который я переписал с чисто структурными операциями как

With[{q = 80},
 Coefficient[Times @@ 
 Table[Plus @@ Table[x^(j*Prime[i]), {j, 0, Floor[q/Prime[i]]}],
        {i, 1, PrimePi[q]}], x, q]] // Timing
(*
-> {8.36357, 10003}
*)

(это просто формирует список терминов, а затем умножает их, поэтому символический анализ не выполняется).

Простое построение полинома происходит мгновенно, но в нем есть несколько тысяч терминов, поэтому, вероятно, происходит то, что Coefficient тратит много времени, чтобы убедиться, что у него правильный коэффициент. На самом деле вы можете решить эту проблему, используя Expand полином. Таким образом:

 With[{q = 80}, Coefficient[Expand[\!\(
 \*UnderoverscriptBox[\(\[Product]\), \(i = 1\), \(PrimePi[q]\)]\((
 \*UnderoverscriptBox[\(\[Sum]\), \(j = 0\), \(\[LeftFloor]
 \*FractionBox[\(q\), \(Prime[i]\)]\[RightFloor]\)]
 \*SuperscriptBox[\(x\), \(j*Prime[i]\)])\)\)], x, q]] // Timing
 (*
 -> {0.240862, 10003}
 *)

и это также работает для моего метода.

Итак, чтобы подвести итог, просто вставьте Expand перед выражением и перед тем, как взять коэффициент.

Answer 2

Я думаю, что причина в том, что оригинальный код медленный, в том, что Coefficient создан для работы даже с очень большими выражениями - такими, которые не поместились бы в памяти, если наивно расширены.

Вот оригинальный полином:

poly[q_, x_] := Product[Sum[ x^(j*Prime[i]), 
                            {j, 0, Floor[q/Prime[i]]}], {i, 1, PrimePi[q]}]

Посмотрите, как для не слишком больших q, расширение полинома занимает гораздо больше памяти и становится довольно медленным:

In[2]:= Through[{LeafCount, ByteCount}[poly[300, x]]] // Timing
        Through[{LeafCount, ByteCount}[Expand@poly[300, x]]] // Timing
Out[2]= { 0.01, { 1859,   55864}}
Out[3]= {25.27, {77368, 3175840}}

Теперь давайте определим коэффициент тремя различными способами и рассчитаем их время

coeff[q_]    := Module[{x}, Coefficient[poly[q, x], x, q]]
exCoeff[q_]  := Module[{x}, Coefficient[Expand@poly[q, x], x, q]]
serCoeff[q_] := Module[{x}, SeriesCoefficient[poly[q, x], {x, 0, q}]]

In[7]:= Table[   coeff[q],{q,1,30}]//Timing
        Table[ exCoeff[q],{q,1,30}]//Timing
        Table[serCoeff[q],{q,1,30}]//Timing
Out[7]= {0.37,{0,1,1,1,2,2,3,3,4,5,6,7,9,10,12,14,17,19,23,26,30,35,40,46,52,60,67,77,87,98}}
Out[8]= {0.12,{0,1,1,1,2,2,3,3,4,5,6,7,9,10,12,14,17,19,23,26,30,35,40,46,52,60,67,77,87,98}}
Out[9]= {0.06,{0,1,1,1,2,2,3,3,4,5,6,7,9,10,12,14,17,19,23,26,30,35,40,46,52,60,67,77,87,98}}

In[10]:=    coeff[100]//Timing
          exCoeff[100]//Timing
         serCoeff[100]//Timing
Out[10]= {56.28,40899}
Out[11]= { 0.84,40899}
Out[12]= { 0.06,40899}

Так что SeriesCoefficient - определенно путь. Если конечно ты не немного лучше в комбинаторике, чем я, и вы знаете следующие формулы простого разбиения ( OEIS )

In[13]:= CoefficientList[Series[1/Product[1-x^Prime[i],{i,1,30}],{x,0,30}],x]
Out[13]= {1,0,1,1,1,2,2,3,3,4,5,6,7,9,10,12,14,17,19,23,26,30,35,40,46,52,60,67,77,87,98}

In[14]:= f[n_]:=Length@IntegerPartitions[n,All,Prime@Range@PrimePi@n]; Array[f,30]
Out[14]= {0,1,1,1,2,2,3,3,4,5,6,7,9,10,12,14,17,19,23,26,30,35,40,46,52,60,67,77,87,98}