12.07.18 15:15, David Mertz пише:

On Thu, Jul 12, 2018, 7:56 AM Serhiy Storchaka mailto:storchaka@gmail.com> wrote:

* isprime(n) Tests if n is a prime number.

How would you test this? The Miller-Rabin Primality test? For large numbers, iterating through candidate prime divisors is pricey.

* primes() Returns an iterator of prime numbers: 2, 3, 5, 7, 11, 13,...

How would you implements this? Sieve of Eratosthenes is really fun to show students as a Python generator function. But the cached primes DO grow unboundedly as you utilize the generator. Wheel factorization as first pass? Do you cached the first N primes so the each instance of iterator can provide initial elements much faster? What is N?

I didn't mean any concrete implementation. Sure there are enough efficient and simple. I sometimes need a sequence of prime numbers for solving toy problems, and quickly write something like: def primes(n): return (i for i in range(2, n) if all(i % j for j in primes(int(sqrt(i)) + 1))) Having such feature in the stdlib would be very convenient. Any reasonable implementation is enough efficient for my purposes.

12.07.18 17:30, Paul Moore пише:

I'm -1 on having inefficient versions like your primes(n) implementation above in the stdlib. People will use them not realising they may not be suitable for anything other than toy problems.

This is just the shortest version that I write every time when I need primes. For the stdlib implementation we don't have such strong limit on the size.

Your divide_and_round might be nice, but there are multiple common ways of rounding 0.5 (up, down, towards 0, away from 0, to even, to odd, ...) Your implementation does round to even, but that's not necessarily the obvious one (schools often teach round up or down, so use of your implementation in education could be confusing).

This is the one that is useful in many applications and can't be trivially written as expression in Python, but is useful in many applications (date and time, Fraction, Decimal). Divide and rounding down: n//m. Divide and rounding up: (n+m-1)//m or -((-n)//m).

12.07.18 19:14, Steven D'Aprano пише:

Sorry, I don't understand why you say that divide and round up can't be trivially written in Python. Don't you give two implementations yourself?

Sorry, perhaps I was unclear, but I said the opposite.

On Thu, Jul 12, 2018, 9:31 AM Serhiy Storchaka wrote:

12.07.18 15:15, David Mertz пише:

...

On Thu, Jul 12, 2018, 7:56 AM Serhiy Storchaka I didn't mean any concrete implementation. Sure there are enough efficient and simple. I sometimes need a sequence of prime numbers for solving toy problems, and quickly write something like:

def primes(n): return (i for i in range(2, n) if all(i % j for j in primes(int(sqrt(i)) + 1)))

Having such feature in the stdlib would be very convenient. Any reasonable implementation is enough efficient for my purposes.

That's the point I was getting at. "For your purposes" isn't general enough to include in the standard library. There are quite a few algorithms and efficiency trade-offs to decide on. It's not clear to me that any choice is sufficiently "one size fits all" to include. I'm not saying there definitely IS NOT a reasonable compromise approach to these things, but I'd have to judge concrete suggestions... Or better still, people with better knowledge of number theory than mine should make those judgements.

On Fri, Jul 13, 2018 at 4:11 AM, Steven D'Aprano wrote:

There is no reason why primality testing can't be deterministic up to 2**64, and probabilistic with a ludicrously small chance of false positives beyond that. The implementation I use can be expected to fail on average once every 18 thousand years if you did nothing but test primes every millisecond of the day, 24 hours a day. That's good enough for most purposes :-)

What about false negatives? Guaranteed none? The failure mode of the function should, IMO, be a defined and documented aspect of it. ChrisA

I take it back... Tim's (really Will Ness's) version is *always* faster and more memory friendly. I'm still trying to understand it though :-). On Thu, Jul 12, 2018 at 6:05 PM David Mertz wrote:

On Thu, Jul 12, 2018 at 2:22 PM Chris Angelico wrote:

...

On Fri, Jul 13, 2018 at 4:11 AM, Steven D'Aprano wrote:

...

There is no reason why primality testing can't be deterministic up to 2**64, and probabilistic with a ludicrously small chance of false positives beyond that. The implementation I use can be expected to fail on average once every 18 thousand years if you did nothing but test primes every millisecond of the day, 24 hours a day. That's good enough for most purposes :-)

What about false negatives? Guaranteed none? The failure mode of the function should, IMO, be a defined and documented aspect of it.

Miller-Rabin or other pseudo-primality tests do not produce false negatives IIUC.

I'm more concerned in the space/time tradeoff in the primes() generator. I like this implementation for teaching purposes, but I'm well aware that it grows in memory usage rather quickly (the one Tim Peters posted is probably a better balance; but it depends on how many of these you create and how long you run them).

from math import sqrt, ceil

def up_to(seq, lim): for n in seq: if n < lim: yield n else: break

def sieve_generator(): "Pretty good Sieve; skip the even numbers, stop at sqrt(candidate)" yield 2 candidate = 3 found = [] while True: lim = int(ceil(sqrt(candidate))) if all(candidate % prime != 0 for prime in up_to(found, lim)): yield candidate found.append(candidate) candidate += 2

So then how do you implement isprime(). One naive way is to compare it against elements of sieve_generator() until we are equal or larger than the test element. But that's not super fast. A pseudo-primality test is much faster (except for in the first few hundred thousand primes, maybe).

-- Keeping medicines from the bloodstreams of the sick; food from the bellies of the hungry; books from the hands of the uneducated; technology from the underdeveloped; and putting advocates of freedom in prisons. Intellectual property is to the 21st century what the slave trade was to the 16th.

-- Keeping medicines from the bloodstreams of the sick; food from the bellies of the hungry; books from the hands of the uneducated; technology from the underdeveloped; and putting advocates of freedom in prisons. Intellectual property is to the 21st century what the slave trade was to the 16th.

Oops... yeah. I had fixed the "up to and including" previously, but somehow I copied the wrong version to the thread. On Thu, Jul 12, 2018 at 6:36 PM Tim Peters wrote:

[David Mertz]

...

Miller-Rabin or other pseudo-primality tests do not produce false negatives IIUC.

That's true if they're properly implemented ;-) If they say False, the input is certainly composite (but there isn't a clue as to what a factor may be); if the say True, the input is "probably" a prime.

...

I'm more concerned in the space/time tradeoff in the primes() generator. I like this implementation for teaching purposes, but I'm well aware that it grows in memory usage rather quickly (the one Tim Peters posted is probably a better balance; but it depends on how many of these you create and how long you run them).

As you noted later:

I take it back... Tim's (really Will Ness's) version is *always* faster

...

and more

memory friendly.

The original version of this came from an ActiveState recipe that many on c.l.py contributed to (including me) years ago, and that's where all the speed came from. There is, for example, no division or square roots. Will's contribution was the unexpected way to slash memory use, via the generator calling itself recursively. Elegant and effective!

...

I'm still trying to understand it though :-).

Picture doing the Sieve of Eratosthenes by hand, crossing out multiples of "the next" prime you find. That's basically what it's doing, except doing the "crossing out" part incrementally, using a dict to remember, for each prime p, "the next" multiple of p you haven't yet gotten to.

...

from math import sqrt, ceil

def up_to(seq, lim): for n in seq: if n < lim: yield n else: break

def sieve_generator(): "Pretty good Sieve; skip the even numbers, stop at sqrt(candidate)" yield 2 candidate = 3 found = [] while True: lim = int(ceil(sqrt(candidate))) if all(candidate % prime != 0 for prime in up_to(found, lim)): yield candidate found.append(candidate) candidate += 2

Note that this generates squares of odd primes too (9, 25, ...). `up_to` needs to be more like `up_to_and_including`.

-- Keeping medicines from the bloodstreams of the sick; food from the bellies of the hungry; books from the hands of the uneducated; technology from the underdeveloped; and putting advocates of freedom in prisons. Intellectual property is to the 21st century what the slave trade was to the 16th.

On Thu, Jul 12, 2018 at 05:35:54PM -0500, Tim Peters wrote:

[David Mertz]

...

Miller-Rabin or other pseudo-primality tests do not produce false negatives IIUC.

That's true if they're properly implemented ;-) If they say False, the input is certainly composite (but there isn't a clue as to what a factor may be); if the say True, the input is "probably" a prime.

I'm sure Tim knows this, but that's only sometimes true. Since I had no formal education on primality testing and had to pick this up in dribs and drabs over many years, I was under the impression for the longest time that Miller-Rabin was always probabilitistic, but that's not quite right. Without going into the gory details, it turns out that for any N, you can do a brute-force primality test by checking against *every* candidate up to some maximum value. (Which is how trial division works, only Miller-Rabin tests are more expensive than a single division.) Such an exhaustive check is not practical, hence the need to fall back on randomly choosing candidates. However, that's in the most general case. At least for small enough N, there are rather small sets of candidates which give a completely deterministic and conclusive result. For N <= 2**64, it never requires more than 12 Miller-Rabin tests to give a conclusive answer, and for some N, as few as a single test is enough. To be clear, these are specially chosen tests, not random tests. So in a good implementation, for N up to 2**64, there is never a need for a "probably prime" result. It either is, or isn't prime, and there's no question which. In principle, there could be similar small sets of conclusive tests for larger N too, but (as far as I know) nobody has discovered them yet. Hence we fall back on choosing Miller-Rabin tests at random. The chances of an arbitrary composite number passing k such tests is on average 1/4**k, so we can make the average probability of failure as small as we like, by doing more tests. With a dozen or twenty tests, the probability of a false positive (a composite number wrongly reported as prime) is of the same order of magnitude as that of a passing cosmic ray flipping a bit in memory and causing the wrong answer to appear. -- Steve

On Fri, Jul 13, 2018 at 10:40 AM, Steven D'Aprano wrote:

On Fri, Jul 13, 2018 at 04:20:55AM +1000, Chris Angelico wrote:

...

On Fri, Jul 13, 2018 at 4:11 AM, Steven D'Aprano wrote:

...

There is no reason why primality testing can't be deterministic up to 2**64, and probabilistic with a ludicrously small chance of false positives beyond that. The implementation I use can be expected to fail on average once every 18 thousand years if you did nothing but test primes every millisecond of the day, 24 hours a day. That's good enough for most purposes :-)

What about false negatives? Guaranteed none? The failure mode of the function should, IMO, be a defined and documented aspect of it.

If a non-deterministic primality tester such as (but not limited to) Miller-Rabin says a number is composite, it is definitely composite (ignoring bugs in the implementation, or hardware failures). If it says it is prime, in the worst case it is merely "probably prime", with an error rate proportional to 1/4**k where we can choose the k.

(The bigger the k the more work may be done in the worst case.)

You can say that about algorithms easily enough. My point is that this ought to be a constraint on the function - implementations may choose other algorithms, but they MUST follow one pattern or the other, meaning that a Python script can depend on it without knowing the implementation. Like guaranteeing that list.sort() is stable without stipulating the actual sort algo used. ChrisA

On Fri, Jul 13, 2018 at 11:48 AM, Steven D'Aprano wrote:

On Fri, Jul 13, 2018 at 10:56:21AM +1000, Chris Angelico wrote:

...

You can say that about algorithms easily enough. My point is that this ought to be a constraint on the function - implementations may choose other algorithms, but they MUST follow one pattern or the other, meaning that a Python script can depend on it without knowing the implementation. Like guaranteeing that list.sort() is stable without stipulating the actual sort algo used.

I cannot imagine an algorithm that wasn't totally brain-dead (like "flip a coin") which could wrongly report a prime number as composite. Maybe I'm not imaginative enough :-)

Haha. Okay. I'm not familiar with every possible primality test, so I had no idea that they ALL have the same failure mode :) ChrisA

On Fri, Jul 13, 2018 at 2:58 PM, Tim Peters wrote:

[Chris Angelico, on "probable prime" testers]

...

You can say that about algorithms easily enough. My point is that this ought to be a constraint on the function - implementations may choose other algorithms, but they MUST follow one pattern or the other, meaning that a Python script can depend on it without knowing the implementation. Like guaranteeing that list.sort() is stable without stipulating the actual sort algo used.

I agree! Don't worry about it - if this module gets added, I'll make sure of it. My own "probable prime" function starts like so:

def pp(n, k=10): """ Return True iff n is a probable prime, using k Miller-Rabin tests.

When False, n is definitely composite. """

In the context of algorithmic number theory, "no false negatives" is implied by that context, but it should be spelled out anyway. A "probable prime", by definition, is an integer that satisfies some property that _all_ primes satisfy, but that some composites may satisfy too. The variation among probable prime algorithms is in which specific property(ies) they test for.

For example:

def pp1(n): return True

meets the promise, but is useless ;-)

def pp2(n): return n % 3 != 0 or n == 3

is slightly less useless.

But

def pp3(n): return n % 3 != 0

would be incorrect, because pp3(3) returns False - `n % 3 != 0` is not a property that all primes satisfy.

Thank you. Great explanation. Much appreciated! ChrisA

On Thu, Jul 12, 2018 at 5:20 AM Daniel Moisset wrote:

(I don't have a good name): something telling that an integer can be represented exactly as a float

One might also ask that question of e.g. decimal.Decimal, fractions.Fraction, so maybe it's better as a method or somewhere else. (Is this any different than float(x).as_integer_ratio() == (x.numerator, x.denominator) for ints/fractions?) -- Devin

About the name, why not intmath ? And why not using sage ? Or create a package on pypi ? Le jeu. 12 juil. 2018 à 15:11, Serhiy Storchaka a écrit :

12.07.18 14:55, Serhiy Storchaka пише:

...

What are your thoughts about adding a new imath module for integer mathematics? It could contain the following functions:

I forgot yet one useful function:

def divide_and_round(a, b): q, r = divmod(a, b) r *= 2 greater_than_half = r > b if b > 0 else r < b if greater_than_half or r == b and q % 2 == 1: q += 1 return q

(This is a copy of the implementation from datetime.py, there are yet few implementations in Python and C in other files.)

_______________________________________________ Python-ideas mailing list Python-ideas@python.org https://mail.python.org/mailman/listinfo/python-ideas Code of Conduct: http://python.org/psf/codeofconduct/

On 12 July 2018 at 23:41, Serhiy Storchaka wrote:

12.07.18 16:15, Robert Vanden Eynde пише:

...

About the name, why not intmath ?

Because cmath. But if most core developers prefer intmath, I have no objections.

My initial reaction from just the subject title was "But we already have cmath, why would we need imath?", based on the fact that mathematicians write complex numbers as "X + Yi", rather than the "X + Yj" that Python borrowed from electrical engineering (where "i" already had a different meaning as the symbol for AC current). Calling the proposed module "intmath" instead would be clearer to me (and I agree with the rationale that as the number of int-specific functions increases, separating them out from the more float-centric math module makes sense). Beyond that, I think the analogy with the statistics module is a good one: 1. There are a number of integer-specific algorithms where naive implementations are going to be slower than they need to be, subtly incorrect in some cases, or both. With a standard library module, we can provide a robust test suite, and pointers towards higher performance alternatives for folks that need them. 2. For educational purposes, being able to introduce the use cases for a capability before delving into the explanation of how that capability works can be quite a powerful technique (the standard implementations can also provide a cross-check on the accuracy of student implementations). And in the intmath case, there are likely to be more opportunities to delete private implementations from other standard library modules in favour of the public intmath functions. Cheers, Nick. -- Nick Coghlan | ncoghlan@gmail.com | Brisbane, Australia

On Thu, Jul 12, 2018 at 11:05 AM, Steven D'Aprano wrote:

I'm not sure that we need a specific module for integer-valued maths. I think a more important change would be to allow math functions to be written in Python, not just C:

- move the current math library to a private _math library;

- add math.py, which calls "from _math import *".

FWIW, Guido rejected that idea when we added math.isclose() -- Victor had even written it up already. Not that people never change their mind... -CHB -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

Hi In https://bugs.python.org/issue40028 https://bugs.python.org/issue40028?, Ross Rhodes suggested adding to the standard library a function that finds the prime factors of a non-negative integer. So far as I know, any general algorithm for doing this requires a list of prime numbers. To this issue I've just added https://bugs.python.org/issue40028?#msg364757 which reads: A pre-computed table of primes might be better. Of course, how long should

the table be. There's an infinity of primes.

Consider

...

...

...

2**32 4294967296

This number is approximately 4 * (10**9). According to https://en.wikipedia.org/wiki/Prime_number_theorem, there are 50,847,534 primes less than 10**9. So, very roughly, there are 200,000,000 primes less than 2**32.

Thus, storing a list of all these prime numbers as 32 bit unsigned integers would occupy about

...

...

...

200_000_000 / (1024**3) * 4 0.7450580596923828 or in other words 3/4 gigabytes on disk.

A binary search into this list, using as starting point the expected location provided by the prime number theorem, might very well require on average less than two block reads into the file that holds the prime number list on disk. And if someone needs to find primes of this size, they've probably got a spare gigabyte or two.

I'm naturally inclined to this approach because by mathematical research involves spending gigahertz days computing tables. I then use the tables to examine hypotheses. See https://arxiv.org/abs/1011.4269. This involves subsets of the vertices of the 5-dimensional cube. There are of course 2**32 such subsets.

-- Jonathan

On Mon, Mar 23, 2020 at 5:13 AM Neil Girdhar wrote:

I mean:

def binom(n, *ks): # Check that there is at least one ki, and that their sum is less than n, and that they are all nonnegative. # Returns n! / (prod(ki! for ki in ks) * (n-sum(ks))!)

This would still work for binom(n, k), but would also work for the mulinomial case.

Thanks for pulling this up again. I actually would have very much liked to have an efficient and accurate binom(n,k) function available - everything I could find was either floating-point (with the limitations thereof) or too inefficient to be used for ridiculously large values of n and/or k. +1 on moving forward with adding an imath module. ChrisA