1 | /* |
2 | * Copyright (c) 2008-2020 Stefan Krah. All rights reserved. |
3 | * |
4 | * Redistribution and use in source and binary forms, with or without |
5 | * modification, are permitted provided that the following conditions |
6 | * are met: |
7 | * |
8 | * 1. Redistributions of source code must retain the above copyright |
9 | * notice, this list of conditions and the following disclaimer. |
10 | * |
11 | * 2. Redistributions in binary form must reproduce the above copyright |
12 | * notice, this list of conditions and the following disclaimer in the |
13 | * documentation and/or other materials provided with the distribution. |
14 | * |
15 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |
16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
17 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
18 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
19 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
20 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
21 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
22 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
23 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
24 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
25 | * SUCH DAMAGE. |
26 | */ |
27 | |
28 | |
29 | #ifndef LIBMPDEC_BASEARITH_H_ |
30 | #define LIBMPDEC_BASEARITH_H_ |
31 | |
32 | |
33 | #include "mpdecimal.h" |
34 | #include "typearith.h" |
35 | |
36 | |
37 | /* Internal header file: all symbols have local scope in the DSO */ |
38 | MPD_PRAGMA(MPD_HIDE_SYMBOLS_START) |
39 | |
40 | |
41 | mpd_uint_t _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
42 | mpd_size_t m, mpd_size_t n); |
43 | void _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n); |
44 | mpd_uint_t _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v); |
45 | mpd_uint_t _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, |
46 | mpd_uint_t b); |
47 | mpd_uint_t _mpd_baseincr(mpd_uint_t *u, mpd_size_t n); |
48 | void _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
49 | mpd_size_t m, mpd_size_t n); |
50 | void _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n); |
51 | void _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v, |
52 | mpd_size_t m, mpd_size_t n); |
53 | void _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
54 | mpd_uint_t v); |
55 | mpd_uint_t _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
56 | mpd_uint_t v); |
57 | mpd_uint_t _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
58 | mpd_uint_t v, mpd_uint_t b); |
59 | mpd_uint_t _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
60 | mpd_uint_t v); |
61 | mpd_uint_t _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, |
62 | mpd_uint_t v, mpd_uint_t b); |
63 | int _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r, const mpd_uint_t *uconst, |
64 | const mpd_uint_t *vconst, mpd_size_t nplusm, mpd_size_t n); |
65 | void _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, |
66 | mpd_size_t m, mpd_size_t shift); |
67 | mpd_uint_t _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen, |
68 | mpd_size_t shift); |
69 | |
70 | |
71 | |
72 | #ifdef CONFIG_64 |
73 | extern const mpd_uint_t mprime_rdx; |
74 | |
75 | /* |
76 | * Algorithm from: Division by Invariant Integers using Multiplication, |
77 | * T. Granlund and P. L. Montgomery, Proceedings of the SIGPLAN '94 |
78 | * Conference on Programming Language Design and Implementation. |
79 | * |
80 | * http://gmplib.org/~tege/divcnst-pldi94.pdf |
81 | * |
82 | * Variables from the paper and their translations (See section 8): |
83 | * |
84 | * N := 64 |
85 | * d := MPD_RADIX |
86 | * l := 64 |
87 | * m' := floor((2**(64+64) - 1)/MPD_RADIX) - 2**64 |
88 | * |
89 | * Since N-l == 0: |
90 | * |
91 | * dnorm := d |
92 | * n2 := hi |
93 | * n10 := lo |
94 | * |
95 | * ACL2 proof: mpd-div-words-r-correct |
96 | */ |
97 | static inline void |
98 | _mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo) |
99 | { |
100 | mpd_uint_t n_adj, h, l, t; |
101 | mpd_uint_t n1_neg; |
102 | |
103 | /* n1_neg = if lo >= 2**63 then MPD_UINT_MAX else 0 */ |
104 | n1_neg = (lo & (1ULL<<63)) ? MPD_UINT_MAX : 0; |
105 | /* n_adj = if lo >= 2**63 then lo+MPD_RADIX else lo */ |
106 | n_adj = lo + (n1_neg & MPD_RADIX); |
107 | |
108 | /* (h, l) = if lo >= 2**63 then m'*(hi+1) else m'*hi */ |
109 | _mpd_mul_words(&h, &l, mprime_rdx, hi-n1_neg); |
110 | l = l + n_adj; |
111 | if (l < n_adj) h++; |
112 | t = h + hi; |
113 | /* At this point t == qest, with q == qest or q == qest+1: |
114 | * 1) 0 <= 2**64*hi + lo - qest*MPD_RADIX < 2*MPD_RADIX |
115 | */ |
116 | |
117 | /* t = 2**64-1 - qest = 2**64 - (qest+1) */ |
118 | t = MPD_UINT_MAX - t; |
119 | |
120 | /* (h, l) = 2**64*MPD_RADIX - (qest+1)*MPD_RADIX */ |
121 | _mpd_mul_words(&h, &l, t, MPD_RADIX); |
122 | l = l + lo; |
123 | if (l < lo) h++; |
124 | h += hi; |
125 | h -= MPD_RADIX; |
126 | /* (h, l) = 2**64*hi + lo - (qest+1)*MPD_RADIX (mod 2**128) |
127 | * Case q == qest+1: |
128 | * a) h == 0, l == r |
129 | * b) q := h - t == qest+1 |
130 | * c) r := l |
131 | * Case q == qest: |
132 | * a) h == MPD_UINT_MAX, l == 2**64-(MPD_RADIX-r) |
133 | * b) q := h - t == qest |
134 | * c) r := l + MPD_RADIX = r |
135 | */ |
136 | |
137 | *q = (h - t); |
138 | *r = l + (MPD_RADIX & h); |
139 | } |
140 | #else |
141 | static inline void |
142 | _mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo) |
143 | { |
144 | _mpd_div_words(q, r, hi, lo, MPD_RADIX); |
145 | } |
146 | #endif |
147 | |
148 | |
149 | /* Multiply two single base MPD_RADIX words, store result in array w[2]. */ |
150 | static inline void |
151 | _mpd_singlemul(mpd_uint_t w[2], mpd_uint_t u, mpd_uint_t v) |
152 | { |
153 | mpd_uint_t hi, lo; |
154 | |
155 | _mpd_mul_words(&hi, &lo, u, v); |
156 | _mpd_div_words_r(&w[1], &w[0], hi, lo); |
157 | } |
158 | |
159 | /* Multiply u (len 2) and v (len m, 1 <= m <= 2). */ |
160 | static inline void |
161 | _mpd_mul_2_le2(mpd_uint_t w[4], mpd_uint_t u[2], mpd_uint_t v[2], mpd_ssize_t m) |
162 | { |
163 | mpd_uint_t hi, lo; |
164 | |
165 | _mpd_mul_words(&hi, &lo, u[0], v[0]); |
166 | _mpd_div_words_r(&w[1], &w[0], hi, lo); |
167 | |
168 | _mpd_mul_words(&hi, &lo, u[1], v[0]); |
169 | lo = w[1] + lo; |
170 | if (lo < w[1]) hi++; |
171 | _mpd_div_words_r(&w[2], &w[1], hi, lo); |
172 | if (m == 1) return; |
173 | |
174 | _mpd_mul_words(&hi, &lo, u[0], v[1]); |
175 | lo = w[1] + lo; |
176 | if (lo < w[1]) hi++; |
177 | _mpd_div_words_r(&w[3], &w[1], hi, lo); |
178 | |
179 | _mpd_mul_words(&hi, &lo, u[1], v[1]); |
180 | lo = w[2] + lo; |
181 | if (lo < w[2]) hi++; |
182 | lo = w[3] + lo; |
183 | if (lo < w[3]) hi++; |
184 | _mpd_div_words_r(&w[3], &w[2], hi, lo); |
185 | } |
186 | |
187 | |
188 | /* |
189 | * Test if all words from data[len-1] to data[0] are zero. If len is 0, nothing |
190 | * is tested and the coefficient is regarded as "all zero". |
191 | */ |
192 | static inline int |
193 | _mpd_isallzero(const mpd_uint_t *data, mpd_ssize_t len) |
194 | { |
195 | while (--len >= 0) { |
196 | if (data[len] != 0) return 0; |
197 | } |
198 | return 1; |
199 | } |
200 | |
201 | /* |
202 | * Test if all full words from data[len-1] to data[0] are MPD_RADIX-1 |
203 | * (all nines). Return true if len == 0. |
204 | */ |
205 | static inline int |
206 | _mpd_isallnine(const mpd_uint_t *data, mpd_ssize_t len) |
207 | { |
208 | while (--len >= 0) { |
209 | if (data[len] != MPD_RADIX-1) return 0; |
210 | } |
211 | return 1; |
212 | } |
213 | |
214 | |
215 | MPD_PRAGMA(MPD_HIDE_SYMBOLS_END) /* restore previous scope rules */ |
216 | |
217 | |
218 | #endif /* LIBMPDEC_BASEARITH_H_ */ |
219 | |