From 5176c62fb6e3f295598e71309a06874dece3c85f Mon Sep 17 00:00:00 2001 From: primo-ppcg Date: Thu, 20 Jul 2023 21:52:01 +0700 Subject: [PATCH 1/4] add functions for modular arithmetic --- project.janet | 4 ++ spork/math.janet | 3 +- src/cmath.c | 184 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 190 insertions(+), 1 deletion(-) create mode 100644 src/cmath.c diff --git a/project.janet b/project.janet index 0a22077..8283813 100644 --- a/project.janet +++ b/project.janet @@ -52,3 +52,7 @@ :source @["src/zip.c" "deps/miniz/miniz.c"] :defines @{"_LARGEFILE64_SOURCE" true} :headers @["deps/miniz/miniz.h"]) + +(declare-native + :name "spork/cmath" + :source @["src/cmath.c"]) diff --git a/spork/math.janet b/spork/math.janet index 03699cb..fb6d30b 100644 --- a/spork/math.janet +++ b/spork/math.janet @@ -1,4 +1,5 @@ (use ./misc) +(import ./cmath :prefix "" :export true) # Statistics @@ -309,7 +310,7 @@ (defn z-score ``` - Gets the standard score for number `x` from mean `m` + Gets the standard score for number `x` from mean `m` and standard deviation `d`. ``` [x m d] diff --git a/src/cmath.c b/src/cmath.c new file mode 100644 index 0000000..ad42f90 --- /dev/null +++ b/src/cmath.c @@ -0,0 +1,184 @@ +/* +* Copyright (c) 2023 Calvin Rose +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to +* deal in the Software without restriction, including without limitation the +* rights to use, copy, modify, merge, publish, distribute, sublicense, and/or +* sell copies of the Software, and to permit persons to whom the Software is +* furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included in +* all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +* IN THE SOFTWARE. +*/ + +#include + +int64_t _mod_impl(int64_t a, int64_t m) { + if (m == 0) return a; + + int64_t x = a % m; + return (((a ^ m) < 0) && (x != 0)) ? x + m : x; +} + +int64_t _jacobi_impl(int64_t a, int64_t m) { + int64_t res = 2; + a = _mod_impl(a, m); + while (a != 0) { + int64_t l = a & -a; + a /= l; + res ^= (a & m) ^ ((l % 3) & (m ^ (m >> 1))); + int64_t tmpa = a; + a = _mod_impl(m, a); + m = tmpa; + } + return m == 1 ? (res & 2) - 1 : 0; +} + +int64_t _invmod_impl(int64_t a, int64_t m) { + int64_t x = 0; + int64_t u = 1; + int64_t n = (m < 0) ? -m : m; + a = _mod_impl(a, n); + while (a != 0) { + int64_t tmpx = x; + x = u; + u = tmpx - (n / a) * u; + int64_t tmpa = a; + a = _mod_impl(n, a); + n = tmpa; + } + return _mod_impl(x, m); +} + +#if defined(__SIZEOF_INT128__) + +int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { + if (m == 0) return a * b; + + int64_t x = (((signed __int128)a) * b) % m; + return (((a ^ m) < 0) && (x != 0)) ? x + m : x; +} + +#elif defined(_WIN64) && defined(_MSC_VER) + +#include +int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { + if (m == 0) return a * b; + + int64_t r; + int64_t s = _mul128(a, b, &r); + (void)_div128(r, s, m, &r); + return r; +} + +#else + +int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { + int64_t res = 0; + a = _mod_impl(a, m); + b = _mod_impl(b, m); + if (a < b) { + int64_t tmpa = a; + a = b; + b = tmpa; + } + if (b < 0) b -= m; + while (b > 0) { + if ((b & 1) == 1) + res += (((m - res - a) ^ m) < 0) ? a - m : a; + a += (((m - a - a) ^ m) < 0) ? a - m : a; + b >>= 1; + } + return res; +} + +#endif + +int64_t _modpow_impl(int64_t a, int64_t b, int64_t m) { + int64_t res = 1; + if (b < 0) { + a = _invmod_impl(a, m); + b = -b; + } + while (b > 0) { + if ((b & 1) == 1) + res = _mulmod_impl(res, a, m); + a = _mulmod_impl(a, a, m); + b >>= 1; + } + return res; +} + +Janet wrap_result(int64_t a, Janet m) { + if (!janet_checktype(m, JANET_ABSTRACT)) + return janet_wrap_number(a); + const JanetAbstractType *at = janet_abstract_type(janet_unwrap_abstract(m)); + int64_t *box = janet_abstract(at, sizeof(int64_t)); + *box = a; + return janet_wrap_abstract(box); +} + +JANET_FN(cfun_cmath_jacobi, + "(math/jacobi a m)", + "Computes the Jacobi Symbol (a|m).") { + janet_fixarity(argc, 2); + int64_t a = janet_getinteger64(argv, 0); + int64_t m = janet_getinteger64(argv, 1); + + return janet_wrap_number(_jacobi_impl(a, m)); +} + +JANET_FN(cfun_cmath_invmod, + "(math/invmod a m)", + "Modular multiplicative inverse of `a` mod `m`. " + "Both arguments must be integer. The return value has the same type as `m`.") { + janet_fixarity(argc, 2); + int64_t a = janet_getinteger64(argv, 0); + int64_t m = janet_getinteger64(argv, 1); + + return wrap_result(_invmod_impl(a, m), argv[1]); +} + +JANET_FN(cfun_cmath_mulmod, + "(math/mulmod a b m)", + "Modular multiplication of `a` and `b` mod `m`. " + "All arguments must be integer. The return value has the same type as `m`.") { + janet_fixarity(argc, 3); + int64_t a = janet_getinteger64(argv, 0); + int64_t b = janet_getinteger64(argv, 1); + int64_t m = janet_getinteger64(argv, 2); + + return wrap_result(_mulmod_impl(a, b, m), argv[2]); +} + +JANET_FN(cfun_cmath_powmod, + "(math/powmod a b m)", + "Modular exponentiation of `a` to the power of `b` mod `m`. " + "All arguments must be integer. The return value has the same type as `m`.") { + janet_fixarity(argc, 3); + int64_t a = janet_getinteger64(argv, 0); + int64_t b = janet_getinteger64(argv, 1); + int64_t m = janet_getinteger64(argv, 2); + + return wrap_result(_modpow_impl(a, b, m), argv[2]); +} + +JANET_MODULE_ENTRY(JanetTable *env) { + JanetRegExt cfuns[] = { + JANET_REG("jacobi", cfun_cmath_jacobi), + JANET_REG("invmod", cfun_cmath_invmod), + JANET_REG("mulmod", cfun_cmath_mulmod), + JANET_REG("powmod", cfun_cmath_powmod), + JANET_REG_END + }; + janet_cfuns_ext(env, "cmath", cfuns); +} From 1b3a3c32737127ac403eec418d1a7c9d4237e0d3 Mon Sep 17 00:00:00 2001 From: primo-ppcg Date: Mon, 24 Jul 2023 17:46:27 +0700 Subject: [PATCH 2/4] add functions for primality testing --- spork/math.janet | 129 ++++++++++++++++++++++++++++++++++++++++++++++- src/cmath.c | 6 +-- 2 files changed, 131 insertions(+), 4 deletions(-) diff --git a/spork/math.janet b/spork/math.janet index fb6d30b..03f428b 100644 --- a/spork/math.janet +++ b/spork/math.janet @@ -1,5 +1,5 @@ (use ./misc) -(import ./cmath :prefix "" :export true) +(import spork/cmath :prefix "" :export true) # Statistics @@ -1169,3 +1169,130 @@ (for i 0 r (+= res (* (get-in m [0 i]) (perm (minor m 0 i)))))))) + +# Number Theory + +(def- small-primes + "primes less than 212" + [2 3 5 7 11 13 17 19 23 29 31 37 + 41 43 47 53 59 61 67 71 73 79 83 89 + 97 101 103 107 109 113 127 131 137 139 149 151 + 157 163 167 173 179 181 191 193 197 199 211]) + +(def- soe-indices + "pre-calced sieve of eratosthenes for 2 3 5 7" + [1 11 13 17 19 23 29 31 37 41 43 47 + 53 59 61 67 71 73 79 83 89 97 101 103 + 107 109 113 121 127 131 137 139 143 149 151 157 + 163 167 169 173 179 181 187 191 193 197 199 209]) + +(def- soe-offsets + "distances between sieve values" + [10 2 4 2 4 6 2 6 4 2 4 6 + 6 2 6 4 2 6 4 6 8 4 2 4 + 2 4 8 6 4 6 2 4 6 2 6 6 + 4 2 4 6 2 6 4 2 4 2 10 2]) + +(def- soe-table + "tabulated offsets mod 105" + [0 10 2 0 4 0 0 0 8 0 0 2 0 4 0 + 0 6 2 0 4 0 0 4 6 0 0 6 0 0 2 + 0 6 2 0 4 0 0 4 6 0 0 2 0 4 2 + 0 6 6 0 0 0 0 6 6 0 0 0 0 4 2 + 0 6 2 0 4 0 0 4 6 0 0 2 0 6 2 + 0 6 0 0 4 0 0 4 6 0 0 2 0 4 8 + 0 0 2 0 10 0 0 4 0 0 0 2 0 4 2]) + +(def- mr-bases-32 + "bases for miller-rabin determinism to 2^32" + [2 7 61]) + +(def- mr-bases-64 + "bases for miller-rabin determinism to 2^64" + [2 325 9375 28178 450775 9780504 1795265022]) + +(def- prime-prod + # (* ;(take 13 small-primes)) + 304250263527210) + +(defn miller-rabin-prp? + ``Performs a Miller-Rabin probable prime test on `n` against all given bases. + If no bases are given, sufficient bases are used to ensure that the check + is deterministic for all `n` less than 2^63.`` + [n & bases] + (def ps + (if (empty? bases) + (if (int? n) mr-bases-32 mr-bases-64) + bases)) + (var d (- n 1)) + (var s 0) + (while (zero? (mod d 2)) + (++ s) + (set d (/ d 2))) + (label result + (each p ps + (var x (powmod p d n)) + (when (compare< 1 x (- n 1)) + (for r 1 s + (set x (mulmod x x n)) + (if-not (compare< 1 x (- n 1)) (break))) + (if-not (compare= x (- n 1)) + (return result false)))) + true)) + +(defn prime? + ``A primality test, deterministic for all `n` less than 2^63.`` + [n] + (cond + (compare<= n 211) (do + (def m (binary-search n small-primes compare<)) + (compare= n (in small-primes m))) + (= 0 (jacobi n prime-prod)) false + (miller-rabin-prp? n))) + +(defn next-prime + "Returns the next prime number strictly larger than `n`." + [n] + (label result + (def x (+ n 1)) + (if (compare<= x 2) (return result 2)) + (when (compare<= x 211) + (def m (binary-search x small-primes compare<)) + (return result (in small-primes m))) + (def j (mod x 210)) + (def m (binary-search j soe-indices compare<)) + (var i (+ x (- (in soe-indices m) j))) + (def offs [;(slice soe-offsets m) ;(slice soe-offsets 0 m)]) + (forever + (each o offs + (if (and (not= 0 (jacobi i prime-prod)) (miller-rabin-prp? i)) + (return result i)) + (+= i o))))) + +(defn primes + "A boundless prime number generator." + [] + (coro + (each n small-primes (yield n)) + (def sieve @{221 13 253 11}) + (def pg (drop 6 (primes))) + (var p (resume pg)) + (var q (* p p)) + (var n 211) + (forever + (each offset soe-offsets + (+= n offset) + (def stp (get sieve n)) + (cond + stp (do + (put sieve n nil) + (var nxt (/ n stp)) + (+= nxt (in soe-table (mod nxt 105))) + (while (get sieve (* nxt stp)) + (+= nxt (in soe-table (mod nxt 105)))) + (put sieve (* nxt stp) stp)) + (< n q) (yield n) + (do + (put sieve (+ q (* p (in soe-table (mod p 105)))) p) + (set p (resume pg)) + (set q (* p p)))))))) diff --git a/src/cmath.c b/src/cmath.c index ad42f90..96beb40 100644 --- a/src/cmath.c +++ b/src/cmath.c @@ -65,7 +65,7 @@ int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { if (m == 0) return a * b; int64_t x = (((signed __int128)a) * b) % m; - return (((a ^ m) < 0) && (x != 0)) ? x + m : x; + return (((a ^ b ^ m) < 0) && (x != 0)) ? x + m : x; } #elif defined(_WIN64) && defined(_MSC_VER) @@ -73,11 +73,11 @@ int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { #include int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { if (m == 0) return a * b; - + int64_t r; int64_t s = _mul128(a, b, &r); (void)_div128(r, s, m, &r); - return r; + return (((a ^ b ^ m) < 0) && (r != 0)) ? r + m : r; } #else From b828418e642d609ff0b765f9cb7d5c17a8e0653d Mon Sep 17 00:00:00 2001 From: primo-ppcg Date: Thu, 27 Jul 2023 17:11:58 +0700 Subject: [PATCH 3/4] cleanup, add test cases Do not export `miller-rabin-prp?`. --- spork/math.janet | 39 +++++----- test/suite-math.janet | 176 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 195 insertions(+), 20 deletions(-) diff --git a/spork/math.janet b/spork/math.janet index 03f428b..9442121 100644 --- a/spork/math.janet +++ b/spork/math.janet @@ -1212,10 +1212,10 @@ [2 325 9375 28178 450775 9780504 1795265022]) (def- prime-prod - # (* ;(take 13 small-primes)) - 304250263527210) + # (* ;(tuple/slice small-primes 1 14)) + 6541380665835015) -(defn miller-rabin-prp? +(defn- miller-rabin-prp? ``Performs a Miller-Rabin probable prime test on `n` against all given bases. If no bases are given, sufficient bases are used to ensure that the check is deterministic for all `n` less than 2^63.`` @@ -1228,7 +1228,7 @@ (var s 0) (while (zero? (mod d 2)) (++ s) - (set d (/ d 2))) + (/= d 2)) (label result (each p ps (var x (powmod p d n)) @@ -1241,33 +1241,32 @@ true)) (defn prime? - ``A primality test, deterministic for all `n` less than 2^63.`` + "A primality test, deterministic for all `n` less than 2^63." [n] (cond (compare<= n 211) (do (def m (binary-search n small-primes compare<)) (compare= n (in small-primes m))) + (zero? (mod n 2)) false (= 0 (jacobi n prime-prod)) false (miller-rabin-prp? n))) (defn next-prime "Returns the next prime number strictly larger than `n`." [n] - (label result - (def x (+ n 1)) - (if (compare<= x 2) (return result 2)) - (when (compare<= x 211) - (def m (binary-search x small-primes compare<)) - (return result (in small-primes m))) - (def j (mod x 210)) - (def m (binary-search j soe-indices compare<)) - (var i (+ x (- (in soe-indices m) j))) - (def offs [;(slice soe-offsets m) ;(slice soe-offsets 0 m)]) - (forever - (each o offs - (if (and (not= 0 (jacobi i prime-prod)) (miller-rabin-prp? i)) - (return result i)) - (+= i o))))) + (var x (+ n 1)) + (if (compare<= x 211) + (in small-primes (binary-search x small-primes compare<)) + (label result + (def i (mod x 210)) + (def m (binary-search i soe-indices compare<)) + (+= x (- (in soe-indices m) i)) + (def offs [;(tuple/slice soe-offsets m) ;(tuple/slice soe-offsets 0 m)]) + (forever + (each o offs + (if (and (not= 0 (jacobi x prime-prod)) (miller-rabin-prp? x)) + (return result x)) + (+= x o)))))) (defn primes "A boundless prime number generator." diff --git a/test/suite-math.janet b/test/suite-math.janet index c87bfce..ecd5876 100644 --- a/test/suite-math.janet +++ b/test/suite-math.janet @@ -421,4 +421,180 @@ @[-3 3 -1]])) "permanent") +(def test-primes + [100123456789 + 107928278317 + 113131311401 + 125411328001 + 137438691329 + 150614187107 + 168409140869 + 191198888863 + 232911191713 + 260389232731 + 313473008141 + 344980016453 + 411379717319 + 526858348381 + 581485876661 + 657835997711 + 701234567897 + 777775777777 + 977973373171 + 1000008000001 + 1110110110111 + 1146890986411 + 1344409044431 + 1534139560403 + 1686910196861 + 1918960000861 + 2748779069441 + 3484606209571 + 4615392897979 + 5599297703999 + 6987191424553 + 7777775552353 + 9977770001777 + 10861196119801 + 12345678987431 + 15154262241479 + 18826507658281 + 24738041398529 + 29202114663949 + 35557777777987 + 50000010000001 + 60818091990661 + 74596893730427 + 88929267773197 + 100033000330001 + 106111115118119 + 123123454321321 + 139239739439839 + 176860696068671 + 233444000011111 + 303160419086407 + 333555577577777 + 485398038695407 + 727777887889889 + 799333555511111 + 973369606963379 + 1003229774283941 + 1235711131175321 + 1510553619999637 + 1889080110806881 + 2357353373727757 + 3391382115599173 + 5111111111111119 + 7005264275346131 + 9007199254740847]) + +(def pseudoprimes + [1373653 + 1530787 + 1987021 + 2284453 + 3116107 + 5173601 + 6787327 + 11541307 + 13694761 + 15978007 + 16070429 + 16879501 + 25326001 + 27509653 + 27664033 + 28527049 + 54029741 + 61832377 + 66096253 + 74927161 + 80375707 + 101649241 + 161304001 + 960946321 + 1157839381 + 3215031751 + 3697278427 + 5764643587 + 6770862367 + 14386156093 + 15579919981 + 18459366157 + 19887974881 + 21276028621 + 27716349961 + 29118033181 + 37131467521 + 41752650241 + 42550716781 + 43536545821 + 118670087467 + 307768373641 + 315962312077 + 354864744877 + 457453568161 + 528929554561 + 546348519181 + 602248359169 + 1362242655901 + 1871186716981 + 2152302898747 + 2273312197621 + 2366338900801 + 3343433905957 + 3461715915661 + 3474749660383 + 3477707481751 + 4341937413061 + 4777422165601 + 5537838510751]) + +(defn trial-division + [n] + (cond + (< n 2) false + (= n 2) true + (= 0 (mod n 2)) false + (label result + (var p 3) + (while (<= (* p p) n) + (if (= 0 (mod n p)) (return result false)) + (+= p 2)) + true))) + +(def- prime-prod + # product of primes from 3 to 43 + 6541380665835015) + +(def pg (take 10000 (primes))) +(assert (all prime? pg) "small primes") +(assert (= 104729 (last pg)) "10000th prime") + +(assert (all = (map next-prime pg) (drop 1 pg)) "next-prime") + +(assert + (all |(= (trial-division $) (prime? $)) (range 1 100000)) + "prime? vs trial division") + +(assert + (all |(= (trial-division $) (not= 0 (jacobi $ prime-prod))) (range 49 2209 2)) + "jacobi vs trial division") + +(assert (all prime? test-primes) "test primes") +(when-let [s64 int/s64 + u64 int/u64] + (assert (all prime? (map s64 test-primes)) "test primes s64") + (assert (all prime? (map u64 test-primes)) "test primes u64")) + +(assert (not (some prime? pseudoprimes)) "pseudoprimes") + +(each p test-primes + (assert + (all |(= 1 (mulmod $ (invmod $ p) p)) (range 1 1000)) + (string "invmod " p)) + (assert + (all |(= 1 (mulmod (powmod $ $ p) (powmod $ (- $) p) p)) (range 1 1000)) + (string "powmod " p))) + (end-suite) From 2587b2f9a0a1de54f3d2368ac90b5c98f8919449 Mon Sep 17 00:00:00 2001 From: primo-ppcg Date: Wed, 9 Aug 2023 18:20:40 +0700 Subject: [PATCH 4/4] return `math/nan` if no modular inverse exists --- spork/math.janet | 6 +++--- src/cmath.c | 49 +++++++++++++++++++++++++++---------------- test/suite-math.janet | 2 +- 3 files changed, 35 insertions(+), 22 deletions(-) diff --git a/spork/math.janet b/spork/math.janet index 9442121..2e98c79 100644 --- a/spork/math.janet +++ b/spork/math.janet @@ -1212,7 +1212,7 @@ [2 325 9375 28178 450775 9780504 1795265022]) (def- prime-prod - # (* ;(tuple/slice small-primes 1 14)) + "product of primes from 3 to 43" 6541380665835015) (defn- miller-rabin-prp? @@ -1264,8 +1264,8 @@ (def offs [;(tuple/slice soe-offsets m) ;(tuple/slice soe-offsets 0 m)]) (forever (each o offs - (if (and (not= 0 (jacobi x prime-prod)) (miller-rabin-prp? x)) - (return result x)) + (if (not= 0 (jacobi x prime-prod)) + (if (miller-rabin-prp? x) (return result x))) (+= x o)))))) (defn primes diff --git a/src/cmath.c b/src/cmath.c index 96beb40..31c774c 100644 --- a/src/cmath.c +++ b/src/cmath.c @@ -40,7 +40,7 @@ int64_t _jacobi_impl(int64_t a, int64_t m) { a = _mod_impl(m, a); m = tmpa; } - return m == 1 ? (res & 2) - 1 : 0; + return (m == 1) ? (res & 2) - 1 : 0; } int64_t _invmod_impl(int64_t a, int64_t m) { @@ -56,7 +56,7 @@ int64_t _invmod_impl(int64_t a, int64_t m) { a = _mod_impl(n, a); n = tmpa; } - return _mod_impl(x, m); + return (n == 1) ? _mod_impl(x, m) : 0; } #if defined(__SIZEOF_INT128__) @@ -103,24 +103,18 @@ int64_t _mulmod_impl(int64_t a, int64_t b, int64_t m) { #endif -int64_t _modpow_impl(int64_t a, int64_t b, int64_t m) { - int64_t res = 1; - if (b < 0) { - a = _invmod_impl(a, m); - b = -b; - } - while (b > 0) { - if ((b & 1) == 1) - res = _mulmod_impl(res, a, m); - a = _mulmod_impl(a, a, m); - b >>= 1; - } - return res; +Janet wrap_nan() { +#ifdef NAN + return janet_wrap_number(NAN); +#else + return janet_wrap_number(0.0 / 0.0); +#endif } Janet wrap_result(int64_t a, Janet m) { if (!janet_checktype(m, JANET_ABSTRACT)) return janet_wrap_number(a); + const JanetAbstractType *at = janet_abstract_type(janet_unwrap_abstract(m)); int64_t *box = janet_abstract(at, sizeof(int64_t)); *box = a; @@ -140,12 +134,17 @@ JANET_FN(cfun_cmath_jacobi, JANET_FN(cfun_cmath_invmod, "(math/invmod a m)", "Modular multiplicative inverse of `a` mod `m`. " - "Both arguments must be integer. The return value has the same type as `m`.") { + "Both arguments must be integer. The return value has the same type as `m`. " + "If no inverse exists, returns `math/nan` instead.") { janet_fixarity(argc, 2); int64_t a = janet_getinteger64(argv, 0); int64_t m = janet_getinteger64(argv, 1); - return wrap_result(_invmod_impl(a, m), argv[1]); + int64_t res = _invmod_impl(a, m); + if (res == 0) + return wrap_nan(); + + return wrap_result(res, argv[1]); } JANET_FN(cfun_cmath_mulmod, @@ -169,7 +168,21 @@ JANET_FN(cfun_cmath_powmod, int64_t b = janet_getinteger64(argv, 1); int64_t m = janet_getinteger64(argv, 2); - return wrap_result(_modpow_impl(a, b, m), argv[2]); + int64_t res = 1; + if (b < 0) { + a = _invmod_impl(a, m); + if (a == 0) + return wrap_nan(); + b = -b; + } + while (b > 0) { + if ((b & 1) == 1) + res = _mulmod_impl(res, a, m); + a = _mulmod_impl(a, a, m); + b >>= 1; + } + + return wrap_result(res, argv[2]); } JANET_MODULE_ENTRY(JanetTable *env) { diff --git a/test/suite-math.janet b/test/suite-math.janet index ecd5876..d4e4e2f 100644 --- a/test/suite-math.janet +++ b/test/suite-math.janet @@ -564,7 +564,7 @@ true))) (def- prime-prod - # product of primes from 3 to 43 + "product of primes from 3 to 43" 6541380665835015) (def pg (take 10000 (primes)))