Opened 8 years ago

Closed 4 years ago

# Change coprime_integers(n) in Integer class to use sieving instead of naive gcd

Reported by: Owned by: spice Simon Spicer minor sage-8.1 basic arithmetic coprime, integer was Simon Spicer, David Roe Frédéric Chapoton, Yann Laigle-Chapuy N/A 482909f 482909f3df4ce42865420b4693135c16483195ad

### Description

As written currently, the coprime_integers(n) method in the Integer class uses gcd naively:

```    def coprime_integers(self, m):
...
v = []
for n in range(1,m):
if self.gcd(n) == 1:
v.append(Integer(n))
return v
```

This method should be rewritten to use sieving to speed things up.

### comment:1 Changed 8 years ago by spice

• Status changed from new to needs_review

### comment:2 Changed 8 years ago by was

• Status changed from needs_review to needs_work

It can't work because this makes no sense:

```for p in self.prime_divisors:
```

### Changed 8 years ago by spice

Included missing parentheses. Doctests now properly pass.

### comment:3 Changed 8 years ago by spice

• Status changed from needs_work to needs_review

### comment:4 Changed 8 years ago by was

• Reviewers set to William Stein, Hao Chen
• Status changed from needs_review to needs_work

Hey, it's not 2011 anymore, as much as I like that number

```    - Simon Spicer (2011-11-12): Rewritten to use sieving methods
```

### Changed 8 years ago by spice

Fixes date typo, now uses naive method for small inputs.

### comment:5 Changed 8 years ago by spice

New patch replaces previous patch. The method now uses naive gcd if both self and m are < 1000, since this seems to be faster.

Some timings. Old code:

```sage: %timeit 100.coprime_integers(100)
10000 loops, best of 3: 44.5 us per loop
sage: %timeit 101.coprime_integers(101)
10000 loops, best of 3: 52.6 us per loop
sage: %timeit 1000.coprime_integers(1000)
1000 loops, best of 3: 529 us per loop
sage: %timeit (10^6+12).coprime_integers(10^6+12)
1 loops, best of 3: 603 ms per loop
sage: %timeit (next_prime(10^6)).coprime_integers(10^6)
1 loops, best of 3: 795 ms per loop
sage: %timeit (10^8+512).coprime_integers(10^5)
10 loops, best of 3: 55.9 ms per loop
sage: %timeit (10^8+512).coprime_integers(10^7)
1 loops, best of 3: 5.61 s per loop
```

New code (without naive method for small inputs):

```sage: %timeit 100.coprime_integers(100)
10000 loops, best of 3: 145 us per loop
sage: %timeit 101.coprime_integers(101)
10000 loops, best of 3: 139 us per loop
sage: %timeit 1000.coprime_integers(1000)
1000 loops, best of 3: 510 us per loop
sage: %timeit (10^6+12).coprime_integers(10^6+12)
1 loops, best of 3: 373 ms per loop
sage: %timeit (next_prime(10^6)).coprime_integers(10^6)
1 loops, best of 3: 383 ms per loop
sage: %timeit (10^8+512).coprime_integers(10^5)
10 loops, best of 3: 39.8 ms per loop
sage: %timeit (10^8+512).coprime_integers(10^7)
1 loops, best of 3: 3.98 s per loop
#MemoryError for m=10^8
```

The new code should be faster when self has few prime factors, but both old and new code scale with m, the bound to which integers coprime to self are returned, since we need to fit that many sage integers in memory.

### comment:6 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:7 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.2 to sage-6.3

### comment:8 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.3 to sage-6.4

### comment:9 Changed 4 years ago by chapoton

• Branch set to public/15404
• Commit set to f52e29cc41ca5fcba60c7a6409f684d280bfd56e
• Status changed from needs_work to needs_review

New commits:

 ​f52e29c `trac 15404 sieve for coprime_integers`

### comment:10 Changed 4 years ago by chapoton

• Milestone changed from sage-6.4 to sage-8.1

### comment:11 Changed 4 years ago by chapoton

• Status changed from needs_review to needs_work

### comment:12 Changed 4 years ago by git

• Commit changed from f52e29cc41ca5fcba60c7a6409f684d280bfd56e to 21081f288e38a3eacfeff00370a5c78b60493d7f

Branch pushed to git repo; I updated commit sha1. New commits:

 ​21081f2 `trac 15404 typo`

### comment:13 Changed 4 years ago by chapoton

• Status changed from needs_work to needs_review

### comment:14 Changed 4 years ago by git

• Commit changed from 21081f288e38a3eacfeff00370a5c78b60493d7f to 6b01cb8d749419f463e91d65618961663464021a

Branch pushed to git repo; I updated commit sha1. New commits:

 ​6b01cb8 `New implementation for coprime_integers`

### comment:15 Changed 4 years ago by roed

• Authors changed from Simon Spicer to Simon Spicer, David Roe

I've written a new implementation that fixes some problems with the old and speeds it up:

• It now handles the case that `self = 0` or negative and that `m` is too large or negative gracefully
• It uses `trial_division()` rather than `prime_divisors()`, which will be faster and deals with the case that `self` has prime factors larger than `m` better.
• It uses bit fiddling rather than a list of booleans.

### comment:16 Changed 4 years ago by git

• Commit changed from 6b01cb8d749419f463e91d65618961663464021a to d0b1014e22a964ac8b2dc06989a0e7fdd64f3ad2

Branch pushed to git repo; I updated commit sha1. New commits:

 ​d0b1014 `Fix bugs in coprime_integers`

### comment:17 Changed 4 years ago by chapoton

• Reviewers changed from William Stein, Hao Chen to Frédéric Chapoton
• Status changed from needs_review to positive_review

ok

### comment:18 follow-up: ↓ 19 Changed 4 years ago by ylchapuy

Just my 2 cents: `mpz_divexact_ui` is expected to be faster than `mpz_fdiv_q`, and `mpz_sgn` instead of comparing to zero. Also as you have `plong` available (why not unsigned?), it's probably better to use `mpz_cmp_si` and `mpz_divisible_ui_p`.

All this is nitpicking though.

### comment:19 in reply to: ↑ 18 Changed 4 years ago by ylchapuy

and `mpz_sgn` instead of comparing to zero.

Oups, forget this was, I misread :)

### comment:20 Changed 4 years ago by git

• Commit changed from d0b1014e22a964ac8b2dc06989a0e7fdd64f3ad2 to 5c9fa776a79efc782ed4fe2d4bbced68fbfb582e
• Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

 ​5c9fa77 `Changing some gmp functions in coprime_integers and adjusting description of algorithm`

### comment:21 Changed 4 years ago by roed

Thanks. I've adjusted `plong` and `ilong` to unsigned longs, adjusted the functions and described the algorithm correctly.

### comment:22 Changed 4 years ago by ylchapuy

Am I missing something (it's late here) or is this never happening: `mpz_tstbit(sieve.value, plong)` ? If `p == m` it can't be composite because it's smallest factor would have been removed before.

### comment:23 Changed 4 years ago by git

• Commit changed from 5c9fa776a79efc782ed4fe2d4bbced68fbfb582e to bbeb57c7a0128b1910f55018c9ce2da52df135a2

Branch pushed to git repo; I updated commit sha1. New commits:

 ​bbeb57c `Remove extraneous tests in coprime_integers implementation`

### comment:24 Changed 4 years ago by roed

Thanks. I removed that and another test, which were left over from before I added the last break statement.

### comment:25 follow-up: ↓ 27 Changed 4 years ago by ylchapuy

• Reviewers changed from Frédéric Chapoton to Frédéric Chapoton, Yann Laigle-Chapuy
• Status changed from needs_review to positive_review

Just curious: what's the rationale for the added argument `algorithm = None`? Are you planning to add another implementation?

Otherwise LGTM.

### comment:26 Changed 4 years ago by git

• Commit changed from bbeb57c7a0128b1910f55018c9ce2da52df135a2 to 224db94a29347649fba38805f8b9f40794d8a99a
• Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

 ​224db94 `Remove algorithm keyword`

### comment:27 in reply to: ↑ 25 Changed 4 years ago by roed

• Status changed from needs_review to positive_review

Just curious: what's the rationale for the added argument `algorithm = None`? Are you planning to add another implementation?

That was left over from when I was still experimenting with other implementations.

Otherwise LGTM.

I've deleted it and set the ticket back to positive review. Thanks!

### comment:28 Changed 4 years ago by fbissey

Hum, gives me time out in `sage/groups/additive_abelian/qmodnz.py` the backtrace points to stuff here

```Stack backtrace
---------------
No symbol table info available.
#1  0x00007f83b5e85b53 in print_enhanced_backtrace () from /usr/lib64/python2.7/site-packages/cysignals/signals.so
No symbol table info available.
#2  0x00007f83b5e85c4a in sigdie () from /usr/lib64/python2.7/site-packages/cysignals/signals.so
No symbol table info available.
#3  0x00007f83b5e88247 in cysigs_signal_handler () from /usr/lib64/python2.7/site-packages/cysignals/signals.so
No symbol table info available.
#4  <signal handler called>
No symbol table info available.
#5  0x00007f83b52f0b6a in __gmpz_divexact_ui () from /usr/lib64/libgmp.so.10
No symbol table info available.
#6  0x00007f83aa3ef26e in __pyx_pw_4sage_5rings_7integer_7Integer_131coprime_integers () from /usr/lib64/python2.7/site-packages/sage/rings/integer.so
No symbol table info available.
```

### comment:29 Changed 4 years ago by chapoton

• Status changed from positive_review to needs_work

### comment:30 follow-up: ↓ 31 Changed 4 years ago by ylchapuy

Hum,

```p = slf.trial_division(mlong, mpz_get_si(p.value)+1)
ilong = plong = mpz_get_ui(p.value)
```

this might be problematic if `p` is set to `slf` and doesn't fit in a long.

If moreover `plong` end up being `1`, the `while mpz_divisible_ui_p(slf.value, plong):` loop will be infinite.

### comment:31 in reply to: ↑ 30 ; follow-up: ↓ 32 Changed 4 years ago by roed

Hum,

```p = slf.trial_division(mlong, mpz_get_si(p.value)+1)
ilong = plong = mpz_get_ui(p.value)
```

this might be problematic if `p` is set to `slf` and doesn't fit in a long.

That can't happen since `trial_division` only goes up to the bound `mlong`.

If moreover `plong` end up being `1`, the `while mpz_divisible_ui_p(slf.value, plong):` loop will be infinite.

That can't happen since `trial_division` starts at `p+1`, which is at least 2.

Last edited 4 years ago by roed (previous) (diff)

### comment:32 in reply to: ↑ 31 Changed 4 years ago by roed

Hum,

```p = slf.trial_division(mlong, mpz_get_si(p.value)+1)
ilong = plong = mpz_get_ui(p.value)
```

this might be problematic if `p` is set to `slf` and doesn't fit in a long.

That can't happen since `trial_division` only goes up to the bound `mlong`.

Oh, I was wrong: in this case `p` will equal `slf`. I'll fix it.

### comment:33 Changed 4 years ago by git

• Commit changed from 224db94a29347649fba38805f8b9f40794d8a99a to 9fc7b4dba99945e66986775fc1325dd9512d612d

Branch pushed to git repo; I updated commit sha1. New commits:

 ​9fc7b4d `Add sig_on/off in coprime_integers and fix broken comparison if p doesn't fit in long`

### comment:34 Changed 4 years ago by roed

• Status changed from needs_work to needs_review

### comment:35 Changed 4 years ago by ylchapuy

• Status changed from needs_review to needs_work

Do not compile.

`mpz_cmp_si(p, mlong)` should be `mpz_cmp_si(p.value, mlong)`

### comment:36 Changed 4 years ago by git

• Commit changed from 9fc7b4dba99945e66986775fc1325dd9512d612d to 4f47b563bf6e18666ab78e89a940fa1ff3934d88

Branch pushed to git repo; I updated commit sha1. New commits:

 ​fd1c03a `Merge commit '9fc7b4dba9' into t/15404/public/15404` ​4f47b56 `Fix typo in coprime_integers`

### comment:37 Changed 4 years ago by roed

• Status changed from needs_work to needs_review

Oops. Fixed.

### comment:38 Changed 4 years ago by ylchapuy

• Status changed from needs_review to needs_work

patchbot failures

### comment:39 Changed 4 years ago by ylchapuy

I think `self = 1` leads to an infinite loop.

### comment:40 Changed 4 years ago by git

• Commit changed from 4f47b563bf6e18666ab78e89a940fa1ff3934d88 to 482909f3df4ce42865420b4693135c16483195ad

Branch pushed to git repo; I updated commit sha1. New commits:

 ​482909f `Fix two bugs in coprime_integers and add test`

### comment:41 Changed 4 years ago by roed

• Status changed from needs_work to needs_review

You're right. I've fixed that (plus a problem that would show up if `self` were negative) and added a test.

Tests pass....

### comment:43 Changed 4 years ago by ylchapuy

• Status changed from needs_review to positive_review

The tests are a useful addition I think :) LGTM, bots problems seems unrelated

### comment:44 Changed 4 years ago by vbraun

• Branch changed from public/15404 to 482909f3df4ce42865420b4693135c16483195ad
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.