Legendre’s Symbol

May 1, 2012

Here are our new versions of the jacobi and legendre symbols, which combine the rules given above into a single function:

(define (jacobi a m)
  (if (not (integer? a)) (error 'jacobi "must be integer")
    (if (not (and (integer? m) (positive? m) (odd? m)))
        (error 'jacobi "modulus must be odd positive integer")
        (let loop1 ((a (modulo a m)) (m m) (t 1))
          (if (zero? a) (if (= m 1) t 0)
            (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
              (let loop2 ((a a) (t t))
                (if (even? a) (loop2 (/ a 2) (* t z))
                  (loop1 (modulo m a) a
                         (if (and (= (modulo a 4) 3)
                                  (= (modulo m 4) 3))
                             (- t) t))))))))))

Although we’ve rearranged things, this code faithfully executes the six rules; it is interesting to trace from the rules to the code. The time complexity is about the same as taking a gcd, so it’s very fast. The legendre symbol is shown below:

(define (legendre a m)
  (if (prime? m) (jacobi a m)
    (error 'legendre "modulus must be prime")))

And here are some tests; the first two failed on the old version of jacobi:

> (jacobi -1 7)
-1
> (jacobi 3 3)
0
> (jacobi 127 703)
-1

You can run the program at http://programmingpraxis.codepad.org/6uVnD0RR.

About these ads

Pages: 1 2

13 Responses to “Legendre’s Symbol”

  1. ardnew said

    I first tried translating this to a completely iterative algorithm, but it was becoming a mess.. I also couldn’t think of a better way to reduce the rule 3 logic.

    use strict;
    use warnings;
    
    sub reduce
    {
      my ($p, $a, $m) = (1, @_);
    
      if ($a > 2)
      {       
        # rule 4
        $p *=  reduce(2, $m) * reduce($a / 2, $m) unless $a & 1;
        
        # rule 6
        $p *= -reduce(-$m, $a) if $a % 4 == 3 and $m % 4 == 3;
      }
    
      # rules 1 and 2
      $p *= $a unless $a >> 1;
    
      # rule 3
      if (2 == $a)
      {
        my $mod8 = $m % 8;
        $p *= -1 if 3 == $mod8 or 5 == $mod8;
        $p *=  1 if 1 == $mod8 or 7 == $mod8;
      }
    
      return $p;
    }
    
    sub legendre($)
    {
      my ($a, $m) = @{(shift)};
      
      die "modulus must be prime positive integer$/"
        unless int($m) == $m and $m > 0 and $m & 1; 
        # need to also verify $m is prime, but meh
      
      # rule 5
      my $n = $a % $m;
      
      print "($a / $m) = ".reduce($n, $m).$/;
    }
    
    legendre(\@ARGV);
    
  2. programmingpraxis said

    Ardnew: the usual iterative solution looks like this:

    a = a mod m;
    t = 1;
    while (a != 0) {
        while (a even) {
            a = a / 2;
            if (m mod 8 in {3,5}) t = -t;
        }
        (a,m) = (m,a) // swap variables
        if (a === m === 3 (mod 4)) t = -t
        a = a mod m;
    }
    if (m == 1) return t;
    return 0;

  3. ardnew said

    very nice! always seems so simple once you see the answer…

    also, there’s a bug in my code that doesn’t reapply rule 5 when rule 6 is evaluated. I’m not familiar enough with the math to prove this behavior is incorrect though

  4. programmingpraxis said

    In a private email, Mike pointed out that I missed the else clause in Rule 6. It is now fixed.

  5. Mike said

    An iterative Python version:

    def jacobi(n, m):
        j = 1
    
        # rule 5
        n %= m
        
        while n:
            # rules 3 and 4
            t = 0
            while not n & 1:
                n /= 2
                t += 1
            if t&1 and m%8 in (3, 5):
                j = -j
    
            # rule 6
            if (n % 4 == m % 4 == 3):
                j = -j
    
            # rules 5 and 6
            n, m = m % n, n
    
        return j if m == 1 else 0
    
    
  6. ardnew said

    Mike, a couple comments:

    a) I love that 3-way equality condition in rule 6. Never seen that before :)
    b) Can you explain the behavior of the expression: “if t&1 and m%8 in (3, 5)”

    And Mr. programmingpraxis (i assume that is your real name), a question:
    Is rule 5 even a necessary reduction? I’m beginning to think its there merely for computational optimization

  7. ardnew said

    Oops, I was caught up in that equality expression, disregard comment b), Mike.

  8. Mike said

    @ardnew,

    The ability to chain relational operators is a nice feature of Python syntax. It compiles to a series of 2-way relational tests ‘and’ed together. “if n%4 == m%4 == 3:” compiles like “if n%4 == m%4 and m%4 == 3:”, but the m%4 is only evaluated once. It is handy when you want to check if a value is between two limits, like “if lowerbound <= x = y <= z:" in a previous exercise to determine that y was no larger than either x or z. You can also chain any number together: "a c == d …”

    Rule 5 is the only way to reduce ‘a’ when it’s not even.

  9. Mike said

    Oops, wordpress ate some of my relational operators as html tags. Should be:

    The ability to chain relational operators is a nice feature of Python syntax. It compiles to a series of 2-way relational tests ‘and’ed together. “if n%4 == m%4 == 3:” compiles like “if n%4 == m%4 and m%4 == 3:”, but the m%4 is only evaluated once. It is handy when you want to check if a value is between two limits, like “if lowerbound <= x < upperbound:". I used something like "if x >= y <= z:" in a previous exercise to determine that y was no larger than either x or z. You can also chain any number together: "a < b > c == d …”

    Rule 5 is the only way to reduce ‘a’ when it’s not even.

  10. Graham said

    Since Mike already did an iterative Python version, here’s a straight translation of the rules to Haskell:

    jacobi :: Integer -> Integer -> Integer
    jacobi a m
        | even m        = error "m must be odd."
        | a == 0        = 0
        | a == 1        = 1
        | a == 2        = if (m `mod` 8) `elem` [3, 5] then -1 else 1
        | even a        = jacobi 2 m * jacobi (a `div` 2) m
        | a >= m        = jacobi (a `mod` m) m
        | otherwise     = let j = jacobi m a in
                            if (a `mod` 4 == 3) && (m `mod` 4 == 3) then -j else j
    

    Not the most efficient solution, probably. Happy to see another math exercise, though! (even though it should be titled “Jacobi’s Symbol,” since we don’t know m is prime)

  11. Mike said

    As a preprocessing step, I think rule 5 needs to be applied if a (7 / -1) by rule 6 (the only rule that could apply
    (7 / -1) => (0 / -1) by rule 5
    (0 / -1) => 0 by rule 1

    However, (-1 / 7) = -1 according to the programmingpraxis solution.

    If rule 5 is applied first, then:
    (-1 / 7) => (6 / 7) by rule 5
    (6 / 7) => (2 / 7) * (3 / 7) => (3 / 7) => -(7 / 3) => -(1 / 3) => -1 by rules 4, 3, 6, 5, and 2, respectively

  12. Mike said

    I did it again. (Maybe I need an option to preview a post)

    The previous post should be:

    As a preprocessing step, I think rule 5 needs to be applied if a < 0.

    For example:
    (-1 / 7) => (7 / -1) by rule 6 (the only rule that could apply)
    (7 / -1) => (0 / -1) by rule 5
    (0 / -1) => 0 by rule 1

    However, (-1 / 7) = -1 according to the programmingpraxis solution.

    If rule 5 is applied first, then:
    (-1 / 7) => (6 / 7) by rule 5
    (6 / 7) => (2 / 7) * (3 / 7) => (3 / 7) => -(7 / 3) => -(1 / 3) => -1; by rules 4, 3, 6, 5, and 2, respectively.

  13. programmingpraxis said

    The jacobi symbol for (-1 7) is -1. Wolfram Alpha agrees.

    In arithmetic mod m, the only numbers allowed are 0 to m-1 inclusive. There are no negative numbers. So any time you see a negative number you have to add m repeatedly until it is in the range 0 to m-1 inclusive. So yes, Rule 5 applies whenever a<0. I'll edit the rule.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 632 other followers

%d bloggers like this: