math/mpx-mul4-x86-sse2.S: Maintain a local copy of the counter.
[catacomb] / math / mpx.c
CommitLineData
d03ab969 1/* -*-c-*-
2 *
d03ab969 3 * Low-level multiprecision arithmetic
4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
d03ab969 9 *
10 * This file is part of Catacomb.
11 *
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
45c0fd36 16 *
d03ab969 17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
45c0fd36 21 *
d03ab969 22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25 * MA 02111-1307, USA.
26 */
27
d03ab969 28/*----- Header files ------------------------------------------------------*/
29
444083ae
MW
30#include "config.h"
31
c8a2f9ef 32#include <assert.h>
d03ab969 33#include <stdio.h>
34#include <stdlib.h>
35#include <string.h>
36
37#include <mLib/bits.h>
23bbea75 38#include <mLib/macros.h>
d03ab969 39
444083ae 40#include "dispatch.h"
d03ab969 41#include "mptypes.h"
42#include "mpx.h"
75263f25 43#include "bitops.h"
d03ab969 44
45/*----- Loading and storing -----------------------------------------------*/
46
0c9ebe47
MW
47/* --- These are all variations on a theme --- *
48 *
49 * Essentially we want to feed bits into a shift register, @ibits@ bits at a
50 * time, and extract them @obits@ bits at a time whenever there are enough.
51 * Of course, @i@ and @o@ will, in general, be different sizes, and we don't
52 * necessarily know which is larger.
53 *
54 * During an operation, we have a shift register @w@ and a most-recent input
55 * @t@. Together, these hold @bits@ significant bits of input. We arrange
56 * that @bits < ibits + obits <= 2*MPW_BITS@, so we can get away with using
57 * an @mpw@ for both of these quantitities.
58 */
59
60/* --- @MPX_GETBITS@ --- *
61 *
62 * Arguments: @ibits@ = width of input units, in bits
63 * @obits@ = width of output units, in bits
64 * @iavail@ = condition expression: is input data available?
65 * @getbits@ = function or macro: set argument to next input
66 *
67 * Use: Read an input unit into @t@ and update the necessary
68 * variables.
69 *
70 * It is assumed on entry that @bits < obits@. On exit, we have
71 * @bits < ibits + obits@, and @t@ is live.
72 */
73
74#define MPX_GETBITS(ibits, obits, iavail, getbits) do { \
75 if (!iavail) goto flush; \
76 if (bits >= ibits) w |= t << (bits - ibits); \
77 getbits(t); \
78 bits += ibits; \
79} while (0)
80
81/* --- @MPX_PUTBITS@ --- *
82 *
83 * Arguments: @ibits@ = width of input units, in bits
84 * @obits@ = width of output units, in bits
85 * @oavail@ = condition expression: is output space available?
86 * @putbits@ = function or macro: write its argument to output
87 *
88 * Use: Emit an output unit, and update the necessary variables. If
89 * the output buffer is full, then force an immediate return.
90 *
91 * We assume that @bits < ibits + obits@, and that @t@ is only
92 * relevant if @bits >= ibits@. (The @MPX_GETBITS@ macro
93 * ensures that this is true.)
94 */
95
96#define SHRW(w, b) ((b) < MPW_BITS ? (w) >> (b) : 0)
97
98#define MPX_PUTBITS(ibits, obits, oavail, putbits) do { \
99 if (!oavail) return; \
100 if (bits < ibits) { \
101 putbits(w); \
102 bits -= obits; \
103 w = SHRW(w, obits); \
104 } else { \
105 putbits(w | (t << (bits - ibits))); \
106 bits -= obits; \
107 if (bits >= ibits) w = SHRW(w, obits) | (t << (bits - ibits)); \
108 else w = SHRW(w, obits) | (t >> (ibits - bits)); \
109 t = 0; \
110 } \
111} while (0)
112
113/* --- @MPX_LOADSTORE@ --- *
114 *
115 * Arguments: @name@ = name of function to create, without @mpx_@ prefix
116 * @wconst@ = qualifiers for @mpw *@ arguments
117 * @oconst@ = qualifiers for octet pointers
118 * @decls@ = additional declarations needed
119 * @ibits@ = width of input units, in bits
120 * @iavail@ = condition expression: is input data available?
121 * @getbits@ = function or macro: set argument to next input
122 * @obits@ = width of output units, in bits
123 * @oavail@ = condition expression: is output space available?
124 * @putbits@ = function or macro: write its argument to output
125 * @clear@ = statements to clear remainder of output
126 *
127 * Use: Generates a function to convert between a sequence of
128 * multiprecision words and a vector of octets.
129 *
130 * The arguments @ibits@, @iavail@ and @getbits@ are passed on
131 * to @MPX_GETBITS@; similarly, @obits@, @oavail@, and @putbits@
132 * are passed on to @MPX_PUTBITS@.
133 *
134 * The following variables are in scope: @v@ and @vl are the
135 * current base and limit of the word vector; @p@ and @q@ are
136 * the base and limit of the octet vector; @w@ and @t@ form the
137 * shift register used during the conversion (see commentary
138 * above); and @bits@ tracks the number of live bits in the
139 * shift register.
140 */
141
142#define MPX_LOADSTORE(name, wconst, oconst, decls, \
143 ibits, iavail, getbits, obits, oavail, putbits, \
144 clear) \
145 \
146void mpx_##name(wconst mpw *v, wconst mpw *vl, \
147 oconst void *pp, size_t sz) \
148{ \
149 mpw t = 0, w = 0; \
150 oconst octet *p = pp, *q = p + sz; \
151 int bits = 0; \
152 decls \
153 \
154 for (;;) { \
155 while (bits < obits) MPX_GETBITS(ibits, obits, iavail, getbits); \
156 while (bits >= obits) MPX_PUTBITS(ibits, obits, oavail, putbits); \
157 } \
158 \
159flush: \
160 while (bits > 0) MPX_PUTBITS(ibits, obits, oavail, putbits); \
161 clear; \
162}
163
164#define EMPTY
165
166/* --- Macros for @getbits@ and @putbits@ --- */
167
168#define GETMPW(t) do { t = *v++; } while (0)
169#define PUTMPW(x) do { *v++ = MPW(x); } while (0)
170
171#define GETOCTETI(t) do { t = *p++; } while (0)
172#define PUTOCTETD(x) do { *--q = U8(x); } while (0)
173
174#define PUTOCTETI(x) do { *p++ = U8(x); } while (0)
175#define GETOCTETD(t) do { t = *--q; } while (0)
176
177/* --- Machinery for two's complement I/O --- */
178
179#define DECL_2CN \
180 unsigned c = 1;
181
182#define GETMPW_2CN(t) do { \
183 t = MPW(~*v++ + c); \
184 c = c && !t; \
185} while (0)
186
187#define PUTMPW_2CN(t) do { \
188 mpw _t = MPW(~(t) + c); \
189 c = c && !_t; \
190 *v++ = _t; \
191} while (0)
192
193#define FLUSHW_2CN do { \
194 if (c) MPX_ONE(v, vl); \
195 else MPX_ZERO(v, vl); \
196} while (0)
197
198#define FLUSHO_2CN do { \
199 memset(p, c ? 0xff : 0, q - p); \
200} while (0)
201
d03ab969 202/* --- @mpx_storel@ --- *
203 *
204 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
c8a2f9ef 205 * @void *pp@ = pointer to octet array
d03ab969 206 * @size_t sz@ = size of octet array
207 *
208 * Returns: ---
209 *
210 * Use: Stores an MP in an octet array, least significant octet
211 * first. High-end octets are silently discarded if there
212 * isn't enough space for them.
213 */
214
0c9ebe47
MW
215MPX_LOADSTORE(storel, const, EMPTY, EMPTY,
216 MPW_BITS, (v < vl), GETMPW,
217 8, (p < q), PUTOCTETI,
218 { memset(p, 0, q - p); })
d03ab969 219
220/* --- @mpx_loadl@ --- *
221 *
222 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
c8a2f9ef 223 * @const void *pp@ = pointer to octet array
d03ab969 224 * @size_t sz@ = size of octet array
225 *
226 * Returns: ---
227 *
228 * Use: Loads an MP in an octet array, least significant octet
229 * first. High-end octets are ignored if there isn't enough
230 * space for them.
231 */
232
0c9ebe47
MW
233MPX_LOADSTORE(loadl, EMPTY, const, EMPTY,
234 8, (p < q), GETOCTETI,
235 MPW_BITS, (v < vl), PUTMPW,
236 { MPX_ZERO(v, vl); })
237
d03ab969 238
239/* --- @mpx_storeb@ --- *
240 *
241 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
c8a2f9ef 242 * @void *pp@ = pointer to octet array
d03ab969 243 * @size_t sz@ = size of octet array
244 *
245 * Returns: ---
246 *
247 * Use: Stores an MP in an octet array, most significant octet
248 * first. High-end octets are silently discarded if there
249 * isn't enough space for them.
250 */
251
0c9ebe47
MW
252MPX_LOADSTORE(storeb, const, EMPTY, EMPTY,
253 MPW_BITS, (v < vl), GETMPW,
254 8, (p < q), PUTOCTETD,
255 { memset(p, 0, q - p); })
d03ab969 256
257/* --- @mpx_loadb@ --- *
258 *
259 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
c8a2f9ef 260 * @const void *pp@ = pointer to octet array
d03ab969 261 * @size_t sz@ = size of octet array
262 *
263 * Returns: ---
264 *
265 * Use: Loads an MP in an octet array, most significant octet
266 * first. High-end octets are ignored if there isn't enough
267 * space for them.
268 */
269
0c9ebe47
MW
270MPX_LOADSTORE(loadb, EMPTY, const, EMPTY,
271 8, (p < q), GETOCTETD,
272 MPW_BITS, (v < vl), PUTMPW,
273 { MPX_ZERO(v, vl); })
d03ab969 274
f09e814a 275/* --- @mpx_storel2cn@ --- *
276 *
277 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
278 * @void *pp@ = pointer to octet array
279 * @size_t sz@ = size of octet array
280 *
281 * Returns: ---
282 *
283 * Use: Stores a negative MP in an octet array, least significant
284 * octet first, as two's complement. High-end octets are
285 * silently discarded if there isn't enough space for them.
286 * This obviously makes the output bad.
287 */
288
0c9ebe47
MW
289MPX_LOADSTORE(storel2cn, const, EMPTY, DECL_2CN,
290 MPW_BITS, (v < vl), GETMPW_2CN,
291 8, (p < q), PUTOCTETI,
292 { FLUSHO_2CN; })
f09e814a 293
294/* --- @mpx_loadl2cn@ --- *
295 *
296 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
297 * @const void *pp@ = pointer to octet array
298 * @size_t sz@ = size of octet array
299 *
300 * Returns: ---
301 *
302 * Use: Loads a negative MP in an octet array, least significant
303 * octet first, as two's complement. High-end octets are
304 * ignored if there isn't enough space for them. This probably
305 * means you made the wrong choice coming here.
306 */
307
0c9ebe47
MW
308MPX_LOADSTORE(loadl2cn, EMPTY, const, DECL_2CN,
309 8, (p < q), GETOCTETI,
310 MPW_BITS, (v < vl), PUTMPW_2CN,
311 { FLUSHW_2CN; })
f09e814a 312
313/* --- @mpx_storeb2cn@ --- *
314 *
315 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
316 * @void *pp@ = pointer to octet array
317 * @size_t sz@ = size of octet array
318 *
319 * Returns: ---
320 *
321 * Use: Stores a negative MP in an octet array, most significant
322 * octet first, as two's complement. High-end octets are
323 * silently discarded if there isn't enough space for them,
324 * which probably isn't what you meant.
325 */
326
0c9ebe47
MW
327MPX_LOADSTORE(storeb2cn, const, EMPTY, DECL_2CN,
328 MPW_BITS, (v < vl), GETMPW_2CN,
329 8, (p < q), PUTOCTETD,
330 { FLUSHO_2CN; })
f09e814a 331
332/* --- @mpx_loadb2cn@ --- *
333 *
334 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
335 * @const void *pp@ = pointer to octet array
336 * @size_t sz@ = size of octet array
337 *
338 * Returns: ---
339 *
340 * Use: Loads a negative MP in an octet array, most significant octet
341 * first as two's complement. High-end octets are ignored if
342 * there isn't enough space for them. This probably means you
343 * chose this function wrongly.
344 */
345
0c9ebe47
MW
346MPX_LOADSTORE(loadb2cn, EMPTY, const, DECL_2CN,
347 8, (p < q), GETOCTETD,
348 MPW_BITS, (v < vl), PUTMPW_2CN,
349 { FLUSHW_2CN; })
f09e814a 350
d03ab969 351/*----- Logical shifting --------------------------------------------------*/
352
5ee480b5 353/* --- @MPX_SHIFT1@ --- *
d03ab969 354 *
5ee480b5
MW
355 * Arguments: @init@ = initial accumulator value
356 * @out@ = expression to store in each output word
357 * @next@ = expression for next accumulator value
d03ab969 358 *
5ee480b5
MW
359 * Use: Performs a single-position shift. The input is scanned
360 * right-to-left. In the expressions @out@ and @next@, the
361 * accumulator is available in @w@ and the current input word is
362 * in @t@.
d03ab969 363 *
5ee480b5
MW
364 * This macro is intended to be used in the @shift1@ argument of
365 * @MPX_SHIFTOP@, and expects variables describing the operation
366 * to be set up accordingly.
d03ab969 367 */
368
5ee480b5
MW
369#define MPX_SHIFT1(init, out, next) do { \
370 mpw t, w = (init); \
371 while (av < avl) { \
372 if (dv >= dvl) break; \
373 t = MPW(*av++); \
374 *dv++ = (out); \
375 w = (next); \
376 } \
377 if (dv < dvl) { *dv++ = MPW(w); MPX_ZERO(dv, dvl); } \
378} while (0)
379
380/* --- @MPX_SHIFTW@ --- *
381 *
382 * Arguments: @max@ = the maximum shift (in words) which is nontrivial
383 * @clear@ = function (or macro) to clear low-order output words
384 * @copy@ = statement to copy words from input to output
385 *
386 * Use: Performs a shift by a whole number of words. If the shift
387 * amount is @max@ or more words, then the destination is
388 * @clear@ed entirely; otherwise, @copy@ is executed.
389 *
390 * This macro is intended to be used in the @shiftw@ argument of
391 * @MPX_SHIFTOP@, and expects variables describing the operation
392 * to be set up accordingly.
393 */
d03ab969 394
5ee480b5
MW
395#define MPX_SHIFTW(max, clear, copy) do { \
396 if (nw >= (max)) clear(dv, dvl); \
397 else copy \
398} while (0)
d03ab969 399
5ee480b5
MW
400/* --- @MPX_SHIFTOP@ --- *
401 *
402 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
403 * @shift1@ = statement to shift by a single bit
404 * @shiftw@ = statement to shift by a whole number of words
405 * @shift@ = statement to perform a general shift
406 *
407 * Use: Emits a shift operation. The input is @av@..@avl@; the
408 * output is @dv@..@dvl@; and the shift amount (in bits) is
409 * @n@. In @shiftw@ and @shift@, @nw@ and @nb@ are set up such
410 * that @n = nw*MPW_BITS + nb@ and @nb < MPW_BITS@.
411 */
d03ab969 412
5ee480b5
MW
413#define MPX_SHIFTOP(name, shift1, shiftw, shift) \
414 \
415void mpx_##name(mpw *dv, mpw *dvl, \
416 const mpw *av, const mpw *avl, \
417 size_t n) \
418{ \
419 \
420 if (n == 0) \
421 MPX_COPY(dv, dvl, av, avl); \
422 else if (n == 1) \
423 do shift1 while (0); \
424 else { \
425 size_t nw = n/MPW_BITS; \
426 unsigned nb = n%MPW_BITS; \
427 if (!nb) do shiftw while (0); \
428 else do shift while (0); \
429 } \
430}
d03ab969 431
5ee480b5
MW
432/* --- @MPX_SHIFT_LEFT@ --- *
433 *
434 * Arguments: @name@ = name of function to define (without `@mpx_@' prefix)
435 * @init1@ = initializer for single-bit shift accumulator
436 * @clear@ = function (or macro) to clear low-order output words
437 * @flush@ = expression for low-order nontrivial output word
438 *
439 * Use: Emits a left-shift operation. This expands to a call on
440 * @MPX_SHIFTOP@, but implements the complicated @shift@
441 * statement.
442 *
443 * The @init1@ argument is as for @MPX_SHIFT1@, and @clear@ is
444 * as for @MPX_SHIFTW@ (though is used elsewhere). In a general
445 * shift, @nw@ whole low-order output words are set using
446 * @clear@; high-order words are zeroed; and the remaining words
447 * set with a left-to-right pass across the input; at the end of
448 * the operation, the least significant output word above those
449 * @clear@ed is set using @flush@, which may use the accumulator
450 * @w@ = @av[0] << nb@.
451 */
d03ab969 452
5ee480b5
MW
453#define MPX_SHIFT_LEFT(name, init1, clear, flush) \
454MPX_SHIFTOP(name, { \
455 MPX_SHIFT1(init1, \
456 w | (t << 1), \
457 t >> (MPW_BITS - 1)); \
458}, { \
459 MPX_SHIFTW(dvl - dv, clear, { \
460 MPX_COPY(dv + nw, dvl, av, avl); \
461 clear(dv, dv + nw); \
462 }); \
463}, { \
464 size_t nr = MPW_BITS - nb; \
465 size_t dvn = dvl - dv; \
466 size_t avn = avl - av; \
467 mpw w; \
468 \
469 if (dvn <= nw) { \
470 clear(dv, dvl); \
471 break; \
472 } \
473 \
474 if (dvn <= avn + nw) { \
475 avl = av + dvn - nw; \
476 w = *--avl << nb; \
477 } else { \
478 size_t off = avn + nw + 1; \
479 MPX_ZERO(dv + off, dvl); \
480 dvl = dv + off; \
481 w = 0; \
482 } \
483 \
484 while (avl > av) { \
485 mpw t = *--avl; \
486 *--dvl = MPW(w | (t >> nr)); \
487 w = t << nb; \
488 } \
489 \
490 *--dvl = MPW(flush); \
491 clear(dv, dvl); \
492})
c8a2f9ef 493
5ee480b5
MW
494/* --- @mpx_lsl@ --- *
495 *
496 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
497 * @const mpw *av, *avl@ = source vector base and limit
498 * @size_t n@ = number of bit positions to shift by
499 *
500 * Returns: ---
501 *
502 * Use: Performs a logical shift left operation on an integer.
503 */
d03ab969 504
5ee480b5 505MPX_SHIFT_LEFT(lsl, 0, MPX_ZERO, w)
d03ab969 506
81578196 507/* --- @mpx_lslc@ --- *
508 *
509 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
510 * @const mpw *av, *avl@ = source vector base and limit
511 * @size_t n@ = number of bit positions to shift by
512 *
513 * Returns: ---
514 *
515 * Use: Performs a logical shift left operation on an integer, only
516 * it fills in the bits with ones instead of zeroes.
517 */
518
5ee480b5 519MPX_SHIFT_LEFT(lslc, 1, MPX_ONE, w | (MPW_MAX >> nr))
81578196 520
d03ab969 521/* --- @mpx_lsr@ --- *
522 *
523 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
524 * @const mpw *av, *avl@ = source vector base and limit
525 * @size_t n@ = number of bit positions to shift by
526 *
527 * Returns: ---
528 *
529 * Use: Performs a logical shift right operation on an integer.
530 */
531
5ee480b5
MW
532MPX_SHIFTOP(lsr, {
533 MPX_SHIFT1(av < avl ? *av++ >> 1 : 0,
534 w | (t << (MPW_BITS - 1)),
535 t >> 1);
536}, {
537 MPX_SHIFTW(avl - av, MPX_ZERO,
538 { MPX_COPY(dv, dvl, av + nw, avl); });
539}, {
540 size_t nr = MPW_BITS - nb;
541 mpw w;
542
543 av += nw;
544 w = av < avl ? *av++ : 0;
545 while (av < avl) {
546 mpw t;
547 if (dv >= dvl) goto done;
548 t = *av++;
549 *dv++ = MPW((w >> nb) | (t << nr));
550 w = t;
d03ab969 551 }
5ee480b5
MW
552 if (dv < dvl) {
553 *dv++ = MPW(w >> nb);
554 MPX_ZERO(dv, dvl);
d03ab969 555 }
d03ab969 556done:;
5ee480b5 557})
d03ab969 558
0f32e0f8 559/*----- Bitwise operations ------------------------------------------------*/
560
f09e814a 561/* --- @mpx_bitop@ --- *
0f32e0f8 562 *
563 * Arguments: @mpw *dv, *dvl@ = destination vector
564 * @const mpw *av, *avl@ = first source vector
565 * @const mpw *bv, *bvl@ = second source vector
566 *
567 * Returns: ---
568 *
f09e814a 569 * Use; Provides the dyadic boolean functions.
0f32e0f8 570 */
571
f09e814a 572#define MPX_BITBINOP(string) \
0f32e0f8 573 \
f09e814a 574void mpx_bit##string(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, \
575 const mpw *bv, const mpw *bvl) \
0f32e0f8 576{ \
577 MPX_SHRINK(av, avl); \
578 MPX_SHRINK(bv, bvl); \
579 \
580 while (dv < dvl) { \
581 mpw a, b; \
582 a = (av < avl) ? *av++ : 0; \
583 b = (bv < bvl) ? *bv++ : 0; \
75263f25 584 *dv++ = B##string(a, b); \
23bbea75 585 IGNORE(a); IGNORE(b); \
0f32e0f8 586 } \
587}
588
f09e814a 589MPX_DOBIN(MPX_BITBINOP)
0f32e0f8 590
591void mpx_not(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
592{
593 MPX_SHRINK(av, avl);
594
595 while (dv < dvl) {
596 mpw a;
597 a = (av < avl) ? *av++ : 0;
598 *dv++ = ~a;
599 }
600}
601
d03ab969 602/*----- Unsigned arithmetic -----------------------------------------------*/
603
f45a00c6 604/* --- @mpx_2c@ --- *
605 *
606 * Arguments: @mpw *dv, *dvl@ = destination vector
607 * @const mpw *v, *vl@ = source vector
608 *
609 * Returns: ---
610 *
611 * Use: Calculates the two's complement of @v@.
612 */
613
614void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
615{
616 mpw c = 0;
617 while (dv < dvl && v < vl)
618 *dv++ = c = MPW(~*v++);
619 if (dv < dvl) {
620 if (c > MPW_MAX / 2)
621 c = MPW(~0);
622 while (dv < dvl)
623 *dv++ = c;
624 }
625 MPX_UADDN(dv, dvl, 1);
626}
627
1a05a8ef 628/* --- @mpx_ueq@ --- *
629 *
630 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
631 * @const mpw *bv, *bvl@ = second argument vector base and limit
632 *
633 * Returns: Nonzero if the two vectors are equal.
634 *
635 * Use: Performs an unsigned integer test for equality.
636 */
637
638int mpx_ueq(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
639{
640 MPX_SHRINK(av, avl);
641 MPX_SHRINK(bv, bvl);
642 if (avl - av != bvl - bv)
643 return (0);
644 while (av < avl) {
645 if (*av++ != *bv++)
646 return (0);
647 }
648 return (1);
649}
650
d03ab969 651/* --- @mpx_ucmp@ --- *
652 *
653 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
654 * @const mpw *bv, *bvl@ = second argument vector base and limit
655 *
656 * Returns: Less than, equal to, or greater than zero depending on
657 * whether @a@ is less than, equal to or greater than @b@,
658 * respectively.
659 *
660 * Use: Performs an unsigned integer comparison.
661 */
662
663int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
664{
665 MPX_SHRINK(av, avl);
666 MPX_SHRINK(bv, bvl);
667
668 if (avl - av > bvl - bv)
669 return (+1);
670 else if (avl - av < bvl - bv)
671 return (-1);
672 else while (avl > av) {
673 mpw a = *--avl, b = *--bvl;
674 if (a > b)
675 return (+1);
676 else if (a < b)
677 return (-1);
678 }
679 return (0);
680}
1a05a8ef 681
d03ab969 682/* --- @mpx_uadd@ --- *
683 *
684 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
685 * @const mpw *av, *avl@ = first addend vector base and limit
686 * @const mpw *bv, *bvl@ = second addend vector base and limit
687 *
688 * Returns: ---
689 *
690 * Use: Performs unsigned integer addition. If the result overflows
691 * the destination vector, high-order bits are discarded. This
692 * means that two's complement addition happens more or less for
693 * free, although that's more a side-effect than anything else.
694 * The result vector may be equal to either or both source
695 * vectors, but may not otherwise overlap them.
696 */
697
698void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
699 const mpw *bv, const mpw *bvl)
700{
701 mpw c = 0;
702
703 while (av < avl || bv < bvl) {
704 mpw a, b;
705 mpd x;
706 if (dv >= dvl)
707 return;
708 a = (av < avl) ? *av++ : 0;
709 b = (bv < bvl) ? *bv++ : 0;
710 x = (mpd)a + (mpd)b + c;
711 *dv++ = MPW(x);
712 c = x >> MPW_BITS;
713 }
714 if (dv < dvl) {
715 *dv++ = c;
716 MPX_ZERO(dv, dvl);
717 }
718}
719
dd517851 720/* --- @mpx_uaddn@ --- *
721 *
722 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
723 * @mpw n@ = other addend
724 *
725 * Returns: ---
726 *
727 * Use: Adds a small integer to a multiprecision number.
728 */
729
730void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
731
f46efa79 732/* --- @mpx_uaddnlsl@ --- *
733 *
734 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
735 * @mpw a@ = second argument
736 * @unsigned o@ = offset in bits
737 *
738 * Returns: ---
739 *
740 * Use: Computes %$d + 2^o a$%. If the result overflows then
741 * high-order bits are discarded, as usual. We must have
742 * @0 < o < MPW_BITS@.
743 */
744
745void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
746{
747 mpd x = (mpd)a << o;
748
749 while (x && dv < dvl) {
750 x += *dv;
751 *dv++ = MPW(x);
752 x >>= MPW_BITS;
753 }
754}
755
d03ab969 756/* --- @mpx_usub@ --- *
757 *
758 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
759 * @const mpw *av, *avl@ = first argument vector base and limit
760 * @const mpw *bv, *bvl@ = second argument vector base and limit
761 *
762 * Returns: ---
763 *
764 * Use: Performs unsigned integer subtraction. If the result
765 * overflows the destination vector, high-order bits are
766 * discarded. This means that two's complement subtraction
767 * happens more or less for free, althuogh that's more a side-
768 * effect than anything else. The result vector may be equal to
769 * either or both source vectors, but may not otherwise overlap
770 * them.
771 */
772
773void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
774 const mpw *bv, const mpw *bvl)
775{
776 mpw c = 0;
777
778 while (av < avl || bv < bvl) {
779 mpw a, b;
780 mpd x;
781 if (dv >= dvl)
782 return;
783 a = (av < avl) ? *av++ : 0;
784 b = (bv < bvl) ? *bv++ : 0;
c8a2f9ef 785 x = (mpd)a - (mpd)b - c;
d03ab969 786 *dv++ = MPW(x);
c8a2f9ef 787 if (x >> MPW_BITS)
788 c = 1;
789 else
790 c = 0;
d03ab969 791 }
c8a2f9ef 792 if (c)
793 c = MPW_MAX;
d03ab969 794 while (dv < dvl)
c8a2f9ef 795 *dv++ = c;
d03ab969 796}
797
dd517851 798/* --- @mpx_usubn@ --- *
799 *
800 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
801 * @n@ = subtrahend
802 *
803 * Returns: ---
804 *
805 * Use: Subtracts a small integer from a multiprecision number.
806 */
807
808void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
809
f46efa79 810/* --- @mpx_uaddnlsl@ --- *
811 *
812 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
813 * @mpw a@ = second argument
814 * @unsigned o@ = offset in bits
815 *
816 * Returns: ---
817 *
818 * Use: Computes %$d + 2^o a$%. If the result overflows then
819 * high-order bits are discarded, as usual. We must have
820 * @0 < o < MPW_BITS@.
821 */
822
823void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
824{
825 mpw b = a >> (MPW_BITS - o);
826 a <<= o;
827
828 if (dv < dvl) {
c29970a7 829 mpd x = (mpd)*dv - MPW(a);
f46efa79 830 *dv++ = MPW(x);
831 if (x >> MPW_BITS)
832 b++;
833 MPX_USUBN(dv, dvl, b);
834 }
835}
836
d03ab969 837/* --- @mpx_umul@ --- *
838 *
839 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
840 * @const mpw *av, *avl@ = multiplicand vector base and limit
841 * @const mpw *bv, *bvl@ = multiplier vector base and limit
842 *
843 * Returns: ---
844 *
845 * Use: Performs unsigned integer multiplication. If the result
846 * overflows the desination vector, high-order bits are
847 * discarded. The result vector may not overlap the argument
848 * vectors in any way.
849 */
850
444083ae
MW
851CPU_DISPATCH(EMPTY, (void), void, mpx_umul,
852 (mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
853 const mpw *bv, const mpw *bvl),
854 (dv, dvl, av, avl, bv, bvl), pick_umul, simple_umul);
855
856static void simple_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
857 const mpw *bv, const mpw *bvl)
d03ab969 858{
859 /* --- This is probably worthwhile on a multiply --- */
860
861 MPX_SHRINK(av, avl);
862 MPX_SHRINK(bv, bvl);
863
864 /* --- Deal with a multiply by zero --- */
45c0fd36 865
d03ab969 866 if (bv == bvl) {
c8a2f9ef 867 MPX_ZERO(dv, dvl);
d03ab969 868 return;
869 }
870
871 /* --- Do the initial multiply and initialize the accumulator --- */
872
873 MPX_UMULN(dv, dvl, av, avl, *bv++);
874
875 /* --- Do the remaining multiply/accumulates --- */
876
c8a2f9ef 877 while (dv < dvl && bv < bvl) {
d03ab969 878 mpw m = *bv++;
c8a2f9ef 879 mpw c = 0;
d03ab969 880 const mpw *avv = av;
881 mpw *dvv = ++dv;
882
883 while (avv < avl) {
884 mpd x;
885 if (dvv >= dvl)
886 goto next;
c8a2f9ef 887 x = (mpd)*dvv + (mpd)m * (mpd)*avv++ + c;
888 *dvv++ = MPW(x);
d03ab969 889 c = x >> MPW_BITS;
890 }
c8a2f9ef 891 MPX_UADDN(dvv, dvl, c);
d03ab969 892 next:;
893 }
894}
895
444083ae
MW
896#define MAYBE_UMUL4(impl) \
897 extern void mpx_umul4_##impl(mpw */*dv*/, \
898 const mpw */*av*/, const mpw */*avl*/, \
899 const mpw */*bv*/, const mpw */*bvl*/); \
900 static void maybe_umul4_##impl(mpw *dv, mpw *dvl, \
901 const mpw *av, const mpw *avl, \
902 const mpw *bv, const mpw *bvl) \
903 { \
904 size_t an = avl - av, bn = bvl - bv, dn = dvl - dv; \
905 if (!an || an%4 != 0 || !bn || bn%4 != 0 || dn < an + bn) \
906 simple_umul(dv, dvl, av, avl, bv, bvl); \
907 else { \
908 mpx_umul4_##impl(dv, av, avl, bv, bvl); \
909 MPX_ZERO(dv + an + bn, dvl); \
910 } \
911 }
912
913#if CPUFAM_X86
914 MAYBE_UMUL4(x86_sse2)
915#endif
916
917static mpx_umul__functype *pick_umul(void)
918{
919#if CPUFAM_X86
920 DISPATCH_PICK_COND(mpx_umul, maybe_umul4_x86_sse2,
921 cpu_feature_p(CPUFEAT_X86_SSE2));
922#endif
923 DISPATCH_PICK_FALLBACK(mpx_umul, simple_umul);
924}
925
dd517851 926/* --- @mpx_umuln@ --- *
927 *
928 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
929 * @const mpw *av, *avl@ = multiplicand vector base and limit
930 * @mpw m@ = multiplier
931 *
932 * Returns: ---
933 *
934 * Use: Multiplies a multiprecision integer by a single-word value.
935 * The destination and source may be equal. The destination
936 * is completely cleared after use.
937 */
938
939void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
106b481c 940 { MPX_UMULN(dv, dvl, av, avl, m); }
dd517851 941
942/* --- @mpx_umlan@ --- *
943 *
944 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
945 * @const mpw *av, *avl@ = multiplicand vector base and limit
946 * @mpw m@ = multiplier
947 *
948 * Returns: ---
949 *
950 * Use: Multiplies a multiprecision integer by a single-word value
951 * and adds the result to an accumulator.
952 */
953
954void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m)
106b481c 955 { MPX_UMLAN(dv, dvl, av, avl, m); }
dd517851 956
c8a2f9ef 957/* --- @mpx_usqr@ --- *
958 *
959 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
960 * @const mpw *av, *av@ = source vector base and limit
961 *
962 * Returns: ---
963 *
964 * Use: Performs unsigned integer squaring. The result vector must
965 * not overlap the source vector in any way.
966 */
967
968void mpx_usqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
969{
970 MPX_ZERO(dv, dvl);
971
972 /* --- Main loop --- */
973
974 while (av < avl) {
975 const mpw *avv = av;
976 mpw *dvv = dv;
977 mpw a = *av;
978 mpd c;
979
980 /* --- Stop if I've run out of destination --- */
981
982 if (dvv >= dvl)
983 break;
984
985 /* --- Work out the square at this point in the proceedings --- */
986
987 {
c8a2f9ef 988 mpd x = (mpd)a * (mpd)a + *dvv;
989 *dvv++ = MPW(x);
990 c = MPW(x >> MPW_BITS);
991 }
992
993 /* --- Now fix up the rest of the vector upwards --- */
994
995 avv++;
996 while (dvv < dvl && avv < avl) {
c8a2f9ef 997 mpd x = (mpd)a * (mpd)*avv++;
998 mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
999 c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
1000 *dvv++ = MPW(y);
1001 }
1002 while (dvv < dvl && c) {
1003 mpd x = c + *dvv;
1004 *dvv++ = MPW(x);
1005 c = x >> MPW_BITS;
1006 }
1007
1008 /* --- Get ready for the next round --- */
1009
1010 av++;
1011 dv += 2;
1012 }
1013}
1014
d03ab969 1015/* --- @mpx_udiv@ --- *
1016 *
1017 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
1018 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
1019 * @const mpw *dv, *dvl@ = divisor vector base and limit
c8a2f9ef 1020 * @mpw *sv, *svl@ = scratch workspace
d03ab969 1021 *
1022 * Returns: ---
1023 *
1024 * Use: Performs unsigned integer division. If the result overflows
1025 * the quotient vector, high-order bits are discarded. (Clearly
1026 * the remainder vector can't overflow.) The various vectors
1027 * may not overlap in any way. Yes, I know it's a bit odd
1028 * requiring the dividend to be in the result position but it
1029 * does make some sense really. The remainder must have
c8a2f9ef 1030 * headroom for at least two extra words. The scratch space
f45a00c6 1031 * must be at least one word larger than the divisor.
d03ab969 1032 */
1033
1034void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
c8a2f9ef 1035 const mpw *dv, const mpw *dvl,
1036 mpw *sv, mpw *svl)
d03ab969 1037{
d03ab969 1038 unsigned norm = 0;
1039 size_t scale;
1040 mpw d, dd;
1041
1042 /* --- Initialize the quotient --- */
1043
1044 MPX_ZERO(qv, qvl);
1045
c8a2f9ef 1046 /* --- Perform some sanity checks --- */
1047
1048 MPX_SHRINK(dv, dvl);
1049 assert(((void)"division by zero in mpx_udiv", dv < dvl));
1050
d03ab969 1051 /* --- Normalize the divisor --- *
1052 *
1053 * The algorithm requires that the divisor be at least two digits long.
1054 * This is easy to fix.
1055 */
1056
c8a2f9ef 1057 {
1058 unsigned b;
d03ab969 1059
c8a2f9ef 1060 d = dvl[-1];
c29970a7 1061 for (b = MPW_P2; b; b >>= 1) {
34e4f738 1062 if (d <= (MPW_MAX >> b)) {
c8a2f9ef 1063 d <<= b;
1064 norm += b;
1065 }
1066 }
1067 if (dv + 1 == dvl)
1068 norm += MPW_BITS;
d03ab969 1069 }
d03ab969 1070
1071 /* --- Normalize the dividend/remainder to match --- */
1072
c8a2f9ef 1073 if (norm) {
c8a2f9ef 1074 mpx_lsl(rv, rvl, rv, rvl, norm);
f45a00c6 1075 mpx_lsl(sv, svl, dv, dvl, norm);
c8a2f9ef 1076 dv = sv;
f45a00c6 1077 dvl = svl;
c8a2f9ef 1078 MPX_SHRINK(dv, dvl);
1079 }
1080
d03ab969 1081 MPX_SHRINK(rv, rvl);
c8a2f9ef 1082 d = dvl[-1];
1083 dd = dvl[-2];
d03ab969 1084
1085 /* --- Work out the relative scales --- */
1086
1087 {
1088 size_t rvn = rvl - rv;
c8a2f9ef 1089 size_t dvn = dvl - dv;
d03ab969 1090
1091 /* --- If the divisor is clearly larger, notice this --- */
1092
1093 if (dvn > rvn) {
1094 mpx_lsr(rv, rvl, rv, rvl, norm);
1095 return;
1096 }
1097
1098 scale = rvn - dvn;
1099 }
1100
1101 /* --- Calculate the most significant quotient digit --- *
1102 *
1103 * Because the divisor has its top bit set, this can only happen once. The
1104 * pointer arithmetic is a little contorted, to make sure that the
1105 * behaviour is defined.
1106 */
1107
1108 if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
1109 mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1110 if (qvl - qv > scale)
1111 qv[scale] = 1;
1112 }
1113
1114 /* --- Now for the main loop --- */
1115
1116 {
c8a2f9ef 1117 mpw *rvv = rvl - 2;
d03ab969 1118
1119 while (scale) {
c8a2f9ef 1120 mpw q;
1121 mpd rh;
d03ab969 1122
1123 /* --- Get an estimate for the next quotient digit --- */
1124
c8a2f9ef 1125 mpw r = rvv[1];
1126 mpw rr = rvv[0];
1127 mpw rrr = *--rvv;
1128
1129 scale--;
1130 rh = ((mpd)r << MPW_BITS) | rr;
d03ab969 1131 if (r == d)
1132 q = MPW_MAX;
c8a2f9ef 1133 else
1134 q = MPW(rh / d);
d03ab969 1135
1136 /* --- Refine the estimate --- */
1137
1138 {
1139 mpd yh = (mpd)d * q;
ce76ff16 1140 mpd yy = (mpd)dd * q;
1141 mpw yl;
c8a2f9ef 1142
ce76ff16 1143 if (yy > MPW_MAX)
1144 yh += yy >> MPW_BITS;
1145 yl = MPW(yy);
c8a2f9ef 1146
1147 while (yh > rh || (yh == rh && yl > rrr)) {
1148 q--;
1149 yh -= d;
ce76ff16 1150 if (yl < dd)
1151 yh--;
99b30c23 1152 yl = MPW(yl - dd);
c8a2f9ef 1153 }
1154 }
1155
1156 /* --- Remove a chunk from the dividend --- */
1157
1158 {
1159 mpw *svv;
1160 const mpw *dvv;
f45a00c6 1161 mpw mc = 0, sc = 0;
c8a2f9ef 1162
f45a00c6 1163 /* --- Calculate the size of the chunk --- *
1164 *
1165 * This does the whole job of calculating @r >> scale - qd@.
1166 */
c8a2f9ef 1167
f45a00c6 1168 for (svv = rv + scale, dvv = dv;
1169 dvv < dvl && svv < rvl;
1170 svv++, dvv++) {
1171 mpd x = (mpd)*dvv * (mpd)q + mc;
1172 mc = x >> MPW_BITS;
1173 x = (mpd)*svv - MPW(x) - sc;
c8a2f9ef 1174 *svv = MPW(x);
f45a00c6 1175 if (x >> MPW_BITS)
1176 sc = 1;
1177 else
1178 sc = 0;
1179 }
1180
1181 if (svv < rvl) {
1182 mpd x = (mpd)*svv - mc - sc;
1183 *svv++ = MPW(x);
1184 if (x >> MPW_BITS)
1185 sc = MPW_MAX;
1186 else
1187 sc = 0;
1188 while (svv < rvl)
1189 *svv++ = sc;
c8a2f9ef 1190 }
c8a2f9ef 1191
f45a00c6 1192 /* --- Fix if the quotient was too large --- *
c8a2f9ef 1193 *
f45a00c6 1194 * This doesn't seem to happen very often.
c8a2f9ef 1195 */
1196
c8a2f9ef 1197 if (rvl[-1] > MPW_MAX / 2) {
1198 mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
1199 q--;
1200 }
1201 }
1202
1203 /* --- Done for another iteration --- */
1204
1205 if (qvl - qv > scale)
1206 qv[scale] = q;
1207 r = rr;
1208 rr = rrr;
1209 }
1210 }
1211
1212 /* --- Now fiddle with unnormalizing and things --- */
1213
1214 mpx_lsr(rv, rvl, rv, rvl, norm);
d03ab969 1215}
1216
698bd937 1217/* --- @mpx_udivn@ --- *
1218 *
1219 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
1220 * dividend)
1221 * @const mpw *rv, *rvl@ = dividend
1222 * @mpw d@ = single-precision divisor
1223 *
1224 * Returns: Remainder after divison.
1225 *
1226 * Use: Performs a single-precision division operation.
1227 */
1228
1229mpw mpx_udivn(mpw *qv, mpw *qvl, const mpw *rv, const mpw *rvl, mpw d)
1230{
1231 size_t i;
1232 size_t ql = qvl - qv;
1233 mpd r = 0;
1234
1235 i = rvl - rv;
1236 while (i > 0) {
1237 i--;
1238 r = (r << MPW_BITS) | rv[i];
1239 if (i < ql)
1240 qv[i] = r / d;
1241 r %= d;
1242 }
1243 return (MPW(r));
1244}
1245
42684bdb 1246/*----- Test rig ----------------------------------------------------------*/
1247
1248#ifdef TEST_RIG
1249
1250#include <mLib/alloc.h>
1251#include <mLib/dstr.h>
1252#include <mLib/quis.h>
1253#include <mLib/testrig.h>
1254
1255#include "mpscan.h"
1256
1257#define ALLOC(v, vl, sz) do { \
1258 size_t _sz = (sz); \
1259 mpw *_vv = xmalloc(MPWS(_sz)); \
1260 mpw *_vvl = _vv + _sz; \
444083ae 1261 memset(_vv, 0xa5, MPWS(_sz)); \
42684bdb 1262 (v) = _vv; \
1263 (vl) = _vvl; \
1264} while (0)
1265
1266#define LOAD(v, vl, d) do { \
1267 const dstr *_d = (d); \
1268 mpw *_v, *_vl; \
1269 ALLOC(_v, _vl, MPW_RQ(_d->len)); \
1270 mpx_loadb(_v, _vl, _d->buf, _d->len); \
1271 (v) = _v; \
1272 (vl) = _vl; \
1273} while (0)
1274
1275#define MAX(x, y) ((x) > (y) ? (x) : (y))
45c0fd36 1276
42684bdb 1277static void dumpbits(const char *msg, const void *pp, size_t sz)
1278{
1279 const octet *p = pp;
1280 fputs(msg, stderr);
1281 for (; sz; sz--)
1282 fprintf(stderr, " %02x", *p++);
1283 fputc('\n', stderr);
1284}
1285
1286static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
1287{
1288 fputs(msg, stderr);
1289 MPX_SHRINK(v, vl);
1290 while (v < vl)
1291 fprintf(stderr, " %08lx", (unsigned long)*--vl);
1292 fputc('\n', stderr);
1293}
1294
1295static int chkscan(const mpw *v, const mpw *vl,
1296 const void *pp, size_t sz, int step)
1297{
1298 mpscan mps;
1299 const octet *p = pp;
1300 unsigned bit = 0;
1301 int ok = 1;
1302
1303 mpscan_initx(&mps, v, vl);
1304 while (sz) {
1305 unsigned x = *p;
1306 int i;
1307 p += step;
1308 for (i = 0; i < 8 && MPSCAN_STEP(&mps); i++) {
1309 if (MPSCAN_BIT(&mps) != (x & 1)) {
1310 fprintf(stderr,
1311 "\n*** error, step %i, bit %u, expected %u, found %u\n",
1312 step, bit, x & 1, MPSCAN_BIT(&mps));
1313 ok = 0;
1314 }
1315 x >>= 1;
1316 bit++;
1317 }
1318 sz--;
1319 }
1320
1321 return (ok);
1322}
1323
1324static int loadstore(dstr *v)
1325{
1326 dstr d = DSTR_INIT;
1327 size_t sz = MPW_RQ(v->len) * 2, diff;
1328 mpw *m, *ml;
1329 int ok = 1;
1330
1331 dstr_ensure(&d, v->len);
1332 m = xmalloc(MPWS(sz));
1333
1334 for (diff = 0; diff < sz; diff += 5) {
1335 size_t oct;
1336
1337 ml = m + sz - diff;
1338
1339 mpx_loadl(m, ml, v->buf, v->len);
1340 if (!chkscan(m, ml, v->buf, v->len, +1))
1341 ok = 0;
1342 MPX_OCTETS(oct, m, ml);
1343 mpx_storel(m, ml, d.buf, d.sz);
1344 if (memcmp(d.buf, v->buf, oct) != 0) {
1345 dumpbits("\n*** storel failed", d.buf, d.sz);
1346 ok = 0;
1347 }
1348
1349 mpx_loadb(m, ml, v->buf, v->len);
1350 if (!chkscan(m, ml, v->buf + v->len - 1, v->len, -1))
1351 ok = 0;
1352 MPX_OCTETS(oct, m, ml);
1353 mpx_storeb(m, ml, d.buf, d.sz);
1354 if (memcmp(d.buf + d.sz - oct, v->buf + v->len - oct, oct) != 0) {
1355 dumpbits("\n*** storeb failed", d.buf, d.sz);
1356 ok = 0;
1357 }
1358 }
1359
1360 if (!ok)
1361 dumpbits("input data", v->buf, v->len);
1362
12ed8a1f 1363 xfree(m);
42684bdb 1364 dstr_destroy(&d);
1365 return (ok);
1366}
1367
f09e814a 1368static int twocl(dstr *v)
1369{
1370 dstr d = DSTR_INIT;
1371 mpw *m, *ml;
1372 size_t sz;
1373 int ok = 1;
1374
1375 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1376 dstr_ensure(&d, sz);
1377
1378 sz = MPW_RQ(sz);
1379 m = xmalloc(MPWS(sz));
1380 ml = m + sz;
1381
1382 mpx_loadl(m, ml, v[0].buf, v[0].len);
1383 mpx_storel2cn(m, ml, d.buf, v[1].len);
1384 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1385 dumpbits("\n*** storel2cn failed", d.buf, v[1].len);
1386 ok = 0;
1387 }
1388
1389 mpx_loadl2cn(m, ml, v[1].buf, v[1].len);
1390 mpx_storel(m, ml, d.buf, v[0].len);
1391 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1392 dumpbits("\n*** loadl2cn failed", d.buf, v[0].len);
1393 ok = 0;
1394 }
1395
1396 if (!ok) {
1397 dumpbits("pos", v[0].buf, v[0].len);
1398 dumpbits("neg", v[1].buf, v[1].len);
1399 }
1400
12ed8a1f 1401 xfree(m);
f09e814a 1402 dstr_destroy(&d);
1403
1404 return (ok);
1405}
1406
1407static int twocb(dstr *v)
1408{
1409 dstr d = DSTR_INIT;
1410 mpw *m, *ml;
1411 size_t sz;
1412 int ok = 1;
1413
1414 sz = v[0].len; if (v[1].len > sz) sz = v[1].len;
1415 dstr_ensure(&d, sz);
1416
1417 sz = MPW_RQ(sz);
1418 m = xmalloc(MPWS(sz));
1419 ml = m + sz;
1420
1421 mpx_loadb(m, ml, v[0].buf, v[0].len);
1422 mpx_storeb2cn(m, ml, d.buf, v[1].len);
1423 if (memcmp(d.buf, v[1].buf, v[1].len)) {
1424 dumpbits("\n*** storeb2cn failed", d.buf, v[1].len);
1425 ok = 0;
1426 }
1427
1428 mpx_loadb2cn(m, ml, v[1].buf, v[1].len);
1429 mpx_storeb(m, ml, d.buf, v[0].len);
1430 if (memcmp(d.buf, v[0].buf, v[0].len)) {
1431 dumpbits("\n*** loadb2cn failed", d.buf, v[0].len);
1432 ok = 0;
1433 }
1434
1435 if (!ok) {
1436 dumpbits("pos", v[0].buf, v[0].len);
1437 dumpbits("neg", v[1].buf, v[1].len);
1438 }
1439
12ed8a1f 1440 xfree(m);
f09e814a 1441 dstr_destroy(&d);
1442
1443 return (ok);
1444}
1445
42684bdb 1446static int lsl(dstr *v)
1447{
1448 mpw *a, *al;
1449 int n = *(int *)v[1].buf;
1450 mpw *c, *cl;
1451 mpw *d, *dl;
1452 int ok = 1;
1453
1454 LOAD(a, al, &v[0]);
1455 LOAD(c, cl, &v[2]);
1456 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1457
1458 mpx_lsl(d, dl, a, al, n);
1a05a8ef 1459 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1460 fprintf(stderr, "\n*** lsl(%i) failed\n", n);
45c0fd36 1461 dumpmp(" a", a, al);
42684bdb 1462 dumpmp("expected", c, cl);
1463 dumpmp(" result", d, dl);
1464 ok = 0;
1465 }
1466
12ed8a1f 1467 xfree(a); xfree(c); xfree(d);
42684bdb 1468 return (ok);
1469}
1470
81578196 1471static int lslc(dstr *v)
1472{
1473 mpw *a, *al;
1474 int n = *(int *)v[1].buf;
1475 mpw *c, *cl;
1476 mpw *d, *dl;
1477 int ok = 1;
1478
1479 LOAD(a, al, &v[0]);
1480 LOAD(c, cl, &v[2]);
1481 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS);
1482
1483 mpx_lslc(d, dl, a, al, n);
1484 if (!mpx_ueq(d, dl, c, cl)) {
1485 fprintf(stderr, "\n*** lslc(%i) failed\n", n);
45c0fd36 1486 dumpmp(" a", a, al);
81578196 1487 dumpmp("expected", c, cl);
1488 dumpmp(" result", d, dl);
1489 ok = 0;
1490 }
1491
12ed8a1f 1492 xfree(a); xfree(c); xfree(d);
81578196 1493 return (ok);
1494}
1495
42684bdb 1496static int lsr(dstr *v)
1497{
1498 mpw *a, *al;
1499 int n = *(int *)v[1].buf;
1500 mpw *c, *cl;
1501 mpw *d, *dl;
1502 int ok = 1;
1503
1504 LOAD(a, al, &v[0]);
1505 LOAD(c, cl, &v[2]);
1506 ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS + 1);
1507
1508 mpx_lsr(d, dl, a, al, n);
1a05a8ef 1509 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1510 fprintf(stderr, "\n*** lsr(%i) failed\n", n);
45c0fd36 1511 dumpmp(" a", a, al);
42684bdb 1512 dumpmp("expected", c, cl);
1513 dumpmp(" result", d, dl);
1514 ok = 0;
1515 }
1516
12ed8a1f 1517 xfree(a); xfree(c); xfree(d);
42684bdb 1518 return (ok);
1519}
1520
1521static int uadd(dstr *v)
1522{
1523 mpw *a, *al;
1524 mpw *b, *bl;
1525 mpw *c, *cl;
1526 mpw *d, *dl;
1527 int ok = 1;
1528
1529 LOAD(a, al, &v[0]);
1530 LOAD(b, bl, &v[1]);
1531 LOAD(c, cl, &v[2]);
1532 ALLOC(d, dl, MAX(al - a, bl - b) + 1);
1533
1534 mpx_uadd(d, dl, a, al, b, bl);
1a05a8ef 1535 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1536 fprintf(stderr, "\n*** uadd failed\n");
45c0fd36
MW
1537 dumpmp(" a", a, al);
1538 dumpmp(" b", b, bl);
42684bdb 1539 dumpmp("expected", c, cl);
1540 dumpmp(" result", d, dl);
1541 ok = 0;
1542 }
1543
12ed8a1f 1544 xfree(a); xfree(b); xfree(c); xfree(d);
42684bdb 1545 return (ok);
1546}
1547
1548static int usub(dstr *v)
1549{
1550 mpw *a, *al;
1551 mpw *b, *bl;
1552 mpw *c, *cl;
1553 mpw *d, *dl;
1554 int ok = 1;
1555
1556 LOAD(a, al, &v[0]);
1557 LOAD(b, bl, &v[1]);
1558 LOAD(c, cl, &v[2]);
1559 ALLOC(d, dl, al - a);
1560
1561 mpx_usub(d, dl, a, al, b, bl);
1a05a8ef 1562 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1563 fprintf(stderr, "\n*** usub failed\n");
45c0fd36
MW
1564 dumpmp(" a", a, al);
1565 dumpmp(" b", b, bl);
42684bdb 1566 dumpmp("expected", c, cl);
1567 dumpmp(" result", d, dl);
1568 ok = 0;
1569 }
1570
12ed8a1f 1571 xfree(a); xfree(b); xfree(c); xfree(d);
42684bdb 1572 return (ok);
1573}
1574
1575static int umul(dstr *v)
1576{
1577 mpw *a, *al;
1578 mpw *b, *bl;
1579 mpw *c, *cl;
1580 mpw *d, *dl;
1581 int ok = 1;
1582
1583 LOAD(a, al, &v[0]);
1584 LOAD(b, bl, &v[1]);
1585 LOAD(c, cl, &v[2]);
1586 ALLOC(d, dl, (al - a) + (bl - b));
1587
1588 mpx_umul(d, dl, a, al, b, bl);
1a05a8ef 1589 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1590 fprintf(stderr, "\n*** umul failed\n");
45c0fd36
MW
1591 dumpmp(" a", a, al);
1592 dumpmp(" b", b, bl);
42684bdb 1593 dumpmp("expected", c, cl);
1594 dumpmp(" result", d, dl);
1595 ok = 0;
1596 }
1597
12ed8a1f 1598 xfree(a); xfree(b); xfree(c); xfree(d);
42684bdb 1599 return (ok);
1600}
1601
1602static int usqr(dstr *v)
1603{
1604 mpw *a, *al;
1605 mpw *c, *cl;
1606 mpw *d, *dl;
1607 int ok = 1;
1608
1609 LOAD(a, al, &v[0]);
1610 LOAD(c, cl, &v[1]);
1611 ALLOC(d, dl, 2 * (al - a));
1612
1613 mpx_usqr(d, dl, a, al);
1a05a8ef 1614 if (!mpx_ueq(d, dl, c, cl)) {
42684bdb 1615 fprintf(stderr, "\n*** usqr failed\n");
45c0fd36 1616 dumpmp(" a", a, al);
42684bdb 1617 dumpmp("expected", c, cl);
1618 dumpmp(" result", d, dl);
1619 ok = 0;
1620 }
1621
12ed8a1f 1622 xfree(a); xfree(c); xfree(d);
42684bdb 1623 return (ok);
1624}
1625
1626static int udiv(dstr *v)
1627{
1628 mpw *a, *al;
1629 mpw *b, *bl;
1630 mpw *q, *ql;
1631 mpw *r, *rl;
1632 mpw *qq, *qql;
1633 mpw *s, *sl;
1634 int ok = 1;
1635
1636 ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
1637 LOAD(b, bl, &v[1]);
1638 LOAD(q, ql, &v[2]);
1639 LOAD(r, rl, &v[3]);
1640 ALLOC(qq, qql, al - a);
1641 ALLOC(s, sl, (bl - b) + 1);
1642
1643 mpx_udiv(qq, qql, a, al, b, bl, s, sl);
1a05a8ef 1644 if (!mpx_ueq(qq, qql, q, ql) ||
1645 !mpx_ueq(a, al, r, rl)) {
42684bdb 1646 fprintf(stderr, "\n*** udiv failed\n");
1647 dumpmp(" divisor", b, bl);
1648 dumpmp("expect r", r, rl);
1649 dumpmp("result r", a, al);
1650 dumpmp("expect q", q, ql);
1651 dumpmp("result q", qq, qql);
1652 ok = 0;
1653 }
1654
12ed8a1f 1655 xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq);
42684bdb 1656 return (ok);
1657}
1658
1659static test_chunk defs[] = {
1660 { "load-store", loadstore, { &type_hex, 0 } },
f09e814a 1661 { "2cl", twocl, { &type_hex, &type_hex, } },
1662 { "2cb", twocb, { &type_hex, &type_hex, } },
42684bdb 1663 { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } },
81578196 1664 { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } },
42684bdb 1665 { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } },
1666 { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } },
1667 { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },
1668 { "umul", umul, { &type_hex, &type_hex, &type_hex, 0 } },
1669 { "usqr", usqr, { &type_hex, &type_hex, 0 } },
1670 { "udiv", udiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
1671 { 0, 0, { 0 } }
1672};
1673
1674int main(int argc, char *argv[])
1675{
0f00dc4c 1676 test_run(argc, argv, defs, SRCDIR"/t/mpx");
42684bdb 1677 return (0);
1678}
1679
42684bdb 1680#endif
1681
d03ab969 1682/*----- That's all, folks -------------------------------------------------*/