Though memory efficient, primes can be much slower than prime_range since it uses a different algorithm and caching strategy.
{{{id=1| time v = list(primes(10^6)) /// Time: CPU 0.50 s, Wall: 0.50 s }}} {{{id=12| time v = prime_range(10^6) /// Time: CPU 0.05 s, Wall: 0.05 s }}} {{{id=16| /// }}}Primality Testing
We can test whether or not a number is prime using the is_prime function. There is also a function is_pseudoprime, which is potentially much, much faster, especially on huge numbers, but in theory could claim a number to be composite even though it is prime (there are no known examples of this, but it surely happens).
{{{id=15| is_prime(2011) /// True }}} {{{id=13| is_prime(2009) /// False }}} {{{id=17| is_pseudoprime(2009) /// False }}}next_prime and next_probable_prime find the next prime (or pseudoprime) after a number.
{{{id=21| n = next_probable_prime(10^300) /// }}} {{{id=19| time is_pseudoprime(n) /// True Time: CPU 0.01 s, Wall: 0.01 s }}}Takes much, much longer, but algorithm is provably correct.
{{{id=20| time is_prime(n) /// True Time: CPU 7.76 s, Wall: 7.77 s }}} {{{id=25| /// }}}Factoring Integers
Use the command factor to factor an integer. It also has a verbose= option, which if you set to 4 or 8 produces lots of output about the factoring algorithms as they are tried. I don't know a good place to read about the format of the verbose= output, except for reading the source code of PARI.
{{{id=35| factor(2012) /// 2^2 * 503 }}} {{{id=36| factor(10^50 + 4) /// 2^2 * 13 * 89 * 21607605877268798617113223854796888504753673293 }}} {{{id=37| factor(10^50 + 4, verbose=8) /// OddPwrs: is 2276944211802351761945170668051973 ...a 3rd, 5th, or 7th power? modulo: resid. (remaining possibilities) 211: 88 (3rd 1, 5th 1, 7th 0) 209: 34 (3rd 0, 5th 1, 7th 0) 61: 30 (3rd 0, 5th 0, 7th 0) OddPwrs: examining 2276944211802351761945170668051973 *** Warning: IFAC: untested integer declared prime. 2276944211802351761945170668051973 Starting APRCL: Choosing t = 840 Solving the triangular system Solving the triangular system Jacobi sums and tables computed Step4: q-values (# = 14): 421 281 211 71 61 43 41 31 29 13 11 7 5 3 Step5: testing conditions lp Step6: testing potential divisors Individual Fermat powerings: 2 : 7 3 : 7 4 : 8 5 : 8 7 : 6 8 : 2 Number of Fermat powerings = 38 Maximal number of nondeterministic steps = 0 2^2 * 13 * 89 * 21607605877268798617113223854796888504753673293 }}}If you multiply together two random primes with 50 digits each say, the resulting number will be too hard for Sage too factor in a reasonable amount of time.
Notice below that often when make our number 2 digits longer, factoring gets twice as time consuming in Sage.
{{{id=39| def hard(k): return next_prime(10^k) * next_prime(10^k + 10^(k-1)) /// }}} {{{id=38| for k in [20..27]: n = hard(k) t = cputime(); f=factor(n) print k, n, cputime(t) /// 20 11000000000000000015590000000000000004407 0.15 21 1100000000000000000151700000000000000002691 0.41 22 110000000000000000000129000000000000000000027 0.5 23 11000000000000000000016970000000000000000004797 0.88 24 1100000000000000000000108700000000000000000000707 0.98 25 110000000000000000000000433000000000000000000000377 1.65 26 11000000000000000000000013070000000000000000000003819 3.97 27 1100000000000000000000000136300000000000000000000002369 5.07 }}} {{{id=41| /// }}} {{{id=30| /// }}}Counting Primes
The function prime_pi computes the number of primes up to x. You can also use the plot command to draw a nice staircase plot of prime_pi.
{{{id=42| prime_pi(10) /// 4 }}} {{{id=45| prime_pi(10.7) /// 4 }}} {{{id=34| prime_pi(100) /// 25 }}} {{{id=43| prime_pi(1000) /// 168 }}} {{{id=44| plot(prime_pi, 1, 100) ///
}}}
Due to the very hard work of Andrew Ohana -- a UW undergraduate math major -- the prime_pi function can compute $\pi(x)$ for surprisingly large values of $\pi$. This illustrates that the function does not just make a list of all primes up to $x$ and see how long the list is, since if it did that, then the example below would cause my computer to run out of memory (since 100 billion is a big number).
{{{id=48| time prime_pi(10^11) /// 4118054813 Time: CPU 1.18 s, Wall: 1.18 s }}}We can now illustrate the Riemann Hypothesis, which asserts that prime_pi(x) and Li(x) are within $\sqrt{x}\log(x)$ of each other.
{{{id=59| x = 10^9 pp = prime_pi(x) print 'pi(x) = %10.1f'%pp print 'Li(x) = %10.1f'%Li(x) print 'x/log(x) = %10.1f'%(x/math.log(x)) print 'sqrt(x)*log(x)= %10.1f'%(math.sqrt(x)*math.log(x)) print '|pi(x)-Li(x)| = %10.1f'%abs(pp - Li(x)) print '|pi(x)-x/l(x)|= %10.1f'%abs(pp - x/math.log(x)) /// pi(x) = 50847534.0 Li(x) = 50849233.9 x/log(x) = 48254942.4 sqrt(x)*log(x)= 655327.2 |pi(x)-Li(x)| = 1699.9 |pi(x)-x/l(x)|= 2592591.6 }}} {{{id=46| x = var('x') B = 10^5 G = plot(prime_pi, 2, B) + plot(Li, 2, B, color='red') + plot(x/log(x), 2, B, color='green') G += plot(lambda x: prime_pi(x) - math.sqrt(x)*math.log(x), 2, B, color='black') G += plot(lambda x: prime_pi(x) + math.sqrt(x)*math.log(x), 2, B, color='black') G ///
}}}
{{{id=58|
///
}}}
{{{id=29|
///
}}}
Evaluating the Riemann Zeta Function
We can use zeta to evaluate the Riemann Zeta function with any complex number as input. Use the N command or coercion to CC (the complex field) to give a numerical answer.
{{{id=33| zeta(2) /// 1/6*pi^2 }}} {{{id=28| zeta(3+I) /// zeta(I + 3) }}} {{{id=27| CC(zeta(3+I)) /// 1.10721440843141 - 0.148290867178175*I }}} {{{id=23| zeta(CC(3+I)) /// 1.10721440843141 - 0.148290867178175*I }}} {{{id=50| N(zeta(3+I)) /// 1.10721440843141 - 0.148290867178175*I }}} {{{id=51| N(zeta(3+I), 100) /// 1.1072144084314091956251002058 - 0.14829086717817534849076412567*I }}} {{{id=55| complex_plot(zeta, (-30,30), (-30,30)) ///
}}}
{{{id=54|
plot(lambda y: abs(zeta(1/2+I*y)), (0,50))
///
}}}
{{{id=52|
G = plot3d(lambda x, y: abs(zeta(x+I*y)), (.2,.7), (0,50), plot_points=100)
G.show(viewer='tachyon', aspect_ratio=[8,.3,1], zoom=1.4)
///
}}}
{{{id=61|
///
}}}