Path: blob/main/Modules/_decimal/libmpdec/literature/umodarith.lisp
12 views
;1; Copyright (c) 2008-2020 Stefan Krah. All rights reserved.2;3; Redistribution and use in source and binary forms, with or without4; modification, are permitted provided that the following conditions5; are met:6;7; 1. Redistributions of source code must retain the above copyright8; notice, this list of conditions and the following disclaimer.9;10; 2. Redistributions in binary form must reproduce the above copyright11; notice, this list of conditions and the following disclaimer in the12; documentation and/or other materials provided with the distribution.13;14; THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND15; ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE16; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE17; ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE18; FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL19; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS20; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)21; HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT22; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY23; OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF24; SUCH DAMAGE.25;262728(in-package "ACL2")2930(include-book "arithmetic/top-with-meta" :dir :system)31(include-book "arithmetic-2/floor-mod/floor-mod" :dir :system)323334;; =====================================================================35;; Proofs for several functions in umodarith.h36;; =====================================================================37383940;; =====================================================================41;; Helper theorems42;; =====================================================================4344(defthm elim-mod-m<x<2*m45(implies (and (<= m x)46(< x (* 2 m))47(rationalp x) (rationalp m))48(equal (mod x m)49(+ x (- m)))))5051(defthm modaux-1a52(implies (and (< x m) (< 0 x) (< 0 m)53(rationalp x) (rationalp m))54(equal (mod (- x) m)55(+ (- x) m))))5657(defthm modaux-1b58(implies (and (< (- x) m) (< x 0) (< 0 m)59(rationalp x) (rationalp m))60(equal (mod x m)61(+ x m)))62:hints (("Goal" :use ((:instance modaux-1a63(x (- x)))))))6465(defthm modaux-1c66(implies (and (< x m) (< 0 x) (< 0 m)67(rationalp x) (rationalp m))68(equal (mod x m)69x)))7071(defthm modaux-2a72(implies (and (< 0 b) (< b m)73(natp x) (natp b) (natp m)74(< (mod (+ b x) m) b))75(equal (mod (+ (- m) b x) m)76(+ (- m) b (mod x m)))))7778(defthm modaux-2b79(implies (and (< 0 b) (< b m)80(natp x) (natp b) (natp m)81(< (mod (+ b x) m) b))82(equal (mod (+ b x) m)83(+ (- m) b (mod x m))))84:hints (("Goal" :use (modaux-2a))))8586(defthm linear-mod-187(implies (and (< x m) (< b m)88(natp x) (natp b)89(rationalp m))90(equal (< x (mod (+ (- b) x) m))91(< x b)))92:hints (("Goal" :use ((:instance modaux-1a93(x (+ b (- x))))))))9495(defthm linear-mod-296(implies (and (< 0 b) (< b m)97(natp x) (natp b)98(natp m))99(equal (< (mod x m)100(mod (+ (- b) x) m))101(< (mod x m) b))))102103(defthm linear-mod-3104(implies (and (< x m) (< b m)105(natp x) (natp b)106(rationalp m))107(equal (<= b (mod (+ b x) m))108(< (+ b x) m)))109:hints (("Goal" :use ((:instance elim-mod-m<x<2*m110(x (+ b x)))))))111112(defthm modaux-2c113(implies (and (< 0 b) (< b m)114(natp x) (natp b) (natp m)115(<= b (mod (+ b x) m)))116(equal (mod (+ b x) m)117(+ b (mod x m))))118:hints (("Subgoal *1/8''" :use (linear-mod-3))))119120(defthmd modaux-2d121(implies (and (< x m) (< 0 x) (< 0 m)122(< (- m) b) (< b 0) (rationalp m)123(<= x (mod (+ b x) m))124(rationalp x) (rationalp b))125(equal (+ (- m) (mod (+ b x) m))126(+ b x)))127:hints (("Goal" :cases ((<= 0 (+ b x))))128("Subgoal 2'" :use ((:instance modaux-1b129(x (+ b x)))))))130131(defthm mod-m-b132(implies (and (< 0 x) (< 0 b) (< 0 m)133(< x b) (< b m)134(natp x) (natp b) (natp m))135(equal (mod (+ (mod (- x) m) b) m)136(mod (- x) b))))137138139;; =====================================================================140;; addmod, submod141;; =====================================================================142143(defun addmod (a b m base)144(let* ((s (mod (+ a b) base))145(s (if (< s a) (mod (- s m) base) s))146(s (if (>= s m) (mod (- s m) base) s)))147s))148149(defthmd addmod-correct150(implies (and (< 0 m) (< m base)151(< a m) (<= b m)152(natp m) (natp base)153(natp a) (natp b))154(equal (addmod a b m base)155(mod (+ a b) m)))156:hints (("Goal" :cases ((<= base (+ a b))))157("Subgoal 2.1'" :use ((:instance elim-mod-m<x<2*m158(x (+ a b)))))))159160(defun submod (a b m base)161(let* ((d (mod (- a b) base))162(d (if (< a d) (mod (+ d m) base) d)))163d))164165(defthmd submod-aux1166(implies (and (< a (mod (+ a (- b)) base))167(< 0 base) (< a base) (<= b base)168(natp base) (natp a) (natp b))169(< a b))170:rule-classes :forward-chaining)171172(defthmd submod-aux2173(implies (and (<= (mod (+ a (- b)) base) a)174(< 0 base) (< a base) (< b base)175(natp base) (natp a) (natp b))176(<= b a))177:rule-classes :forward-chaining)178179(defthmd submod-correct180(implies (and (< 0 m) (< m base)181(< a m) (<= b m)182(natp m) (natp base)183(natp a) (natp b))184(equal (submod a b m base)185(mod (- a b) m)))186:hints (("Goal" :cases ((<= base (+ a b))))187("Subgoal 2.2" :use ((:instance submod-aux1)))188("Subgoal 2.2'''" :cases ((and (< 0 (+ a (- b) m))189(< (+ a (- b) m) m))))190("Subgoal 2.1" :use ((:instance submod-aux2)))191("Subgoal 1.2" :use ((:instance submod-aux1)))192("Subgoal 1.1" :use ((:instance submod-aux2)))))193194195(defun submod-2 (a b m base)196(let* ((d (mod (- a b) base))197(d (if (< a b) (mod (+ d m) base) d)))198d))199200(defthm submod-2-correct201(implies (and (< 0 m) (< m base)202(< a m) (<= b m)203(natp m) (natp base)204(natp a) (natp b))205(equal (submod-2 a b m base)206(mod (- a b) m)))207:hints (("Subgoal 2'" :cases ((and (< 0 (+ a (- b) m))208(< (+ a (- b) m) m))))))209210211;; =========================================================================212;; ext-submod is correct213;; =========================================================================214215; a < 2*m, b < 2*m216(defun ext-submod (a b m base)217(let* ((a (if (>= a m) (- a m) a))218(b (if (>= b m) (- b m) b))219(d (mod (- a b) base))220(d (if (< a b) (mod (+ d m) base) d)))221d))222223; a < 2*m, b < 2*m224(defun ext-submod-2 (a b m base)225(let* ((a (mod a m))226(b (mod b m))227(d (mod (- a b) base))228(d (if (< a b) (mod (+ d m) base) d)))229d))230231(defthmd ext-submod-ext-submod-2-equal232(implies (and (< 0 m) (< m base)233(< a (* 2 m)) (< b (* 2 m))234(natp m) (natp base)235(natp a) (natp b))236(equal (ext-submod a b m base)237(ext-submod-2 a b m base))))238239(defthmd ext-submod-2-correct240(implies (and (< 0 m) (< m base)241(< a (* 2 m)) (< b (* 2 m))242(natp m) (natp base)243(natp a) (natp b))244(equal (ext-submod-2 a b m base)245(mod (- a b) m))))246247248;; =========================================================================249;; dw-reduce is correct250;; =========================================================================251252(defun dw-reduce (hi lo m base)253(let* ((r1 (mod hi m))254(r2 (mod (+ (* r1 base) lo) m)))255r2))256257(defthmd dw-reduce-correct258(implies (and (< 0 m) (< m base)259(< hi base) (< lo base)260(natp m) (natp base)261(natp hi) (natp lo))262(equal (dw-reduce hi lo m base)263(mod (+ (* hi base) lo) m))))264265(defthmd <=-multiply-both-sides-by-z266(implies (and (rationalp x) (rationalp y)267(< 0 z) (rationalp z))268(equal (<= x y)269(<= (* z x) (* z y)))))270271(defthmd dw-reduce-aux1272(implies (and (< 0 m) (< m base)273(natp m) (natp base)274(< lo base) (natp lo)275(< x m) (natp x))276(< (+ lo (* base x)) (* base m)))277:hints (("Goal" :cases ((<= (+ x 1) m)))278("Subgoal 1''" :cases ((<= (* base (+ x 1)) (* base m))))279("subgoal 1.2" :use ((:instance <=-multiply-both-sides-by-z280(x (+ 1 x))281(y m)282(z base))))))283284(defthm dw-reduce-aux2285(implies (and (< x (* base m))286(< 0 m) (< m base)287(natp m) (natp base) (natp x))288(< (floor x m) base)))289290;; This is the necessary condition for using _mpd_div_words().291(defthmd dw-reduce-second-quotient-fits-in-single-word292(implies (and (< 0 m) (< m base)293(< hi base) (< lo base)294(natp m) (natp base)295(natp hi) (natp lo)296(equal r1 (mod hi m)))297(< (floor (+ (* r1 base) lo) m)298base))299:hints (("Goal" :cases ((< r1 m)))300("Subgoal 1''" :cases ((< (+ lo (* base (mod hi m))) (* base m))))301("Subgoal 1.2" :use ((:instance dw-reduce-aux1302(x (mod hi m)))))))303304305;; =========================================================================306;; dw-submod is correct307;; =========================================================================308309(defun dw-submod (a hi lo m base)310(let* ((r (dw-reduce hi lo m base))311(d (mod (- a r) base))312(d (if (< a r) (mod (+ d m) base) d)))313d))314315(defthmd dw-submod-aux1316(implies (and (natp a) (< 0 m) (natp m)317(natp x) (equal r (mod x m)))318(equal (mod (- a x) m)319(mod (- a r) m))))320321(defthmd dw-submod-correct322(implies (and (< 0 m) (< m base)323(natp a) (< a m)324(< hi base) (< lo base)325(natp m) (natp base)326(natp hi) (natp lo))327(equal (dw-submod a hi lo m base)328(mod (- a (+ (* base hi) lo)) m)))329:hints (("Goal" :in-theory (disable dw-reduce)330:use ((:instance dw-submod-aux1331(x (+ lo (* base hi)))332(r (dw-reduce hi lo m base)))333(:instance dw-reduce-correct)))))334335336;; =========================================================================337;; ANSI C arithmetic for uint64_t338;; =========================================================================339340(defun add (a b)341(mod (+ a b)342(expt 2 64)))343344(defun sub (a b)345(mod (- a b)346(expt 2 64)))347348(defun << (w n)349(mod (* w (expt 2 n))350(expt 2 64)))351352(defun >> (w n)353(floor w (expt 2 n)))354355;; join upper and lower half of a double word, yielding a 128 bit number356(defun join (hi lo)357(+ (* (expt 2 64) hi) lo))358359360;; =============================================================================361;; Fast modular reduction362;; =============================================================================363364;; These are the three primes used in the Number Theoretic Transform.365;; A fast modular reduction scheme exists for all of them.366(defmacro p1 ()367(+ (expt 2 64) (- (expt 2 32)) 1))368369(defmacro p2 ()370(+ (expt 2 64) (- (expt 2 34)) 1))371372(defmacro p3 ()373(+ (expt 2 64) (- (expt 2 40)) 1))374375376;; reduce the double word number hi*2**64 + lo (mod p1)377(defun simple-mod-reduce-p1 (hi lo)378(+ (* (expt 2 32) hi) (- hi) lo))379380;; reduce the double word number hi*2**64 + lo (mod p2)381(defun simple-mod-reduce-p2 (hi lo)382(+ (* (expt 2 34) hi) (- hi) lo))383384;; reduce the double word number hi*2**64 + lo (mod p3)385(defun simple-mod-reduce-p3 (hi lo)386(+ (* (expt 2 40) hi) (- hi) lo))387388389; ----------------------------------------------------------390; The modular reductions given above are correct391; ----------------------------------------------------------392393(defthmd congruence-p1-aux394(equal (* (expt 2 64) hi)395(+ (* (p1) hi)396(* (expt 2 32) hi)397(- hi))))398399(defthmd congruence-p2-aux400(equal (* (expt 2 64) hi)401(+ (* (p2) hi)402(* (expt 2 34) hi)403(- hi))))404405(defthmd congruence-p3-aux406(equal (* (expt 2 64) hi)407(+ (* (p3) hi)408(* (expt 2 40) hi)409(- hi))))410411(defthmd mod-augment412(implies (and (rationalp x)413(rationalp y)414(rationalp m))415(equal (mod (+ x y) m)416(mod (+ x (mod y m)) m))))417418(defthmd simple-mod-reduce-p1-congruent419(implies (and (integerp hi)420(integerp lo))421(equal (mod (simple-mod-reduce-p1 hi lo) (p1))422(mod (join hi lo) (p1))))423:hints (("Goal''" :use ((:instance congruence-p1-aux)424(:instance mod-augment425(m (p1))426(x (+ (- hi) lo (* (expt 2 32) hi)))427(y (* (p1) hi)))))))428429(defthmd simple-mod-reduce-p2-congruent430(implies (and (integerp hi)431(integerp lo))432(equal (mod (simple-mod-reduce-p2 hi lo) (p2))433(mod (join hi lo) (p2))))434:hints (("Goal''" :use ((:instance congruence-p2-aux)435(:instance mod-augment436(m (p2))437(x (+ (- hi) lo (* (expt 2 34) hi)))438(y (* (p2) hi)))))))439440(defthmd simple-mod-reduce-p3-congruent441(implies (and (integerp hi)442(integerp lo))443(equal (mod (simple-mod-reduce-p3 hi lo) (p3))444(mod (join hi lo) (p3))))445:hints (("Goal''" :use ((:instance congruence-p3-aux)446(:instance mod-augment447(m (p3))448(x (+ (- hi) lo (* (expt 2 40) hi)))449(y (* (p3) hi)))))))450451452; ---------------------------------------------------------------------453; We need a number less than 2*p, so that we can use the trick from454; elim-mod-m<x<2*m for the final reduction.455; For p1, two modular reductions are sufficient, for p2 and p3 three.456; ---------------------------------------------------------------------457458;; p1: the first reduction is less than 2**96459(defthmd simple-mod-reduce-p1-<-2**96460(implies (and (< hi (expt 2 64))461(< lo (expt 2 64))462(natp hi) (natp lo))463(< (simple-mod-reduce-p1 hi lo)464(expt 2 96))))465466;; p1: the second reduction is less than 2*p1467(defthmd simple-mod-reduce-p1-<-2*p1468(implies (and (< hi (expt 2 64))469(< lo (expt 2 64))470(< (join hi lo) (expt 2 96))471(natp hi) (natp lo))472(< (simple-mod-reduce-p1 hi lo)473(* 2 (p1)))))474475476;; p2: the first reduction is less than 2**98477(defthmd simple-mod-reduce-p2-<-2**98478(implies (and (< hi (expt 2 64))479(< lo (expt 2 64))480(natp hi) (natp lo))481(< (simple-mod-reduce-p2 hi lo)482(expt 2 98))))483484;; p2: the second reduction is less than 2**69485(defthmd simple-mod-reduce-p2-<-2*69486(implies (and (< hi (expt 2 64))487(< lo (expt 2 64))488(< (join hi lo) (expt 2 98))489(natp hi) (natp lo))490(< (simple-mod-reduce-p2 hi lo)491(expt 2 69))))492493;; p3: the third reduction is less than 2*p2494(defthmd simple-mod-reduce-p2-<-2*p2495(implies (and (< hi (expt 2 64))496(< lo (expt 2 64))497(< (join hi lo) (expt 2 69))498(natp hi) (natp lo))499(< (simple-mod-reduce-p2 hi lo)500(* 2 (p2)))))501502503;; p3: the first reduction is less than 2**104504(defthmd simple-mod-reduce-p3-<-2**104505(implies (and (< hi (expt 2 64))506(< lo (expt 2 64))507(natp hi) (natp lo))508(< (simple-mod-reduce-p3 hi lo)509(expt 2 104))))510511;; p3: the second reduction is less than 2**81512(defthmd simple-mod-reduce-p3-<-2**81513(implies (and (< hi (expt 2 64))514(< lo (expt 2 64))515(< (join hi lo) (expt 2 104))516(natp hi) (natp lo))517(< (simple-mod-reduce-p3 hi lo)518(expt 2 81))))519520;; p3: the third reduction is less than 2*p3521(defthmd simple-mod-reduce-p3-<-2*p3522(implies (and (< hi (expt 2 64))523(< lo (expt 2 64))524(< (join hi lo) (expt 2 81))525(natp hi) (natp lo))526(< (simple-mod-reduce-p3 hi lo)527(* 2 (p3)))))528529530; -------------------------------------------------------------------------531; The simple modular reductions, adapted for compiler friendly C532; -------------------------------------------------------------------------533534(defun mod-reduce-p1 (hi lo)535(let* ((y hi)536(x y)537(hi (>> hi 32))538(x (sub lo x))539(hi (if (> x lo) (+ hi -1) hi))540(y (<< y 32))541(lo (add y x))542(hi (if (< lo y) (+ hi 1) hi)))543(+ (* hi (expt 2 64)) lo)))544545(defun mod-reduce-p2 (hi lo)546(let* ((y hi)547(x y)548(hi (>> hi 30))549(x (sub lo x))550(hi (if (> x lo) (+ hi -1) hi))551(y (<< y 34))552(lo (add y x))553(hi (if (< lo y) (+ hi 1) hi)))554(+ (* hi (expt 2 64)) lo)))555556(defun mod-reduce-p3 (hi lo)557(let* ((y hi)558(x y)559(hi (>> hi 24))560(x (sub lo x))561(hi (if (> x lo) (+ hi -1) hi))562(y (<< y 40))563(lo (add y x))564(hi (if (< lo y) (+ hi 1) hi)))565(+ (* hi (expt 2 64)) lo)))566567568; -------------------------------------------------------------------------569; The compiler friendly versions are equal to the simple versions570; -------------------------------------------------------------------------571572(defthm mod-reduce-aux1573(implies (and (<= 0 a) (natp a) (natp m)574(< (- m) b) (<= b 0)575(integerp b)576(< (mod (+ b a) m)577(mod a m)))578(equal (mod (+ b a) m)579(+ b (mod a m))))580:hints (("Subgoal 2" :use ((:instance modaux-1b581(x (+ a b)))))))582583(defthm mod-reduce-aux2584(implies (and (<= 0 a) (natp a) (natp m)585(< b m) (natp b)586(< (mod (+ b a) m)587(mod a m)))588(equal (+ m (mod (+ b a) m))589(+ b (mod a m)))))590591592(defthm mod-reduce-aux3593(implies (and (< 0 a) (natp a) (natp m)594(< (- m) b) (< b 0)595(integerp b)596(<= (mod a m)597(mod (+ b a) m)))598(equal (+ (- m) (mod (+ b a) m))599(+ b (mod a m))))600:hints (("Subgoal 1.2'" :use ((:instance modaux-1b601(x b))))602("Subgoal 1''" :use ((:instance modaux-2d603(x I))))))604605606(defthm mod-reduce-aux4607(implies (and (< 0 a) (natp a) (natp m)608(< b m) (natp b)609(<= (mod a m)610(mod (+ b a) m)))611(equal (mod (+ b a) m)612(+ b (mod a m)))))613614615(defthm mod-reduce-p1==simple-mod-reduce-p1616(implies (and (< hi (expt 2 64))617(< lo (expt 2 64))618(natp hi) (natp lo))619(equal (mod-reduce-p1 hi lo)620(simple-mod-reduce-p1 hi lo)))621:hints (("Goal" :in-theory (disable expt)622:cases ((< 0 hi)))623("Subgoal 1.2.2'" :use ((:instance mod-reduce-aux1624(m (expt 2 64))625(b (+ (- HI) LO))626(a (* (expt 2 32) hi)))))627("Subgoal 1.2.1'" :use ((:instance mod-reduce-aux3628(m (expt 2 64))629(b (+ (- HI) LO))630(a (* (expt 2 32) hi)))))631("Subgoal 1.1.2'" :use ((:instance mod-reduce-aux2632(m (expt 2 64))633(b (+ (- HI) LO))634(a (* (expt 2 32) hi)))))635("Subgoal 1.1.1'" :use ((:instance mod-reduce-aux4636(m (expt 2 64))637(b (+ (- HI) LO))638(a (* (expt 2 32) hi)))))))639640641(defthm mod-reduce-p2==simple-mod-reduce-p2642(implies (and (< hi (expt 2 64))643(< lo (expt 2 64))644(natp hi) (natp lo))645(equal (mod-reduce-p2 hi lo)646(simple-mod-reduce-p2 hi lo)))647:hints (("Goal" :cases ((< 0 hi)))648("Subgoal 1.2.2'" :use ((:instance mod-reduce-aux1649(m (expt 2 64))650(b (+ (- HI) LO))651(a (* (expt 2 34) hi)))))652("Subgoal 1.2.1'" :use ((:instance mod-reduce-aux3653(m (expt 2 64))654(b (+ (- HI) LO))655(a (* (expt 2 34) hi)))))656("Subgoal 1.1.2'" :use ((:instance mod-reduce-aux2657(m (expt 2 64))658(b (+ (- HI) LO))659(a (* (expt 2 34) hi)))))660("Subgoal 1.1.1'" :use ((:instance mod-reduce-aux4661(m (expt 2 64))662(b (+ (- HI) LO))663(a (* (expt 2 34) hi)))))))664665666(defthm mod-reduce-p3==simple-mod-reduce-p3667(implies (and (< hi (expt 2 64))668(< lo (expt 2 64))669(natp hi) (natp lo))670(equal (mod-reduce-p3 hi lo)671(simple-mod-reduce-p3 hi lo)))672:hints (("Goal" :cases ((< 0 hi)))673("Subgoal 1.2.2'" :use ((:instance mod-reduce-aux1674(m (expt 2 64))675(b (+ (- HI) LO))676(a (* (expt 2 40) hi)))))677("Subgoal 1.2.1'" :use ((:instance mod-reduce-aux3678(m (expt 2 64))679(b (+ (- HI) LO))680(a (* (expt 2 40) hi)))))681("Subgoal 1.1.2'" :use ((:instance mod-reduce-aux2682(m (expt 2 64))683(b (+ (- HI) LO))684(a (* (expt 2 40) hi)))))685("Subgoal 1.1.1'" :use ((:instance mod-reduce-aux4686(m (expt 2 64))687(b (+ (- HI) LO))688(a (* (expt 2 40) hi)))))))689690691692693694