Opened 8 years ago
Closed 4 years ago
#15404 closed enhancement (fixed)
Change coprime_integers(n) in Integer class to use sieving instead of naive gcd
Reported by:  spice  Owned by:  Simon Spicer 

Priority:  minor  Milestone:  sage8.1 
Component:  basic arithmetic  Keywords:  coprime, integer 
Cc:  was  Merged in:  
Authors:  Simon Spicer, David Roe  Reviewers:  Frédéric Chapoton, Yann LaigleChapuy 
Report Upstream:  N/A  Work issues:  
Branch:  482909f (Commits, GitHub, GitLab)  Commit:  482909f3df4ce42865420b4693135c16483195ad 
Dependencies:  Stopgaps: 
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.
Attachments (2)
Change History (46)
comment:1 Changed 8 years ago by
 Status changed from new to needs_review
comment:2 Changed 8 years ago by
 Status changed from needs_review to needs_work
comment:3 Changed 8 years ago by
 Status changed from needs_work to needs_review
Parentheses have been added.
comment:4 Changed 8 years ago by
 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 (20111112): Rewritten to use sieving methods
comment:5 Changed 8 years ago by
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
 Milestone changed from sage6.1 to sage6.2
comment:7 Changed 7 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:8 Changed 7 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:9 Changed 4 years ago by
 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
 Milestone changed from sage6.4 to sage8.1
comment:11 Changed 4 years ago by
 Status changed from needs_review to needs_work
comment:12 Changed 4 years ago by
 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
 Status changed from needs_work to needs_review
comment:14 Changed 4 years ago by
 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
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 thatm
is too large or negative gracefully  It uses
trial_division()
rather thanprime_divisors()
, which will be faster and deals with the case thatself
has prime factors larger thanm
better.  It uses bit fiddling rather than a list of booleans.
comment:16 Changed 4 years ago by
 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
 Reviewers changed from William Stein, Hao Chen to Frédéric Chapoton
 Status changed from needs_review to positive_review
ok
comment:18 followup: ↓ 19 Changed 4 years ago by
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
comment:20 Changed 4 years ago by
 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
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
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
 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
Thanks. I removed that and another test, which were left over from before I added the last break statement.
comment:25 followup: ↓ 27 Changed 4 years ago by
 Reviewers changed from Frédéric Chapoton to Frédéric Chapoton, Yann LaigleChapuy
 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
 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
 Status changed from needs_review to positive_review
Replying to ylchapuy:
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
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/sitepackages/cysignals/signals.so No symbol table info available. #2 0x00007f83b5e85c4a in sigdie () from /usr/lib64/python2.7/sitepackages/cysignals/signals.so No symbol table info available. #3 0x00007f83b5e88247 in cysigs_signal_handler () from /usr/lib64/python2.7/sitepackages/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/sitepackages/sage/rings/integer.so No symbol table info available.
comment:29 Changed 4 years ago by
 Status changed from positive_review to needs_work
comment:30 followup: ↓ 31 Changed 4 years ago by
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 ; followup: ↓ 32 Changed 4 years ago by
Replying to 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 toslf
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 being1
, thewhile 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.
comment:32 in reply to: ↑ 31 Changed 4 years ago by
Replying to roed:
Replying to 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 toslf
and doesn't fit in a long.That can't happen since
trial_division
only goes up to the boundmlong
.
Oh, I was wrong: in this case p
will equal slf
. I'll fix it.
comment:33 Changed 4 years ago by
 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
 Status changed from needs_work to needs_review
comment:35 Changed 4 years ago by
 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
 Commit changed from 9fc7b4dba99945e66986775fc1325dd9512d612d to 4f47b563bf6e18666ab78e89a940fa1ff3934d88
comment:39 Changed 4 years ago by
I think self = 1
leads to an infinite loop.
comment:40 Changed 4 years ago by
 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
 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.
comment:42 Changed 4 years ago by
Tests pass....
comment:43 Changed 4 years ago by
 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
 Branch changed from public/15404 to 482909f3df4ce42865420b4693135c16483195ad
 Resolution set to fixed
 Status changed from positive_review to closed
It can't work because this makes no sense: