Line data Source code
1 : /*
2 : * Copyright (c) 2008-2012 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 BASEARITH_H
30 : #define BASEARITH_H
31 :
32 :
33 : #include "mpdecimal.h"
34 : #include <stdio.h>
35 : #include "typearith.h"
36 :
37 :
38 : mpd_uint_t _mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
39 : mpd_size_t m, mpd_size_t n);
40 : void _mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n);
41 : mpd_uint_t _mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v);
42 : mpd_uint_t _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v,
43 : mpd_uint_t b);
44 : mpd_uint_t _mpd_baseincr(mpd_uint_t *u, mpd_size_t n);
45 : void _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
46 : mpd_size_t m, mpd_size_t n);
47 : void _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n);
48 : void _mpd_basemul(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_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
51 : mpd_uint_t v);
52 : mpd_uint_t _mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
53 : mpd_uint_t v);
54 : mpd_uint_t _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
55 : mpd_uint_t v, mpd_uint_t b);
56 : mpd_uint_t _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
57 : mpd_uint_t v);
58 : mpd_uint_t _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
59 : mpd_uint_t v, mpd_uint_t b);
60 : int _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r, const mpd_uint_t *uconst,
61 : const mpd_uint_t *vconst, mpd_size_t nplusm, mpd_size_t n);
62 : void _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n,
63 : mpd_size_t m, mpd_size_t shift);
64 : mpd_uint_t _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,
65 : mpd_size_t shift);
66 :
67 :
68 :
69 : #ifdef CONFIG_64
70 : extern const mpd_uint_t mprime_rdx;
71 :
72 : /*
73 : * Algorithm from: Division by Invariant Integers using Multiplication,
74 : * T. Granlund and P. L. Montgomery, Proceedings of the SIGPLAN '94
75 : * Conference on Programming Language Design and Implementation.
76 : *
77 : * http://gmplib.org/~tege/divcnst-pldi94.pdf
78 : *
79 : * Variables from the paper and their translations (See section 8):
80 : *
81 : * N := 64
82 : * d := MPD_RADIX
83 : * l := 64
84 : * m' := floor((2**(64+64) - 1)/MPD_RADIX) - 2**64
85 : *
86 : * Since N-l == 0:
87 : *
88 : * dnorm := d
89 : * n2 := hi
90 : * n10 := lo
91 : *
92 : * ACL2 proof: mpd-div-words-r-correct
93 : */
94 : static inline void
95 : _mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo)
96 : {
97 : mpd_uint_t n_adj, h, l, t;
98 : mpd_uint_t n1_neg;
99 :
100 : /* n1_neg = if lo >= 2**63 then MPD_UINT_MAX else 0 */
101 : n1_neg = (lo & (1ULL<<63)) ? MPD_UINT_MAX : 0;
102 : /* n_adj = if lo >= 2**63 then lo+MPD_RADIX else lo */
103 : n_adj = lo + (n1_neg & MPD_RADIX);
104 :
105 : /* (h, l) = if lo >= 2**63 then m'*(hi+1) else m'*hi */
106 : _mpd_mul_words(&h, &l, mprime_rdx, hi-n1_neg);
107 : l = l + n_adj;
108 : if (l < n_adj) h++;
109 : t = h + hi;
110 : /* At this point t == qest, with q == qest or q == qest+1:
111 : * 1) 0 <= 2**64*hi + lo - qest*MPD_RADIX < 2*MPD_RADIX
112 : */
113 :
114 : /* t = 2**64-1 - qest = 2**64 - (qest+1) */
115 : t = MPD_UINT_MAX - t;
116 :
117 : /* (h, l) = 2**64*MPD_RADIX - (qest+1)*MPD_RADIX */
118 : _mpd_mul_words(&h, &l, t, MPD_RADIX);
119 : l = l + lo;
120 : if (l < lo) h++;
121 : h += hi;
122 : h -= MPD_RADIX;
123 : /* (h, l) = 2**64*hi + lo - (qest+1)*MPD_RADIX (mod 2**128)
124 : * Case q == qest+1:
125 : * a) h == 0, l == r
126 : * b) q := h - t == qest+1
127 : * c) r := l
128 : * Case q == qest:
129 : * a) h == MPD_UINT_MAX, l == 2**64-(MPD_RADIX-r)
130 : * b) q := h - t == qest
131 : * c) r := l + MPD_RADIX = r
132 : */
133 :
134 : *q = (h - t);
135 : *r = l + (MPD_RADIX & h);
136 : }
137 : #else
138 : static inline void
139 0 : _mpd_div_words_r(mpd_uint_t *q, mpd_uint_t *r, mpd_uint_t hi, mpd_uint_t lo)
140 : {
141 0 : _mpd_div_words(q, r, hi, lo, MPD_RADIX);
142 0 : }
143 : #endif
144 :
145 :
146 : /* Multiply two single base MPD_RADIX words, store result in array w[2]. */
147 : static inline void
148 0 : _mpd_singlemul(mpd_uint_t w[2], mpd_uint_t u, mpd_uint_t v)
149 : {
150 : mpd_uint_t hi, lo;
151 :
152 0 : _mpd_mul_words(&hi, &lo, u, v);
153 0 : _mpd_div_words_r(&w[1], &w[0], hi, lo);
154 0 : }
155 :
156 : /* Multiply u (len 2) and v (len m, 1 <= m <= 2). */
157 : static inline void
158 0 : _mpd_mul_2_le2(mpd_uint_t w[4], mpd_uint_t u[2], mpd_uint_t v[2], mpd_ssize_t m)
159 : {
160 : mpd_uint_t hi, lo;
161 :
162 0 : _mpd_mul_words(&hi, &lo, u[0], v[0]);
163 0 : _mpd_div_words_r(&w[1], &w[0], hi, lo);
164 :
165 0 : _mpd_mul_words(&hi, &lo, u[1], v[0]);
166 0 : lo = w[1] + lo;
167 0 : if (lo < w[1]) hi++;
168 0 : _mpd_div_words_r(&w[2], &w[1], hi, lo);
169 0 : if (m == 1) return;
170 :
171 0 : _mpd_mul_words(&hi, &lo, u[0], v[1]);
172 0 : lo = w[1] + lo;
173 0 : if (lo < w[1]) hi++;
174 0 : _mpd_div_words_r(&w[3], &w[1], hi, lo);
175 :
176 0 : _mpd_mul_words(&hi, &lo, u[1], v[1]);
177 0 : lo = w[2] + lo;
178 0 : if (lo < w[2]) hi++;
179 0 : lo = w[3] + lo;
180 0 : if (lo < w[3]) hi++;
181 0 : _mpd_div_words_r(&w[3], &w[2], hi, lo);
182 : }
183 :
184 :
185 : /*
186 : * Test if all words from data[len-1] to data[0] are zero. If len is 0, nothing
187 : * is tested and the coefficient is regarded as "all zero".
188 : */
189 : static inline int
190 0 : _mpd_isallzero(const mpd_uint_t *data, mpd_ssize_t len)
191 : {
192 0 : while (--len >= 0) {
193 0 : if (data[len] != 0) return 0;
194 : }
195 0 : return 1;
196 : }
197 :
198 : /*
199 : * Test if all full words from data[len-1] to data[0] are MPD_RADIX-1
200 : * (all nines). Return true if len == 0.
201 : */
202 : static inline int
203 0 : _mpd_isallnine(const mpd_uint_t *data, mpd_ssize_t len)
204 : {
205 0 : while (--len >= 0) {
206 0 : if (data[len] != MPD_RADIX-1) return 0;
207 : }
208 0 : return 1;
209 : }
210 :
211 :
212 : #endif /* BASEARITH_H */
213 :
214 :
215 :
|
