# Tag Info

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 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 wrote this answer as I was figuring things out. If you just want the answer, copy the definitions of fundQ, Ast, ksFactors, foo, makeOneIndex and magicJ out of the code blocks below. You should have DirichletCharacter[d, magicJ[d], n] == KroneckerSymbol[d, n] Disclaimer 1: This answer is based on plausible interpretations of things not quite said in the ... 6 I get this result as well in both V8 and V9. Product[n^MoebiusMu[n], {n, 1, Infinity}] (* Out: 1/(4*Pi^2) *) It's a simple fact, though, that an infinite product can converge to a non-zero value only if the general term tends to 1. As MoebiusMu takes each of the values$\pm 1$(as well as zero) infinitly often, this product simply can't converge. We ... 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 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 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, primePowerCountcache, iPrimePowerCount] iPrimePowerCount[n1_, n2_] := Count[Range[n1, n2], _?PrimePowerQ] primePowerCountcache = {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 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 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 Here's a functional way to use the property (the property, which has been removed from the original question, wasN'/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 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 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 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 ... 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. ... 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

Artes's solution is the best, I think. If you just want to treat this as an ordinary Diophantine problem, you can do that with Solve[] (making this approach more or less equivalent to Yves's): {p, q} = {-4, 11}; {r, s} = {16, -1}; {x, y} /. Solve[{(q - s) x - (p - r) y == -Det[{{p, q}, {r, s}}], x > 0, y > 0, Min[p, r] < x < ...

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 ...

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