Tag Info

Hot answers tagged

14

In general Mathematica cannot compute symbolically infinite sums over primes because of the lack of appropriate mathematical tools. However there are infinite products over primes which are basically well understood on the mathematical level. One famous example is the Euler formula for the Riemann zeta function, one of the most beautiful (and mysterious ...


11

You could use FactorInteger to find out whether or not there are exactly two primes building up a number: SemiPrimeQ[n_Integer] := With[{factors = FactorInteger[n]}, Total[factors[[All, 2]]] == 2 ] The rest is easy: Select[Range[50], SemiPrimeQ] (* {4, 6, 9, 10, 14, 15, 21, 22, 25, 26, 33, 34, 35, 38, 39, 46, 49} *) And for those who like inline ...


10

A naive approach would be this: primePower[n_] := Count[ Range @ n, _?PrimePowerQ] This function works well however it might be very inefficient for large n. It takes a bit to evaluate e.g. primePower[10^6] 78734 which is only a little bigger than PrimePi[10^6] 78498 The latter is much more efficient since it uses advanced algorithms for ...


9

Number theory questions are always a huge accumulator for up votes. :) From my experience I can say that the builtin MangoldtLambda function is pretty slow. So let's define a Mangoldt function on our own. The Mangoldt function is defined by: $\Lambda(n) \equiv \left\{ \begin{array}{1 1} ln\ p & \quad \text{if n = $p^k$ for p a prime}\\ 0 & \quad ...


8

Even though Mathematica has a broad range of powerful capabilites (see e.g. this comparison of computer algebra systems) in related fields (number theory, quantifier elimination) it sometimes doesn't appear to be clever enough to prove simple theorems, e.g. this should yield False however we get back the input: Resolve[ Exists[p, p ∈ Integers && ...


7

Perfect numbers: Select[ Range[10^6], Total[Divisors @ #] == 2 # &] {6, 28, 496, 8128} abundant numbers: Select[ Range[10^3], Total[ Most @ Divisors @ #] > # &]//Short { 12, 18, 20, 24, 30, 36, 40, 42, 48, <<228>>, 968, 972, 978, 980, 984, 990, 992, 996, 1000} I used Short to to get only a few since there are: ...


7

Following Daniel's comment, I post my solution: PolynomialCRT[pol_List, mod_List, x_] := Module[{m, q, l}, l = Length[pol]; m = Table[Times @@ Drop[mod, {i}], {i, l}]; q = Table[PolynomialExtendedGCD[m[[i]], mod[[i]], x][[2, 1]], {i, l}]; Simplify[(q pol).m] ] There is no error checking code (for instance, pol and mod have the same length).


6

In order to understand the issue, we should provide the underlying definitions. Mathematica helps in verifying appropriate relations and definitions. The main functional equation relating Riemann's zeta function $\zeta\;$, to Euler's $\Gamma\;$, established in Riemann's famous paper Über die Anzahl der Primzahlen unter einer gegebener Grösse (1859, English ...


6

We can use Reduce to give us: expr = (2^p - (2^2) (3^2))/(3^3); Reduce[{expr == n, n ∈ Integers}, p, Integers] (n | p) ∈ Integers && n >= -1 && p == Log2[9 (4 + 3 n)] Since the only x such that Log2[x] ∈ Integers are powers of two, we must have 9 (4 + 3 n) equal a power of two while n is simultaneously an integer. This clearly ...


6

I've found that NSum[] takes a bit too long here to compute Riemann's prime-counting function, so I've resorted to generating the terms and summing them: With[{n = 8*^3}, Total[MoebiusMu[Range[n]] N[LogIntegral[1000^(1/Range[n])]/Range[n]], Method -> "CompensatedSummation"]] 168.35915686601484 which gives a result close to that of ...


6

A faster approach to finding Perfect numbers using DivisorSigma Select[Range[10^6], DivisorSigma[1, #] == 2 # &] {6, 28, 496, 8128} Here's an even faster approach: Pick[#, MapThread[Equal, {DivisorSigma[1, #], 2 #}], True] &[Range[10^6]] and a little bit faster: Pick[#, DivisorSigma[1, #] - 2 #, 0] &@Range[10^6] For Abundant ...


6

Here's a fancy memoized solution: Clear[primePowerCount, primePowerCount`cache, iPrimePowerCount] iPrimePowerCount[n1_, n2_] := Count[Range[n1, n2], _?PrimePowerQ] primePowerCount`cache = {1}; primePowerCount[1] = 0; primePowerCount[n_?Positive] := Module[{n0, res}, n0 = First@Nearest[primePowerCount`cache, n]; If[n0 < n, res = ...


6

One can implement a general formula: $$ x = a_1 + \cfrac{b}{a_2 + \cfrac{b}{a_3 + \cfrac{b}{\ddots + \cfrac{b}{a_n} } } } $$ continuedFraction[x_, n_, b_: 1] := Sign[b] Reap[NestWhile[b/# - Sow@Floor[b/#] &, Abs[b]/x, # != 0 &, 1, n]][[2, 1]]; If b = 1 it is the regular continued fraction ContinuedFraction[Sqrt[2], 10] ...


5

If I understand correctly: f[n_] := 13 + Ceiling[n, 4]/2; f[Range[20]] {15, 15, 15, 15, 17, 17, 17, 17, 19, 19, 19, 19, 21, 21, 21, 21, 23, 23, 23, 23} More general approach: sample = {15, 15, 15, 15, 17, 17, 17, 17, 19, 19, 19, 19, 21, 21, 21, 21}; linrec = FindLinearRecurrence[sample] {1, 0, 0, 1, -1} f2[n_] := LinearRecurrence[linrec, ...


5

!! WARNING !! Given the above circumstances the answer below is completely wrong and should not be replicated rather than used for guideline. The problem lies in the fact that Assumptions DON'T affect summation indices hence they are disregarded in the code below. You can evaluate the series using the following snippet: Sum[1/(k n (k + n)^2), {k, 1, ∞}, ...


5

Not too hard; all that's needed is a simple application of matrix identities: ColumnHermiteDecomposition[mat_ /; MatrixQ[mat, IntegerQ]] := Transpose /@ HermiteDecomposition[Transpose[mat]] Test: mat = {{1, 2, 3, 2, 2}, {1, 2, 3, 4, 0}, {0, 5, 4, 2, 1}, {3, 2, 4, 0, 2}}; {u, t} = ColumnHermiteDecomposition[mat]; u {{8, 24, 22, ...


5

Without considering the algorithm the only thing that stands out to me is the use of Delete/Partition/Range rather than the native Drop function, as specific native functions are usually faster. A complication of using Drop is that you will need a different exit condition for the loop. That would look something like this: lucky[z_List] := Module[{i = ...


5

A function from the article that cormullion linked is shorter and faster than what I proposed below. Transcribed in terse style: uf[m_, 1] := {{}} uf[1, n_] := {{}} uf[m_, n_?PrimeQ] := If[m < n, {}, {{n}}] uf[m_, n_] := uf[m, n] = Join @@ Table[Prepend[#, d] & /@ uf[d, n/d], {d, Select[Rest@Divisors@n, # <= m &]}] uf[n_] := uf[n, n] ...


5

(This is a math question, not a Mathematica question.) To add to Artes's answer, there is the well-known(!) identity $$\zeta(-n)=\frac{(-1)^n}{n+1}B_{n+1}$$ so you might as well ask why $$\begin{align*} -\frac12\times B_2&=-\frac1{14}\times B_{14}\\ -\frac12\times \frac16&=-\frac1{14}\times \frac76 \end{align*}$$ A justification for the ...


5

Here's a functional way to use the property (the property, which has been removed from the original question, was $N'/D' = N/D + 1/D'D$ or equivalently $N'D-D'N=1$): farey1[n_] := NestWhileList[ With[{num0 = Numerator[#], den0 = Denominator[#]}, First @ Minimize[{num/den, num den0 - num0 den == 1 && 1 <= den <= n && 1 ...


5

Here's a way to exploit the mediant property of the Farey series. To calculate the mediant: med[{a_, b_}] := (Numerator[a] + Numerator[b])/(Denominator[a] + Denominator[b]); Then the Farey series is: farey[n_] := farey[n] = DeleteCases[ Riffle[ farey[n - 1], med /@ Partition[farey[n - 1], 2, 1]], _?(Denominator[#] > n &)]; with initial ...


5

In short NSum cannot handle this sort of sequence. Indeed, strictly the series is not convergent, and some notion of summability/regularization needs to be chosen. Given the nature of MoebiusMu, "Dirichlet" seems appropriate: Sum[MoebiusMu[k], {k, 1, \[Infinity]}, Regularization -> "Dirichlet"] (* -2 *) Here's how one can see NSum is not working ...


5

I can't complete with Artes's mathematical knowledge and approach, but simply as a point of reference, for formulating a brute-force approach it will be more memory efficient to use Sum, though it will still be very slow for large input. Sum[Boole @ PrimePowerQ @ i, {i, 5*^6}] // Timing MaxMemoryUsed[] {46.535, 348940} 15083688


4

Graham, Knuth, and Patashnik in their book Concrete Mathematics (pages 118 and 150) discuss the Farey series. Their recurrence does not require finding Subsets, computing the elements in order starting with $0/1$ and $1/n$. Although very fast, Subsets can use too much memory when very large $n$ are required, as for some PE problems. ...


4

Here's a way to brute-force search for numbers that have the property that the sum of the digits raised to an integer power is equal to the number itself. list = {}; Do[If[Total[IntegerDigits[b^e]] == b, AppendTo[list, {b, e}]], {b, 2, 800}, {e, 2, 100}]; This returns a list of the numbers and powers (here's just the first 50), ordered so that they ...


3

Be careful. You cannot prove an identity through a finite number of tests! However, you can test the identity directly. ClearAll[q] ((-1)^(1/Denominator[q]))^Numerator[q] === (-1)^q True In fact, the two terms have the same underlying representation: ((-1)^(1/Denominator[q]))^Numerator[q] // FullForm (-1)^q // FullForm Power[-1, q] ...


3

Here is a way to get close. ParallelTable[ Sum[1/(k n (k + n)^2), {k, N@Prime[Range[i]]}, {n, N@Prime[Range[i]]}], {i, 10, 2000, 10}] This gives the following output: {0.0445365, 0.0448078, 0.0448455, 0.0448556, 0.0448597, 0.0448617, <<188>>, 0.0448652, 0.0448652, 0.0448652, 0.0448652, 0.0448652, 0.0448652} Which can be plotted ...


3

The most straightforward way of evaluating your sum in Mathematica is as follows: Sum[1/(Prime[k] Prime[n] (Prime[k] + Prime[n])^2), {k, 1, Infinity}, {n, 1, Infinity}] See for instance the last example on the how-to EvaluateInfiniteSumsAndProducts. This returns unevaluated however, so it seems that if there is a closed form for this series then ...


3

You should use FrobeniusNumber instead of FrobeniusSolve, since it serves this purpose The Frobenius number of $a_{1}, ... a_{n}$,.is the largest integer b for which the Frobenius equation $a_{1} x_{1}+ ... a_{n} x_{n} = k$ has no non-negative integer solutions. The $a_{i}$ must be positive integers. FrobeniusNumber[{6, 9, 20}] 43 Nevertheless ...


3

How you might approach this is going to depend on your function. Here's something to get the discussion rolling. Say the function is z[m_] := Sqrt[m]; and you wish to know what values of m in the range 1 to 1000 have integer-valued z[m]. Then you can use: Position[Boole[IntegerQ /@ z[Range[1, 200]]], 1] // Flatten which gives a list of all the numbers ...



Only top voted, non community-wiki answers of a minimum length are eligible