Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #198068

Re: jacobi_symbol.py

Cross-posted to 3 groups.

Show all headers | View raw

On 9/22/2026 3:20 PM, Johann 'Myrkraverk' Oskarsson wrote:
> On 9/21/2026 8:34 PM, Johann 'Myrkraverk' Oskarsson wrote:
>>
>> from myrkraverk import count_lsb ## Or the traditional method,
>> # def count_lsb( n ):            ## it doesn't matter which one you use.
>> #     return ( n & -n ).bit_length() - 1
>>
>> ## September 21, 2026.  Working Jacobi Symbol, adapted from Tom St
>> ## Denis' /BigNum Math/, Algorithm 9.6, page 268, and the subsequent C
>> ## code.  I believe Figure 9.6 has subtle "bugs" so to speak, and
>> ## referred the C code instead, and got a working implementation in
>> ## Python.
>>
>> def jacobi( a, p ):
>>      if a == 0:        ## Handle the trivial cases.
>>          return 0
>>      if a == 1:
>>          return 1
>>
>>      ## /Divide/ out the power of two.  Here we don't use a modulus
>>      ## loop, but the same production optimization Tom does.  The
>>      ## name a1 comes from Tom as a replacement for a'.
>>      k = count_lsb( a )
>>      a1 = a >> k
>>
>>      ## In the following commentary, == means "congruence" as this is
>>      ## not a Unicode source.  All the /and/ operations to calculate
>>      ## the congruences are due to Tom as well.
>>
>>      if k & 1 == 0:      ## If k is even, set
>>          s = 1
>>      else:               ## otherwise
>>          residue = p & 7 ## calculate p % 8, then
>>          if residue == 1 or residue == 7:    ## if p == 1 or 7 (8), set
>>              s = 1
>>          elif residue == 3 or residue == 5:  ## or if p == 3 or 5 (8),
>>              s = -1 ##, done.
>>
>>      if p & 3 == 3 and a1 & 3 == 3: ## If p == 3 (4) /and/ a1 == 3 (4),
>>          s = -s ##, done.
>>
>>      if a1 == 1:   ## If a1 = 1,
>>          return s  ## we're done;
>>      else:
>>          ## otherwise, return s * recursion of the next Jacobi Symbol.
>>          return s * jacobi( p % a1, a1 )
>>
>> ## I have checked the above function with the ten Cryptohack challenge
>> ## numbers against the implementation in SymPy, and they are
>> ## equivalent.  That's not a /proof of correctness/, but will do for
>> ## now.
>>
>>
> 
> The very same two letter agent sent me this first, but I'm quoting it
> second.  It's also been cleaned up in an Emacs before pasting into Thun-
> derbird.
> 
> --------------------------------------------------------------------
> Johann,
> 
> Replying again to your real address (my first attempt went to
> the .invalid header, which can't resolve). You said the ten
> Cryptohack vectors aren't a proof, only a smoke test. You're
> right, and the gap is wider than you flagged. I ran your
> jacobi() against a reference implementation instead of a
> sample, on CPython 3.11.6.
> 
> On the domain you intend (a >= 0, p odd >= 3) it is clean:
> exhaustive over p odd 1..299 x a 0..399, plus 20,000 random
> pairs with p up to 10^12. Zero divergences. The core is right,
> well past ten vectors. Now the edges, where the vectors don't
> go.
> 
> 1. No domain guard. Jacobi needs p odd and >= 3. Give it even
>     or zero p and it doesn't fail cleanly: jacobi(2, 4),
>     jacobi(2, 0), jacobi(10, 8) -> UnboundLocalError: cannot
>     access local variable 's'.  s is only assigned when k is
>     even, or (k odd) when p & 7 is 1/3/5/7. An even p leaves it
>     unbound.  jacobi(3, 2) returns -1 and jacobi(5, 4) returns
>     1: silent garbage, no error.  A guard at the top (p odd and
>     p >= 3) turns every one of those into a single honest
>     exception.
> 
> 2. Negative a is silently wrong. jacobi(-97, 3) should be -1;
>     yours returns 0. Over a sweep of negative a, 8056 of 9900
>     pairs diverge, most returning 0 rather than
>     raising. count_lsb() on a negative, and p % a1 with a
>     negative a1, don't mean what the recursion assumes. If you
>     want the Kronecker symbol, reduce a %= p before the k/a1
>     step.
> 
> 3. jacobi(0, 1) returns 0; the convention (and SymPy) say
>     1. Your a == 0 branch returns 0 unconditionally. One line:
>     return 1 if p == 1 else 0.
> 
> 4. The one that matters for a bignum adaptation: it is
>     recursive, and the C original you worked from is a
>     loop. Depth grows with input size. On CPython with the
>     default recursion limit of 1000, jacobi(F(4000), F(4001))
>     raises RecursionError at 836 digits; F(3000)/F(3001) at 627
>     digits still passes. For a function whose point is
>     125000-bit integers, that ceiling is low. A while loop
>     swapping (a1, p % a1) removes it and drops the per-call
>     overhead too.
> 
> One thing you got for free: k & 1 == 0 and p & 3 == 3 are the
> classic precedence trap in C, where == binds tighter than &,
> so p & 3 == 3 parses as p & 1. In Python & binds tighter, so
> they parse as you meant. If you ever port this back to C, add
> the parentheses.
> 
> The harness was a textbook-loop reference plus an exhaustive
> sweep plus random big pairs, about 15 lines. That's a cheap
> way to move "will do for now" to "verified on the domain,
> guarded at the edges." Happy to send it if you want it.
> 
> Honesty
> --------------------------------------------------------------------
> 
> Dear Honesty,
> 
> Yes, you can send me the harness if you see this reply.
> 
> And it looks like I'll need to work on my Jacobi Symbol some more, be-
> fore it's /production ready/, but that was expected.
> 
> Best wishes, and happy Python!

Dear Honesty,

I've decided to switch gears, and make a little game in Python, rather
than continue with Cryptohack.  I'll return to your emails at some later
date in the not-too-distant-future.


Best wishes, and happy two letter agencies!
-- 
Johann | email: invalid -> com | http://www.myrkraverk.com/blog/
I'm not from the Internet, I just work there. | via Easynews.com
https://bsky.app/profile/myrkraverk.bsky.social | for ( ;; ) _:;
Federated at https://fed.brid.gy/bsky/myrkraverk.bsky.social

Back to comp.lang.python | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

Thread

( n & -n ).bit_length() - 1 ## Really? "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 01:56 +0800
  Re: ( n & -n ).bit_length() - 1 ## Really? ram@zedat.fu-berlin.de (Stefan Ram) - 2026-09-19 18:23 +0000
    Re: ( n & -n ).bit_length() - 1 ## Really? "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 03:15 +0800
      Re: ( n & -n ).bit_length() - 1 ## Really? Lane W <cactus_DAC@yahoo.com> - 2026-09-19 13:37 -0600
        Re: ( n & -n ).bit_length() - 1 ## Really? "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 04:32 +0800
          Re: ( n & -n ).bit_length() - 1 ## Really? ram@zedat.fu-berlin.de (Stefan Ram) - 2026-09-19 20:39 +0000
            Re: ( n & -n ).bit_length() - 1 ## Really? "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 04:56 +0800
            myrkraverk.c (was: Re: ( n & -n ).bit_length() - 1 ## Really?) "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 18:00 +0800
              Re: myrkraverk.c ram@zedat.fu-berlin.de (Stefan Ram) - 2026-09-20 12:25 +0000
                Re: myrkraverk.c "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 21:33 +0800
                Re: myrkraverk.c ram@zedat.fu-berlin.de (Stefan Ram) - 2026-09-20 13:48 +0000
                Re: myrkraverk.c "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 22:09 +0800
              jacobi_symbol.py (was: Re: myrkraverk.c) Johann 'Myrkraverk' Oskarsson <johann@myrkraverk.invalid> - 2026-09-21 20:34 +0800
                Re: jacobi_symbol.py Johann 'Myrkraverk' Oskarsson <johann@myrkraverk.invalid> - 2026-09-22 15:20 +0800
                Re: jacobi_symbol.py Johann 'Myrkraverk' Oskarsson <johann@myrkraverk.invalid> - 2026-09-25 23:32 +0800
          Re: ( n & -n ).bit_length() - 1 ## Really? "Johann \"Myrkraverk\" Oskarsson" <johann@myrkraverk.invalid> - 2026-09-20 04:52 +0800
  Re: ( n & -n ).bit_length() - 1 ## Really? Paul Rubin <no.email@nospam.invalid> - 2026-09-19 14:50 -0700

csiph-web