• jacobi_symbol.py (was: Re: myrkraverk.c)

    From Johann 'Myrkraverk' Oskarsson@johann@myrkraverk.invalid to comp.lang.python,sci.math,sci.crypt on Mon Sep 21 20:34:52 2026
    From Newsgroup: sci.crypt


    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.
    --
    Johann | email: invalid -> com | http://www.myrkraverk.com/blog/
    I'm not from the Internet, I just work there. | via XS News https://bsky.app/profile/myrkraverk.bsky.social | for ( ;; ) _:;
    --- Synchronet 3.22a-Linux NewsLink 1.2
  • From Johann 'Myrkraverk' Oskarsson@johann@myrkraverk.invalid to comp.lang.python,sci.math,sci.crypt on Tue Sep 22 15:20:46 2026
    From Newsgroup: sci.crypt

    On 9/21/2026 8:34 PM, Johann 'Myrkraverk' Oskarsson wrote:

    from myrkraverk import count_lsb ## Or the traditional method,
    # def count_lsb( n ):-a-a-a-a-a-a-a-a-a-a-a ## it doesn't matter which one you use.
    #-a-a-a-a return ( n & -n ).bit_length() - 1

    ## September 21, 2026.-a Working Jacobi Symbol, adapted from Tom St
    ## Denis' /BigNum Math/, Algorithm 9.6, page 268, and the subsequent C
    ## code.-a 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 ):
    -a-a-a if a == 0:-a-a-a-a-a-a-a ## Handle the trivial cases.
    -a-a-a-a-a-a-a return 0
    -a-a-a if a == 1:
    -a-a-a-a-a-a-a return 1

    -a-a-a ## /Divide/ out the power of two.-a Here we don't use a modulus
    -a-a-a ## loop, but the same production optimization Tom does.-a The
    -a-a-a ## name a1 comes from Tom as a replacement for a'.
    -a-a-a k = count_lsb( a )
    -a-a-a a1 = a >> k

    -a-a-a ## In the following commentary, == means "congruence" as this is
    -a-a-a ## not a Unicode source.-a All the /and/ operations to calculate
    -a-a-a ## the congruences are due to Tom as well.

    -a-a-a if k & 1 == 0:-a-a-a-a-a ## If k is even, set
    -a-a-a-a-a-a-a s = 1
    -a-a-a else:-a-a-a-a-a-a-a-a-a-a-a-a-a-a ## otherwise
    -a-a-a-a-a-a-a residue = p & 7 ## calculate p % 8, then
    -a-a-a-a-a-a-a if residue == 1 or residue == 7:-a-a-a ## if p == 1 or 7 (8), set
    -a-a-a-a-a-a-a-a-a-a-a s = 1
    -a-a-a-a-a-a-a elif residue == 3 or residue == 5:-a ## or if p == 3 or 5 (8),
    -a-a-a-a-a-a-a-a-a-a-a s = -1 ##, done.

    -a-a-a if p & 3 == 3 and a1 & 3 == 3: ## If p == 3 (4) /and/ a1 == 3 (4),
    -a-a-a-a-a-a-a s = -s ##, done.

    -a-a-a if a1 == 1:-a-a ## If a1 = 1,
    -a-a-a-a-a-a-a return s-a ## we're done;
    -a-a-a else:
    -a-a-a-a-a-a-a ## otherwise, return s * recursion of the next Jacobi Symbol.
    -a-a-a-a-a-a-a 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.-a 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!
    --
    Johann | email: invalid -> com | http://www.myrkraverk.com/blog/
    I'm not from the Internet, I just work there. | via XS News https://bsky.app/profile/myrkraverk.bsky.social | for ( ;; ) _:;
    --- Synchronet 3.22a-Linux NewsLink 1.2
  • From Johann 'Myrkraverk' Oskarsson@johann@myrkraverk.invalid to comp.lang.python,sci.math,sci.crypt on Fri Sep 25 23:32:42 2026
    From Newsgroup: sci.crypt

    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 ):-a-a-a-a-a-a-a-a-a-a-a ## it doesn't matter which one you use.
    #-a-a-a-a return ( n & -n ).bit_length() - 1

    ## September 21, 2026.-a Working Jacobi Symbol, adapted from Tom St
    ## Denis' /BigNum Math/, Algorithm 9.6, page 268, and the subsequent C
    ## code.-a 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 ):
    -a-a-a-a if a == 0:-a-a-a-a-a-a-a ## Handle the trivial cases.
    -a-a-a-a-a-a-a-a return 0
    -a-a-a-a if a == 1:
    -a-a-a-a-a-a-a-a return 1

    -a-a-a-a ## /Divide/ out the power of two.-a Here we don't use a modulus
    -a-a-a-a ## loop, but the same production optimization Tom does.-a The
    -a-a-a-a ## name a1 comes from Tom as a replacement for a'.
    -a-a-a-a k = count_lsb( a )
    -a-a-a-a a1 = a >> k

    -a-a-a-a ## In the following commentary, == means "congruence" as this is
    -a-a-a-a ## not a Unicode source.-a All the /and/ operations to calculate
    -a-a-a-a ## the congruences are due to Tom as well.

    -a-a-a-a if k & 1 == 0:-a-a-a-a-a ## If k is even, set
    -a-a-a-a-a-a-a-a s = 1
    -a-a-a-a else:-a-a-a-a-a-a-a-a-a-a-a-a-a-a ## otherwise
    -a-a-a-a-a-a-a-a residue = p & 7 ## calculate p % 8, then
    -a-a-a-a-a-a-a-a if residue == 1 or residue == 7:-a-a-a ## if p == 1 or 7 (8), set
    -a-a-a-a-a-a-a-a-a-a-a-a s = 1
    -a-a-a-a-a-a-a-a elif residue == 3 or residue == 5:-a ## or if p == 3 or 5 (8),
    -a-a-a-a-a-a-a-a-a-a-a-a s = -1 ##, done.

    -a-a-a-a if p & 3 == 3 and a1 & 3 == 3: ## If p == 3 (4) /and/ a1 == 3 (4), >> -a-a-a-a-a-a-a-a s = -s ##, done.

    -a-a-a-a if a1 == 1:-a-a ## If a1 = 1,
    -a-a-a-a-a-a-a-a return s-a ## we're done;
    -a-a-a-a else:
    -a-a-a-a-a-a-a-a ## otherwise, return s * recursion of the next Jacobi Symbol.
    -a-a-a-a-a-a-a-a 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.-a 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.-a 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
    -a-a or zero p and it doesn't fail cleanly: jacobi(2, 4),
    -a-a jacobi(2, 0), jacobi(10, 8) -> UnboundLocalError: cannot
    -a-a access local variable 's'.-a s is only assigned when k is
    -a-a even, or (k odd) when p & 7 is 1/3/5/7. An even p leaves it
    -a-a unbound.-a jacobi(3, 2) returns -1 and jacobi(5, 4) returns
    -a-a 1: silent garbage, no error.-a A guard at the top (p odd and
    -a-a p >= 3) turns every one of those into a single honest
    -a-a exception.

    2. Negative a is silently wrong. jacobi(-97, 3) should be -1;
    -a-a yours returns 0. Over a sweep of negative a, 8056 of 9900
    -a-a pairs diverge, most returning 0 rather than
    -a-a raising. count_lsb() on a negative, and p % a1 with a
    -a-a negative a1, don't mean what the recursion assumes. If you
    -a-a want the Kronecker symbol, reduce a %= p before the k/a1
    -a-a step.

    3. jacobi(0, 1) returns 0; the convention (and SymPy) say
    -a-a 1. Your a == 0 branch returns 0 unconditionally. One line:
    -a-a return 1 if p == 1 else 0.

    4. The one that matters for a bignum adaptation: it is
    -a-a recursive, and the C original you worked from is a
    -a-a loop. Depth grows with input size. On CPython with the
    -a-a default recursion limit of 1000, jacobi(F(4000), F(4001))
    -a-a raises RecursionError at 836 digits; F(3000)/F(3001) at 627
    -a-a digits still passes. For a function whose point is
    -a-a 125000-bit integers, that ceiling is low. A while loop
    -a-a swapping (a1, p % a1) removes it and drops the per-call
    -a-a 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
    --- Synchronet 3.22a-Linux NewsLink 1.2