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: sage-8.1
Component: basic arithmetic Keywords: coprime, integer
Cc: was Merged in:
Authors: Simon Spicer, David Roe Reviewers: Frédéric Chapoton, Yann Laigle-Chapuy
Report Upstream: N/A Work issues:
Branch: 482909f (Commits, GitHub, GitLab) Commit: 482909f3df4ce42865420b4693135c16483195ad
Dependencies: Stopgaps:

Status badges

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)

trac_15404.patch (1.9 KB) - added by spice 8 years ago.
Included missing parentheses. Doctests now properly pass.
trac_15404.2.patch (2.1 KB) - added by spice 8 years ago.
Fixes date typo, now uses naive method for small inputs.

Download all attachments as: .zip

Change History (46)

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

Parentheses have been added.

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:

f52e29ctrac 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:

21081f2trac 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:

6b01cb8New 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:

d0b1014Fix 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: 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

Replying to 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:

5c9fa77Changing 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:

bbeb57cRemove 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: 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:

224db94Remove algorithm keyword

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

  • 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 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: 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: Changed 4 years ago by 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 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

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

9fc7b4dAdd 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:

fd1c03aMerge commit '9fc7b4dba9' into t/15404/public/15404
4f47b56Fix 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:

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

comment:42 Changed 4 years ago by roed

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.